I am trying to model a phase field model (allen cahn) being advected by a viscoelastic oldroyd-b model. Here is the entire code. Could you let me know the mistake i am making?
from fenics import *
import numpy as np
#mesh and function space
mesh = UnitSquareMesh(64, 64)
V = VectorFunctionSpace(mesh, 'P', 2)#for velocity u
Q = FunctionSpace(mesh, 'P', 1) #for pressure p
S = TensorFunctionSpace(mesh, 'P', 1) #for stress tensor tau
P = FunctionSpace(mesh, 'P', 1) #for phase field phi
#Define current solution functions
W = FunctionSpace(mesh, MixedElement([V.ufl_element(), Q.ufl_element(), S.ufl_element()]))
w = Function(W) # Current solution (u, p, tau)
u, p, tau = split(w)
phi = Function(P) # Current solution (phi)
# Define previous solution functions
u_n = Function(V)
p_n = Function(Q)
tau_n = Function(S)
phi_n = Function(P)
# Define test functions
v, q, chi = TestFunctions(W)
psi = TestFunction(P)
#define material parameters
gamma = 0.1 # Adjusted gamma to control the phase field evolution
def eta_s(phi):
return 1.0 + 0.5 * (1.0 + phi) # Tumor has higher eta_s, phi=1
def eta_p(phi):
return 1.0 + 0.5 * (1.0 - phi) # Non-tumor has higher eta_p, phi=-1
def lambda_(phi):
return 1.0 + 0.5 * (1.0 - phi) # Non-Tumor has higher lambda_, phi=1
def rho(phi):
return 1.0 + 0.5 * (1.0 - phi) # Non-Tumor has higher rho, phi=1
def alpha(phi):
return 0.5 * (1 + phi) # Tumor has higher alpha, phi=1
# Body force (proportional to velocity)
f_body = -alpha(phi)*u # Corrected the sign of the body force
# Define the stress tensor
D = sym(grad(u))
sigma = -p * Identity(len(u)) + 2 * eta_s(phi) * D + tau
# Variational form for polymeric stress evolution
F_tau = (inner((tau - tau_n) / dt, chi) * dx
+ inner(dot(u, nabla_grad(tau)), chi) * dx
- inner(dot(tau, grad(u)), chi) * dx
- inner(dot(grad(u).T, tau), chi) * dx
- (2 * eta_p(phi) / lambda_(phi)) * inner(D, chi) * dx
+ (1 / lambda_(phi)) * inner(tau, chi) * dx)
# Variational form for momentum balance (mixed formulation)
F_u = (rho(phi) * inner((u - u_n) / dt, v) * dx
+ rho(phi) * inner(dot(u, nabla_grad(u)), v) * dx
+ inner(sigma, grad(v)) * dx # Stress tensor applied to gradient of test function
- dot(f_body, v) * dx) # Body force term
# variational form for the Allen-Cahn phase field
epsilon = 0.01 # Interface thickness parameter
F_phi = (inner((phi - phi_n) / dt, psi) * dx
+ inner(dot(u, grad(phi)), psi) * dx
+ gamma * epsilon**2 * inner(grad(phi), grad(psi)) * dx
- gamma * inner(phi**3 - phi, psi) * dx)
# Define the Jacobians for the nonlinear solvers
J_tau = derivative(F_tau, w.sub(2))
J_u = derivative(F_u, w)
J_phi = derivative(F_phi, phi)
# Define solvers with adjusted parameters
newton_solver_parameters = {
"nonlinear_solver": "newton",
"newton_solver": {
"relative_tolerance": 1e-6,
"absolute_tolerance": 1e-10,
"maximum_iterations": 100,
"relaxation_parameter": 1.0,
}
}
# Define solvers
tau_problem = NonlinearVariationalProblem(F_tau, w.sub(2), J=J_tau)
tau_solver = NonlinearVariationalSolver(tau_problem)
tau_solver.parameters.update(newton_solver_parameters)
u_problem = NonlinearVariationalProblem(F_u, w, J=J_u)
u_solver = NonlinearVariationalSolver(u_problem)
u_solver.parameters.update(newton_solver_parameters)
phi_problem = NonlinearVariationalProblem(F_phi, phi, J=J_phi)
phi_solver = NonlinearVariationalSolver(phi_problem)
phi_solver.parameters.update(newton_solver_parameters)
# Initialize the previous time step solutions
u_n.assign(Constant((0.1, 0.1))) # Small initial velocity
p_n.assign(Constant(0.0)) # Initial pressure
tau_n.assign(Constant(((0.01, 0.01), (0.01, 0.01)))) # Small initial stress
# Initialize the phase field with the generalized tanh function
class InitialCondition(UserExpression):
def eval(self, values, x):
R0 = 0.2 # Critical radius
epsilon_ = 0.01 # Controls the width of the transition
r = np.sqrt((x[0] - 0.5)**2 + (x[1] - 0.5)**2)
values[0] = np.tanh((R0 - r) / epsilon_)
def value_shape(self):
return ()
phi.interpolate(InitialCondition(degree=1))
phi_n.assign(phi)
# Parameters
dt = 0.005 # Smaller time step to improve convergence
num_steps = 200 # Increased number of steps due to reduced dt
t = 0
phi_values = []
for n in range(num_steps):
t += dt
print(n,t)
tau_solver.solve()
u_solver.solve()
phi_solver.solve()
phi_values.append(phi.copy(deepcopy=True))
u_sol, p_sol, tau_sol = w.split()
u_n.assign(u_sol.copy(deepcopy=True))
p_n.assign(p_sol.copy(deepcopy=True))
tau_n.assign(tau_sol.copy(deepcopy=True))
phi_n.assign(phi.copy(deepcopy=True))
The error i got is this: ---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
in <cell line: 2>()
5 print(n,t)
6
----> 7 tau_solver.solve()
8
9
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
*** https://fenicsproject.discourse.group/
*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.
*** -------------------------------------------------------------------------
*** Error: Unable to access vector of degrees of freedom.
*** Reason: Cannot access a non-const vector from a subfunction.
*** Where: This error was encountered inside Function.cpp.
*** Process: 0
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: 1c52e837eb54cc34627f90bde254be4aa8a2ae17
*** -------------------------------------------------------------------------