I'm trying to solve a system of non-linear reaction-diffusion-advection PDEs numerically using NDSolve.
The two equations are (taken from here):
where $H(x)$ is the Heaviside function, and $P(x,t)$ and $X(x,t)$ are the functions I am solving for. My conditions for flux are $\frac{\partial{P}}{\partial{x}}\bigg|_{x=0} \leq \frac{q}{Av}$ and $\frac{\partial{P}}{\partial{x}}\bigg|_{x=0} \leq 0$. For my code, I just use the equality condition.
The boundary conditions at $-\infty$ and $\infty$ are $X(-\infty) = S$ where $S$ is just $0.01$, $P(-\infty) = 0$, $X(\infty) = S-\frac{K_2 q}{K_1 \alpha}$ and $P(\infty) = \frac{q}{K_1 A}(1+\frac{k}{X(\infty)})$(in the paper I am referring to they are solving for steady states of the system analytically, but I don't see why these conditions can't be true for all values of t).
Now, since I am solving the PDEs between $-\infty<x<\infty$, and I have flux conditions at the origin, I thought I would use NDSolve twice, once between $-\infty<x<0$, and once between $0<x<\infty$, where my Dirichlet boundary conditions for the origin in the second case just become the values that I found in the first case.
The parameters are given by
A = 2100;
v = 4.32 * 10^4;
k1 = 8.27;
k2 = 44.10;
k = 0.007;
alpha = 16.50;
S = 0.01;
dp = 3456*1000;
dx = 3456*1000;
L = 325000;
q = 0.06;
Before I solve it in Mathematica, I non-dimensionalize my equations as follows: take $y=\frac{x}{L}$, where $L$ is the length of the river, and $\tau=K_{1}t$, and divide throughout by $A S$. Thus my equations become:
\begin{align*} K_1\dfrac{\partial{(P)}}{\partial{\tau}} &= D_P L^2\dfrac{\partial^2{(P)}}{\partial{y^2}}-vL\dfrac{\partial{(P)}}{\partial{y}}-K_1\dfrac{X}{X+\frac{k}{S}}P + \frac{f(y,\tau)}{A}\\ K_{1}\dfrac{\partial{(X)}}{\partial{\tau}} &= D_X L^2\dfrac{\partial^2{(X)}}{\partial{y^2}}-vL\dfrac{\partial{(X)}}{\partial{y}}-K_2\dfrac{X}{X+\frac{k}{S}}P +\frac{\alpha}{A}(1-X)\\ \end{align*}
where note that I should really have a $P'(x,t)$ where $P'(x,t)=\dfrac{P(x,t)}{S}$.
Here's my code for finding a solution:
A = 2100;
v = 4.32*10^4;
k1 = 8.27;
k2 = 44.1;
k = 0.007;
alpha = 16.5;
S = 0.01;
dp = 3.456*10^6;
dx = 3.456*10^6;
L = 325000;
q = 0.06;
xmin = 100;
xmax = 100;
Tmax = 50;
PollutantAddition[x_, t_] := (q/A)*HeavisideTheta[x];
eqn11 =
k1*D[P[x, t], t] == dp*L*L*D[P[x, t], x, x] - v*L*D[P[x, t], x] - k1*(Q[x, t]/(Q[x, t] + 0.7))*P[x, t] + PollutantAddition[x, t];
eqn21 =
k1*D[Q[x, t], t] == dx*L*L*D[Q[x, t], x, x] - v*L*D[Q[x, t], x] - k2*(Q[x, t]/(Q[x, t] + 0.7))*P[x, t] + (alpha/A)*(1 - Q[x, t]);
ic2 = Q[x, 0] == 1;
ic1 = P[x, 0] == 0;
{P2, X2} =
{P, Q} /.
NDSolve[
{eqn11, eqn21, ic2, ic1, P[-xmin, t] == 0, Q[0, t] == 1,
Derivative[1, 0][Q][0, t] == 0, Derivative[1, 0][P][0, t] == q/(A*v)},
{P, Q}, {x, -xmin, 0}, {t, 0, Tmax}][[1]];
Plot3D[P2[x, t], {x, -xmin, 0}, {t, 0, Tmax},
PlotRange -> All,
PlotLabel -> "Pollutant concentration for constant boundary conditions",
AxesLabel -> {x, t}]
Plot3D[X2[x, t], {x, -xmin, 0}, {t, 0, Tmax},
PlotRange -> All,
PlotLabel -> "Oxygen concentration for constant boundary conditions",
AxesLabel -> {x, t}]
ContourPlot[P2[x, t], {x, -xmin, 0}, {t, 0, Tmax},
PlotRange -> All, PlotLegends -> Automatic, FrameLabel -> {x, t}]
ContourPlot[X2[x, t], {x, -xmin, 0}, {t, 0, Tmax},
PlotRange -> All, PlotLegends -> Automatic, FrameLabel -> {x, t}, Contours -> 10]
And these are the errors I get:
NDSolve::mxsst: Using maximum number of grid points 10000 allowed by the MaxPoints or MinStepSize options for independent variable x..
NDSolve::ndcf: Repeated convergence test failure at t == 3.9054643319313938`; unable to continue.
NDSolve::eerr: Warning: scaled local spatial error estimate of 1.3210326900129446
*^11 at t = 3.9054643319313938in the direction of independent variable x is much greater than the prescribed error tolerance. Grid spacing with 10001 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options.
I still get a weird looking solution with highly negative values (note that I am solving for concentrations, so they really shouldn't be negative). I try to solve the MaxStepSize issue by putting MaxStepSize -> 0.5, after which this is what the Contour Plot of P(x,t) and X(x,t) look like:
This is really not the behavior I would have expected. P(x,t) is constant with time over the length of the river, and X(x,t) is huge ($~10^{26}$) when I wouldn't expect it to go over $0.02$. Also, the oxygen concentration is increasing with time for small x, which again is counterintuitive since I have a term adding a constant amount of pollution (the term with the Heaviside function) for all $x>0$.
Now I clearly don't have the right values for the first part of my solution, but I also am not sure how I'm going to solve for the second part. What boundary values and initial values should I be using?
{P3, X3} =
{P, Q} /.
NDSolve[
{eqn11 == NeumannValue[(D[P2[x, 0], x] /. x -> 0 ), x == 0],
eqn21 == NeumannValue[(D[X2[x, 0], x] /. x -> 0 ), x == 0],
P[0, t] == P2[0, t], Q[0, t] == X2[0, t],
Q[xmax, t] == S - (k2 q)/(k1 alpha),
P[xmax, t] == q/(k1 A) (1 + k/xinf)},
{P, Q}, {x, 0, xmax}, {t, 0, Tmax}][[1]] ;
These are the errors I'm getting:
NDSolve::femnlmdor: The maximum derivative order of the nonlinear PDE coefficients for the Finite Element Method is larger than 1. It may help to rewrite the PDE in inactive form.
ReplaceAll::reps: {(8.27 (P^(0,1))[x,t]==0.0000285714 HeavisideTheta[x]-(8.27 P[x,t] Q[x,t])/(0.7 +Q[<<2>>])-1.404*10^10 (P^(1,0))[x,t]+3.6504*10^17 (P^(2,0))[x,t])==NeumannValue[0.,x==0],(8.27 (Q^(0,1))[x,t]==0.00785714 (1-Q[<<2>>])-(44.1 P[x,t] Q[x,t])/(0.7 +Q[<<2>>])-1.404*10^10 (Q^(1,0))[x,t]+3.6504*10^17 (Q^(2,0))[x,t])==NeumannValue[0.,x==0],P[0,t]==<<1>>,Q[0,t]==<<1>>,Q[100,t]==-0.00939101,P[100,t]==8.7962*10^-7} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.
Set::shape: Lists {P3,X3} and {P,Q}/. {(8.27 (P^(0,1))[x,t]==0.0000285714 HeavisideTheta[x]-(8.27 P[x,t] Q[x,t])/Plus[<<2>>]-1.404*10^10 (P^(1,0))[x,t]+3.6504*10^17 (P^(2,0))[x,t])==NeumannValue[0.,x==0],(8.27 (Q^(0,1))[x,t]==0.00785714 (1+Times[<<2>>])-(44.1 P[x,t] Q[x,t])/Plus[<<2>>]-1.404*10^10 (Q^(1,0))[x,t]+3.6504*10^17 (Q^(2,0))[x,t])==NeumannValue[0.,x==0],P[0,t]==InterpolatingFunction[{{-100.,0.},{0.,4.10469}},{5,5,1,{10001,161},{6,4},0,0,0,0,Automatic,{},{},False},{{-100.,-99.99,-99.98,-99.97,-99.96,-99.95,-99.94,-99.93,-99.92,-99.91,-99.9,-99.89,-99.88,-99.87,-99.86,-99.85,-99.84,-99.83,-99.82,-99.81,-99.8,-99.79,-99.78,-99.77,-99.76,-99.75,-99.74,-99.73,-99.72,-99.71,-99.7,-99.69,-99.68,-99.67,-99.66,-99.65,-99.64,-99.63,-99.62,-99.61,-99.6,-99.59,-99.58,-99.57,-99.56,-99.55,-99.54,-99.53,-99.52,-99.51,<<9951>>},{0.,6.99322*10^-17,1.39864*10^-16,2.09797*10^-16,<<44>>,0.0226525,0.0270858,<<111>>}},<<1>>,{Automatic,Automatic}][0,t],<<1>>==<<1>>,Q[100,t]==-0.00939101,P[100,t]==8.7962*10^-7} are not the same shape.
About the inactive form bit - I am reading about it in the documentation and will try and figure out how to write down my equations.
Does anyone have ideas on what I'm doing wrong, at least with the first part? Any help is appreciated.
