Hi Bhavesh,
Thank you for your answer. Here is the exact code of the entire solver. Everything is carefully commented, but the relevant part is the variational formulation, which is where I want to be calling diffu(r+R) instead of diffu(r) (which I am doing now since Python won’t let me insert R as an argument of the function).
Could you perhaps indicate where and how I should make the change in the code?
Running the code should give you two solution plots, one for c and one for R. For example, try calling the agentonpores(T, num_steps, r, a, b, R_0, c_0, beta, gamma, bc_bottom) function with T = 3, num_steps = 100, a = b = 1, r = 0.1, beta = gamma = 1, R_0 = 0.1, c_0 = 1.
from __future__ import print_function
from fenics import * # imports the key classes UnitSquareMesh, FunctionSpace, Function etc.
import numpy as np # we'll need this, always - most matrix-y operations use numpy
import sympy as sym
from matplotlib import pyplot as plt # go-to plotter, MATLAB style
from scipy.interpolate import interp1d
def agentonpores(T, num_steps, r, a, b, R_0, c_0, beta, gamma, bc_bottom):
dt = T / num_steps # time step size (predetermined by the number of timesteps)
degree = 1 # the degree that determines the kind of elements we are using for the function space
def vol(R,r): # Area of cleanser occupied region.
return 1 - np.pi *(r+R)**2
def leng(R,r): # Local length of agent-cleanser interface.
return 2*np.pi*(r+R)
def diffu(r): # We define a function which outputs the diffusivity for a given radius
R1 = np.linspace(0,0.47,48) # Radii 0, 0.1, 0.2, ... , 0.47
R2 = np.linspace(0.471,0.5,30) # Smaller steps for radii: 0.471, 0.472, ... , 0.5
R = np.concatenate((R1, R2)) # All 78 arguments for calculations placed into one array.
# Effective diffusivities for the particular radii of pores.
D = np.array([1.0000, 0.9997, 0.9988, 0.9972, 0.9950, 0.9922, 0.9888, 0.9849, 0.9803, 0.9752,
0.9696, 0.9634, 0.9567, 0.9496, 0.9420, 0.9340, 0.9256, 0.9168, 0.9076, 0.8982,
0.8884, 0.8783, 0.8680, 0.8574, 0.8467, 0.8357, 0.8246, 0.8133, 0.8019, 0.7904,
0.7787, 0.7668, 0.7548, 0.7427, 0.7303, 0.7177, 0.7048, 0.6915, 0.6776, 0.6631,
0.6476, 0.6310, 0.6127, 0.5921, 0.5684, 0.5401, 0.5051, 0.4596, 0.4543, 0.4487,
0.4430, 0.4370, 0.4309, 0.4244, 0.4178, 0.4109, 0.4036, 0.3961, 0.3883, 0.3800,
0.3714, 0.3624, 0.3528, 0.3428, 0.3322, 0.3209, 0.3089, 0.2962, 0.2824, 0.2677,
0.2516, 0.2340, 0.2145, 0.1925, 0.1671, 0.1365, 0.0962, 0.0000])
# We use cubic interpolation in order to smooth out the resulting curve
f_cubic = interp1d(R, D, 'cubic')
return float(f_cubic(r))
# Create mesh and define function space
nx = ny = 100 # Specifying the resolution of our triangular mesh
mesh = RectangleMesh(Point(0.0, 0.0), Point(a, b), nx, ny, "crossed")
# Uniform, finite element mesh over the unit square.
# The mesh consists of cells, which in 2D are triangles with straight lines.
# Defining a MixedElement as the product space of two simple finite elements.
# Afterwards, using the mixed element to define the function space.
# V = VectorFunctionSpace(mesh, 'P', 2, degree)
P1 = FiniteElement('P', triangle, degree)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)
# 'P' means that we are using the standard Lagrange family of elements.
# The third argument specifies the degree of the finite element.
# Patrick says that this is a sensible choice and should do the job, so I will stick with it.
# Define boundary condition - there is no cleanser initially at t = 0!
# The first condition corresponds to no cleanser, the second to the initial radius.
u_0 = Expression((str(c_0), str(R_0)), degree=1)
tol = 1e-14 # a tiny tolerance that helps us locate the suitable parts of the boundary.
class BoundaryX0(SubDomain): # defining the boundary at x = 0
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0, tol)
class BoundaryXa(SubDomain): # defining the boundary at x = a
def inside(self, x, on_boundary):
return on_boundary and near(x[0], a, tol)
class BoundaryY0(SubDomain): # defining the boundary at y = 0
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0, tol)
class BoundaryYb(SubDomain): # defining the boundary at y = b
def inside(self, x, on_boundary):
return on_boundary and near(x[1], b, tol)
# We are working with facets, so from the FEniCS documentation we need -1 of the dimension of the mesh topology.
# If instead we wanted to work with vertices, the dimension would be -2.
boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim()-1, 0)
# This was previously done with a FacetFunction, which has now deprecated within FEniCS.
boundary_markers.set_all(9999)
bx0 = BoundaryX0()
bxa = BoundaryXa()
by0 = BoundaryY0()
byb = BoundaryYb()
bx0.mark(boundary_markers, 0) # Assigning a "0" marker to x = 0
bxa.mark(boundary_markers, 1) # Assigning a "1" marker to x = 1
by0.mark(boundary_markers, 2) # Assigning a "2" marker to y = 0
byb.mark(boundary_markers, 3) # Assigning a "3" marker to y = 1
# In this dictionary, we specify the kinds of boundary conditions we are dealing with.
# There are three possibilities: Dirichlet, Neumann and Robin.
# These conditions refer to the cleanser concentration, not to the radius of agent.
# For each kind of condition, to the right we specify the function that this condition is assigned.
boundary_conditions_cleanser = {0: {'Neumann': Expression('0', degree=1)}, # insulated left end
1: {'Neumann': Expression('0', degree=1)}, # insulated right end
2: {'Dirichlet': Expression('1', degree=1)}, # applied cleanser concentration on surface
3: {bc_bottom : Expression('0', degree=1)}} # insulated maximum depth / no cleanser
bcs = [] # We first start with an empty array of boundary conditions, which we then fill.
# Next, we fill this array only with Dirichlet conditions, since these are the only ones which are taken
# into account explicitly, as opposed to Neumann conditions within the variational formulation itself.
# We apply the boundary conditions to V.sub(0) - the cleanser (indexing in Python starts from zero)
# This next for loop skips over anything that is not Dirichlet
for i in boundary_conditions_cleanser:
if 'Dirichlet' in boundary_conditions_cleanser[i]:
bc = DirichletBC(V.sub(0), boundary_conditions_cleanser[i]['Dirichlet'],
boundary_markers, i)
bcs.append(bc)
# Define initial value for the solution vector, using the initial condition we set above.
u_n = interpolate(u_0, V)
#c_n = project(u_D, V), depending on whether we want to interpolate or project the boundary data.
u = interpolate(u_0, V) # A good initial guess = initial value! Newton is local, so we start close.
v_1, v_2 = TestFunctions(V) # Test functions for cleanser PDE and radius PDE, respectively.
c, R = split(u) # The components of the solution vector (cleanser first, then radius).
c_n, R_n = split(u_n) # The components of the discretised iterations, as above.
c_= 0.5*(c+c_n) # The "c tilde" for the implicit midpoint method
# Variational formulation
F = (c-c_n)/dt*v_1*vol(R,r)*dx + dot(grad(c_), grad(v_1))*vol(R,r)*diffu(r)*dx \
+ v_1*leng(R,r)*beta*c_*(1+gamma*c_)*dx + v_2*(R-R_n)/dt*dx \
+ v_2*beta*gamma*c_*dx
# language used to express this syntax is UFL (Unified Form Language) and is integral to FEniCS.
# Time-stepping
t = 0
for n in range(num_steps):
# Update current time
t += dt # add the timestep to the current time to move forward
u_0.t = t # update the time
# Compute solution
solve(F == 0, u, bcs) # Apply Newton-Kantorovich to nonlinear problem, with the specified boundary conditions
# Plot solution at the end of the simulation
plt.figure(0)
simulationplot1 = plot(c, title = "Cleanser concentration at end of simulation", mode = 'color')
plt.colorbar(simulationplot1)
plt.xlabel("Horizontal distance $(y)$")
plt.ylabel("Depth within porous medium $(x)$")
plt.figure(1)
simulationplot2 = plot(R, title = "Radius of agent coating at end of simulation", mode = 'color')
plt.colorbar(simulationplot2)
plt.xlabel("Horizontal distance $(y)$")
plt.ylabel("Depth within porous medium $(x)$")
