Physicist here, new to SE and fairly new to MMA. I'm looking to solve a Brownian motion problem with inertia (m ≠ 0) and a velocity dependent diffusion coefficient, which in dimensionless units and with no drag, boils down to a 2+1 dimensional PDE, the Fokker-Planck equation for the probability distribution at position x and velocity v: $\partial_t P= -v \partial_x + \partial_{vv}[D(v)P]$ . The initial condition is $P(x,v,0)=\delta(x)\delta(v-v_0)$ and the boundary condition is that the probability current vanish at infinity.
For constant D=1 the exact solution is a simple Gaussian in x and v (with coefficients that are powers of t). My goal is to identify, for certain polynomial D(v) ≠ 1, the long-distance and long-time functional form of P(x,v,t), e.g., power-law vs. exponential vs. Gaussian. (I.e., less interested in accuracy and more interested in scaling behavior.) The best change-of-variables approach to achieve that might be the subject of another post, or private discussion if someone wishes to engage. But even for just the warm-up problem with D=1 I keep running into stiff system crashes in the limit of large x and t, regardless of choice of variables or minor tweaks to NDSolve parameters (precision goal, integration order, etc.).
One obvious trigger of a stiffness-crash is when the mesh spacing (mcm below) exceeds than the width of the initial value peak used to approximate the initial-value delta function (sig below). But a small mesh spacing, which is fixed as x_max and v_max grow, creates impractically long computation times.
A custom mesh with higher density near the origin (for $v_0$=0) and lower density elsewhere might help. Here's the catch: NDSolve seems to do OK solving the PDE above on a strictly rectangular grid, but poorly on a triangular mesh. And the built in approach to mesh refinement only produces triangle element grids, e.g.:
mrfCirc := Function[{v, a}, If[(Mean[v][[1]])^2 + (Mean[v][[2]])^2 <= r^2, a>aInt, a>aExt]]
mrfSq := Function[{v, a}, If[Abs[Mean[v][[1]]] < r && Abs[Mean[v][[2]]] < r, a>aInt, a>aExt]]
r=4; aInt=.2; aExt=4;
{mesh = ToElementMesh[Rectangle[{-10, -10}, {10, 10}], MaxCellMeasure -> 4,
MeshRefinementFunction -> #], mesh["Wireframe"]} & /@ {mrfCirc, mrfSq}
My question is what's the easiest way, if any, to create a strictly rectangular mesh with higher density of rectangles centered around some point (x0,v0) and lower density further away? Analogous to above, but with rectangles rather than triangles. (I tried ToQuadMesh from the FEMAddOns plugin but it generates an irregular quad element mesh, no better for this purpose than the triangle mesh).
I have other questions but this seems like a good place to start and see what comes up.
Here's the NDsolve code for log(P) (of the many change of variables I've tried, eg Sinh transform, Tanh transform, I've had the most success simply solving for log(P) and retaining x,v).
fpe := D[f[x,v,t],t] + v D[f[x,v,t],x] - (D[f[x,v,t],v])^2 - D[f[x,v,t],v,v] == 0 ; (* f=log[P] *)
ic := f[x,v,0] == -2 Log[sig] - Log[2 \[Pi]] - (x^2 + (v-v0)^2)/(2 sig^2) ; (* log[normal dist], width=sig *)
ndsOpt := {AccuracyGoal->Infinity, PrecisionGoal->pg, Method->{"TimeIntegration"->{"IDA"},
"PDEDiscretization"->{"MethodOfLines", "SpatialDiscretization"->{"FiniteElement",
"MeshOptions"->{"MaxCellMeasure"->mcm} }}}} ;
sig = 1; xmax = vmax = 15; mcm = 0.6; pg = 8; v0 = 0; tmax = 10; (* example *)
{Timing[fsol = NDSolveValue[{fpe, ic}, f, {x,-xmax,xmax}, {v,-vmax,vmax}, {t,0,tmax},
ndsOpt // Evaluate]][[1]], fsol["ElementMesh"]}
tlist = Flatten[{0.5, Range[10]}];
Plot[Table[fsol[x, 0, t], {t, tlist}] // Evaluate, {x, -xmax, xmax}, PlotRange->#, PlotLegends->tlist,
GridLines->Automatic, Frame->True, FrameTicks->{All, All}] & /@ {All, {0, -10}}
