NDSolve – Refine Spatial Mesh with TimeIntegration Method

Question

I am trying to solve the following system of non-linear PDEs with 1 time + 2 space dimensions:

 kp = 0.504; kt = 1.31; kd = 2.8*10^-7;

 equ1 = D[Ini[t, y, z], t] == -kd*Iabs[y, z]*Ini[t, y, z];
 equ2 = D[Rad[t, y, z], t] == 2*kd*Iabs[y, z]*Ini[t, y, z] - 2*kt*(Rad[t, y, z])^2;
 equ3 = D[DB[t, y, z], t] == -kp*Rad[t, y, z]*DB[t, y, z];
 sys = {equ1, equ2, equ3}; 
 iniconditions = {DB[0, y, z] == 12.89*10^3, Rad[0, y, z] == 0, Ini[0, y, z] == 66.2};

The “perturbation” Iabs[x,y] of the system is symmetric:

Here is the code to generate it:

wo = 5.1*10^-7; zR = 2.49*10^-6; I1 = 9.4*10^6; Eps = 387.6; zo = 200*10^-6;
w[z_] = wo*Sqrt[1 + (z/zR)^2];
Iabs[y_, z_] = I1*(wo/w[z])*Exp[(-2*y^2)/w[z]^2]*Exp[-Eps*(-z + zo)]; 

And I am therefore expecting a symmetric result. NDSolve is able to provide a numeric solution to the system of equation using the "TimeIntegration" Method. However the resulting interpolations are asymmetric:

 sol = NDSolve[Join[sys, iniconditions], {DB, Rad, Ini}, {t, 0, 1}, {y, -100*wo, 100*wo}, {z, -100*10^-6, 200*10^-6}]
 DB[t_, y_, z_] = DB[t, y, z] /. sol;
 DensityPlot[(1 - DB[0.1, y, z]/(12.89*10^3))*100, {y, -100*wo,100*wo}, {z, -100*10^-6, 100*10^-6}, AspectRatio -> Automatic, PlotLegends -> Automatic, PlotRange -> All,ColorFunction -> Hue, PlotPoints -> 100];

Following this reference I tried to increase the number of plotted points with PlotPoints, it only smoothed out the result. Therefore I am also trying to refine the spatial mesh over the yz domain.

I have not found an option to do it with the TimeIntegration Method and using a rectangle domain to solve the system:

Omega = Rectangle[{-100*wo, -100*10^-6}, {100*wo, 200*10^-6}];
sol = NDSolve[Join[sys, iniconditions], {DB, Rad, Ini}, {t, 0, 1}, Element[{y,z},Omega], Method -> {"TimeIntegration" -> {"Adams"}}]

Yields the following error messages

NDSolve::femnonlinear: Nonlinear coefficients are not supported in this version of NDSolve

My question is then:

• Is there a way to controllably semi-discretize the problem with a more refined mesh of the yz spatial domain while maintaining the TimeIntegration method?

---

Edit: Using ToElementMesh I manually defined the spatial mesh over the yz-domain and passed it to NDSolve:

Needs["NDSolve`FEM`"]
Omega = ImplicitRegion[True, {{y, -100*wo, 100*wo}, {z,-100*10^-6, 200*10^-6}}];
mesh = ToElementMesh[Omega, MaxCellMeasure -> wo/10 000];
mesh["Wireframe"]

With MaxCellMeasure -> wo/10 000, the mesh is a bit coarse but I also tried with MaxCellMeasure-> wo/100 000 unsuccessfully:

When I pass this mesh to NDSolve, I obtain the same error message as when using ImplicitRegion

sol = NDSolve[Join[sys, iniconditions], {DB, Rad, Ini}, {t, 0, 1}, Element[{y, z}, mesh], Method -> {"TimeIntegration" -> {"Adams"}}];
NDSolve::femnonlinear: Nonlinear coefficients are not supported in this version of NDSolve. >>

So I don't understand why NDSolve is able to solve the system by itself using the TimeIntegration Method but whenever I try to improve the mesh quality it doesn't work.

Answer 1

Thanks to @MichaelE2's comment, I tried the following:

sol = NDSolve[Join[sys, iniconditions], {DB, Rad, Ini}, {t, 0, 1}, {y,-100*wo, 100*wo}, {z, -100*10^-6, 200*10^-6},Method -> {"PDEDiscretization" ->{"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 400}}}]

Which gives a correct symmetric solution to the system of PDEs:

DB[t_, y_, z_] = DB[t, y, z] /. sol;
DensityPlot[(1 - DB[0.3, y, z]/(12.89*10^3))*100, {y, -40*wo,40*wo}, {z,-74*10^-6, 74*10^-6}, AspectRatio -> Automatic, PlotLegends -> Automatic, PlotRange -> All, ColorFunction -> Hue, PlotPoints -> 100]

The trick was here to put "MinPoints" to at least 200 otherwise NDSolve was returning:

Power::infy: Infinite expression 1/0. encountered. >>
Encountered non-numerical value for a derivative at t == 0.`. >>