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


282
views
-2
11 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 11 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 10 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 11 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 11 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 11 months ago by Hamilton  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »