Hi, I have bumped into a problem when trying to solve a linear problem using solve() - well two different solve() routines. Just for the sake of clarity I will present the whole code (so you can run it).
Here I have the UFL formulation of problem (the most important are the last two lines defining the forms "a" and "L")
# Define function spaces (P2-P1)
V = VectorElement('DG',triangle, 2)
Q = FiniteElement('DG',triangle, 1)
W = V*Q
# Define trial and test functions
(u,p) = TrialFunctions(W)
(v,q) = TestFunctions(W)
# Define coefficients
dt = Constant(triangle)
u1 = Coefficient(V)
f = Coefficient(V)
sigma = 10.0
nu = 0.01
# Define bilinear and linear forms in variational formulation
n = FacetNormal(triangle)
F = 0.5*nu*(grad(u)+grad(u).T) # Symmetric part of gradient
def a_(u,v):
M = inner(F,grad(v))*dx - inner(avg(F)*n('+'),jump(v))*dS
return M
def J(u,v):
M = sigma*inner(jump(u),jump(v))*dS
return M
def b(p,v):
M = -p*div(v)*dx + avg(p)*dot(jump(v),n('+'))*dS
return M
def c(u,u_,v):
M = -0.5*inner(grad(v)*u,u_)*dx
return M
def _L(v):
M = inner(f,v)*dx
return M
a = (1/dt)*inner(u, v)*dx + a_(u,v) + J(u,v) + b(p,v) + b(q,u) + c(u,u1,v)
L = (1/dt)*inner(u1, v)*dx + _L(v)
From this file I generate the header file (with ffc -l dolfin NavierStokes2d.ufl command). This file is then included into the main.cpp (I present here just the core):
#include "dolfin.h"
#include "NavierStokes2d.h"
#define RIGHT 5.0
#define TOP 1.0
using namespace dolfin;
// Define noslip domain
class NoslipDomain : public SubDomain
{ ... };
// Define inflow domain
class InflowDomain : public SubDomain
{ ... };
// Define outflow domain
class OutflowDomain : public SubDomain
{ ... };
// Define velocity boundary value at inflow
class InflowVelocity : public Expression
{ ... };
// Define noslip condition for velocity
class NoSlipCondition : public Expression
{ ... };
class InitialVelocity : public Expression
{ ... };
int main()
{
parameters["ghost_mode"] = "shared_facet";
RectangleMesh mesh(Point(0.0,0.0),Point(RIGHT,TOP),(int)10*RIGHT,(int)10*TOP,"crossed");
NavierStokes2d::FunctionSpace W(mesh);
SubSpace V = SubSpace(W,0);
SubSpace Q = SubSpace(W,1);
// Set parameter values
double dt = 1;
double T = 5;
// Define values for boundary conditions
InflowVelocity u_in;
NoSlipCondition u_noslip;
// Define subdomains for boundary conditions
NoslipDomain noslip_domain;
InflowDomain inflow_domain;
OutflowDomain outflow_domain;
// Define boundary conditions
DirichletBC noslip(V, u_noslip, noslip_domain,"pointwise");
DirichletBC inflow(V, u_in, inflow_domain,"pointwise");
std::vector<const DirichletBC*> bcu = {{&noslip,&inflow}};
// Create functions
Function w(W);
// Create coefficients
Constant k(dt);
Constant f(0, 0);
InitialVelocity u0;
// Create forms
NavierStokes2d::BilinearForm a(W, W);
NavierStokes2d::LinearForm L(W);
// Set coefficients
a.dt = k; L.dt = k;
a.u1 = w[0]; L.u1 = w[0]; L.f = f;
Matrix A;
Vector b;
// Use amg preconditioner if available
const std::string prec(has_krylov_solver_preconditioner("amg") ? "amg" : "default");
double t = dt;
while (t < T + DOLFIN_EPS)
{
// Compute
begin("Computing the solution");
if (true) {
assemble(b, L);
assemble(A, a);
for(size_t i=0;i<bcu.size();i++)
bcu[i]->apply(A,b);
solve(A, *w.vector(), b, "gmres", prec);
} else {
solve(a == L,w,bcu);
}
end();
t += dt;
}
plot(w[1], "Pressure", "color");
plot(w[0], "Velocity", "color");
interactive();
return 0;
}
Again, ignore most of the code - just at the end look at the "Compute part". Here is the procedure of assembling the matrix "A" and rhs vector "b". application of dirichlet boundary conditions and calling of the solve() routine for linear system Ax=b. Just under this is the alternative solution with different solve() routine.
The problem is that the two results of these two different solve() approaches are different and I cannot see why.
The "direct linear algebra" approach results in a solution with zero pressure (it looks like the pressure is somehow ignored in the computation).
The "equation form" approach computes the pressure and even though the results are bullshit they are much closer to reality (the bad results I hope are caused by wrong numerical flux approximation of convective term - I have problem with proper upwinding formulation, but that is a different story).
/////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
Well, if you would like to help me even with additional minor problems you can choose from the following questions of mine:
q1) Can I use the definition of V and Q in main.cpp? Because I declare mixedFunctionSpace W by the line NavierStokes2d::FunctionSpace W(mesh) but I cannot see how is secured that it gives me the W space and not any of the other two spaces that are created.
q2) I found all the dolfin header files in usr/include/dolfin. That helps a lot. But I would like to find even the source files .cpp ... I tried the "find" command in command line but it didn't yield a match.
q3) I defined a subclass "InitialVelocity" of Expression class. I defined u0 as InitialVelocity and found in Function.h that there is an operator "=" to assign Expression to a function. Nevertheless, when I try u = u0 (where Function& u = w[0] where Function w(W)) I get the error
Unable to interpolate function into function space.
*** Reason: Wrong size of vector.
There might be many more questions but I guess its enough for now. Thanks ahead for any help or suggestions.
