Hi
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.
Thanks
Hatef
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//Noslip function
class VelNoslip : public Function
{
public:
void eval(real* values, const real* x) const
{
values[0]=0.0 ;
values[1]=0.0 ;
}
};
////////////////////////////////////////////////////////////////////////////////////
//-Noslip BC subdomain
class Noslip_boundary_end : public SubDomain
{
bool inside(const real* x, bool on_boundary) const
{
return
(x[0]>1.0-DOLFIN_EPS) && on_boundary;
}
};
class Noslip_boundary_side : public SubDomain
{
bool inside(const real* x, bool on_boundary) const
{
return
(x[0]>DOLFIN_EPS && x[0]<1.0-DOLFIN_EPS) && on_boundary;
}
};
// define subdomain for inflow
class Inflow_boundary : public SubDomain
{
bool inside(const real* x, bool on_boundary) const
{
return
x[0]< DOLFIN_EPS && on_boundary;
}
};
////////////////////////////////////////////////////////////////////////////////////////////
int main()
{
//Read Mesh and Subdomain marker and define function space
UnitSquare mesh(32,32);
//Define fuction space
NSWEFunctionSpace V(mesh);
SubSpace velocity(V,0);
SubSpace pressure(V,1);
MeshFunction <unsigned int> sub_domains(mesh, mesh.topology().dim() -1);
// Mark all facets as sub domain 3
sub_domains = 3;
// Mark end Noslip subdomain
Noslip_boundary_end noslip_boundary_end;
noslip_boundary_end.mark(sub_domains,0);
// mark side noslip subdomain
Noslip_boundary_side noslip_boundary_side;
noslip_boundary_side.mark(sub_domains,1);
//mark Inflow subdomain
Inflow_boundary inflow_boundary;
inflow_boundary.mark(sub_domains,2);
//pseudo time
real t = 0.0;
//solution vector
Function v(V);
Function U(velocity);
Function P(pressure);
Function Un(velocity);
Un.vector().zero();
Function Pn(pressure);
Pn.vector().zero();
Constant f(2,0.0);
VelNoslip velnoslip;
Constant inflow(0.1);
//set up Forms
NSWEBilinearForm a(V,V) ;
a.U0=U;
NSWELinearForm L(V);
L.Un=Un; L.Pn=Pn ; L.F=f; L.U0=U; L.P0=P;
//set up boundary condition
DirichletBC bc0(velocity, velnoslip,sub_domains,0);
DirichletBC bc1(velocity, velnoslip,sub_domains,1);
DirichletBC bc2(velocity, velnoslip,sub_domains,2);
DirichletBC bc3(pressure, inflow,sub_domains,2);
//collect all the BCs
Array<BoundaryCondition*> bcs(&bc0,&bc1,&bc2,&bc3);
//parameters for time stepping
real T = 0.005;
real k = 0.001;
File file_2("Pressure.pvd");
File file_1("Velocity.pvd");
//time stepping
Progress prg("Time-stepping");
while ( t < T )
{
VariationalProblem nlpde(a,L,bcs,true);
nlpde.solve(v);
U = v[0];
P = v[1];
// Save the solution to file
file_2 << P;
file_1 << U;
// Move to the next interval
Un.vector()=U.vector() ;
Pn.vector()=P.vector() ;
prg = t / T;
t += k;
}
}
_______________________________________________ DOLFIN-dev mailing list [email protected] http://www.fenics.org/mailman/listinfo/dolfin-dev
