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
Attachments: