Convergence Order Kovasznay Flow

variational formulation
1

Hi everyone. I’am runing Steady NS using the Kovasznay’s solution, but I have problems with convergence order. I should get a rate of order 1 but I get anything but what I should obtain.

from dolfin import *
import numpy
import math


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

#parameters
Re= 30; k=1; num = 4; 
Lx = 1.5; Ly = 2.0; #d = 0.1
lam = int(Re/2-sqrt(Re*Re/4+4*pi*pi))
nu = 1.0/Re
h = []; E = []  # mesh size and errors

########################################################################################

def sigma(u, p):
    return nu*grad(u) - p*Identity(len(u))

########################################################################################

def compute(dk):
    d = 0.2/pow(2,dk)
    mesh = RectangleMesh(Point(-0.5, 0.0), Point(Lx, Ly), int(Lx/d), int(Ly/d),"crossed")
    sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
    hh = mesh.hmax()
    h.append(hh) 
    # define function spaces
    V = VectorElement('P', mesh.ufl_cell(), k+1) #Velocity
    Q = FiniteElement('P', mesh.ufl_cell(), k) #Presure
    R = FiniteElement('R', mesh.ufl_cell(), k-1) #Lagrange Multiplier
    W = FunctionSpace(mesh, MixedElement([V,Q,R]))


    #exact solution
    u1 = '1-exp(lam*x[0])*cos(2*pi*x[1])'
    u2 = '(lam/(2*pi))*exp(lam*x[0])*sin(2*pi*x[1])'
    u_ex = Expression((u1,u2), domain=mesh, lam = lam, degree=3)
    p1 = '0.5*exp(lam*x[0])-(1/(8*lam))*(exp(3*lam)-exp(-lam))'
    p_ex = Expression(p1, domain=mesh, lam = lam, degree=3)
    #f1 = '(pow(lam,2)-4*pow(pi,2))*exp(lam*x[0])*cos(2*pi*x[1])-lam*exp(2*lam*x[0])'  
    #f2 = '(2*pi*lam-pow(lam,3)/(2*pi))*exp(lam*x[0])*sin(2*pi*x[1])'    
    #f = Expression((f1,f2), domain=mesh, lam = lam, degree=3)     

    bc =  DirichletBC(W.sub(0), u_ex, "on_boundary")

    # test and trial functions
    w     = Function(W)
    (u,p,r) = (as_vector((w[0], w[1])), w[2], w[3])
    (v,q,s) = TestFunctions(W)
    n = FacetNormal(mesh)
    ds = Measure('ds', domain = mesh, subdomain_data = sub_domains)
       
    #Bilinear FOrm
    f = -nu*div(grad(u_ex))+grad(u_ex)*u_ex+grad(p_ex)
    F =   inner(grad(u)*u, v)*dx \
        + nu*inner(grad(u), grad(v))*dx \
        - div(v)*p*dx + p*s*dx +q*r*dx\
        - q*div(u)*dx -inner(f,v)*dx\
        - dot(sigma(u,p)*n,v)*ds
        #+ dot(p*n,v)*ds- dot(nu*nabla_grad(u)*n, v)*ds
        
  

    #Solver
    dw = TrialFunction(W)
    dF = derivative(F, w, dw)
    nsproblem = NonlinearVariationalProblem(F, w, bc, dF)
    solver = NonlinearVariationalSolver(nsproblem)
    prm = solver.parameters
    prm['nonlinear_solver'] = 'snes'
    prm['snes_solver']['line_search'] = 'bt'
    prm['snes_solver']['linear_solver'] = 'mumps'
    prm['snes_solver']['report'] = False
    solver.solve()
    (u,p,r) = w.split()
    u.rename("velocity","velocity"); p.rename("pressure","pressure");
    File("kovasznay/velocity.pvd", "compressed") << u
    File("kovasznay/pressure.pvd", "compressed") << p

    p_i = interpolate(p_ex, W.sub(1).collapse())
    u_i = interpolate(u_ex, W.sub(0).collapse())
    File("kovasznay/velocity_ex.pvd", "compressed") << u_i
    File("kovasznay/pressure_ex.pvd", "compressed") << p_i

    # Compute errors
    #p_i = interpolate(p_ex, W.sub(1).collapse())
    origin = (0.0, 0.0)
    value = p(origin) 
    p_err = sqrt(assemble((p_ex - p)*(p_ex - p)*dx))
    errors = {'p_L2_err': p_err}
    return errors

########################################################################################




for dk in range(num):
  E.append(compute(dk))  # list of dicts
  

# Compute convergence rates
from math import log as ln          
error_types = E[0].keys()
for error_type in sorted(error_types):
  l = len(E)
  print('\nError norm based on', error_type)
  for i in range(0, l):
    if i==0:
      print('h=%.5f Error=%.10f Convergence rate= -- ' % (h[i], E[i][error_type]))
    else:
      r = ln(E[i][error_type]/E[i-1][error_type])/ln(h[i]/h[i-1])
      print('h=%.5f Error=%.10f Convergence rate=%.4f' % (h[i], E[i][error_type], r))

can anyone give me feedback please, because I cannot seeing the problem. Thanks.

2

The problem here is that you are enforcing the constraint that \int_\Omega p = 0 in the discrete solution, while \int_\Omega p_{\text{ex}} \neq 0. This can be resolved by computing

p_avg = Constant(assemble(p_ex*dx)/assemble(1.0*dx(domain=mesh)))

and modifying the constraint equation term in the variational form to

(p-p_avg)*s*dx
3

Ok, I have an improve but still I don’t the correct order which is 1. I don’t understand what happen.

error-convergence

4

Don’t you expect a rate of two here? That would be the expected convergence rate for the P_2–P_1 Taylor–Hood element when looking at the velocity in H^1 and the pressure in L^2.

2 Likes
5

Thanks David. Your help was very useful to me.