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};

  Reynolds[pp_, t_ /;t == 0, xn_?NumericQ, yn_?NumericQ] := 0;
  (*Uncomment below, it works;*)
  Reynolds[pp_, t_, xA_?MatrixQ, yA_?MatrixQ] := 0;
   (*Some integral like below, I use 0 here, not affect the error*)
    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}];*)
  PrintTemporary@Dynamic@{foo, tmax, Clock[Infinity]};
  (*broken down NDSolve call*)
  Internal`InheritedBlock[{MapThread}, {state} = 
    NDSolve`ProcessEquations[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;
   NDSolve`Iterate[state, {0, tmax}];
   sol = NDSolve`ProcessSolutions[state]]] // AbsoluteTiming

But the MMA has error like

NDSolve`Iterate::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 NDSolve`Iterate[] ! 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!

  • $\begingroup$ Could you explain your problem since it is not clear from your code what do you try to solve? $\endgroup$ May 25, 2021 at 22:20
  • $\begingroup$ @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. $\endgroup$
    – Yan QH
    May 25, 2021 at 23:50
  • $\begingroup$ Could you give a link to a paper or book with model description? $\endgroup$ May 26, 2021 at 10:10
  • $\begingroup$ @AlexTrounev Thank you for your kind help. The problem is solved now. I checked the state@"NumericalFunction" in NDSolve`StateData["<"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. $\endgroup$
    – Yan QH
    May 26, 2021 at 11:31
  • $\begingroup$ Could you add this piece of code to your question and show one working example? $\endgroup$ May 26, 2021 at 20:52


Your Answer

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