Hello everyone,
I am working with a code that numerically solves a system of parabolic pde.
The proposed code seems to disregard the bc I am trying to impose with some difficulty. In fact if for the first 2 variables, phi and mu I don’t have to enter anything since they are subject to Neumann homogeneous bc, as far as n is concerned I need to enter Dirichlet conditions and I have it as in underlined lines.
from dolfin import *
import numpy as np
import sys
# Compiler options
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
parameters["form_compiler"]["quadrature_degree"] = 4
mesh= Mesh("brain.xml")
#--------------------------------- Tensors D and T ----------------------------------------#
# Code for C++ evaluation of DTI and construction of tensor D
defineMatrix_code_D = """
#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>
namespace py = pybind11;
#include <dolfin/function/Expression.h>
#include <dolfin/mesh/MeshFunction.h>
class Components_DT_D : public dolfin::Expression
{
public:
// Create expression with 6 components
Components_DT_D() : dolfin::Expression(6) {}
// Function for evaluating expression on each cell
void eval(Eigen::Ref<Eigen::VectorXd> values, Eigen::Ref<const Eigen::VectorXd> x, const ufc::cell& cell) const override
{
const uint topDim = cell.topological_dimension;
const uint cell_index = cell.index;
values[0] = (*d11)[cell_index];
values[1] = (*d12)[cell_index];
values[2] = (*d13)[cell_index];
values[3] = (*d22)[cell_index];
values[4] = (*d23)[cell_index];
values[5] = (*d33)[cell_index];
}
// The data stored in mesh functions
std::shared_ptr<dolfin::MeshFunction<double> > d11;
std::shared_ptr<dolfin::MeshFunction<double> > d12;
std::shared_ptr<dolfin::MeshFunction<double> > d13;
std::shared_ptr<dolfin::MeshFunction<double> > d22;
std::shared_ptr<dolfin::MeshFunction<double> > d23;
std::shared_ptr<dolfin::MeshFunction<double> > d33;
};
PYBIND11_MODULE(SIGNATURE, m)
{
py::class_<Components_DT_D, std::shared_ptr<Components_DT_D>, dolfin::Expression>
(m, "Components_DT_D")
.def(py::init<>())
.def_readwrite("d11", &Components_DT_D::d11)
.def_readwrite("d12", &Components_DT_D::d12)
.def_readwrite("d13", &Components_DT_D::d13)
.def_readwrite("d22", &Components_DT_D::d22)
.def_readwrite("d23", &Components_DT_D::d23)
.def_readwrite("d33", &Components_DT_D::d33);
}
"""
# Define DT components expression and matrix
d11 = MeshFunction("double",mesh,"d11S.xml.gz")
d22 = MeshFunction("double",mesh,"d22S.xml.gz")
d33 = MeshFunction("double",mesh,"d33S.xml.gz")
d12 = MeshFunction("double",mesh,"d12S.xml.gz")
d13 = MeshFunction("double",mesh,"d13S.xml.gz")
d23 = MeshFunction("double",mesh,"d23S.xml.gz")
d11_file_pvd = File("d11S.pvd","compressed")
d22_file_pvd = File("d22S.pvd","compressed")
d33_file_pvd = File("d33S.pvd","compressed")
d12_file_pvd = File("d12S.pvd","compressed")
d13_file_pvd = File("d13S.pvd","compressed")
d23_file_pvd = File("d23S.pvd","compressed")
d11_file_pvd << d11
d22_file_pvd << d22
d33_file_pvd << d33
d12_file_pvd << d12
d13_file_pvd << d13
d23_file_pvd << d23
d = CompiledExpression(compile_cpp_code(defineMatrix_code_D).Components_DT_D(), d11 = d11, d12 = d12, d13 = d13, d22 = d22, d23 = d23, d33 = d33, degree=2)
# Diffusion tensor
D0 = as_matrix([ [d[0], d[1], d[2]], [d[1], d[3], d[4]], [d[2], d[4], d[5]] ])
# Code for C++ evaluation of DTI and construction of tensor T
defineMatrix_code_T = """
#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>
namespace py = pybind11;
#include <dolfin/function/Expression.h>
#include <dolfin/mesh/MeshFunction.h>
class Components_DT_T : public dolfin::Expression
{
public:
// Create expression with 6 components
Components_DT_T() : Expression(6) {}
// Function for evaluating expression on each cell
void eval(Eigen::Ref<Eigen::VectorXd> values, Eigen::Ref<const Eigen::VectorXd> x, const ufc::cell& cell) const override
{
const uint topDim = cell.topological_dimension;
const uint cell_index = cell.index;
values[0] = (*t11)[cell_index];
values[1] = (*t12)[cell_index];
values[2] = (*t13)[cell_index];
values[3] = (*t22)[cell_index];
values[4] = (*t23)[cell_index];
values[5] = (*t33)[cell_index];
}
// The data stored in mesh functions
std::shared_ptr<dolfin::MeshFunction<double> > t11;
std::shared_ptr<dolfin::MeshFunction<double> > t12;
std::shared_ptr<dolfin::MeshFunction<double> > t13;
std::shared_ptr<dolfin::MeshFunction<double> > t22;
std::shared_ptr<dolfin::MeshFunction<double> > t23;
std::shared_ptr<dolfin::MeshFunction<double> > t33;
};
PYBIND11_MODULE(SIGNATURE, m)
{
py::class_<Components_DT_T, std::shared_ptr<Components_DT_T>, dolfin::Expression>
(m, "Components_DT_T")
.def(py::init<>())
.def_readwrite("t11", &Components_DT_T::t11)
.def_readwrite("t12", &Components_DT_T::t12)
.def_readwrite("t13", &Components_DT_T::t13)
.def_readwrite("t22", &Components_DT_T::t22)
.def_readwrite("t23", &Components_DT_T::t23)
.def_readwrite("t33", &Components_DT_T::t33);
}
"""
# Define DT components expression and matrix
t11 = MeshFunction("double",mesh,"t11S.xml.gz")
t22 = MeshFunction("double",mesh,"t22S.xml.gz")
t33 = MeshFunction("double",mesh,"t33S.xml.gz")
t12 = MeshFunction("double",mesh,"t12S.xml.gz")
t13 = MeshFunction("double",mesh,"t13S.xml.gz")
t23 = MeshFunction("double",mesh,"t23S.xml.gz")
t11_file_pvd = File("t11S.pvd","compressed")
t22_file_pvd = File("t22S.pvd","compressed")
t33_file_pvd = File("t33S.pvd","compressed")
t12_file_pvd = File("t12S.pvd","compressed")
t13_file_pvd = File("t13S.pvd","compressed")
t23_file_pvd = File("t23S.pvd","compressed")
t11_file_pvd << t11
t22_file_pvd << t22
t33_file_pvd << t33
t12_file_pvd << t12
t13_file_pvd << t13
t23_file_pvd << t23
tmat = CompiledExpression(compile_cpp_code(defineMatrix_code_T).Components_DT_T(), t11 = t11, t12 = t12, t13 = t13, t22 = t22, t23 = t23, t33 = t33, degree=2)
# Tensor of preferential directions
mat_T = as_matrix([ [tmat[0], tmat[1], tmat[2]], [tmat[1], tmat[3], tmat[4]], [tmat[2], tmat[4], tmat[5]]])
# ---------------------- Define physical and numerical parameters --------------------------- #
# Physical parameters
nu = Constant(1/345600)
M_0 = Constant(390625*8.64*pow(10,10))
epsilon = Constant(0.00079)
kappa = Constant(622.7)
delta = Constant(0.15)
D_n = Constant(pow(10,-10)/8.64)
delta_n = Constant(15552/86400)
S_n = Constant(25/216)
print ('Evolution of the tumor up to the third cycle of chemo')
Day_after_surgery = int(input('How many days after the surgery did the treatment start? '))
T_tot = 131 + Day_after_surgery
# Numerical parameters (time loop)
#dt = 0.1
deltaT = (T_tot/301)*86400
dt1 = Constant(deltaT)
dx = Measure("dx", domain=mesh)
# ---------------------- Define finite elements and function spaces --------------------------- #
# Finite Element function spaces
V = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, V)
# FE-space to determine the initial condition of phi
element1 = MixedElement([V,V,V])
ME1 = FunctionSpace(mesh, element1)
# Trial functions at time n+1
solution1 = Function(ME1)
phi, mu, n = split(solution1)
# Trial functions at time n
solution1_old = Function(ME1)
phi_old, mu_old, n_old = split(solution1_old)
# Trial functions for derivatives
d_solution = TrialFunction(ME1)
d_phi, d_mu, d_n = split(d_solution)
# Test functions
varphi, vv, q = TestFunctions(ME1)
# ---------------------- Cahn-Hilliard Class ----------------------------------------------- #
# Interface-class for the Newton solver
class CahnHilliardEquation(NonlinearProblem):
def __init__(self, a, L, bc_n):
#def __init__(self, a, L):
NonlinearProblem.__init__(self)
self.L = L
self.a = a
self.bc_n = bc_n
self.reset_sparsity = True
def F(self, b, x):
assemble(self.L, tensor=b)
self.bc.apply(b, x)
def J(self, A, x):
assemble(self.a, tensor=A)
self.bc.apply(A)
self.reset_sparsity = False
# ---------------------- Initial Conditions ----------------------------------------------- #
# Define analytically the initial condition for phi
phi_initial = TrialFunction(ME)
psi = TestFunction(ME)
# function to test the GBM code
test1 = '2*exp(-pow((x[0]-195.56295),2)-pow((x[1]-313.18395),2)-pow((x[2]-15.66619),2))-1'
phi_initial = Expression(test1, degree = 2)
# Initial data
phi0 = interpolate(phi_initial, ME1.sub(0).collapse())
mu0 = interpolate(Expression(("0.0"), degree=1), ME1.sub(1).collapse())
n0 = interpolate(Expression(("1.0"), degree=1), ME1.sub(2).collapse())
assign(solution1_old, [phi0, mu0, n0])
assign(solution1, [phi0, mu0, n0])
**#Boundary condition for nutrient n**
**n0_bc = Constant(1.0)**
**bc_n = DirichletBC(ME1.sub(2), n0_bc, mesh)**
# Define nutrient initial condition by solving the stationary equation
n0_new = TrialFunction(ME)
w = TestFunction(ME)
a_n = D_n * inner(grad(n0_new), grad(w)) * dx \
+ (1/3) * S_n * n0_new * (2-phi_old) * w * dx \
+ delta_n * n0_new * (1/2) * (1+phi_old) * w * dx
L_n = (1/3) * S_n * (2-phi_old) * w * dx
# Compute solution
print('Computing nutrient initial condition.')
n0_new = Function(ME)
problem_n = LinearVariationalProblem(a_n, L_n, n0_new)
solver_n = LinearVariationalSolver(problem_n)
prm = solver_n.parameters
prm['linear_solver'] = 'mumps'
solver_n.solve()
assign(solution1_old.sub(2), n0_new)
# Saving initial conditions on a VTK file
vtkfile_phi = File("output/condition_phi.pvd")
vtkfile_mu = File("output/condition_mu.pvd")
vtkfile_n = File("output/condition_n.pvd")
vtkfile_phi << solution1_old.sub(0)
vtkfile_mu << solution1_old.sub(1)
vtkfile_n << solution1_old.sub(2)
# ---------------------- Time Loop ----------------------------------------------- #
volume = []
time_v = []
# Initializing the variables to write on the VTK file
vtkfile1 = File("output/phi_evolution.pvd")
vtkfile2 = File("output/nutrient_evolution.pvd")
t1 = 0.0
k_c=5
k_r=10
for time in range(301):
t1 += deltaT
# Compute the tumor volume at each time iteration
phi_c = 0.5 * (1 + phi_old)
volume_c = assemble(phi_c * dx)
volume.append(volume_c)
time_v.append(t1)
np.savetxt("output/tumor_volume", volume) # save the current tumor volume on a txt
np.savetxt("output/time", time_v) # save the correspondent timestep on a txt
F0 = phi * varphi * dx - phi_old * varphi * dx \
+ deltaT * (1/M_0) * phi_old * pow((1-phi_old),2) * inner(mat_T * grad(mu), grad(varphi)) * dx \
- deltaT * nu * (n - delta) * (1/2) * (1 + phi_old) * varphi * dx \
+ deltaT * (k_r + k_c) * phi_old * varphi * dx \
F1 = mu * vv * dx - pow(epsilon, 2) * inner(grad(phi), grad(vv)) * dx \
- kappa * pow(phi, 3) * vv * dx + kappa * phi_old * vv * dx
F2 = n * q * dx - n_old * q * dx + deltaT * D_n * inner(D0*grad(n), grad(q)) * dx \
- deltaT * S_n * (1 - n) * (1/3) * (2 - phi_old) * q * dx \
+ deltaT * delta_n * n * (1/2) * (1 + phi_old) * q * dx
# system assembly
FF = F0 + F1 + F2
# Compiler options
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
#parameters["form_compiler"]["quadrature_degree"] = 2
# Nonlinear solver: Newton's method
JJ = derivative(FF, solution1, d_solution)
problem = CahnHilliardEquation(JJ, FF, [bc_n])
solver = NewtonSolver()
prm = solver.parameters
prm["linear_solver"] = "mumps"
prm["convergence_criterion"] = "incremental"
prm["relative_tolerance"] = 1e-6
# solution for CH + nutrient
print('Solving CH + nutrient')
solver.solve(problem, solution1.vector())
#writing the solution on vtk (every 10 timesteps)
if (time % 10 == 0):
vtkfile1 << (solution1.split()[0], t1)
vtkfile2 << (solution1.split()[2], t1)
# assignment
solution1_old.vector()[:] = solution1.vector()
but I get the following error:
File "/Users/lorenzomarta/Desktop/Script simulazioni/Test17/GBM_Besta_17.py", line 290, in <module>
bc_n = DirichletBC(ME1.sub(2), n0_bc, mesh)
File "/Applications/anaconda3/envs/fenicsproject/lib/python3.10/site-packages/dolfin/fem/dirichletbc.py", line 116, in __init__
raise RuntimeError("Invalid argument")
RuntimeError: Invalid argument