3

I want to step-monitor at each time step of a spatio-temporal NDSolve (using the finite element method) the maximum bending of the function sampled over the curent mesh. For this I first need access to all u[xi,t] over the mesh {xi}. (Note here that I let NDSolve define the mesh)

reg = Line[{{0}, {1}}];
shape = D[0.125 Erf[(x - 0.5)/0.125], x];
op = D[u[t, x], {x, 2}] - D[u[t, x], {t, 2}];
ics = {u[0, x] == shape, Derivative[1, 0][u][0, x] == 0};
sol = NDSolveValue[{op == 0, ics}, u, {t, 0, 2}, {x} ∈ reg, 
      StepMonitor :> Print["xis,uis=" <> ToString["???"]]]

1) Is it possible?

2) How to access (from inside NDSolve) the mesh that has been created?

3) How to access the ui's function at each timestep (replacement of "???")?

Thank you for your help.

Denis

user21
  • 39,710
  • 8
  • 110
  • 167
user3650925
  • 333
  • 1
  • 7
  • Try sol["ElementMesh"]. For 3) I suggest that you take a look at NDSolveValue documentation for option StepMonitor. – Pinti Feb 18 '19 at 16:56
  • Thanks Pinti, but question 2 is "How to access (from inside NDSolve) the mesh that has been created? And for question 3, it is not in the doc. – user3650925 Feb 18 '19 at 17:15
  • 1
    Yes, I apologize, I misread 2). Then you can create mesh upfront with ToElementMesh[Line[{{0}, {1}}]] and pass that to NDSolveValue. For 3) you can try something like StepMonitor :> Print["t=", t, " u=", u[t, x]]. – Pinti Feb 18 '19 at 17:38
  • I have already tried with u[t, 0] because x is nothing, whereas t is the current time at each timestep. But I don't know why, even u is not known from inside NDSolve for this spatio-temporal case... – user3650925 Feb 18 '19 at 17:48

1 Answers1

5

I'm not sure how to get the mesh out of NDSolve but you can use the subprocesses described in Components and Data Structures to reimplement NDSolve with a hook to the element mesh. Note that below is the idea of an implementation that works on the OP's example. I'm not sure what will happen if it is tried on a problem in which there is no FEM ElementMesh created. 

ndsolve`mesh;  (* hook for access to ElementMesh *)
ndsolve[stuff___] := Module[{state, s},
   state = NDSolve`ProcessEquations[stuff];
   ndsolve`mesh = First[state]["FiniteElementData"]["FEMMethodData"]["ElementMesh"];
   Block[{s = #},
      NDSolve`Iterate[s, NDSolve`SolutionDataComponent[s@"VariableRanges", "T"]];
      NDSolve`ProcessSolutions[s]
      ] & /@ state
   ];

foo = {};
ndsolve[{op == 0, ics}, u, {t, 0, 2}, {x} ∈ reg,
 StepMonitor :> (foo = {foo, 
     Plot[Evaluate@ElementMeshInterpolation[{ndsolve`mesh}, u[t, x]][x2],
      {x2} ∈ reg, PlotRange -> {-0.02, 1.3}, 
      PlotLabel -> Row[{HoldForm[t], " = ", t}]]})]

ListAnimate@Flatten[foo]

enter image description here

The use of linked lists to accumulate results (in foo above) and its efficiency is discussed here and here.

Michael E2
  • 235,386
  • 17
  • 334
  • 747