Hi,
I am solving time dependent fluid flow poblem in dolfin (2019.01) and I am struggling between interpolation and PETSc error while working with mixed elements. If I manage to get rid of the interpolation error then PETSc is unhappy while fixing the PETSc causes the interpolation become impossible. Additionally, I need to solve it using nonlinear solver instead of the linear algorithms illustrated in tutorials. I will sincerely appreciate help from the experts.
Below is an MWE:
import gmsh, meshio, dolfin
import matplotlib.pyplot as plt
import tqdm.autonotebook
import numpy as np
from scipy.optimize import fsolve
# Print log messages only from the root process in parallel
dolfin.parameters["std_out_all_processes"] = False;
# Enable extrapolation when mapping fields between different meshes
# https://fenicsproject.discourse.group/t/zero-extrapolation-of-a-function-to-a-greater-domain/1826
dolfin.parameters['allow_extrapolation'] = True
def generate_gmsh_mesh():
gmsh.initialize()
gmsh.option.setNumber("General.Terminal",0)
gdim = 2
model_rank = 0
occ = gmsh.model.occ
fov = (1, 1)
rectangle = occ.addRectangle(0,0,0, fov[0], fov[1])
occ.synchronize()
domains = gmsh.model.getEntities(gdim)
for j, dom in enumerate(domains):
gmsh.model.addPhysicalGroup(gdim, [dom[1]], tag=j+1)
all_bnd = gmsh.model.getEntities(gdim - 1)
for bnd in all_bnd:
gmsh.model.addPhysicalGroup(gdim - 1, [bnd[1]], tag=bnd[1])
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.refine()
gmsh.write("dolfin_ss_mesh.msh")
# close gmsh
gmsh.finalize()
# obtain dolfin mesh
msh = meshio.read("dolfin_ss_mesh.msh")
############### new block for saving and loading mesh with markers
# https://fenicsproject.discourse.group/t/accessing-and-marking-imported-boundaries/5753/8
# https://fenicsproject.discourse.group/t/transitioning-from-mesh-xml-to-mesh-xdmf-from-dolfin-convert-to-meshio/412/8
def create_mesh(mesh, cell_type, prune_z=True):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
mesh_points = mesh.points[:, :2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=mesh_points,
cells={cell_type: cells},
cell_data={"name_to_read": [cell_data]})
return out_mesh
triangle_mesh = create_mesh(msh, "triangle",prune_z=True)
line_mesh = create_mesh(msh, "line",prune_z=True)
meshio.write("dolfin_ss_mesh.xdmf", triangle_mesh)
meshio.write("dolfin_ss_mesh_mf.xdmf", line_mesh)
# now read it from xdmf format
mesh = dolfin.Mesh()
with dolfin.XDMFFile("dolfin_ss_mesh.xdmf") as infile:
infile.read(mesh)
mvc = dolfin.MeshValueCollection("size_t", mesh, dim=1)
with dolfin.XDMFFile("dolfin_ss_mesh_mf.xdmf") as infile:
infile.read(mvc, "name_to_read")
mesh_bnd = dolfin.cpp.mesh.MeshFunctionSizet(mesh, mvc)
return mesh, mesh_bnd
dt, T = 0.01, 1
num_steps = dt/T
Um = 1
L = 1
rhof, mu = 1, 1
u_in = dolfin.Expression(("4*{}*x[1]*({}-x[1])/{}".format(Um, L, L**2), "0.0"), degree=2)
u_noslip = dolfin.Constant((0.0, 0.0))
p_out = dolfin.Constant(0.0)
# initialize mesh
mesh, mesh_bnd = generate_gmsh_mesh()
# Define coefficients
k = dolfin.Constant(dt)
f = dolfin.Constant((0, 0))
inner = dolfin.inner
dx = dolfin.dx
grad = dolfin.grad
div = dolfin.div
V_el = dolfin.VectorElement("CG", mesh.ufl_cell(), 2)
Q_el = dolfin.FiniteElement("CG", mesh.ufl_cell(), 1)
W = dolfin.FunctionSpace(mesh, V_el*Q_el)
up0 = dolfin.Function(W)
u0 = up0.split()[0] # get a Function
# Time-stepping
progress = tqdm.autonotebook.tqdm(desc="Solving PDE", total=num_steps)
t = dt # time for the first iteration
t_idx = 1
while t < T + dolfin.DOLFIN_EPS:
progress.update(1)
# generate new mesh given the position of the cylinder
mesh, mesh_bnd = generate_gmsh_mesh() # in MWE, the mesh does not change, but later it will
# define the physics given the new mesh
# Define function spaces (P2-P1)
V_el = dolfin.VectorElement("CG", mesh.ufl_cell(), 2)
Q_el = dolfin.FiniteElement("CG", mesh.ufl_cell(), 1)
W = dolfin.FunctionSpace(mesh, V_el*Q_el)
ut, pt = dolfin.TrialFunctions(W)
up = dolfin.Function(W)
upn = dolfin.Function(W)
u, p = up.split(deepcopy=True)
un, pn = upn.split(deepcopy=True)
u.interpolate(u0) # set the starting point of the nonlinear solver to the old solution
un.interpolate(u0) # interpolate old solution to the new mesh
F = rhof/k*inner(u - un, ut)*dx + rhof*inner(grad(u)*u, ut)*dx + mu*inner(grad(u), grad(ut))*dx - p*div(ut)*dx
+ rhof*pt*div(u)*dx
x = dolfin.SpatialCoordinate(mesh)
# boundary conditions
bc_bot = dolfin.DirichletBC(W.sub(0), u_noslip, mesh_bnd, 1)
bc_outlet = dolfin.DirichletBC(W.sub(1), p_out, mesh_bnd, 2)
bc_top = dolfin.DirichletBC(W.sub(0), u_noslip, mesh_bnd, 3)
bc_inlet = dolfin.DirichletBC(W.sub(0), u_in, mesh_bnd, 4)
bcs = [bc_top, bc_bot, bc_inlet, bc_outlet]
# https://fenicsproject.discourse.group/t/attributeerror-listtensor-has-no-attribute-cpp-object/4307
dolfin.solve(F==0, up, bcs=bcs,
solver_parameters={"newton_solver":{"relative_tolerance": 1e-6,
# "linear_solver": "mumps",
"maximum_iterations" : 100,
}},
form_compiler_parameters={"cpp_optimize": True, "representation": "uflacs"} )
# save the current distribution of velocity and pressure
ufile << u
pfile << p
# Prepare for the next time step
up0 = dolfin.Function(W) # Function for the current mesh
u0 = up0.split(deepcopy=True)[0]
u0.assign(u) # store the current solution for the next iteration
t += dt
t_idx += 1
