Hi!
For my application, I want to solve a phase field based on information from a mechanical problem. I have this information at the quadrature, and so far, I’ve been projecting this to the nodes, solving the phase field, and then projecting the phase field variable back to the quadrature points.
Now I would like to clean this up, and satisfy my equations directly at the quadrature points. First of all I am unsure whether this conceptually makes sense, given a general FEM formulation. Secondly, in the MWE below, I get the phase field to run in a “DG” 0 space, however the result in paraview appears to be noise rather than a phase field solution.
MWE: (switch lines 20-21 and 31-36 to switch between the two options)
import numpy as np
import dolfinx.fem.petsc
from dolfinx import mesh, fem
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from dolfinx.io import XDMFFile
import ufl
print(f"dolfinx version: {dolfinx.__version__}")
dt = 2e-1
epsilon = 1e-2
mesh_domain = mesh.create_unit_square(MPI.COMM_WORLD, 50, 50, mesh.CellType.triangle)
"""
Function spaces
"""
# V = fem.functionspace(mesh_domain, ("CG", 1)) # works as expected
V = fem.functionspace(mesh_domain, ("DG", 0)) # gives noisy outcome
phi = fem.Function(V, name="phase")
phi_prev = fem.Function(V, name="phase_old")
v = ufl.TestFunction(V)
"""
Phase field formulation
"""
# Phase-field energy
# Corresponding to CG 1
# F = (phi - phi_prev) / dt * v * ufl.dx + epsilon**2 * ufl.inner(ufl.grad(phi), ufl.grad(v)) * ufl.dx + 2 * phi * ((1-phi)**2 - phi * (1-phi)) * v * ufl.dx
# Corresponding to DG 0
dx = ufl.Measure("dx", domain=mesh_domain)
F = (phi - phi_prev) / dt * v * dx + epsilon**2 * ufl.inner(ufl.grad(phi), ufl.grad(v)) * dx + 2 * phi * ((1-phi)**2 - phi * (1-phi)) * v * dx # ufl.dx no longer is a valid integration domain. New "dx" is.
"""
Initial conditions
"""
rng = np.random.default_rng(42)
phi.interpolate(lambda x: 0.5 + 0.5 * (0.5 - rng.random(x.shape[1]))) # random
"""
Solver setup: Create nonlinear problem and Newton solver
"""
problem = NonlinearProblem(F, phi)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
"""
Solving
"""
file = XDMFFile(MPI.COMM_WORLD, "phasefield_demo.xdmf", "w")
file.write_mesh(mesh_domain)
t = 0.0
timesteps = 20
file.write_function(phi, 0)
for i in range(timesteps):
phi_prev.x.array[:] = phi.x.array
t += dt
r = solver.solve(phi)
print(f"Step {int(t / dt)}: num iterations: {r[0]}")
# Process & save solution
file.write_function(phi, i+1)
file.close()
CG 1 outcome:
DG 0 outcome:
Any help with solving this would be much appreciated!