Hey @dokken,
The solution for legacy FEniCS after @bleyerj recommendations worked well and the updated code looks like below:
from fenics import *
from fenics_adjoint import *
import numpy as np
import ufl as ufl
from mpi4py import MPI
mu_fluid = 0.0089
rho_fluid = 1.0
L = 1.0
u_inlet = 0.1
p_outlet = 0.0
Re = rho_fluid*u_inlet*L/mu_fluid
lx, ly, lz = 2.0, 1.0, 1.0
nelx, nely, nelz = 20, 10, 10
mesh = BoxMesh.create(MPI.COMM_WORLD, [Point(0.0, 0.0, 0.0), Point(lx, ly, lz)],
[nelx, nely, nelz],
CellType.Type.hexahedron)
V = VectorElement("CG", mesh.ufl_cell(), 2)
Q = FiniteElement("CG", mesh.ufl_cell(), 1)
solution_space = MixedElement([V,Q])
W = FunctionSpace(mesh, solution_space)
parameters["form_compiler"]["quadrature_degree"] = 2
bcs = []
fluid_boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
fluid_boundary_parts.set_all(3)
no_slip_boundary_right = CompiledSubDomain("near(x[1],0.0)")
no_slip_boundary_left = CompiledSubDomain("near(x[1],ly)", ly = ly)
no_slip_boundary_bottom = CompiledSubDomain("near(x[2],0.0)")
no_slip_boundary_top = CompiledSubDomain("near(x[2],lz)", lz=lz)
inlet_front = CompiledSubDomain("on_boundary && near(x[0],0.0)")
outlet_back = CompiledSubDomain("on_boundary && near(x[0],lx)", lx=lx)
no_slip_boundary_right.mark(fluid_boundary_parts, 2)
no_slip_boundary_left.mark(fluid_boundary_parts, 2)
no_slip_boundary_bottom.mark(fluid_boundary_parts, 2)
no_slip_boundary_top.mark(fluid_boundary_parts, 2)
inlet_front.mark(fluid_boundary_parts, 0)
outlet_back.mark(fluid_boundary_parts, 1)
bcu_no_slip_right = DirichletBC(W.sub(0), Constant((0.0, 0.0, 0.0)), no_slip_boundary_right)
bcu_no_slip_left = DirichletBC(W.sub(0), Constant((0.0, 0.0, 0.0)), no_slip_boundary_left)
bcu_no_slip_bottom = DirichletBC(W.sub(0), Constant((0.0, 0.0, 0.0)), no_slip_boundary_bottom)
bcu_no_slip_top = DirichletBC(W.sub(0), Constant((0.0, 0.0, 0.0)), no_slip_boundary_top)
inflow_profile = Expression(('u_inlet','0','0'), u_inlet=u_inlet, degree=2)
bcu_inlet = DirichletBC(W.sub(0), inflow_profile, inlet_front)
bcu_back_pressure = DirichletBC(W.sub(1), Constant(0), outlet_back)
bcs = [bcu_no_slip_right, bcu_no_slip_left,\
bcu_no_slip_bottom, bcu_no_slip_top,\
bcu_inlet, bcu_back_pressure]
w = Function(W)
dw = TrialFunction(W)
(u, p) = split(w)
(du, dp) = TrialFunctions(W)
(v, q) = TestFunctions(W)
F_navier_stokes = (
rho_fluid*dot(dot(u,nabla_grad(u)),v)*dx
+ mu_fluid*inner(grad(u),grad(v))*dx
- inner(p,div(v))*dx
+ q*div(u)*dx
)
J_navier_stokes = derivative(F_navier_stokes, w, dw)
problem_navier_stokes = NonlinearVariationalProblem(F_navier_stokes, w, bcs, J_navier_stokes)
solver_navier_stokes = NonlinearVariationalSolver(problem_navier_stokes)
solver_navier_stokes.parameters["newton_solver"]["absolute_tolerance"] = 1E-4
solver_navier_stokes.parameters["newton_solver"]["relative_tolerance"] = 1E-3
solver_navier_stokes.parameters["newton_solver"]["maximum_iterations"] = 1000
solver_navier_stokes.parameters["newton_solver"]["linear_solver"] = "mumps"
solver_navier_stokes.solve()
dGamma = Measure("ds")(subdomain_data=fluid_boundary_parts)
inlet_area_local = L**2
pressure_drop = (1/inlet_area_local)*p*dGamma(0)
dpressure_drop_dU = derivative(pressure_drop, p, dp)
dpressure_drop_dU_vector = assemble(dpressure_drop_dU)
bcs_new = []
inflow_profile = Expression(('0','0','0'), u_inlet=u_inlet, degree=2)
bcu_inlet_new = DirichletBC(W.sub(0), inflow_profile, inlet_front)
bcs_new = [bcu_no_slip_right, bcu_no_slip_left,\
bcu_no_slip_bottom, bcu_no_slip_top,\
bcu_inlet_new, bcu_back_pressure]
dResidual_dU = J_navier_stokes
A, b = assemble_system(dResidual_dU, -dpressure_drop_dU, bcs_new)
solution = Function(W)
linear_solver = KrylovSolver('default', 'ilu')
linear_solver.parameters["relative_tolerance"] = 1E-5
linear_solver.solve(A, solution.vector(), b)
print("Norm of solution: ", norm(solution.vector(), 'l2'), flush=True)
print("Solution[23]: ", solution.vector()[23], flush=True)
where the norm of solution and value at element 23 of dolfin mesh are:
Norm of solution: 2.55711615777551e-35
Solution[23]: -9.427552233559459e-38
In the translated code to dolfinx syntax, the solution of the element of dolfinx mesh located at the same position of element 23 (numbering is different between dolfin and dolfinx) is completely different. The code in dolfinx looks like below. I’m using dolfinx 0.6.0:
from mpi4py import MPI
from dolfinx.fem import (
Function, FunctionSpace, form, locate_dofs_topological, dirichletbc,
Constant, assemble_scalar
)
from dolfinx.fem.petsc import (
apply_lifting, assemble_matrix, assemble_vector, create_vector,
set_bc, NonlinearProblem
)
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.mesh import (
CellType, create_box, locate_entities_boundary, meshtags
)
from ufl import (
dx, grad, inner, derivative, TrialFunction, TrialFunctions, Measure, dot,
adjoint, VectorElement, FiniteElement, MixedElement, TestFunctions,
CellDiameter, nabla_grad, div, split
)
from petsc4py.PETSc import ScalarType
from dolfinx.io import XDMFFile
import numpy as np
from dolfinx import log
from petsc4py import PETSc
mu_fluid = 0.0089
rho_fluid = 1.0
L = 1.0
radius = L
u_inlet = 0.1
p_outlet = 0.0
alpha_solid = 1.0e6
alpha_fluid = 0.0
lx, ly, lz = 2.0, 1.0, 1.0
nelx, nely, nelz = 20, 10, 10
domain = create_box(MPI.COMM_WORLD, [(0., 0., 0.), (lx, ly, lz)], [nelx, nely, nelz], CellType.hexahedron)
V = VectorElement("CG", domain.ufl_cell(), 2)
Q = FiniteElement("CG", domain.ufl_cell(), 1)
solution_space = MixedElement([V,Q])
W = FunctionSpace(domain, solution_space)
D = FunctionSpace(domain, ("DG", 0))
volume_fractions = Function(D)
tolerance = 1e-6
boundary_conditions_navier_stokes = []
def no_slip_boundary_right(x):
return (x[1] < tolerance)
def no_slip_boundary_left(x):
return (x[1] > (ly - tolerance))
def no_slip_boundary_bottom(x):
return (x[2] < tolerance)
def no_slip_boundary_top(x):
return (x[2] > (lz - tolerance))
def inlet_front(x):
return (x[0] < tolerance)
def outlet_back(x):
return (x[0] > (lx - tolerance))
def inlet_velocity(x):
values = np.zeros((3, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = u_inlet
return values
def inlet_velocity_homogeneous(x):
values = np.zeros((3, x.shape[1]), dtype=PETSc.ScalarType)
return values
facets_no_slip_boundary_right = locate_entities_boundary(domain, domain.topology.dim - 1, no_slip_boundary_right)
facets_no_slip_boundary_left = locate_entities_boundary(domain, domain.topology.dim - 1, no_slip_boundary_left)
facets_no_slip_boundary_bottom = locate_entities_boundary(domain, domain.topology.dim - 1, no_slip_boundary_bottom)
facets_no_slip_boundary_top = locate_entities_boundary(domain, domain.topology.dim - 1, no_slip_boundary_top)
facets_inlet_front = locate_entities_boundary(domain, domain.topology.dim - 1, inlet_front)
facets_outlet_back = locate_entities_boundary(domain, domain.topology.dim - 1, outlet_back)
facet_indices = [facets_no_slip_boundary_right, facets_no_slip_boundary_left,\
facets_no_slip_boundary_bottom, facets_no_slip_boundary_top,\
facets_inlet_front, facets_outlet_back]
facet_markers = []
facet_markers.append(np.full_like(facet_indices[0], 2))
facet_markers.append(np.full_like(facet_indices[1], 2))
facet_markers.append(np.full_like(facet_indices[2], 2))
facet_markers.append(np.full_like(facet_indices[3], 2))
facet_markers.append(np.full_like(facet_indices[4], 0))
facet_markers.append(np.full_like(facet_indices[5], 1))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tags = meshtags(domain, domain.topology.dim - 1, facet_indices[sorted_facets], facet_markers[sorted_facets])
B, B_to_W = W.sub(0).collapse()
u_no_slip = Function(B)
u_no_slip.x.set(0)
bcu_no_slip_right = dirichletbc(u_no_slip, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_no_slip_boundary_right), W.sub(0))
bcu_no_slip_left = dirichletbc(u_no_slip, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_no_slip_boundary_left), W.sub(0))
bcu_no_slip_bottom = dirichletbc(u_no_slip, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_no_slip_boundary_bottom), W.sub(0))
bcu_no_slip_top = dirichletbc(u_no_slip, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_no_slip_boundary_top), W.sub(0))
u_inlet_profile = Function(B)
u_inlet_profile.interpolate(inlet_velocity)
bcu_inlet = dirichletbc(u_inlet_profile, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_inlet_front), W.sub(0))
G, G_to_W = W.sub(1).collapse()
zero_pressure = Function(G)
zero_pressure.x.set(0)
bcu_back_pressure = dirichletbc(zero_pressure, locate_dofs_topological((W.sub(1), G), domain.topology.dim - 1, facets_outlet_back), W.sub(1))
boundary_conditions_navier_stokes = [bcu_no_slip_right, bcu_no_slip_left,\
bcu_no_slip_bottom, bcu_no_slip_top,\
bcu_inlet, bcu_back_pressure]
def alpha(volume_fractions):
return alpha_solid + (alpha_fluid - alpha_solid)*volume_fractions
volume_fractions.x.set(1.0)
w = Function(W)
dw = TrialFunction(W)
(u, p) = split(w)
(du, dp) = TrialFunctions(W)
(v, q) = TestFunctions(W)
dx = Measure("dx", domain=domain, metadata={"quadrature_degree": 2})
F_navier_stokes = (
rho_fluid*dot(dot(u,nabla_grad(u)),v)*dx
+ mu_fluid*inner(grad(u),grad(v))*dx
- inner(p,div(v))*dx
+ alpha(volume_fractions)*dot(v,u)*dx
+ q*div(u)*dx
)
J_navier_stokes = derivative(F_navier_stokes, w, dw)
problem_navier_stokes = NonlinearProblem(F_navier_stokes, w,
boundary_conditions_navier_stokes,
J_navier_stokes)
solver_navier_stokes = NewtonSolver(MPI.COMM_WORLD, problem_navier_stokes)
solver_navier_stokes.convergence_criteria = "incremental"
solver_navier_stokes.report = True
solver_navier_stokes.atol = 1E-4
solver_navier_stokes.rtol = 1E-3
solver_navier_stokes.max_it = 1000
krylov_subspace = solver_navier_stokes.krylov_solver
options = PETSc.Options()
option_prefix = krylov_subspace.getOptionsPrefix()
options[f"{option_prefix}kps_type"] = "preonly"
krylov_subspace.setFromOptions()
newton_solver_iterations, converged = solver_navier_stokes.solve(w)
assert (converged)
w.x.scatter_forward()
dGamma = Measure("ds", domain=domain, subdomain_data=facet_tags, metadata=None)
inlet_area = L**2
pressure_drop = (1/inlet_area)*p*dGamma(0)
dpressure_drop_dp = derivative(pressure_drop, p, dp)
dpressure_drop_dp_form = form(dpressure_drop_dp)
bcs_adjoint = []
u_inlet_profile = Function(B)
u_inlet_profile.interpolate(inlet_velocity_homogeneous)
bcu_inlet_new = dirichletbc(u_inlet_profile, locate_dofs_topological((W.sub(0), B), domain.topology.dim - 1, facets_inlet_front), W.sub(0))
bcs_adjoint = [bcu_no_slip_right, bcu_no_slip_left,\
bcu_no_slip_bottom, bcu_no_slip_top,\
bcu_inlet_new, bcu_back_pressure]
solution = Function(W)
dresidual_dU = J_navier_stokes
dresidual_dU_form = form(dresidual_dU)
dresidual_dU_matrix = assemble_matrix(dresidual_dU_form, bcs_adjoint)
dresidual_dU_matrix.assemble()
dpressure_drop_dp_vector = create_vector(dpressure_drop_dp_form)
assemble_vector(dpressure_drop_dp_vector, dpressure_drop_dp_form)
apply_lifting(dpressure_drop_dp_vector, [dresidual_dU_form], [bcs_adjoint])
dpressure_drop_dp_vector.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(dpressure_drop_dp_vector, bcs_adjoint)
options = PETSc.Options()
linear_solver = PETSc.KSP().create(MPI.COMM_WORLD)
linear_solver.setType("preonly")
linear_solver.setTolerances(rtol=1e-8)
linear_solver.getPC().setType("ilu")
options = PETSc.Options()
linear_solver.setFromOptions()
linear_solver.setInitialGuessNonzero(False)
linear_solver.setOperators(dresidual_dU_matrix)
linear_solver.solve(dpressure_drop_dp_vector, solution.vector)
solution.x.scatter_forward()
print("Norm of solution: ", np.linalg.norm(solution.x.array), flush=True)
print("solution[21]: ", solution.x.array[21])
where the norm of solution and value obtained at element 21 of dolfinx mesh are:
Norm of solution: 237.3371158413562
solution[21]: -9.345168135762837
I’ve already checked the solution in paraview, they are the same. I also checked the boundary conditions application and couldn’t find any difference. Actually, I printed the locations where I apply the boundary conditions to confirm they were being properly selected. For me it seems a matter of syntax, i.e., on how to work with mixed elements in dolfinx. What do you think? Thank you.