How to define a PDE with coefficients that depend on the independent variables?

[legofernando1999]

Hi,

I'm trying to solve the PDE:

This is what I have done so far:
`def _2D_C_pypde_Solver(self, dt, M, N, c_old, Da_A, Da_P, Da_D, Ceq, Pe, L, R):
from pde import CylindricalSymGrid, ScalarField, PDEBase
grid = CylindricalSymGrid(radius=1., bounds_z=(0., 1.), shape=(M, N))
state = ScalarField(grid=grid, data=c_old, dtype=np.float64)
bc_z_bottom = {'value': 0.}
bc_z_top = {'value': np.nan} # dummy value for non-existent boundary, #351
bc_r_left = {'derivative': 0.}
bc_r_right = {'type': 'mixed', 'value': -Da_D}
bc_r = [bc_r_left, bc_r_right]
bc_z = [bc_z_bottom, bc_z_top]

    class TransportPDE(PDEBase):

        def __init__(self, bc):
            """ 
            Initialize boundary conditions.
            """
            self.bc = bc

        def evolution_rate(self, state, t):
            dC_dz = state.gradient(self.bc)['z']
            nabla2_rC = state.laplace(self.bc)['r']
            C2 = state**2

            return -2*(1 - r)*dC_dz + L/(R*Pe) * nabla2_rC - Da_A*C2 \
                    + Da_P*(1 - Ceq)*np.exp(-Da_P*z)
        
    eq = TransportPDE([bc_r, bc_z])
    return eq.solve(state, t_range=1, dt=dt)`

But I know this won't work because I haven't defined r nor z in the evolution_rate function, I don't know how to do that. Also, how can I tell the solver that I want to use backward differences in the z direction (because there's only one boundary condition)? Finally, how can I let the solver know that the variables R, Pe and Ceq depend on t and z (these are vectors), while Da_A, Da_P and Da_D depend on t, r, and z (these are matrices)?

Thank you for your time.

[david-zwicker]

You're asking multiple questions, and I try to answer at least a few of them:

• You can access the values of the (discretized) coordinates along each axis using grid.axes_coords or the coordinates of each of the discretized cells using grid.cell_coords. The grid instance is available from state.grid in evolution_rate.
• The backward differences in the z-direction are not implemented, although it is in principle possible. If you absolutely need this, you could implement it yourself by changing the make_laplace function of the cylindrical grid. The respective implementation of the spherical coordinates shows how different implementations of the operators can be done. If you do this, it would be great if you create a pull request, so this can be added to the package.
• I don't understand your last question. Looking at the equation, all the quantities you list there look like scalar fields to me.