HI experts! I'm trying to use a simple time step forward loop for computing a PDE for geophysical fluids. The problem is that [Stigma]np1 = [Stigma]n + dt*equaz is not recognized as an interpolating function. It should be fixed as a given function in NDSolve, but it isn't.

Thanks a lot, Renzo

NON_LINEAR Wind-Driven
(Stommel&Munk)


\[CapitalOmega] = Rectangle[{0, 0}, {a, b}] /. {a -> 1, b -> 1};
RegionPlot[\[CapitalOmega], AspectRatio -> Automatic]

Needs["NDSolve`FEM`"];
mesh = ToElementMesh[\[CapitalOmega], "MaxCellMeasure" -> 1/1000, 
  "MeshElementType" -> TriangleElement]
mesh["Wireframe"]

Tau = Pi;
dt = 0.1;
T = 2;
\[Epsilon] = 0.0*3.5*10^-4;
a = 0.1^2;
s = 0.025;
jacob[x, y] = jacobi[\[Psi], \[Stigma]]
pesce = -\[Stigma][x, y] + \!\(
\*SubscriptBox[\(\[Del]\), \({x, y}\)] . \(
\*SubscriptBox[\(\[Del]\), \({x, y}\)]\[Psi][x, y]\)\)
psiz = {\[Psi][x, y], \[Stigma][x, y]};

QG = \[Alpha]*jacobi[\[Psi], \[Stigma]] + D[\[Psi][x, y], x] + 
  Tau*Sin[Pi*y] - \[Epsilon]*\!\(
\*SubscriptBox[\(\[Del]\), \({x, y}\)] . \(
\*SubscriptBox[\(\[Del]\), \({x, y}\)]\[Stigma][x, y]\)\) + s*\[Stigma][x, y]
NSOperator = {QG, pesce} // Flatten
bcs = {DirichletCondition[\[Psi][x, y] == 0., True], 
   DirichletCondition[\[Stigma][x, y] == 0., True]};
bcpsi = {DirichletCondition[\[Psi][x, y] == 0, True]};
bczeta = {DirichletCondition[\[Stigma][x, y] == 0, True]};





StokesOperator = NSOperator /. \[Alpha] -> 0.;
pde = StokesOperator == {0, 0};
{\[Psi]0, \[Stigma]0} = 
  NDSolveValue[{pde, bcs}, {\[Psi], \[Stigma]}, {x, y} \[Element] mesh, 
   Method -> {"FiniteElement", 
     "InterpolationOrder" -> {\[Psi] -> 2, \[Stigma] -> 2}, 
     "IntegrationOrder" -> 5}];
stokes = ContourPlot[\[Psi]0[x, y], {x, 0, 1}, {y, 0, 1}, 
  PlotLabel -> Style["Linear flow ", Black, FontFamily -> "Times", Bold]]
\[Stigma]n = \[Stigma]0;
\[Psi]n = \[Psi]0;
t = 0.;
For[j = 1, j < T, j++,
 t = t + dt;
 equaz = QG /. {\[Psi] -> \[Psi]n, \[Stigma][x, y] -> \[Stigma]n[x, y], \!\(
\*SubscriptBox[\(\[Del]\), \({x, y}\)] . \(
\*SubscriptBox[\(\[Del]\), \({x, y}\)]\[Stigma][x, y]\)\) -> \!\(
\*SubscriptBox[\(\[Del]\), \({x, y}\)] . \(
\*SubscriptBox[\(\[Del]\), \({x, y}\)]\[Stigma]n[x, y]\)\), 
    jacobi[\[Psi], \[Stigma]] -> jacobi[\[Psi]n, \[Stigma]n], 
    D[\[Psi][x, y], x] -> D[\[Psi]n[x, y], x], \[Alpha] -> a};
 \[Stigma]np1 = \[Stigma]n + dt*equaz
     pde = -\[Stigma]np1[x, y] + \!\(
\*SubscriptBox[\(\[Del]\), \({x, y}\)] . \(
\*SubscriptBox[\(\[Del]\), \({x, y}\)]\[Psi][x, y]\)\) == 0.;
 \[Psi]np1 = 
  NDSolveValue[{pde, bcpsi}, \[Psi], {x, y} \[Element] mesh, 
   Method -> {"FiniteElement", "InterpolationOrder" -> {\[Psi] -> 2}, 
     "IntegrationOrder" -> 5}];
 \[Stigma]n = \[Stigma]np1;
 \[Psi]n = \[Psi]np1;]

ContourPlot[\[Psi]n[x, y], {x, 0, 1}, {y, 0, 1}, 
 PlotLabel -> "NL-flow:stream_function "]
ContourPlot[\[Stigma]np1[x, y], {x, 0, 1}, {y, 0, 1}, 
 PlotLabel -> "NL-flow:Vorticity "]

Evaluate[\[Stigma]n + dt*equaz]












\.10