Hi Hatef, It sounds as if you are interested in the Navier-Stokes equations, or similar? If so, you may want to check out the Unicorn solver: http://www.fenics.org/wiki/Unicorn
Direct any Unicorn questions to: [email protected] Best, Johan > <p>Hi</p><p>I installed the new version of dolfin (0.9.0) as you wanted > and am trying to model a coupled set of PDEs which are nonlinear. As I > told in my previous emails, there are currenlty not much documents on > nonlinear problems. so based on what I have seen on the demos I came up > with the following code which does not seem to be working. I would > appreciate if anyone could suggest any material or hint on modeling these > kind of problems in which there are more that one variable. in my case > there is a vector variable (Velocity) and a scalar one (Pressure) and > Nonlinearity is only due to the term (u.grad)u and I am trying to use > Newton's solver for nonlinear iterations.</p><p><br /><br > />Thanks</p><p>Hatef</p><p>//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////<br > />//Noslip function<br /><br />class VelNoslip : public Function<br > /> {<br /> public:<br /><br /> void > eval(real* values, const real* x) const<br /> {<br > /> <br > /> > values[0]=0.0 ;<br > /> > values[1]=0.0 ;<br /> <br > /> }<br /><br > /> };<br /><br > />////////////////////////////////////////////////////////////////////////////////////<br > />//-Noslip BC subdomain<br />class Noslip_boundary_end : public > SubDomain<br /> {<br /> bool inside(const real* x, > bool on_boundary) const<br /> {<br > /> return <br > /> > (x[0]>1.0-DOLFIN_EPS) && on_boundary;<br > /> <br /> }<br /> > };<br /><br />class Noslip_boundary_side : public SubDomain<br /> > {<br /> bool inside(const real* x, bool on_boundary) > const<br /> {<br /> return > <br > /> > (x[0]>DOLFIN_EPS && x[0]<1.0-DOLFIN_EPS) && > on_boundary;<br /> <br /> > }<br /> };<br /><br /><br />// define subdomain for inflow<br > />class Inflow_boundary : public SubDomain<br /> {<br > /> bool inside(const real* x, bool on_boundary) const<br > /> {<br /> return <br > /> > x[0]< DOLFIN_EPS && on_boundary;<br > /> <br /> }<br /> > };<br /><br > />////////////////////////////////////////////////////////////////////////////////////////////<br > />int main()<br />{<br /><br /><br />//Read Mesh and Subdomain marker and > define function space<br />UnitSquare mesh(32,32);</p><p>//Define > fuction space<br /> NSWEFunctionSpace V(mesh);<br /> SubSpace > velocity(V,0);<br /> SubSpace pressure(V,1);</p><p><br /><br > />MeshFunction <unsigned int> sub_domains(mesh, > mesh.topology().dim() -1);<br /><br /><br />// Mark all facets as sub > domain 3<br /> sub_domains = 3;<br />// Mark end Noslip subdomain<br > /> Noslip_boundary_end noslip_boundary_end;<br /> > noslip_boundary_end.mark(sub_domains,0);<br /><br /><br />// mark side > noslip subdomain<br /> Noslip_boundary_side noslip_boundary_side;<br > /> noslip_boundary_side.mark(sub_domains,1);<br /><br /><br />//mark > Inflow subdomain<br /> Inflow_boundary inflow_boundary;<br /> > inflow_boundary.mark(sub_domains,2);<br /> //pseudo time<br /> > real t = 0.0;<br /> <br /> //solution vector<br /> > Function v(V);<br /> Function U(velocity);<br /> Function > P(pressure); <br /> <br /> Function Un(velocity);<br /> > Un.vector().zero();<br /> Function Pn(pressure);<br /> > Pn.vector().zero();<br /> <br /> <br /> Constant > f(2,0.0);<br /> VelNoslip velnoslip;<br /> Constant > inflow(0.1);<br /> <br /> <br /><br /> <br /> <br > /> //set up Forms<br /> NSWEBilinearForm a(V,V) ;<br /> > a.U0=U;<br /> NSWELinearForm L(V);<br /> L.Un=Un; L.Pn=Pn ; > L.F=f; L.U0=U; L.P0=P;<br /> <br /> <br /> //set up > boundary condition<br /> <br /> DirichletBC bc0(velocity, > velnoslip,sub_domains,0);<br /> DirichletBC bc1(velocity, > velnoslip,sub_domains,1);<br /> DirichletBC bc2(velocity, > velnoslip,sub_domains,2);<br /> DirichletBC bc3(pressure, > inflow,sub_domains,2);<br /> //collect all the BCs<br /> > Array<BoundaryCondition*> > bcs(&bc0,&bc1,&bc2,&bc3);<br /> <br /> > //parameters for time stepping<br /><br /> real T = 0.005;<br > /> real k = 0.001;<br /> <br /> <br /> <br > /> File file_2("Pressure.pvd");<br /> File > file_1("Velocity.pvd");<br /><br /> <br /><br /> > //time stepping<br /> Progress prg("Time-stepping");<br > /> while ( t < T )<br /> { <br /> <br /> > VariationalProblem nlpde(a,L,bcs,true);<br /> <br /> > nlpde.solve(v);<br /> <br /> U = v[0];<br /> P = > v[1];</p><p><br /> // Save the solution to file<br /> file_2 > << P;<br /> file_1 << U;<br /> <br /> // > Move to the next interval </p><p> Un.vector()=U.vector() ; <br > /> Pn.vector()=P.vector() ; <br /> prg = t / T;<br /> t > += k;<br /> <br /> }<br /><br />}</p> > _______________________________________________ > DOLFIN-dev mailing list > [email protected] > http://www.fenics.org/mailman/listinfo/dolfin-dev > _______________________________________________ DOLFIN-dev mailing list [email protected] http://www.fenics.org/mailman/listinfo/dolfin-dev
