Hi, all. Could you help me to check my code?


242
views
-2
9 months ago by
Hi, all,
I am implementing Gauge method for Naver-Stokes equations, see
https://web.math.princeton.edu/~weinan/pdf%20files/gauge%20method.pdf in details.
The Gauge method contains three steps, as follows.
  $m^i,\quad u^i,\quad\phi^i$mi, ui, ϕi  are known as the intial values
  • find  $m^{i+1}$mi+1 
            $\frac{m^{i+1}-m^i}{dt}-\mu\Delta m^{i+1}+\left(u^i\cdot\nabla\right)u^i=f,\qquad in\quad\Omega$mi+1midt μΔmi+1+(ui·)ui=ƒ , in Ω  
           with the boundary conditions as follows
              $m^{i+1}=\left(\nabla\phi^i\cdot\tau\right)\tau,\quad\quad on\quad\partial\Omega$mi+1=(ϕi·τ)τ, on ∂Ω 
  • find  $\phi^{i+1}$ϕi+1 
            $\Delta\phi^{i+1}=\nabla\cdot m^{i+1}$Δϕi+1=·mi+1 ,     in $\Omega$Ω 
            with pure Neumann boundary condition
             $\frac{\partial\phi^{i+1}}{\partial n}=0$ϕi+1n =0

  • find  $u^{i+1}$ui+1 
           $u^{i+1}=m^{i+1}-\nabla\phi^{i+1}$ui+1=mi+1ϕi+1 

Could you tell me how to implement the procedure?
Community: FEniCS Project
3
With patience and diligent research.
written 9 months ago by pf4d  
Hi, Evan Cummings
Could you tell me how to implement the boundary condition of  $m^{i+1}$mi+1 for the first step?
I did my best, however, it still does not work.
Best,
Hamilton
written 9 months ago by Hamilton  
######   Gauge method with Dirichlet Boundary conditon   #######
######   ( u \dot n )\dot n +  ( u \dot t )\dot t = u    #######
######         u = (\grad(\phi) \cdot \tau)\tau          #######

from dolfin import *
import time
import numpy as np
parameters['reorder_dofs_serial'] = False
parameters["std_out_all_processes"] = False;

#######==================================#########
THETA    =  0.05
NU       =  0.001
Pt       =  30
t        =  0.0
dt       =  1.0/(Pt*Pt)
no_slip  =  Expression(["0.0","0.0"])
MAXTIME  =  3 

#######==================================#########
Th = UnitSquareMesh(Pt,Pt)
def NoSlip(x, on_boundary):
   tol=1E-14
   return on_boundary and (abs(x[1]) < tol or abs(x[0]) < tol
                     or abs(x[1]-1) < tol or abs(x[0]-1) < tol)
n = FacetNormal(Th)

def right_term(phi):
	n = FacetNormal(Th)
	grad_phi = grad(phi)
	t = as_vector([-n[1],n[0]])
	grad_phi_t = dot(grad_phi, t)
	grad_phi_t_t = as_vector([grad_phi_t*t])
	return grad_phi_t_t
   
Vh_cg = VectorFunctionSpace(Th, "Lagrange", 2)
Vh_dg = VectorFunctionSpace(Th, "DG", 2)
Mh = FunctionSpace(Th, "Lagrange", 1)
grad_phi = Function(Vh_cg)

######################=============================##########################
U_exat = Expression(["pi*pow(sin(2*pi*x[0]),2)*sin(4*pi*x[1])*cos(t)","-pi*sin(4*pi*x[0])*pow(sin(2*pi*x[1]),2)*cos(t)"],THETA=THETA,t=0.0)
p_exat = Expression("-cos(2*pi*x[0])*sin(2*pi*x[1])*cos(t)", t=0.0)
phi_i = project(Expression("0.0"), Mh)

bcs  = DirichletBC(Vh_cg,right_term(phi_i),NoSlip)
# Get U_0
U_0 = project(U_exat, Vh_cg)

f = Expression([" 2*pi*sin(2*pi*x[0])*sin(2*pi*x[1])*cos(t) - pi*pow(sin(2*pi*x[0]),2)*sin(4*pi*x[1])*sin(t) - NU*(8*pow(pi,3)*pow(cos(2*pi*x[0]),2)*sin(4*pi*x[1])*cos(t) - 24*pow(pi,3)*pow(sin(2*pi*x[0]),2)*sin(4*pi*x[1])*cos(t)) + 4*pow(pi,3)*cos(2*pi*x[0])*pow(sin(2*pi*x[0]),3)*pow(sin(4*pi*x[1]),2)*pow(cos(t),2) - 4*pow(pi,3)*cos(4*pi*x[1])*pow(sin(2*pi*x[0]),2)*sin(4*pi*x[0])*pow(sin(2*pi*x[1]),2)*pow(cos(t),2)","NU*(8*pow(pi,3)*pow(cos(2*pi*x[1]),2)*sin(4*pi*x[0])*cos(t) - 24*pow(pi,3)*sin(4*pi*x[0])*pow(sin(2*pi*x[1]),2)*cos(t)) + pi*sin(4*pi*x[0])*pow(sin(2*pi*x[1]),2)*sin(t) - 2*pi*cos(2*pi*x[0])*cos(2*pi*x[1])*cos(t) + 4*pow(pi,3)*cos(2*pi*x[1])*pow(sin(4*pi*x[0]),2)*pow(sin(2*pi*x[1]),3)*pow(cos(t),2) - 4*pow(pi,3)*cos(4*pi*x[0])*pow(sin(2*pi*x[0]),2)*pow(sin(2*pi*x[1]),2)*sin(4*pi*x[1])*pow(cos(t),2)"],THETA=THETA,t=0.0,NU=NU) 

M_i   = Function(Vh_cg)
M_j   = Function(Vh_cg)
phi_i = Function(Mh)
phi_j = Function(Mh)
U_j   = Function(Vh_cg)
U_i   = Function(Vh_cg)

M = TrialFunction(Vh_cg)
V = TestFunction(Vh_cg)
p = TrialFunction(Mh)
q = TestFunction(Mh)

M_i.assign(U_0)
U_i.assign(U_0)
U_j_dg  = Function(Vh_dg)
######==============  steps ===============######
F1 = (1/dt)*inner(M-M_i, V)*dx \
    + (NU)*inner(grad(M), grad(V))*dx \
    + inner(grad(U_i)*U_i, V)*dx \
    - inner(f, V)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Step 2
F2 = inner(grad(p), grad(q))*dx - div(M_j)*q*dx + 1e-12*p*q*dx
a2 = lhs(F2)
L2 = rhs(F2)

#######=============  Solve Process ===========#####
while t < MAXTIME :
    t += dt
    print "t =", t    
    # Step 1
    solve(a1 == L1, M_j, bcs)
    # Step 2
    solve(a2 == L2, phi_j)
    '''
    # Assemble system
    A2 = assemble(a2)
    b2 = assemble(L2)
    solver = KrylovSolver(A2, "gmres")
    null_vec = Vector(phi_j.vector())
    Mh.dofmap().set(null_vec, 1.0)
    null_vec *= 1.0/null_vec.norm("l2")
    null_space = VectorSpaceBasis([null_vec])
    solver.set_nullspace(null_space)
    null_space.orthogonalize(b2)
    solver.solve(phi_j.vector(), b2)
    '''
    # Step 3
    Temp_1 = interpolate(M_j, Vh_dg)
    Temp_2 = project(grad(phi_j), Vh_dg)
    U_j_dg.vector()[:] = Temp_1.vector().array() + Temp_2.vector().array() 
    #solve(a3 == L3, U_j_dg, bcs)
    U_j = project(U_j_dg, Vh_cg)
    #  Update values 
    U_i.assign(U_j)
    M_i.assign(M_j)  
    phi_i.assign(phi_j)
    grad_phi = project(grad(phi_i), Vh_cg)
    plot(U_i, title="U_i")#,interactive =True)
written 9 months ago by Hamilton  
The code above does not work, mainly because the boundary condition is not right. Namely, I do not know how to implement the follow boundary condition
 $m^{i+1}=\left(\nabla\phi^{i+1}\cdot\tau\right)\tau$mi+1=(ϕi+1·τ)τ 
The related topic are can be seen,
https://www.allanswered.com/post/mgml/how-to-solve-poisson-equation-with-tangential-boundary-condition/
written 9 months ago by Hamilton  
Thanks for the comments of Nicolás Barnafi and Adam Janecka. I have try my best to implement it, however, it does not work. Could you help me to check my code, mainly the boundary conditions?
written 9 months ago by Hamilton  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »