Hello everyone,
I am trying to extend the method written in A Newton Method at the Algebraic Level to mixed element in order to solve Navier-Stokes eqeuation. I first solve for the Stokes equation in order to get the initial value to the custom nonlinear solver. However, the residual keep increasing for each iteration, I think it is because I use the wrong commands to update and assign functions inside the iteration loop.
The following is the MWE of what I am trying to do, if anybody has experience on the mixed element, could you please suggest me how to fix this?
from fenics import *
# set_log_level(LogLevel.ERROR)
from petsc4py import PETSc
import numpy as np
#%% classes for subdomains
class Fluid(SubDomain):
def inside(self, x, on_boundary):
return x[1] >= 0.0
class Solid(SubDomain):
def inside(self, x, on_boundary):
return x[1] <= 0.0
# classes for boundaries
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.1)
class LeftF(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0) and between(x[1], (0.0, 0.1))
class RightF(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0) and between(x[1], (0.0, 0.1))
class Interface(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.0)
class LeftS(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0) and between(x[1], (-0.1, 0.0))
class RightS(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0) and between(x[1], (-0.1, 0.0))
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], -0.1)
#%% Varying cases
Re = Constant(5E-1, name="Re") #500
Pr = Constant(0.01, name="Pr") #0.01, 100
kCHT = Constant(20.0, name="kCHT") #1.0, 5.0, 20.0
#%% Parameters for fluid and solid
# Fluid parameters for the variational problem
muF = Constant(1.002, name="muF")
rhoF = Constant(998, name="rhoF")
cdF = Constant(0.5984, name="cdF") #0.5984
cpF = Constant(4181, name="cpF")
alpF = Constant(cdF/rhoF/cpF, name="alpF")
# Solid parameters for the variational problem
muS = Constant(1, name="muS")
rhoS = Constant(2705, name="rhoS")
cdS = Constant(227, name="cdS") #227
cpS = Constant(900, name="cpS")
alpS = Constant(cdS/rhoS/cpS, name="alpS")
#%% meshing and defining the subdomains
mesh = RectangleMesh(Point(0.0, -0.1), Point(1.0, 0.1), 100, 10)
# Define markers for Subdomians, Dirichlet and Neuman bcs
subdomain_marker = MeshFunction("size_t", mesh, mesh.geometry().dim(), 0)
Fluid().mark(subdomain_marker, 1)
Solid().mark(subdomain_marker, 2)
boundary_marker = MeshFunction("size_t", mesh, mesh.geometry().dim()-1, 0)
Top().mark(boundary_marker, 1)
LeftF().mark(boundary_marker, 2)
RightF().mark(boundary_marker, 3)
Interface().mark(boundary_marker, 4) # interface
LeftS().mark(boundary_marker, 5)
RightS().mark(boundary_marker, 6)
Bottom().mark(boundary_marker, 7)
#%% Function Spaces
V = VectorElement("CG", triangle, 2) # Velocity on the fluid domain
P = FiniteElement("CG", triangle, 1) # Pressure on the fluid domain
element = MixedElement([V, P])
W = FunctionSpace(mesh, element)
#%% Boundary conditions
U_noslip = Constant((0.0,0.0), name="U_noslip")
pressure_left = Constant(15.0, name="pressure_left")
pressure_right = Constant(0.0, name="pressure_right")
bc1 = DirichletBC(W.sub(0), U_noslip, boundary_marker, 1)
bc2 = DirichletBC(W.sub(0), U_noslip, boundary_marker, 4)
bcs = [bc1,bc2]
#%% define trial and test functions
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
#%% define the domains for integration
dx = Measure("dx", domain=mesh, subdomain_data=subdomain_marker)
ds = Measure("ds", domain=mesh, subdomain_data=boundary_marker)
#%% Define the nonlinear problem
# unit normal vector
n = FacetNormal(mesh)
# weak form for the Stroke-flow
a = muF*inner(grad(u), grad(v))*dx(1) \
- p*div(v)*dx(1) \
- q*div(u)*dx(1)
l = - pressure_left*inner(n, v)*ds(2)\
- pressure_right*inner(n, v)*ds(3)
# Assembly the lhs and rhs matrices
A = assemble(a, keep_diagonal=True)
L = assemble(l)
# Apply boundary conditions
[bc.apply(A) for bc in bcs]
[bc.apply(L) for bc in bcs]
# solve for the solution
A.ident_zeros()
w_k = Function(W, name="fluid variables")
solver_lin = LUSolver("mumps")
solver_lin.solve(A, w_k.vector(), L)
# Split using deep copy
(u_k, p_k) = w_k.split(True)
u_kVec = u_k.vector().get_local()
p_kVec = p_k.vector().get_local()
# Note that all Dirichlet conditions must be zero for
# the correction function in a Newton-type method
#%% Boundary conditions
bc1N = DirichletBC(W.sub(0), Constant((0.0,0.0)), boundary_marker, 1)
bc2N = DirichletBC(W.sub(0), Constant((0.0,0.0)), boundary_marker, 4)
bcsN = [bc1N,bc2N]
# Define variational problem for the matrix and vector
# in a Newton iteration
dw = TrialFunctions(W) # w = w_k + omega*dw
# weak form for the Navier-Stroke equation
F = rhoF*inner(grad(u_k)*u_k, v)*dx(1) \
+ inner(grad(u_k), grad(v))*dx(1) \
- p_k*div(v)*dx(1) \
- q*div(u_k)*dx(1) \
- pressure_left*inner(n, v)*ds(2)\
- pressure_right*inner(n, v)*ds(3)
# derivative of the weak form
a = derivative(F, w_k)
L = -F
# Newton iteration at the algebraic level
dw = Function(W)
(du, dp) = dw.split(True)
w = Function(W) # u = u_k + omega*du
(u, p) = w.split(True)
omega = 1.0 # relaxation parameter
eps = 1.0
tol = 1.0E-10
iter = 0
maxiter = 25
while eps > tol and iter < maxiter:
iter += 1
print('iteration', iter)
# Assembly the lhs and rhs matrices
A = assemble(a, keep_diagonal=True)
b = assemble(L)
# Apply boundary conditions
[bc.apply(A) for bc in bcsN]
[bc.apply(b) for bc in bcsN]
A.ident_zeros()
solve(A, dw.vector(), b)
eps = np.linalg.norm(dw.vector().get_local(), ord=np.Inf)
print('Norm:', eps)
u.vector()[:] = u_k.vector() + omega*du.vector()
p.vector()[:] = p_k.vector() + omega*dp.vector()
assigner = FunctionAssigner(W, W)
assigner.assign(w_k, w)
convergence = 'convergence after %d Newton iterations at the algebraic level' % iter
if iter >= maxiter:
convergence = 'no ' + convergence
Many thanks in advance,
Phuris