Let me start by saying that I have found similar problems to mine on the NARKIVE FiPy mailing list archive but since the equations won't load, they are not very useful. For example Convection-diffusion problem on a 1D cylindrical grid, or on another mailing list archive Re: FiPy Heat Transfer Solution. In the second linked mail Daniel says:
There are two ways to solve on a cylindrical domain in FiPy. You can either use the standard diffusion equation in Cartesian coordinates (2nd equation below) and with a mesh that is actually cylindrical in shape or you can use the diffusion equation formulated on a cylindrical coordinate system (1st equation below) and use a standard 2D / 1D grid mesh.
And the equations are not there. In this case it is actually fine because I understand the first solution and I want to use that.
I want to solve the following equation on a 1D cylindrical grid (sorry I don't have 10 reputation yet so I cannot post the nice rendered equations):
with boundary conditions:
where rho_core is the left side of the mesh, and rho_edge is the right side of the mesh. rho is the normalized radius, J is the Jacobian:
R is the real radius in meters, so the dimension of the Jacobian is distance. The initial conditions doesn't really matter, but in my code example I will use a numerical Dirac-delta at R=0.8.
I have a working example without(!) the Jacobian, but it's quite long, and it doesn't use FiPy's Viewers so I'll link a gist: https://gist.github.com/leferi99/142b90bb686cdf5116ef5aee425a4736
The main part in question is the following:
import fipy as fp ## finite volume PDE solver
from fipy.tools import numerix ## requirement for FiPy, in practice same as numpy
import copy ## we need the deepcopy() function because some FiPy objects are mutable
import numpy as np
import math
## numeric implementation of Dirac delta function
def delta_func(x, epsilon, coeff):
return ((x < epsilon) & (x > -epsilon)) * \
(coeff * (1 + numerix.cos(numerix.pi * x / epsilon)) / (2 * epsilon))
rho_from = 0.7 ## normalized inner radius
rho_to = 1. ## normalized outer radius
nr = 1000 ## number of mesh cells
dr = (rho_to - rho_from) / nr ## normalized distance between the centers of the mesh cells
duration = 0.001 ## length of examined time evolution in seconds
nt = 1000 ## number of timesteps
dt = duration / nt ## length of one timestep
## 3D array for storing the density with the correspondant normalized radius values
## the density values corresponding to the n-th timestep will be in the n-th line
solution = np.zeros((nt,nr,2))
## loading the normalized radial coordinates into the array
for j in range(nr):
solution[:,j,0] = (j * dr) + (dr / 2) + rho_from
mesh = fp.CylindricalGrid1D(dx=dr, nx=nr) ## 1D mesh based on the normalized radial coordinates
mesh = mesh + (0.7,) ## translation of the mesh to rho=0.7
n = fp.CellVariable(mesh=mesh) ## fipy.CellVariable for the density solution in each timestep
diracLoc = 0.8 ## location of the middle of the Dirac delta
diracCoeff = 1. ## Dirac delta coefficient ("height")
diracPercentage = 2 ## width of Dirac delta (full width from 0 to 0) in percentage of full examined radius
diracWidth = int((nr / 100) * diracPercentage)
## diffusion coefficient
diffCoeff = fp.CellVariable(mesh=mesh, value=100.)
## convection coefficient - must be a vector
convCoeff = fp.CellVariable(mesh=mesh, value=(1000.,))
## applying initial condition - uniform density distribution
n.setValue(1)
## boundary conditions
gradLeft = (0.,) ## density gradient (at the "left side of the radius") - must be a vector
valueRight = 0. ## density value (at the "right end of the radius")
n.faceGrad.constrain(gradLeft, where=mesh.facesLeft) ## applying Neumann boundary condition
n.constrain(valueRight, mesh.facesRight) ## applying Dirichlet boundary condition
convCoeff.setValue(0, where=mesh.x<(R_from + dr)) ## convection coefficient 0 at the inner edge
## the PDE
eq = (fp.TransientTerm() == fp.DiffusionTerm(coeff=diffCoeff)
- fp.ConvectionTerm(coeff=convCoeff))
## Solving the PDE and storing the data
for i in range(nt):
eq.solve(var=n, dt=dt)
solution[i,0:nr,1]=copy.deepcopy(n.value)
My code can solve the following equation with the same boundary conditions as indicated above:
To keep it simple I use spatially independent coefficients with the only exeption on the inner edge, where the convection coefficient is 0, and the diffusion coefficient is almost 0. In the linked code I am using a uniform distribution initial condition.
My first question is why do I get the exact same results when using fipy.Grid1D and fipy.CylindricalGrid1D? I should get different results, right? How should I rewrite my code for it to be able to differentiate between the simple 1D Grid and the 1D Cylindrical Grid?
My actual problem is not with this exact code, I just wanted to simplify my problem, but as indicated in the comments this code doesn't produce the same results with the different Grids. So I will just post a GitHub link to a Jupyter Notebook, which may stop working in the future.
The Jupyter Notebook If you want to run it, the first code cell should be run first and after that only the very last cell is relevant. Ignore the reference images. The line plots show the diffusion and convection coefficients. When I ran the last cell with Grid1D or CylindricalGrid1D I got the same results (I compared the plots very precisely)
Sorry but I just cannot rename all my variables, so I hope that based on my comment, and the changed code above (I changed the comments in the code too) you can understand what I'm trying to do.
My other question is regarding the Jacobian. How can I implement it? I've looked at the only example in the documentation which uses a Jacobian, but that Jacobian is a matrix and also it uses the scipy.optimize.fsolve() function.
