Hi everyone,
Im trying to reproduce an algorithm for shape optimization in stokes flows by phase-field approaches. It is actually a simple algorithm but I dont know why I cant get the right output.
Here is the algorithm:
STEP 1. Update the velocity field u_h\in H_1(D) \cap P_2(T_h) by the following state equation using the Taylor-Hood element,
(\nabla u_h, \nabla w_h) + (\nabla p_h, w_h) + (\alpha(\phi_0) u_h, w_h) = 0 , w_h \in H_0^1(D) \cap P_2(T_h)
(\nabla \cdot u_h, q_h) = 0 , q_h \in L^2(D) \cap P_1(T_h)
Where \alpha(\phi) = \alpha_0(1 - \phi)^2
STEP 2. Set \phi_h^{0} = \phi_0, \lambda = 0 and k = 0.
While k < K do
(1) Update the phase-field variable with the piecewise linear element and obtain \phi_h^{k+1,*},
\frac{(\phi_h^{k+1,*}, \psi_h)}{\Delta t} + \epsilon \eta (\nabla \phi_h^{k+1,*}, \nabla \psi_h) + (\left(\frac{\alpha_0 |u_h|^2}{2} + \tilde{S}\right) \phi_h^{k+1,*}, \psi_h) =
\frac{(\phi_h^{k}, \psi_h)}{\Delta t} + (-\frac{\eta}{\epsilon}f(\phi_h^k) + \alpha_0|u_h|^2 - \lambda,\psi_h) + ((\tilde{S} - \frac{\alpha_0 |u_h|^2}{2})\phi_h^k), \psi_h) , \psi_h \in H_0^1(D) \cap P_1(T_h)
Where f(\phi) = \phi (1 - \phi) (1/2 - \phi)
(2) Restrict the phase-field function in [0, 1] by using a cut-off postprocessing strategy and obtain \phi_h^{k+1}
\phi_h^{k+1}(x) = P(\phi_h^{k+1,*}(x)) := \min(\max(\phi_h^{k+1,*}(x), 0), 1)
For STEP 1,I refer the demo downloaded from Stokes equations using Taylor-Hood elements and the u solved out in STEP1 is correct. But when I try to run the codes I write for STEP2, the output is absolutely wrong.
My codes ( some steps of solving stokes equation are omitted):
......
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)
with XDMFFile(MPI.COMM_WORLD, "out_stokes_1/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_1/pressure.xdmf", "w") as pfile_xdmf:
p.x.scatter_forward()
pfile_xdmf.write_mesh(msh)
pfile_xdmf.write_function(p)
P = functionspace(msh, ("Lagrange", 1))
(phi_next, phi_test) = ufl.TrialFunction(P), ufl.TestFunction(P)
k = 0
lamda = 0
while k < 10:
if k == 0:
velocity_effect_1 = 0.5 * alpha_0 * inner(u, u) + S
velocity_effect_2 = - eta / epsilon * phi_0 * (1.0-phi_0) * (0.5-phi_0) + alpha_0 * inner(u, u) - lamda
velocity_effect_3 = S - 0.5 * alpha_0 * inner(u, u)
a1 = 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
L1 = 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 * inner(u, u)+ S
velocity_effect_2 = - eta / epsilon * phi_next * (1.0-phi_next) * (0.5-phi_next) + alpha_0 * inner(u, u) - lamda
velocity_effect_3 = S - 0.5 * alpha_0 * inner(u, u)
a1 = 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
L1 = 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
problem = LinearProblem(a1, L1, bcs=[], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
phi_next = problem.solve()
phi_next_array = phi_next.x.array
phi_next_array[phi_next_array < 0] = 0
phi_next_array[phi_next_array > 1] = 1
phi_next.x.array[:] = phi_next_array
filename1 = f"TEST2/phase_{k:03d}.xdmf"
with XDMFFile(MPI.COMM_WORLD, filename1, "w") as phifile_xdmf:
phi_next.x.scatter_forward()
phifile_xdmf.write_mesh(msh)
phifile_xdmf.write_function(phi_next)
k = k + 1
It is certain that all the parameters and velocity_effect_1,velocity_effect_2,velocity_effect_3 are correctly set since I have checked it for many times. But the output of \phi^k is wrong even afer the first iteration of STEP2.
Initial \phi^0:
My output \phi^1:
And the phase_field function finally become constant 1 over the mesh after 10 times of iteration:
\phi^{10}
Which is abolutely wrong.
I cant figure out why the output is wrong. It seems that the equation in STEP2 looks like a very simple linear problem but why my output is wrong?