# Problems with Dirichlet conditions of a quantity on a pde system

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"

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<>())
}

"""

# 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<>())
}

"""

# 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)

+ (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

# 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

``````

Why are you sending in the mesh as the third argument?

Are you trying to apply the boundary condition on the whole exterior boundary of the mesh?

If so, you should use the string `"on_boundary"`.

Thank you!
Yes, I have to apply the condition to the whole mesh. I just replaced it and ran the simulation but I have another problem:

``````  File "/Users/lorenzomarta/Desktop/Script simulazioni/Test17/GBM_Besta_17.py", line 408, in <module>
solver.solve(problem, solution1.vector())
File "/Users/lorenzomarta/Desktop/Script simulazioni/Test17/GBM_Besta_17.py", line 246, in F
self.bc.apply(b, x)
AttributeError: 'list' object has no attribute 'apply'

``````

How can I solve it? Before inserting bc this problem did not occur.

This should be

``[bc.apply(b, x) for bc in self.bc]``
1 Like

It seems to be working properly.
Thank you