1D advection-diffusion-reaction PDE system

[Marcus Vinicius Pereira]

Hello everyone!

I've been trying to solve a 1D advection-diffusion-reaction PDE system of gas adsorption and can't really know what's going on. I've followed the approach in [1] using Crank-Nicholson for the diffusion step and FD for the reaction equation.

I've even tried to reproduce the solution for this discussion. The riemann solver I've been using is below:

def conservative_advection_1D(q_l, q_r, aux_l, aux_r, problem_data):

    num_rp = q_l.shape[1]

    wave = np.empty((num_eqn, num_waves, num_rp))
    s = np.empty((num_waves, num_rp))
    amdq = np.zeros((num_eqn, num_rp))
    apdq = np.zeros((num_eqn, num_rp))

    s[0, :] = aux_l[0, :].copy()

    for i in range(num_eqn):

        wave[i, 0, :] = q_r[i, :] - q_l[i, :]

        fmh1 = (s[0, :] >= 0) * s[0, :] * q_l[i, :]
        fmh2 = (s[0, :] < 0) * s[0, :] * q_r[i, :]
        fmh = fmh1 + fmh2

        amdq[i, :] = fmh - s[0, :] * q_l[i, :]
        apdq[i, :] = -fmh

    return wave, s, amdq, apdq

Are there any other approach I can try?

[1] Tyson, R., Stern, L. G., & LeVeque, R. J. (2000). Fractional step methods applied to a chemotaxis model. Journal of mathematical biology, 41, 455-475.

[David Ketcheson]

Hi Marcus,

Your question isn't clear to me -- is your code not doing what you wanted/expected?  In what way?  And what equation are you trying to solve (I think in the previous discussion it was noted that there are different "advection" equations -- conservative or non-conservative).

David

[Marcus Vinicius Pereira]

Dear David,

I'm trying to solve the following conservative equation

dC/dt =  -d(uC)/dx +  D*d2C/dx2 - K*dn/dt

By following a similar approach to the reference paper, the second part (D*d2C/dx2) is solved implicitly by Crank-Nicholson and the third part (dn/dt) by forward difference.

[David Ketcheson]

Okay, thanks for clarifying the equation.  I still don't know what kind of behavior you're getting from your code, and I can't see your code except the Riemann solver, so the help I can offer is limited by the information you've provided.

The Riemann solver definitely seems not correct, since it isn't conservative (the fluctuations don't add up to the flux difference).  Why not use an existing Riemann solver?

David

[Marcus Vinicius Pereira]

Dear David,

After a deep review in the code I've discovered that the velocity calculation in the before_step were wrong and the code works fine by using the 1D riemann solver with variable velocity.

thanks a lot!

Marcus