Solving PDE on 1D cylindrical coordinates with FiPy

Question

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.

Answer 1

[cobbling an answer from the discussion in the comments]

The results are similar between a Grid1D and a CylindricalGrid1D, particularly in the early steps, but they are not the same. They are quite different as the problem evolves.

FiPy doesn't like things outside the divergence, but you should be able to multiply the equation by J and put it in the coefficient of the TransientTerm, e.g.,

or

eq = fp.TransientTerm(J) == fp.DiffusionTerm(coeff=J * diffCoeff) - fp.ConvectionTerm(coef=J * convCoeff)

For the Jacobian, you could create a CellVariable for the real radius in terms of the normalized radius, and then take its gradient:

real_radius = fp.CellVariable(mesh=mesh, value=...)
J = real_radius.grad.dot([[1]])

.grad returns a vector, even in 1D, but the coefficient must be scalar, so take the dot product to get the x component.