# NDSolve – Refine Spatial Mesh with TimeIntegration Method

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["NDSolveFEM"]
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.

• Something is going wrong with the brackets, did you mean I1*(wo/w[z])*Exp[(-2*y^2)/w[z]^2] Exp[-Eps*(-z + zo)]? Aug 15, 2016 at 8:10
• Yes, thanks for the correction. Aug 15, 2016 at 8:16
• What code do you actually use to generate the (asymmetric) plot? The NDSolve[] generates three 3D InterpolationFunctions Aug 15, 2016 at 8:48
• Please chek the edited code above the picture of the solution Aug 15, 2016 at 8:56
• Have you seen ToElementMesh and the related FEM tutorials? Aug 15, 2016 at 11:15

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.. >>
`