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

330
views
-2
13 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?
3
With patience and diligent research.
written 13 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 13 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)
t = as_vector([-n[1],n[0]])

Vh_cg = VectorFunctionSpace(Th, "Lagrange", 2)
Vh_dg = VectorFunctionSpace(Th, "DG", 2)
Mh = FunctionSpace(Th, "Lagrange", 1)

######################=============================##########################
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 \
- inner(f, V)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Step 2
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)
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)
plot(U_i, title="U_i")#,interactive =True)

written 13 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,