Hello,
I am trying to modify this Fenics example for the Navier-Stokes equations. For reasons related to my project, I would like to solve the three steps of the splitting scheme at once, with mixed finite-elements. I know that mixed FE is not needed in this case, but I will need it in a more complex code that this minimal working example that I prepared will be merged in. And this minimal working example reproduces the issue I have.
The minimal working example is example.py, where I merge steps 1 and 2 in a step ‘12’, in which I solve simultaneously for v^\ast and \phi:
from __future__ import print_function
from fenics import *
from mshr import *
import argparse
from mshr import *
import ufl as ufl
parser = argparse.ArgumentParser()
parser.add_argument("input_directory")
parser.add_argument("output_directory")
parser.add_argument("T")
parser.add_argument("N")
args = parser.parse_args()
L = 2.2
h = 0.41
r = 0.05
c_r = [0.2, 0.2]
T = (float)(args.T)
N = (int)(args.N)
dt = T / N
rho = 1.0
eta = 0.001
l_profile_v = Expression(('4.0*1.5*x[1]*(h - x[1]) / pow(h, 2)', '0'), degree=2, h=h)
i, j, k, l = ufl.indices(4)
#read mesh
mesh=Mesh()
with XDMFFile((args.input_directory) + "/triangle_mesh.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 2)
with XDMFFile((args.input_directory) + "/line_mesh.xdmf") as infile:
infile.read(mvc, "name_to_read")
# Define boundaries
boundary = 'on_boundary'
boundary_l = 'near(x[0], 0.0)'
boundary_r = 'near(x[0], 2.2)'
boundary_tb = 'near(x[1], 0) || near(x[1], 0.41)'
boundary_square = 'on_boundary && sqrt(pow(x[0] - 0.2, 2) + pow(x[1] - 0.2, 2)) > 0.1'
boundary_circle = 'on_boundary && sqrt(pow(x[0] - 0.2, 2) + pow(x[1] - 0.2, 2)) < 0.1'
# Define function spaces
P_v_bar = VectorElement( 'P', triangle, 2 )
P_phi = FiniteElement('P', triangle, 1)
element = MixedElement( [P_v_bar, P_phi ] )
Q = FunctionSpace(mesh, element)
Q_v_bar = Q.sub(0).collapse()
Q_phi = Q.sub(1).collapse()
Q_v_n = VectorFunctionSpace(mesh, 'P', 2)
# Create XDMF files for visualization output
xdmffile_v_n = XDMFFile( (args.output_directory) + '/v_n.xdmf' )
xdmffile_sigma_n = XDMFFile( (args.output_directory) + '/sigma_n.xdmf' )
# Define functions
nu_v_bar, nu_phi = TestFunctions(Q)
nu_v_n= TestFunction(Q_v_n)
v_n = Function(Q_v_n)
v_n_1 = Function(Q_v_n)
v_n_2 = Function(Q_v_n)
v_bar_ = Function(Q_v_bar)
sigma_n = Function(Q_phi)
sigma_n_1 = Function(Q_phi)
sigma_n_2 = Function(Q_phi)
phi_ = Function(Q_phi)
psi = Function(Q)
J_psi = TrialFunction(Q)
J_v_n = TrialFunction(Q_v_n)
v_bar, phi = split( psi )
v_bar_0 = Function( Q_v_bar )
sigma_0 = Function(Q_phi)
v_n_0 = Function( Q_v_n )
V = (v_bar + v_n_1) / 2.0
sigma_ast = (sigma_n_1 + sigma_n_2)/2.0
# Define expressions used in variational forms
Deltat = Constant(dt)
rho = Constant(rho)
F1 = ( \
rho * (((v_bar[i] - v_n_1[i]) / Deltat \
+ ( 3.0/2.0*v_n_1[j] - 1.0/2.0*v_n_2[j] ) * ((V[i]).dx(j)) ) * nu_v_bar[i] ) \
+ (sigma_n_1 + sigma_n_2)/2.0 * ((nu_v_bar[i]).dx(i)) \
+ eta * ((V[i]).dx(j)) * ((nu_v_bar[i]).dx(j)) \
) * dx
F2 = ((phi.dx( i )) * (nu_phi.dx( i )) + rho / Deltat * ((v_bar[i]).dx( i )) * nu_phi) * dx
F3 = ( ( ( v_n[i] - v_bar[i] ) + Deltat/rho * (phi.dx(i)) ) * nu_v_n[i] ) * dx
#total functional for the mixed problem: (1,2) + 3
F12 = F1 + F2
# boundary conditions
bc_v_bar_l = DirichletBC(Q.sub(0), l_profile_v, boundary_l)
bc_v_bar_tb = DirichletBC(Q.sub(0), Constant((0, 0)), boundary_tb)
bc_v_bar_circle = DirichletBC(Q.sub(0), Constant((0, 0)), boundary_circle)
bc_phi = DirichletBC(Q.sub(1), Constant(0), boundary_r)
bcs12 = [bc_v_bar_l, bc_v_bar_tb, bc_v_bar_circle, bc_phi]
bcs3 = []
#variational problem
J12 = derivative( F12, psi, J_psi )
problem12 = NonlinearVariationalProblem( F12, psi, bcs12, J12 )
solver12 = NonlinearVariationalSolver( problem12 )
J3 = derivative( F3, v_n, J_v_n )
problem3 = NonlinearVariationalProblem( F3, v_n, bcs3, J3 )
solver3 = NonlinearVariationalSolver( problem3 )
# Time-stepping
t = 0
for n in range(N):
print("\n* n = ", n, "\n")
xdmffile_v_n.write(v_n_1, t)
xdmffile_sigma_n.write(sigma_n_1, t)
t+=dt
solver12.solve()
solver3.solve()
#update previous solution: update v_n_2 and w_n_2
v_n_2.assign(v_n_1)
v_n_1.assign(v_n)
#get the solution and write it to file
v_bar_, phi_ = psi.split(deepcopy=True)
#update previous solution: update sigma
sigma_n.assign(sigma_n_2-2.0*project(phi_, Q_phi))
sigma_n_2.assign(sigma_n_1)
sigma_n_1.assign(sigma_n)
To run it, do
$python3 example.py [path to the mesh files] [path where to store the solution] 1 128, where the mesh files are here.
The solution for v agrees with the code without mixed FE in the Fenics example, but the tension \sigma does not, and it blows up at the left edge, see screenshot.
Why???
Thanks
