Hi, I have been working on Stokes equation for a university project lately and was willing to add a small part moving from a condition where Re<<<1 to a turbulent flow, hence considering the transport term for a scaled problem. I changed the bilinear form taking into account the transport term but I get something weird as warning message when I run the code:
Traceback (most recent call last):
File "/media/sf_Numerics_Scripts/Project/Reynolds/reynolds_100.py", line 122, in <module>
a = form(inner(grad(u), grad(v))*dx + Re*dot(dot(grad(u),u), v)* dx - p*div(v)*dx - q*div(u)*dx)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py", line 166, in form
return _create_form(form)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py", line 161, in _create_form
return _form(form)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/forms.py", line 133, in _form
raise RuntimeError("Expecting to find a Mesh in the form.")
RuntimeError: Expecting to find a Mesh in the form.
Find the complete code hereafter:
# import the modules and functions that the program uses:
# +
import numpy as np
import dolfinx
print(dolfinx.__version__, dolfinx.git_commit_hash)
import petsc4py
print(petsc4py.__version__)
import ufl
from dolfinx import cpp as _cpp
from dolfinx import fem, plot
from dolfinx.fem import (Constant, Function, FunctionSpace, dirichletbc,
extract_function_spaces, form,
locate_dofs_geometrical, locate_dofs_topological)
from dolfinx.io import XDMFFile
from dolfinx.mesh import (CellType, GhostMode, create_rectangle,
locate_entities_boundary)
from ufl import div, dx, grad, inner, SpatialCoordinate, dot, MixedElement, TrialFunctions, TestFunctions
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.io import gmshio
import gmsh
# -
# MESH DEFINITION - REFINED
gmsh.initialize()
gmsh.model.add("refined")
fine_dim = 0.01
coarse_dim= 0.1
gmsh.model.geo.addPoint(0, 0, 0, coarse_dim, 1)
gmsh.model.geo.addPoint(1, 0, 0, coarse_dim, 2)
gmsh.model.geo.addPoint(0, 1, 0, fine_dim, 3)
gmsh.model.geo.addPoint(1, 1, 0, fine_dim, 4)
gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(2, 4, 2)
gmsh.model.geo.addLine(4, 3, 3)
gmsh.model.geo.addLine(3, 1, 4)
gmsh.model.geo.addCurveLoop([1, 2, 3, 4], 1)
gmsh.model.geo.addPlaneSurface([1], 1)
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(1, [1, 2, 4], 5)
gmsh.model.addPhysicalGroup(2, [1], name = "My surface")
gmsh.model.mesh.generate(2)
msh, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, MPI.COMM_SELF, 0, gdim=2)
x = SpatialCoordinate(msh)
# Function to mark x = 0, x = 1 and y = 0
def noslip_boundary(x):
return np.logical_or(np.logical_or(np.isclose(x[0], 0.0),
np.isclose(x[0], 1.0)),
np.isclose(x[1], 0.0))
# Function to mark the lid (y = 1)
def lid(x):
return np.isclose(x[1], 1.0)
# Lid velocity
def lid_velocity_expression(x):
return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1])))
#Definition of the spaces
P2 = ufl.VectorElement("Lagrange", msh.ufl_cell(), 2)
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
mixed = MixedElement([P2, P1])
V, Q = FunctionSpace(msh, P2), FunctionSpace(msh, P1)
# No-slip boundary condition
# where x = 0, x = 1, and y = 0
noslip = np.zeros(msh.geometry.dim, dtype=PETSc.ScalarType)
facets = locate_entities_boundary(msh, 1, noslip_boundary)
bc0 = dirichletbc(noslip, locate_dofs_topological(V, 1, facets), V)
# Driving velocity condition u = (1, 0) on top boundary (y = 1)
lid_velocity = Function(V)
lid_velocity.interpolate(lid_velocity_expression)
facets = locate_entities_boundary(msh, 1, lid)
bc1 = dirichletbc(lid_velocity, locate_dofs_topological(V, 1, facets))
# Collect Dirichlet boundary conditions
bcs = [bc0, bc1]
#### UNCOMMENT IF YOU WANT TO PLOT THE MESH ####
# #PLOT OF THE MESH
# import pyvista
# cells, types, geometry = plot.create_vtk_mesh(Q)
# grid = pyvista.UnstructuredGrid(cells, types, geometry)
# plotter = pyvista.Plotter()
# plotter.add_mesh(grid, show_edges=True)
# plotter.view_xy()
# plotter.show(screenshot='mesh_medium.png')
# We now define the bilinear and linear forms corresponding to the weak
# mixed formulation of the Stokes equations in a blocked structure:
Re= 100;
# Define variational problem
# (u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
# (v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)
u,p = TrialFunctions(mixed)
v,q = TestFunctions(mixed)
f = Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0)))
a = form(inner(grad(u), grad(v))*dx + Re*dot(dot(grad(u),u), v)* dx - p*div(v)*dx - q*div(u)*dx)
L = form([inner(f, v) * dx, inner(Constant(msh, PETSc.ScalarType(0)), q) * dx])
# We will use a block-diagonal preconditioner to solve this problem:
a_p11 = form(inner(p, q) * dx)
a_p = [[a[0][0], None],
[None, a_p11]]
# ### Nested matrix solver
A = fem.petsc.assemble_matrix_nest(a, bcs=bcs)
A.assemble()
# We create a nested matrix `P` to use as the preconditioner.
P11 = fem.petsc.assemble_matrix(a_p11, [])
P = PETSc.Mat().createNest([[A.getNestSubMatrix(0, 0), None], [None, P11]])
P.assemble()
# Next, the right-hand side vector is assembled and then modified to
# account for non-homogeneous Dirichlet boundary conditions:
# +
b = fem.petsc.assemble_vector_nest(L)
# Modify ('lift') the RHS for Dirichlet boundary conditions
fem.petsc.apply_lifting_nest(b, a, bcs=bcs)
# Sum contributions from ghost entries on the owner
for b_sub in b.getNestSubVecs():
b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
# Set Dirichlet boundary condition values in the RHS
bcs0 = fem.bcs_by_block(extract_function_spaces(L), bcs)
fem.petsc.set_bc_nest(b, bcs0)
# -
# The pressure field for this problem is determined only up to a
# constant. We can supply the vector that spans the nullspace and any
# component of the solution in this direction will be eliminated during
# the iterative linear solution process.
# +
# Create nullspace vector
null_vec = fem.petsc.create_vector_nest(L)
# Set velocity part to zero and the pressure part to a non-zero constant
null_vecs = null_vec.getNestSubVecs()
null_vecs[0].set(0.0), null_vecs[1].set(1.0)
# Normalize the vector, create a nullspace object, and attach it to the
# matrix
null_vec.normalize()
nsp = PETSc.NullSpace().create(vectors=[null_vec])
assert nsp.test(A)
A.setNullSpace(nsp)
# -
# Now we create a Krylov Subspace Solver `ksp`. We configure it to use
# the MINRES method, and a block-diagonal preconditioner using PETSc's
# additive fieldsplit type preconditioner:
# +
ksp = PETSc.KSP().create(msh.comm)
ksp.setOperators(A, P)
ksp.setType("minres")
ksp.setTolerances(rtol=1e-9)
ksp.getPC().setType("fieldsplit")
ksp.getPC().setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE)
# Define the matrix blocks in the preconditioner with the velocity and
# pressure matrix index sets
nested_IS = P.getNestISs()
ksp.getPC().setFieldSplitIS(
("u", nested_IS[0][0]),
("p", nested_IS[0][1]))
# Set the preconditioners for each block
ksp_u, ksp_p = ksp.getPC().getFieldSplitSubKSP()
ksp_u.setType("preonly")
ksp_u.getPC().setType("gamg")
ksp_p.setType("preonly")
ksp_p.getPC().setType("jacobi")
# Monitor the convergence of the KSP
ksp.setFromOptions()
#Solve
u, p = Function(V), Function(Q)
x = PETSc.Vec().createNest([_cpp.la.petsc.create_vector_wrap(u.x), _cpp.la.petsc.create_vector_wrap(p.x)])
ksp.solve(b, x)
# Norms of the solution vectors are computed:
ux=u.sub(0)
uy=u.sub(1)
# print(ux.vector().get_local())
norm_u_0 = u.x.norm()
norm_p_0 = p.x.norm()
if MPI.COMM_WORLD.rank == 0:
print("(A) Norm of velocity coefficient vector (nested, iterative): {}".format(norm_u_0))
print("(A) Norm of pressure coefficient vector (nested, iterative): {}".format(norm_p_0))
#plot of the nested matrix solver
################ PLOTTING THE PRESSURE ###################################################
import pyvista
cells, types, geometry = plot.create_vtk_mesh(Q)
grid = pyvista.UnstructuredGrid(cells, types, geometry)
grid.point_data["p"] = p.x.array.real
grid.set_active_scalars("p")
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=False)
plotter.view_xy()
plotter.show(screenshot='pressure.png')
###### SAVE SOLUTIONS #############
with XDMFFile(MPI.COMM_WORLD, "out_reynolds/velocity_r100.xdmf", "w") as ufile_xdmf:
u.x.scatter_forward()
ufile_xdmf.write_mesh(msh)
ufile_xdmf.write_function(u)
with XDMFFile(MPI.COMM_WORLD, "out_reynolds/pressure_r100.xdmf", "w") as pfile_xdmf:
p.x.scatter_forward()
pfile_xdmf.write_mesh(msh)
pfile_xdmf.write_function(p)
Would you help me out with this? Thank you in advance.