Convergence of Stokes with data from velocity of Navier-Stokes

variational formulation
1

Hello everyone. I am solving the Stokes equation where the data f=f(u), being u the P1 interpolated analytic velocity from a Navier-Stokes equation. I have problem with the convergence order. I have tried to solve using a lagrangian, impposing the zero mean for the pressure but I cannot get the order 2 corresponding. Any help will be welcome.

#Code
from dolfin import *
#from matplotlib import pyplot as plt
from math import log as ln
import math
import numpy as np


parameters["allow_extrapolation"] = True
set_log_active(True)


num = 4;  h = []; Ldos = []; logerror = []; 
logh = [];  nn = []; bo = [];
logerror.append(0.0); 
Lx = 1.0; Ly = 1.0

for dk in range(num):
    d = 0.2/pow(2,dk)
    mesh = RectangleMesh(Point(0.0, 0.0), Point(Lx, Ly), int(Lx/d), int(Ly/d))
    hh = mesh.hmax()
    
    #Defining the boundaries
    boundary = 'near(x[0], 0.0) || near(x[0], 1.0) || near(x[1], 0.0) || near(x[1], 1.0)'
# Define function spaces
    k = 2
    H = VectorElement('P', mesh.ufl_cell(), k) #Velocity
    Q = FiniteElement('P', mesh.ufl_cell(), k-1) #Presure
    #UU = VectorElement('P', mesh.ufl_cell(), k)
    UU = MixedElement([Q, Q])
    R = FunctionSpace(mesh, UU*Q)
    W = FunctionSpace(mesh, H*Q)
    
    #Phisical Constans
    mu = Constant(0.035)
    rho = Constant(1.2)
    nu = 1 #mu/rho
    hK = CellDiameter(mesh)
    hh = mesh.hmax()
    h.append(hh)
    nn.append(W.dim())  

   #Initial Data from velocity of NS
    a = '-256*pow(x[0],2)*pow(x[0]-1,2)*x[1]*(x[1]-1)*(2*x[1]-1)'
    b = '256*pow(x[1],2)*pow(x[1]-1,2)*x[0]*(x[0]-1)*(2*x[0]-1)'
    u_ex = Expression((a,b), domain=mesh, degree=7)
    u = interpolate(u_ex, R.sub(0).collapse())
  
    p_ex = Expression('150*(x[0]-0.5)*(x[1]-0.5)', domain=mesh, degree=2)
    p_in = interpolate(p_ex,W.sub(1).collapse())
    p_avg = Constant(assemble(p_ex*dx)/assemble(1.0*dx(domain=mesh)))   

  # Inflow boundary condition for velocity
    boundP = ('0.0', '0.0')
    bound = Expression(boundP, degree=1)
    bcs = DirichletBC(W.sub(0), bound, boundary)

    # Define variational problem
    (w, q) = TrialFunctions(W)
    (v, r) = TestFunctions(W)

    f = -nu*div(grad(u_ex))+grad(u_ex)*u_ex+grad(p_ex)
    L = (inner(grad(w), grad(v)) - div(v)*(q-p_avg) + r*div(w))*dx
    R = inner(grad(u)*u, v)*dx - nu*inner(grad(u), grad(v))*dx\
        +inner(f,v)*dx 

    # Compute solution
    U = Function(W)
    solver=solve(L == R, U, bcs)

   # save solutions
    (w, q) = U.split(True); w.rename("velocity","velocity"); q.rename("pressure","pressure");
    q_avg = Constant(assemble(q*dx)/assemble(1.0*dx(domain=mesh)))

   # Explicit computation of L2 norm 
    error = (p_ex  -q )**2*dx # 
    errorL2 = sqrt(abs(assemble(error)))    
    Ldos.append(errorL2)
 
   #Logarithmic rates
    if(dk>0):
      logerror.append(ln(Ldos[dk]/Ldos[dk-1])/(ln(float(h[dk])/h[dk-1])))

# ********  Generating table of convergences **** #

print( 'h  &  ||q-q_h||_0 &  & Log error )
print('===============================================================')
for i in range(num):
    print('%4.4g & %4.4g & %4.4g  '% (h[i], Ldos[i], logerror[i]))
2

Hi,
Please format your code properly by encapsulating it with ```, and make sure it is executable by others by copy-paste. Also remove all lines of code not required for the question, as specified in: Read before posting: How do I get my question answered? (This includes unnecessary imports, lines that are deactivated by commenting etc.)

To set the pressure average to zero you can also just solve without a Lagrange multipler, and subtract the average (as you are computing as q_avg) after solving.

Another way of enforcing a pressure condition is to apply a pointwise Dirichlet boundary condition of the pressure, where you apply the exact pressure solution at one point. See for instance: Pointsource on specific points python.

3

Thanks. I will try the last idea, because the second on I have already tried and doesn’t work because It isn’t a tradition Stokes problem.

4

Well, you are not supposed to subtract the average until after you have solved the problem. I.e. Solve stokes with standard force, then do q-=q_avg