Possible error in assembling with Lagrange multipliers on boundaries.

Question

Hallo
I'm making some tests on the simple problem of the flow around a cylinder with zero tangential stress condition using Newton iterations. I'm imposing null normal velocity on the cylinder and assigned inflow velocity using Lagrange multipliers, as in the code below.

Ppressure = FiniteElement("Lagrange",triangle,1)
Pvelocity = VectorElement("Lagrange",triangle,2)
Plegrange = VectorElement("Lagrange",triangle,2)

V = MixedElement([Ppressure,Pvelocity,Plagrange])

n = FacetNormal(triangle)

(p, u, l) = TrialFunctions(V)
(q, v, m) = TestFunctions(V)
Re = Constant(triangle)
U0 = Coefficient(Pvelocity)
Uinflow = Coefficient(Pvelocity)

A = (1/Re*inner(grad(v), grad(u))+1/Re*inner(grad(v), transpose(grad(u))) - div(v)*p -    q*div(u))*dx  \
+inner(v,grad(U0)*u)*dx + inner(v,grad(u)*U0)*dx  -dot(v,l)*ds(3) -dot(m, u)*ds(3) \
-dot(m[0], dot(u,n))*ds(2)- l[0]*dot(v,n)*ds(2)

L = + inner(v,grad(U0)*U0)*dx-dot(m, Uinflow)*ds(3)

forms = [A,L]

With the elements listed in the code I get the right solution.
If I try to use MINI elements, so P1-P1bubble, leaving P2 for the multiplier I get nonsense solutions using UMFPACK and this error using MUMPS, that is singularity in the first line.

[0]PETSC ERROR: --------------------- Error Message ------------------------------------
[0]PETSC ERROR: Error in external library!
[0]PETSC ERROR: Error reported by MUMPS in numerical factorization phase: INFO(1)=-10, INFO(2)=0

If I use P1 elements for Lagrange multipliers too the solution works but the results have obviously poor accuracy.

So P1-P2-P2 works, P1-P1bubble-P1 works, P1-P1bubble-P2 gives error.

My question is if this error is related to the assembling procedure with this combination of elements. And in case if there is an enrichment that affects also the facets, to improve the accuracy even using P1 Lagrange multipliers.

Thanks for any suggestion.
Stefano

Comments

Are you sure that the combination P1-P1bubble-P2 is uniformly inf-sup stable? And what happens to the DOFs inside the domain? Can you please provide a working example?

---

I prepared an example using the demos. Thank you very much.

Answer 1

I'm working in c++. I prepared an example using parts of the demos. You can take the mesh and CMakeList files from the demo documented/subdomains.
The UFL file, named LinearizedNS.ufl is

####### Navier Stokes formulation with Lagrange multipliers on boundary.
Ppressure = FiniteElement("Lagrange",triangle,1)

#### Velocity P2, Lagrange P2.
####Works.#####
Pvelocity = VectorElement("Lagrange",triangle,2)
Plagrange = VectorElement("Lagrange",triangle,2)

#### Velocity P1bubble, Lagrange P1.
#####Works with poor accuracy(oscillations in presure field).#####
#Pbubble = VectorElement("Bubble",triangle,3)
#P1 = VectorElement("Lagrange",triangle,1)
#Pvelocity = P1+Pbubble
#Plagrange = VectorElement("Lagrange",triangle,1)

#### Velocity P1bubble, Lagrange P2.             
#####Does not work.######
#Pbubble = VectorElement("Bubble",triangle,3)
#P1 = VectorElement("Lagrange",triangle,1)
#Pvelocity = P1+Pbubble
#Plagrange = VectorElement("Lagrange",triangle,2)

V = MixedElement([Ppressure,Pvelocity,Plagrange])

n = FacetNormal(triangle)

(p, u, l) = TrialFunctions(V)
(q, v, m) = TestFunctions(V)

Re = Constant(triangle)

U = Coefficient(Pvelocity)
U0 = Coefficient(Pvelocity)

Uinflow = Coefficient(Pvelocity)

A = (1/Re*inner(grad(v), grad(u))+1/Re*inner(grad(v), transpose(grad(u))) - div(v)*p -  q*div(u))*dx  \
+inner(v,grad(U0)*u)*dx + inner(v,grad(u)*U0)*dx  -dot(v,l)*ds(1) -dot(m, u)*ds(1) \
-dot(m[0], dot(u,n))*ds(2)- l[0]*dot(v,n)*ds(2)

L = + inner(v,grad(U0)*U0)*dx-dot(m, Uinflow)*ds(1)

increment = dot((U-U0),(U-U0))*dx

forms=[A, L, increment]

Just uncomment the set of elements you want of the three implemented.
The main.cpp program.

#include <dolfin.h>
#include "LinearizedNS.h"

using namespace dolfin;

// Inflow condition.
class PrescribedValue : public Expression
{
     public: PrescribedValue() : Expression(2) {}
     void eval(Array<double>& values, const Array<double>& x) const{
             values[0]= -1.0;
             values[1]= 0.0;
      }
};

// Sub domain for no-slip (everything except inflow and outflow)
class Noslip : public SubDomain
{
      bool inside(const Array<double>& x, bool on_boundary) const
      {
          return on_boundary;
       }
};

// Sub domain for inflow (right)
class Inflow : public SubDomain
{
      bool inside(const Array<double>& x, bool on_boundary) const
      {
        return x[0] > 1.0 - DOLFIN_EPS && on_boundary;
      }
};

// Sub domain for outflow (left)
class Outflow : public SubDomain
{
      bool inside(const Array<double>& x, bool on_boundary) const
      {
       return x[0] < DOLFIN_EPS && on_boundary;
       }
};

int main()
{
// Read mesh.
Mesh mesh("./dolfin_fine.xml.gz");

// Define paramenters.
double tollerance = 0.0000001;
double max = 10;
double step;
Constant Re(50.0);

// Define subdomains.
Noslip noslip;
Inflow inflow;
Outflow outflow;

// Mark subdomains for boundary conditions.
MeshFunction<std::size_t> sub_domains(mesh, mesh.topology().dim() - 1);

sub_domains.set_all(0);
noslip.mark(sub_domains,2);
inflow.mark(sub_domains,1);
outflow.mark(sub_domains,0);

// Define finite element space.
LinearizedNS::FunctionSpace V(mesh);

// Define function for solution at two iterations.
Function U(V);
Function U0(V);

//// Define forms.
// Bilinear form.
LinearizedNS::Form_A A(V,V);
// Linear form.
LinearizedNS::Form_L L(V);
// Increment form for stop criteria.
LinearizedNS::Form_increment increment(mesh);

// Define matrix and vector.
PETScMatrix K;
PETScVector b;

// Inflow condition.
PrescribedValue vinflow;

//// First step. Stokes.
// Coefficients assigned to forms.
A.Re = Re;
A.U0 = U[1];
A.ds = sub_domains;
L.U0 =U[1];
L.ds = sub_domains;
L.Uinflow = vinflow;

// Assembly of matrix and right hand sise term.
assemble(K, A);
assemble(b, L);

// De singularization of matrix.
K.ident_zeros();

// Solution.
solve(K,*U.vector(),b,"lu" );

// Calculation of step.
increment.U = U[1];
increment.U0 = U0[1];
step = assemble(increment);

// Update the solution.
U0 = U;
printf("(U-U0)L2= %e\n", step);

//// Newton iterations.
int iterations = 1;
while (step > tollerance and iterations < max) {

A.U0 = U0[1];
L.U0 =U0[1];

assemble(K, A);
assemble(b, L);

K.ident_zeros();

solve(K,*U.vector(),b);

increment.U = U[1];
increment.U0 = U0[1];
step = assemble(increment);
printf("iterations= %i\n", iterations);
printf("(U-U0)L2= %e\n", step);

U0 = U;

iterations = iterations +1;
};

// Separate functions for plot.
Function u = U[1];
Function p = U[0];

// Save solution in VTK format
File ufile_pvd("velocity.pvd");
ufile_pvd << u;
File pfile_pvd("pressure.pvd");
pfile_pvd << p;


// Plot solution
plot(u);
plot(p);
interactive();

return 1;
};

This combination should be stable. Anyway I don't get spurious solutions, but the matrix itself seems wrong. More with the less stable P1-P1bubble-P1 the solution works.

The internal degrees of freedom are put to zero with the ident_zeros() function.

If there is a simpler way to upload examples and codes, I don't know it, sorry.
Thank you very much.