Im trying to reproduce an algorithm for shape optimization in stokes flow by phase-field approaches，and the details of this algorithm ：
Where α(φ) = α_0(1 − φ)2 and f (φ) = φ(1 − φ)(1/2 − φ).
Some definitions on the finite-dimensional spaces:
To reproduce the algorithm,I wrote the codes below:
code 1：update the velociaty field
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import adios4dolfinx as adx
import pathlib
import ufl
from basix.ufl import element, mixed_element
from dolfinx import fem, la, default_scalar_type
from dolfinx.fem import (Constant, Function, dirichletbc,
extract_function_spaces, form, functionspace,
locate_dofs_topological)
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block
from dolfinx.io import XDMFFile,VTXWriter
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary, locate_entities
from ufl import div, dx, grad, inner, ds
# We create a {py:class}`Mesh <dolfinx.mesh.Mesh>`, define functions for
# locating geometrically subsets of the boundary, and define a function
# for the velocity on the lid:
# +
# Create mesh
msh = create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])],
[32, 32], CellType.triangle)
# Function to mark y = 0, y = 1
def noslip_boundary(x):
return np.logical_or(np.isclose(x[1], 0.0), np.isclose(x[1], 1.0))
# Function to mark the lid (x = 0)
def lid(x):
return np.isclose(x[0], 0.0)
# Lid velocity (y(1 - y)/0.25, 0)
def lid_velocity_expression(x):
y = x[1]
velocity_x = y * (1 - y) / 0.25
velocity_y = np.zeros_like(y)
return np.stack((velocity_x, velocity_y))
P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
P1 = element("Lagrange", msh.basix_cell(), 1)
V, Q = functionspace(msh, P2), functionspace(msh, P1)
# No-slip condition on boundaries where y = 0, and y = 1
noslip = np.zeros(msh.geometry.dim, dtype=PETSc.ScalarType) # type: ignore
facets = locate_entities_boundary(msh, 1, noslip_boundary)
bc0 = dirichletbc(noslip, locate_dofs_topological(V, 1, facets), V)
# Driving (lid) velocity condition on left boundary (x = 0)
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]
# -
#alpha_0
alpha_0=Constant(msh, PETSc.ScalarType(0.2))
#define phi_0
W = functionspace(msh, ("DG", 0))
def condition_0(x):
return np.logical_or(x[0] <= 0.25, np.logical_and(x[1] >= 1/3, x[1] <= 2/3))
def condition_1(x):
return np.logical_or(np.logical_and(x[0] > 0.25,x[1] < 1/3),np.logical_and(x[0] > 0.25,x[1] > 2/3))
phi_0 = Function(W)
cells_0 = locate_entities(msh, msh.topology.dim, condition_0)
cells_1 = locate_entities(msh, msh.topology.dim, condition_1)
phi_0.x.array[cells_0] = np.full_like(cells_0, 1, dtype=default_scalar_type)
phi_0.x.array[cells_1] = np.full_like(cells_1, 0, dtype=default_scalar_type)
# The bilinear and linear forms for the Stokes equations are defined
# using a a blocked structure:
# +
# Define variational problem
(u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
(v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)
f = Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0)))
a = form([[inner(grad(u), grad(v)) * dx + inner(alpha_0*(Constant(msh, PETSc.ScalarType(1.0))-phi_0)*(Constant(msh, PETSc.ScalarType(1.0))-phi_0)*u, v)*dx, inner(p, div(v)) * dx],
[inner(div(u), q) * dx, None]])
L = form([inner(f, v) * dx, inner(Constant(msh, PETSc.ScalarType(0)), q) * dx])
# -
# A block-diagonal preconditioner will be used with the iterative
# solvers for this problem:
a_p11 = form(inner(p, q) * dx)
a_p = [[a[0][0], None],
[None, a_p11]]
# ### Nested matrix solver
#
# We assemble the bilinear form into a nested matrix `A`, and call the
# `assemble()` method to communicate shared entries in parallel. Rows
# and columns in `A` that correspond to degrees-of-freedom with
# Dirichlet boundary conditions wil be zeroed by the assembler, and a
# value of 1 will be set on the diagonal for these rows.
def nested_iterative_solver():
"""Solve the Stokes problem using nest matrices and an iterative solver."""
# Assemble nested matrix operators
A = fem.petsc.assemble_matrix_nest(a, bcs=bcs)
A.assemble()
# Create a nested matrix P to use as the preconditioner. The
# top-left block of P is shared with the top-left block of A. The
# bottom-right diagonal entry is assembled from the form a_p11:
P11 = fem.petsc.assemble_matrix(a_p11, [])
P = PETSc.Mat().createNest([[A.getNestSubMatrix(0, 0), None], [None, P11]])
P.assemble()
A00 = A.getNestSubMatrix(0, 0)
A00.setOption(PETSc.Mat.Option.SPD, True)
P00, P11 = P.getNestSubMatrix(0, 0), P.getNestSubMatrix(1, 1)
P00.setOption(PETSc.Mat.Option.SPD, True)
P11.setOption(PETSc.Mat.Option.SPD, True)
# Assemble right-hand side vector
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 for vector entries that are share across
# parallel processes
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 vector
bcs0 = fem.bcs_by_block(extract_function_spaces(L), bcs)
fem.petsc.set_bc_nest(b, bcs0)
# The pressure field is determined only up to a constant. We supply
# a vector that spans the nullspace to the solver, and any component
# of the solution in this direction will be eliminated during the
# solution process.
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 that spans the nullspace, create a nullspace
# object, and attach it to the matrix
null_vec.normalize()
# Create a MINRES Krylov solver and a block-diagonal preconditioner
# using PETSc's additive fieldsplit 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. For the top-left
# Laplace-type operator we use algebraic multigrid. For the
# lower-right block we use a Jacobi preconditioner. By default, GAMG
# will infer the correct near-nullspace from the matrix block size.
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")
# Create finite element {py:class}`Function <dolfinx.fem.Function>`s
# for the velocity (on the space `V`) and for the pressure (on the
# space `Q`). The vectors for `u` and `p` are combined to form a
# nested vector and the system is solved.
u, p = Function(V), Function(Q)
x = PETSc.Vec().createNest([la.create_petsc_vector_wrap(u.x),
la.create_petsc_vector_wrap(p.x)])
ksp.solve(b, x)
#save mesh and u with adio4dolfinx
checkpoint_file = pathlib.Path("function_checkpoint.bp")
adx.write_mesh(msh, checkpoint_file, engine="BP4")
adx.write_function(u, checkpoint_file, engine="BP4")
# Save solution to file in XDMF format for visualization, e.g. with
# ParaView. Before writing to file, ghost values are updated using
# `scatter_forward`.
with XDMFFile(MPI.COMM_WORLD, "out_stokes/velocity.xdmf", "w") as ufile_xdmf:
u.x.scatter_forward()
P1 = element("Lagrange", msh.basix_cell(), 1, shape=(msh.geometry.dim,))
u1 = Function(functionspace(msh, P1))
u1.interpolate(u)
ufile_xdmf.write_mesh(msh)
ufile_xdmf.write_function(u1)
with XDMFFile(MPI.COMM_WORLD, "out_stokes/pressure.xdmf", "w") as pfile_xdmf:
p.x.scatter_forward()
pfile_xdmf.write_mesh(msh)
pfile_xdmf.write_function(p)
# Compute norms of the solution vectors
norm_u = u.x.norm()
norm_p = p.x.norm()
if MPI.COMM_WORLD.rank == 0:
print(f"(A) Norm of velocity coefficient vector (nested, iterative): {norm_u}")
print(f"(A) Norm of pressure coefficient vector (nested, iterative): {norm_p}")
return norm_u, norm_p
# Solve using PETSc MatNest
norm_u_0, norm_p_0 = nested_iterative_solver()
code 2：update the phase-field functions and variables：
from mpi4py import MPI
from petsc4py import PETSc
import ufl
import adios4dolfinx
import numpy as np
import dolfinx
import adios4dolfinx as adx
import pathlib
from pathlib import Path
from dolfinx import fem, io, default_scalar_type
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, LinearProblem
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary, locate_entities
from basix.ufl import element
from dolfinx.fem import functionspace, Function ,Constant, assemble_scalar, form
from ufl import dot, div, dx, grad, inner, ds
# Read msh
checkpoint_file = Path("function_checkpoint.bp")
msh = adios4dolfinx.read_mesh(MPI.COMM_WORLD, checkpoint_file, engine='BP4', ghost_mode=dolfinx.mesh.GhostMode.shared_facet)
# Read u
P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
U = functionspace(msh, P2)
u = Function(U)
adios4dolfinx.read_function(u, checkpoint_file, engine="BP4")
P1 = element("Lagrange", msh.basix_cell(), 1)
V = functionspace(msh, P1)
#define phi_0
W = functionspace(msh, ("DG", 0))
def condition_0(x):
return np.logical_or(x[0] <= 0.25, np.logical_and(x[1] >= 1/3, x[1] <= 2/3))
def condition_1(x):
return np.logical_or(np.logical_and(x[0] > 0.25,x[1] < 1/3),np.logical_and(x[0] > 0.25,x[1] > 2/3))
phi_0 = Function(W)
cells_0 = locate_entities(msh, msh.topology.dim, condition_0)
cells_1 = locate_entities(msh, msh.topology.dim, condition_1)
phi_0.x.array[cells_0] = np.full_like(cells_0, 1, dtype=default_scalar_type)
phi_0.x.array[cells_1] = np.full_like(cells_1, 0, dtype=default_scalar_type)
with io.XDMFFile(MPI.COMM_WORLD, "output.xdmf", "w") as file:
file.write_mesh(msh)
file.write_function(phi_0)
# Assign to the parameter
alpha_0= 0.2
epsilon = 0.01
eta = 0.01
t = 1.0
beta_0 = 1.0
S = 0.5
beta = 0.5
V_0 = 1.0
#Start the inner iteration
k = 0
while k <= 10:
# Assign to λ
lamda = 0
if k == 0:
lamda_next = lamda
else:
lamda_next = lamda_next
# Define variational problem
phi = ufl.TrialFunction(V)
phi_test = ufl.TestFunction(V)
if k == 0:
velocity_effect_1 = 0.5 * alpha_0 * dot(u, u) + S
velocity_effect_2 = - eta / epsilon * phi_0 * (1.0-phi_0) * (0.5-phi_0) + alpha_0 * dot(u, u) - lamda
velocity_effect_3 = S - 0.5 * alpha_0 * dot(u, u)
a = 1.0 / t * (inner(phi, phi_test) * dx)+ epsilon * eta * (inner(grad(phi), grad(phi_test)) * dx) + inner(velocity_effect_1 * phi, phi_test) * dx
L = 1.0 / t * (inner(phi_0, phi_test) * dx) + inner(velocity_effect_2, phi_test) * dx + inner(velocity_effect_3 * phi_0, phi_test) * dx
else:
velocity_effect_1 = alpha_0 * 0.5 * dot(u, u)+ S
velocity_effect_2 = - eta / epsilon * phi_next * (1.0-phi_next) * (0.5-phi_next) + alpha_0 * dot(u, u) - lamda
velocity_effect_3 = S - 0.5 * alpha_0 * dot(u, u)
a = 1.0 / t * (inner(phi, phi_test) * dx)+ epsilon * eta * (inner(grad(phi), grad(phi_test)) * dx) + inner(velocity_effect_1 * phi, phi_test) * dx
L = 1.0 / t * (inner(phi_next, phi_test) * dx) + inner(velocity_effect_2, phi_test) * dx + inner(velocity_effect_3 * phi_next, phi_test) * dx
#Solve out the problem
problem = LinearProblem(a, L, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
phi_next = problem.solve()
#Cauculate ∫_D φ_{h} ^{n + 1, k + 1, *} dx
if k == 0:
integral_form_1 = form(phi_0 * dx)
integral_value_1 = assemble_scalar(integral_form_1)
else:
integral_form_1 = form(phi_next * dx)
integral_value_1 = assemble_scalar(integral_form_1)
#Restrict the phase-field function in [0, 1] by using a cut-off postprocessing strategy
phi_array = phi_next.x.array
phi_array[phi_array < 0] = 0
phi_array[phi_array > 1] = 1
phi_next.x.array[:] = phi_array
#visualize phi_next
filename = f"out_phase/phase_{k:03d}.xdmf"
with io.XDMFFile(MPI.COMM_WORLD, filename, "w") as phifile_xdmf:
phi_next.x.scatter_forward()
phifile_xdmf.write_mesh(msh)
phifile_xdmf.write_function(phi_next)
#update the parameter λ
def update_lamda():
#Cauculate ∫_D φ_{h} ^{n + 1, k + 1} dx
integral_form_2 = form(phi_next * dx)
integral_value_2 = assemble_scalar(integral_form_2)
#Cauculate Jv(φ_{h} ^{n + 1, k + 1,*}) and Jv(φ_{h} ^{n + 1, k + 1})
Jv_1 = integral_value_1 - beta * V_0
Jv_2 = integral_value_2 - beta * V_0
delta_Jv = Jv_1 - Jv_2
#Cauculate β_0^*
beta_0_star = lamda_next * (Jv_2 - Jv_1) / (Jv_2 * Jv_2)
#update λ
if delta_Jv >= 0:
lamda_new = lamda_next - beta_0 * Jv_2
else:
lamda_new = lamda_next - beta_0_star * Jv_2
print(lamda_new)
return lamda_new
lamda_next = update_lamda()
k = k + 1
The output is supposed to be:
But the output of my code is:
visuallized phase-field function phi:
You can see difference between my output and the expected output. One of the differences is that the value of the expected phase-field function φ at the boundary y = 0 and y = 1 is zero, but the φ I solve out didnt appear that.
I had checked my code many times but coudnt figure out The error part in the my code.