NDSolve to solve 2D space 1D time partial integro-differential equations

I have a set of pdes in 2D space and 1D time, one of them is related to the integral in space by Reynolds function. The equations are similar to $$\frac{\partial \phi}{\partial t}=-\frac{\partial \phi_Z}{\partial x}\phi$$ $$\frac{\partial\phi_Z}{\partial t}=-\int_0^1 \frac{\partial\phi}{\partial x} \frac{\partial \phi}{\partial y} dy$$ $$Reynolds\equiv -\int_0^1 \frac{\partial\phi}{\partial x} \frac{\partial \phi}{\partial y} dy$$ To calculate the Reynolds stress (averaging in y direction), I used the method from Michael E2's answer on question 175080, which breaks the NDSolve[] apart. That method works well in 1D space, so I modified it like below:

db = 0.01; rho = 0.005; tmax = 0.8;
eqns = {D[ph[t, x, y], t] == -D[PZ[t, x, y], x]*D[ph[t, x, y], y],
D[PZ[t, x, y], t] == Reynolds[ph[t, x, y], t, x, y]};
incs = {ph[0, x, y] == E^(-((x - 0.5)^2 + (y - 0.5)^2)/0.01), PZ[0, x, y] == 0};
bcs = {ph[t, 0, y] == ph[t, 1, y], ph[t, x, 0] == ph[t, x, 1],
PZ[t, 0, y] == PZ[t, 1, y], PZ[t, x, 0] == PZ[t, x, 1]};
sys = {eqns, incs, bcs};

Block[{Reynolds},

Reynolds[pp_, t_ /;t == 0, xn_?NumericQ, yn_?NumericQ] := 0;
(*Uncomment below, it works;*)
(*Reynolds[pp_,t_,xn_?NumericQ,yn_?NumericQ]:=0;*)

Reynolds[pp_, t_, xA_?MatrixQ, yA_?MatrixQ] := 0;
(*Some integral like below, I use 0 here, not affect the error*)
(*ParallelTable[
Module[{INTfunc = Interpolation[Transpose@{Transpose@{xA, yA}, pp}]},
-NIntegrate[Derivative[1, 0][INTfunc][xx, yi]*Derivative[0, 1][INTfunc][xx, yi]
, {yi, 0, 1}]], {xx, xA}, {yy, yA}];*)

Clear[foo];
PrintTemporary@Dynamic@{foo, tmax, Clock[Infinity]};
(*broken down NDSolve call*)

InternalInheritedBlock[{MapThread}, {state} =
NDSolveProcessEquations[sys, {ph, PZ}, {x, 0, 1}, {y, 0, 1}, {t, 0, tmax},
InterpolationOrder -> All, StepMonitor :> (foo = t),
Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" ->
{"TensorProductGrid", "MinPoints" -> 10}}}];
MapThread[f_, data_, 1] /; ! FreeQ[f, Reynolds] := f @@ data;
NDSolveIterate[state, {0, tmax}];
sol = NDSolveProcessSolutions[state]]] // AbsoluteTiming


But the MMA has error like

NDSolveIterate::nlnum: The function value {0.,0....0.,<<238>>} is not a list of numbers with dimensions {288} at {t,ph[t,x,y],PZ[t,x,y]} = {8.9003*10^-308,{<<1>>},...}.


I believe the problem comes from the definition of

Reynolds[pp_, t_, xA_?MatrixQ, yA_?MatrixQ] := 0;


This function is not used by the NDSolveIterate[] ! Here I want to interpolate pp[xA,yA], so I can do the integral numerically. But seems the argument passed to Reynolds function is Numeric instead of Matrix/Array or Vector. VectorQ is tried:

Reynolds[pp_, t_, xA_?VectorQ, yA_?VectorQ] := 0;


But only NumericQ works.

Reynolds[pp_, t_, xA_?NumericQ, yA_?NumericQ] := 0;


Look like it is not like the 1D space case, where x is a Vector and passed to function. I can't do interpolation with only one numeric value.

Is there any solution for this? For example, force passing Matrix/Array to Reynolds function here? Or other suggestion? Thank you very much!

• Could you explain your problem since it is not clear from your code what do you try to solve? May 25, 2021 at 22:20
• @AlexTrounev I updated the question with descriptions on equations. My question is mainly about how to calculate the averaging in y in Reynolds stress. Thanks. May 25, 2021 at 23:50
• Could you give a link to a paper or book with model description? May 26, 2021 at 10:10
• @AlexTrounev Thank you for your kind help. The problem is solved now. I checked the state@"NumericalFunction" in NDSolveStateData["<"0.">"] , and I found that the MapThread used in StateData for 2D-space is MapThread[Function[{x, y, ph}, Reynolds[ph, t, x, y]], {x, y, ph}, 2]}]]. Therefore, I only need to change  MapThread[f_, data_, 1] /; ! FreeQ[f, Reynolds] := f @@ data; into MapThread[f_, data_, 2]...`. Then the Matrix will be passed to the Reynolds function I defined. This problem is related to zonal flow in tokamak, the zonal flow is driven by Reynolds' stress. May 26, 2021 at 11:31
• Could you add this piece of code to your question and show one working example? May 26, 2021 at 20:52