4
$\begingroup$

I am attempting to solve the nonlinear Richards equation for the dependent variable $\Psi(t,x,y)$. The equation is as follows:

$\frac{\partial\theta(\Psi)}{\partial t}=\nabla. (k(\Psi)\nabla(\Psi + y) )$

where $\theta(\Psi)$ is a function representing the soil water content and $k(\Psi)$ is the hydraulic conductivity of the medium. Both are defined by the equations:

enter image description here

I am using Mathematica 14 for this task, and the code snippet below is what I have written to tackle the equation:

subst = {n -> 2.9, \[Theta]s -> 0.42, \[Theta]r -> 0.026, \[Alpha] -> 
    0.95, ks -> 0.12};
\[Theta][\[CapitalPsi]_] := 
  If[\[CapitalPsi] <= 
     0, \[Theta]r + (\[Theta]s - \[Theta]r) (1/(1 + (-\[Alpha]  \
\[CapitalPsi])^n))^((n - 1)/n), \[Theta]s] /. subst;
k[\[CapitalPsi]_] := 
  If[\[CapitalPsi] <= 0, 
    ks  \[Theta][\[CapitalPsi]]^(1/
        2) (1 - ((1 - \[Theta][\[CapitalPsi]]^(n - 1)/n))^(n - 1)/
         n)^2, ks] /. subst;

LogLinearPlot[\[Theta][-\[CapitalPsi]], {\[CapitalPsi], 0.1, 30000}, 
 PlotRange -> All]
LogLinearPlot[k[-\[CapitalPsi]], {\[CapitalPsi], 0, 30}, 
 PlotRange -> All]

<< NDSolve`FEM`

Clear[x, t]

m = ToElementMesh[FullRegion[2], {{0, 10}, {0, 10}}, 
  MaxCellMeasure -> 0.1, "MeshOrder" -> 2]

m["Wireframe"]

op = (D[\[Theta][u[t, x, y]], t] - 
     Div[k[u[t, x, y]]  Grad[u[t, x, y], {x, y}], {x, 
       y}]) - (NeumannValue[-10., y == 10 && x <= 2] + 
     NeumannValue[10., x == 10 && y <= 2]);

ic = u[0, x, y] == 0

sol = NDSolveValue[{op == 0, ic}, u, {x, y} \[Element] m, {t, 0, 3}];

However, Mathematica returns an error message:

CoefficientArrays::ivar

Reference paper:

https://link.springer.com/article/10.1007/s10596-016-9566-3

EDIT

@user21 and @Alex Trounev thank you both for your insights and guidance. Following your suggestions, I have solved the problem with Neumann boundary conditions:

subst = {n -> 2.06, \[Theta]s -> 0.396, \[Theta]r -> 0.131, \[Alpha] ->
     0.423, ks -> 4.96 10^-2};
\[Theta][\[CapitalPsi]_] = 
  If[\[CapitalPsi] <= 
     0, \[Theta]r + (\[Theta]s - \[Theta]r)  (1/(1 + (-\[Alpha] \
\[CapitalPsi])^n))^((n - 1)/n), \[Theta]s] /. subst;
k[\[CapitalPsi]_] = 
  If[\[CapitalPsi] <= 0, 
    ks \[Theta][\[CapitalPsi]]^(1/
        2) (1 - ((1 - \[Theta][\[CapitalPsi]]^(n - 1)/n))^(n - 1)/
         n)^2, ks] /. subst;
<< NDSolve`FEM`
m = ToElementMesh[FullRegion[2], {{0, 2}, {0, 3}}, 
   MaxCellMeasure -> 0.25, "MeshOrder" -> 2];

coeff = (D[\[Theta][\[CapitalPsi][t, x, y]], t] /. 
    Derivative[1, 0, 0][\[CapitalPsi]][t, x, y] -> 1);
op = coeff*D[\[CapitalPsi][t, x, y], t] + 
   Inactive[
     Div][-k[\[CapitalPsi][t, x, 
        y]] . (Inactive[Grad][(\[CapitalPsi][t, x, y]), {x, y}] + 
        Inactive[Grad][y, {x, y}]), {x, y}];

T = 3/16.;
steps = 100;
\[CapitalDelta]t = T/steps;
\[CapitalDelta]td = T/3;

(*Solving with a flux boundary condition and a continuous negative \
pressure in all the domain*)
sol = Monitor[
  NDSolveValue[{op == 
     NeumannValue[-1, y == 3 && 0 <= x <= 1], \[CapitalPsi][0, x, 
      y] == -1}, \[CapitalPsi], {x, y} \[Element] m, {t, 0, 
    T, \[CapitalDelta]t}, Method -> Automatic, 
   EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])], monitor]

vecplot2[t_] := Block[{vfield, vecplot, streamplot, cplot},
  vfield = {Evaluate[D[(sol[t, x, y]), x]], 
    Evaluate[D[(sol[t, x, y]), y]]};
  streamplot = 
   StreamPlot[{vfield}, {x, y} \[Element] m, StreamPoints -> Fine];
  (*vecplot=VectorPlot[{vfield},{x,y}\[Element]m];*);
  streamplot
  ]
n = 10
delta = T/n
tab = Table[vecplot2[i  T], {i, 0.001, T, delta}]

enter image description here

However, I'm facing a difficulties with Dirichlet boundary conditions, especially when trying to implement the initial condition for Psi[0, x, y] = 1 - y, as outlined in the paper.

How can I apply the initial condition Psi[0, x, y] = 1 - y ?

Here is the code I can't make it work (The same problem of example 2 of the paper):

ic1 = \[CapitalPsi][0, x, y] == 1 - y
     ic2 = 
  DirichletCondition[\[CapitalPsi][t, x, y] == -2 + 
       2.2 t/\[CapitalDelta]td, 
     y == 3 && x <= 1 && t <= \[CapitalDelta]td]
    ic3 = 
   DirichletCondition[\[CapitalPsi][t, x, y] == 0.2, 
      y == 3 && x <= 1 && t > \[CapitalDelta]td]
     ic4 = DirichletCondition[\[CapitalPsi][0, x, y] == 1 - y, True]
      sol = 
     Monitor[NDSolveValue[{op == 0, ic1, ic2, ic3, 
        ic4}, \[CapitalPsi], {x, y} \[Element] m, {t, 0, 
        T, \[CapitalDelta]t}, Method -> Automatic, 
       EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])], 
      monitor]

enter image description here

enter image description here

Improve this question
$\endgroup$
3
  • 2
    $\begingroup$ I have added some code that this will print a better message in a future version (>14.0). Sorry about that. $\endgroup$
    – user21
    Commented Apr 5 at 7:29
  • $\begingroup$ user21 and AlexTrounev thank you for the advice. I've incorporated some adjustments into the code, ensuring the parameters match those described in the example 2 of the paper. This has allowed me to effectively address the issue related to Neumann boundary conditions. However, I'm currently facing challenges in solving the problem when applying Dirichlet boundary conditions. Any guidance on this would be greatly appreciated. $\endgroup$
    – Stratus
    Commented Apr 6 at 16:19
  • $\begingroup$ ic4 does not make any sense. It's an initial condition in a boundary condition. $\endgroup$
    – user21
    Commented Apr 9 at 15:08

2 Answers 2

Reset to default
3
$\begingroup$

The parsing error comes from the time derivative term D[θ[u[t, x, y]], t]. The problem is that the Derivative[1,0,0][u][t,x,y] is hidden inside the If statement and the parser can not get access to it. I'll try to add some code such that the parser fails more gracefully in a future version.

Now, I have made a few other changes. I have replaced the SetDelayed with Set, to get evaluation of the values in the If statements. (See here)

subst = {n -> 2.9, θs -> 0.42, θr -> 0.026, α -> 0.95, ks -> 0.12};
θ[Ψ_] = If[Ψ <= 0, θr + (θs - θr)  (1/(1 + (-α Ψ)^n))^((n - 1)/n), θs] /. subst;
k[Ψ_] = If[Ψ <= 0, ks θ[Ψ]^(1/2) (1 - ((1 - θ[Ψ]^(n - 1)/n))^(n - 1)/n)^2, ks] /. subst;

Made the mesh smaller - less elements, second order elements are the default:

<< NDSolve`FEM`
m = ToElementMesh[FullRegion[2], {{0, 10}, {0, 10}}, MaxCellMeasure -> 1];

Now, we rewrite the time derivative. We just replace Derivative[1, 0, 0][u][t, x, y] with 1 and store the result in the variable coeff

coeff = (D[\[Theta][u[t, x, y]], t] /. 
   Derivative[1, 0, 0][u][t, x, y] -> 1)

The inactive PDE operator is then (with small term rewriting):

op = coeff*D[u[t, x, y], t] + 
   Inactive[
     Div][-k[u[t, x, y]] . Inactive[Grad][u[t, x, y], {x, y}], {x,y}];

Another issue is the initial condition. I suspect you need a slightly negative value to get this started, but this something that would need to come from the physics of the problem.

ic = u[0, x, y] == -1;

sol = Monitor[
  NDSolveValue[{op == 
     NeumannValue[-10., y == 10 && x <= 2] + 
      NeumannValue[10., x == 10 && y <= 2], ic}, 
   u, {x, y} ∈ m, {t, 0, 0.02}, 
   EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])], monitor]

When we time integrate this, we get a warning that the system is stiff:

NDSolveValue::ndsz: At t == 0.008175579484439153`, step size is effectively zero; singularity or stiff system suspected.

To me this indicates that something with the system setup is not quite right, also seen that Alex had to change the coefficients that might be a likely cause. Please also check out the NDSolveValue::ndsz message ref page. Click the blue circle-i at the end of the message.

Here is an plot of how the solution blows up:

Plot3D[sol[0.005, x, y], {x, y} \[Element] m, PlotRange -> All]

enter image description here

Improve this answer
$\endgroup$
3
$\begingroup$

We can solve this problem using implicit difference scheme in time and linear FEM in space described here and here. Practically this is same method as in the paper. We use some simplifications in the functions $\theta, K$ definitions in the form

subst = {n -> 2.9, \[Theta]s -> 0.42, \[Theta]r -> 0.026, \[Alpha] -> 
   0.95, ks -> 0.12}; appro = 
 With[{k = 2.  10^4}, ArcTan[k  #]/Pi + 1/2 &];
\[Theta][\[CapitalPsi]_] := 
  If[\[CapitalPsi] <= 
     0, \[Theta]r + (\[Theta]s - \[Theta]r)  (1/(1 + (real[\[Alpha]   \
\[CapitalPsi]])^n))^((n - 1)/n), \[Theta]s] /. UnitStep -> appro;
k[\[CapitalPsi]_] := 
  Simplify`PWToUnitStep@
    PiecewiseExpand@
     If[\[CapitalPsi] <= 0, 
      ks   \[Theta][\[CapitalPsi]]^(1/
          2)  (1 - ((1 - \[Theta][\[CapitalPsi]]^(n - 1)/n))^(n - 1)/
           n)^2, ks] /. UnitStep -> appro;
real[x_] := Sqrt[x^2];

f[\[CapitalPsi]_] := 
  If[\[CapitalPsi] <= 
    0, (-1 + 
      n)  \[Alpha]  (-\[Theta]r + \[Theta]s)  (real[\[Alpha]  \
\[CapitalPsi]])^(-1 + n)  (1/(
     1 + (real[\[Alpha]  \[CapitalPsi]])^n))^(1 + (-1 + n)/n), 0];

Then we define mesh and inactive form of equation

<< NDSolve`FEM`

m = ToElementMesh[FullRegion[2], {{0, 10}, {0, 10}}, 
  MaxCellMeasure -> 0.1]


op = (f[U[i - 1][x, y]] (u[x, y] - U[i - 1][x, y])/dt - 
    Inactive[Div][
     k[U[i - 1][x, y]]   Inactive[Grad][u[x, y], {x, y}], {x, y}]);

Finally we solve this problem using NDSolve

U[0][x_, y_] := 0; dt = 1/20; 
Do[U[i] = 
    Quiet@NDSolveValue[
      op == (NeumannValue[-10., y == 10 && x <= 2] + 
          NeumannValue[10., x == 10 && y <= 2]) /. subst, 
      u, {x, y} \[Element] m];, {i, 1, 20}] // AbsoluteTiming

It takes about 53 s on my laptop. Visualization

{Table[DensityPlot[Re[U[i][x, y]], {x, y} \[Element] m, 
   ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
   PlotLabel -> Row[{"t = ", N[i  dt]}], FrameLabel -> Automatic], {i,
    1, 20, 3}], 
 ListLinePlot[Table[{i  dt, U[i][5, 5]}, {i, 0, 20}], 
  AxesLabel -> {"t", "u"}]}

Figure 1

From these pictures it follows that solution very quickly comes to stationary solution. Therefore, we don't need large time and can set {t,0,0.02} as follows

U[0][x_, y_] := 0; dt = 1/500; 
Do[U[i] = 
    Quiet@NDSolveValue[
      op == (NeumannValue[-10., y == 10 && x <= 2] + 
          NeumannValue[10., x == 10 && y <= 2]) /. subst, 
      u, {x, y} \[Element] m];, {i, 1, 10}] // AbsoluteTiming

Visualization Figure 2

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.