Error in Time dependent NS equation

1

Hi,

here u=(u1, u2)
u1= 1+ x^2+y^2+2*t
u2= 1+xy+t^2
f=(f1, f2)
f1=-2 +2*u1*x+2*u2*y +2*cost(2y-1)
f2=2*t +y*u1+x*u2+2*cost(2*x-1)
p=(2x-1)(2y-1)cost

from __future__ import print_function
import matplotlib.pyplot as plt
from fenics import *

def  NS(N):
    T=1.0
    num_steps=pow(N,2)
    dt=T/num_steps
    #Creat Mesh and function spaces
    mesh=UnitSquareMesh(4, 4)
    V= VectorFunctionSpace(mesh, 'P', 2)
    Q= FunctionSpace(mesh,'P', 1)
    # Define boundary conditions
    u_D = Expression(('1 + x[0]*x[0] + x[1]*x[1] + 2*t','1 + x[0]*x[1] + pow(t,2)'),
                 degree=2, t=0)
                      
    bc = DirichletBC(V, u_D , "on_boundary")
    bcu=[bc]
    #  Define trial and test functions
    u = TrialFunction(V)
    v = TestFunction(V)
    p = TrialFunction(Q)
    q = TestFunction(Q)


    # Define functions for solutions at previous and current time stepsfrom _future_import print function
    u_n = interpolate(u_D, V)
   #u_n = Function(V)
    u_  = Function(V)
    p  = Function(Q)
    q   = Function(Q)

    # Define expressions used in variational forms

    f = Expression(('-2+2*(1 + x[0]*x[0] + x[1]*x[1] + 2*t)*x[0]+2*(1 + x[0]*x[1] + pow(t,2)*x[1]) +2*cos(t)*(2*x[1]-1)',\
                '2*t +x[1]*(1 + x[0]*x[1] + pow(t,2)) +x[0]*(1 + x[0]*x[1] + pow(t,2))+2*cos(t)*(2*x[0]-1)'), t=0, degree=2)
    k   = Constant(dt)
    # Define variational problem 
    F = inner((u - u_n) / k, v)*dx - inner(f, v)*dx + inner(grad(u),grad(v))*dx + inner(grad(p),v)*dx \
    +inner(grad(u_n)*u, v)*dx\
    + inner(grad(q), u)*dx

    a = lhs(F)
    L = rhs(F)

    # Time-stepping
    t = 0
    for n in range(num_steps):
        # Update current time
     
        t += dt
        u_D.t=t
        f.t=t
        A=assemble(a)
        [bc.apply(A) for bc in bcu]
        b= assemble(L)
        [bc.apply(b) for bc in bcu]
        solve(A, u_.vector(), b)
        u_n.assign(u_)
      
    error_L2 = errornorm(u_D, u_, 'L2')
    return error_L2

for N in [2, 4, 8, 16, 32 ]:
    print("{0:d} error_L2 ={1:.2e}".format(N, NS(N)))

The Output is

2 error_L2 =5.11e-02
4 error_L2 =5.40e-02
8 error_L2 =5.48e-02
16 error_L2 =5.51e-02
32 error_L2 =5.51e-02

The error should decrease at final time t=1, but here it is not decreasing?

Thank You,
Debendra

2

Assuming your manufactured solution is a solution of the system you want to solve: you’re only doing temporal refinement, your spatial approximation error is limited by the coarseness of your mesh. Refining the mesh too should show convergence.

3

Dear @nate,
Here I want to apply dirichlet boundary condition u= u_D on boundry of unitSquaremesh

from __future__ import print_function
import matplotlib.pyplot as plt
from fenics import *

def  NS(N):
    T=1.0
    num_steps=pow(N,3)
    dt=T/num_steps
    #Creat Mesh and function spaces
    mesh=UnitSquareMesh(N, N)
    V= VectorFunctionSpace(mesh, 'P', 2)
    Q= FunctionSpace(mesh,'P', 1)
    # Define boundary conditions
    u_D = Expression(('1 + x[0]*x[0] + x[1]*x[1] + 2*t','1 + x[0]*x[1] + pow(t,2)'),
                 degree=2, t=0)
                      
    bc = DirichletBC(V, u_D , "on_boundary")
    bcu=[bc]
    #  Define trial and test functions
    u = TrialFunction(V)
    v = TestFunction(V)
    p = TrialFunction(Q)
    q = TestFunction(Q)


    # Define functions for solutions at previous and current time stepsfrom _future_import print function
    u_n = interpolate(u_D, V)
   #u_n = Function(V)
    u_  = Function(V)
    p  = Function(Q)
    q   = Function(Q)

    # Define expressions used in variational forms

    f = Expression(('-2+2*(1 + x[0]*x[0] + x[1]*x[1] + 2*t)*x[0]+2*(1 + x[0]*x[1] + pow(t,2)*x[1]) +2*cos(t)*(2*x[1]-1)',\
                '2*t +x[1]*(1 + x[0]*x[1] + pow(t,2)) +x[0]*(1 + x[0]*x[1] + pow(t,2))+2*cos(t)*(2*x[0]-1)'), t=0, degree=2)
    k   = Constant(dt)
    # Define variational problem 
    F = inner((u - u_n) / k, v)*dx - inner(f, v)*dx + inner(grad(u),grad(v))*dx + inner(grad(p),v)*dx \
    +inner(grad(u_n)*u, v)*dx\
    + inner(grad(q), u)*dx

    a = lhs(F)
    L = rhs(F)

    # Time-stepping
    t = 0
    for n in range(num_steps):
        # Update current time
     
        t += dt
        u_D.t=t
        f.t=t
        A=assemble(a)
        [bc.apply(A) for bc in bcu]
        b= assemble(L)
        [bc.apply(b) for bc in bcu]
        solve(A, u_.vector(), b)
        u_n.assign(u_)
      
    error_L2 = errornorm(u_D, u_, 'L2')
    return error_L2

for N in [2, 4, 8, 16 ]:
    print("{0:d} error_L2 ={1:.2e}".format(N, NS(N)))

The out put is

2 error_L2 =5.21e-02
4 error_L2 =5.48e-02
8 error_L2 =5.51e-02
16 error_L2 =5.51e-02

still it is increasing
here i used spatial refinement , Kindly explain elaborately .
Thank You ,
Debendra

4

Looks like you need to carefully debug your code and/or the manufactured problem. Seems like you’re converging to a solution you’re not expecting.

5

Dear @nate,

ns equation
ns equation869Ă—1280 81.5 KB 

The exact solution is u_D= (tsiny,xcost),
 I want to apply on dirichlet boundry condition , where u= u_D on boundry,
p*dx=0 on /omega,
 Could You suggest me,  whether the   weak formulation is wrong ?
What should I do for this problem ? i want to go on this way .
I do not want to use splitting method.

Kindly, see the code,

from __future__ import print_function
import matplotlib.pyplot as plt
from fenics import *

def  NS(N):
    T=1.0
    num_steps=pow(N,3)
    dt=T/num_steps
    #Creat Mesh and function spaces
    mesh=UnitSquareMesh(N, N)
    V= VectorFunctionSpace(mesh, 'P', 2)
    Q= FunctionSpace(mesh,'P', 1)
    # Define boundary conditions
    u_D = Expression(('t*sin(x[1])','x[0]*cos(t)'), degree=2, t=0)
                 
                      
    bc = DirichletBC(V, u_D , "on_boundary")
    bcu=[bc]
    #  Define trial and test functions
    u = TrialFunction(V)
    v = TestFunction(V)
    p = TrialFunction(Q)
    q = TestFunction(Q)


    # Define functions for solutions at previous and current time stepsfrom _future_import print function
    u_n = interpolate(u_D, V)
   #u_n = Function(V)
    u_  = Function(V)
    p  = Function(Q)
    q   = Function(Q)

    # Define expressions used in variational forms

    f = Expression(('sin(x[1]) + x[0]*t*t*cos(t)+ t*sin(x[1])+2*(2*x[1]-1)',\
                   '-x[0]*sin(t)+t*cos(t)*sin(x[1])+ 2*(2*x[0]-1)'), t=0, degree=2)
    k   = Constant(dt)
    # Define variational problem 
    F = inner((u - u_n) / k, v)*dx - inner(f, v)*dx + inner(grad(u),grad(v))*dx + inner(grad(p),v)*dx \
               +inner(grad(u_n)*u, v)*dx\
               + inner(grad(q), u)*dx

    a = lhs(F)
    L = rhs(F)

    # Time-stepping
    t = 0
    for n in range(num_steps):
        # Update current time
     
        t += dt
        u_D.t=t
        f.t=t
        A=assemble(a)
        [bc.apply(A) for bc in bcu]
        b= assemble(L)
        [bc.apply(b) for bc in bcu]
        solve(A, u_.vector(), b)
        u_n.assign(u_)
      
    error_L2 = errornorm(u_D, u_, 'L2')
    return error_L2

for N in [2, 4, 8, 16 ]:
    print("{0:d} error_L2 ={1:.2e}".format(N, NS(N)))

The out put is

2 error_L2 =2.08e-02
4 error_L2 =2.38e-02
8 error_L2 =2.41e-02
16 error_L2 =2.41e-02

Here error is not decreasing, Divergence of exact solution is zero .

6

I do not understand your variational formulation. You do not use test and trial functions for p and q, i.e.

You are also not using a splitting scheme, but you creating separate function spaces for u and p, which does not make sense.

I would strongly suggest you have a look at: https://link.springer.com/content/pdf/10.1007%2F978-3-319-52462-7.pdf (chapter 3.4)

1 Like
7

Hello Debendra
Have you solved your L2 error problem? And if so, how?