First some comments. This same calculation works beautifully in 2D with essentially the same code. I have copied it largely unaltered into 3D as part of a migration to 3D so that I can do more complex problems there next. This particular case has an analytic solution (which I calculate by hand).
Edit Below is some code which captures the essence of what I actually want to do, and with most of the complexity I want to achieve. This runs fine, but lacks precision.
Clear["Global`*"];
(*---------------------------------------------------------------------------\
*)
FindPsi[mesh_] := Module[{\[Psi], zMnSlv = 0, zmxOfst = 0},
Print[Now, " Finding \[Psi]"];
\[Psi] =
NDSolveValue[{Inactive[Laplacian][\[Psi][x, y, z], {x, y, z}] -
Load[x, y, z] == 0,
DirichletCondition[\[Psi][x, y, z] == 0, dz == 0],
DirichletCondition[\[Psi][x, y, z] == BcPsi[x, y, z],
z == rng]}, \[Psi], {x, y, z} \[Element] mesh,
InitialSeeding -> {\[Psi][x, y, dz] == BcPsi[x, y, z]},
Method -> {"FiniteElement"}, PrecisionGoal -> 6];
Print[Now, " Finished finding \[Psi]"];
Return[\[Psi]]
]
(*---------------------------------------------------------------------------\
*)
Uu[fn_] := Module[{uu},
Print["Finding Uu"];
uu = {D[fn[x, y, z], z], -D[fn[x, y, z], x]};
Return[uu];
]
(*---------------------------------------------------------------------------\
*)
UuTh[x_?NumericQ, y_?NumericQ, z_?NumericQ] := Module[{d, u, dz, m},
dz = z - h[x, y];
u = UU0 Log[(dz + z0)/z0];
Return[u]
]
(*---------------------------------------------------------------------------\
*)
Needs["NDSolve`FEM`"];
z0 = 0.0001;
UU0 = 0.10857;
rng = 10.;
NotebookSave[];
h[x_, y_] := 1/(1 + (x - 0.2 y)^2);
dz = z - h[x, y];
BcPsi[x_, y_, z_] := UU0 ((z + z0) Log[(z + z0)/z0] - z);
Load[x_, y_, z_] := UU0/(dz + z0);
solnRegn =
ImplicitRegion[
z >= h[x, y], {{x, -rng, rng}, {y, -rng, rng}, {z, 0, rng}}];
mesh = ToElementMesh[solnRegn, "MaxCellMeasure" -> 0.20,
"MaxBoundaryCellMeasure" -> 0.1]
Print["MeshOrder = ", mesh["MeshOrder"]];
Print[Magnify[mesh["Wireframe"], 1.8]];
psi = FindPsi[mesh];
psi[-3, -2, 2]
uu[x_, y_, z_] = Uu[psi];
uu[-3, -2, 2]
x = -9.;
y = 0;
LogLinearPlot[{uu[x, y, z][[1]], UuTh[x, y, z]}, {z, 0.003, 10.},
PlotRange -> {{0.01, 10}, {0., 1.3}}]
And here is my output:
The blue curve is the output from the FEM calculation. The brown/orange curve is a good approximation to the right curve, calculated analytically. I would like better precision in the numerical solution, without it running out of memory, and without it taking hours to run :-). I don't need 6 digits of precision, but a sensible target would be a relative error of <= 3% along the whole curve, or at least down to z=0.02.
I should add that the function h[x,y] could be any single-valued, continuous surface, which may be inscribed in the box, including possibly an interpolating function.
The steps around z=0.10 would seem to be an artefact of the mesh, while the upturn at z~10 is, I think, due to the bevelled edges of the solution region.
Any advice on improving this would be very gratefully received. Thanks.
Edit Here is another output. Same code, but using "MaxBoundaryCellMeasure" -> 0.015. This uses 991000 mesh elements and takes about 20' on my machine. It is clearly better, but not yet near my desired precision. I don't think I can use more mesh elements without running out of memory.
Further Development
I provide below one attempt to use a mesh refinement function. It does not behave as I expect, but presumably, there is a bug in it. The idea here is that for all vertices within the solution volume for which dz<1, the volume is limited to an amount which is proportional to dz itself. The idea of that is to make the elements get smaller and smaller as we reach the boundary at dz=0, where the BCs apply. First the code:
mrf = Function[{vertices, volume},
Block[{x, y, z, dz}, {x, y, z} = Mean[vertices];
dz = z - h[x, y]; If[dz < 1., volume > 0.002 dz, False]]];
mesh = ToElementMesh[solnRegn, MaxCellMeasure -> 0.25,
MeshRefinementFunction -> mrf]
Print[Magnify[
Show[{Plot3D[h[x, y], {x, -10, 10}, {y, -10, 10},
PlotRange -> {All, {-0.3, 0.3}, {0, 10}}],
Graphics3D[{PointSize[0.004], Point[mesh[[1]]]},
AspectRatio -> Automatic, Axes -> True,
PlotRange -> {All, {-0.3, 0.3}, All}]}]
, 1.8]];
The rest of the code is unchanged from that posted further above. I show here the plot of part of the mesh volume, which has a dot for each vertex. The brown surface is defined by dz=0.
The thing that surprises me is that despite the fact that there are many more mesh vertices in the region around 0.3 <~ dz < 1, there is a marked gap for dz <~ 0.4, containing just a single "surface" of points not the dense cloud I expected. Why the gap?
Finally, the plot of the quantity I am hoping to extract:
Which is very similar to before.
After Adopting Provided Answer
Just to close the loop, I provide here my test plot after adopting the solution provided by @Tim Laska and adapting the step size and ratio.
It is clearly very much improved. I believe the small residual disagreement in the unstepped part of the curve is "physics", ie. a small mismatch such as this is anticipated between the theoretical solution and the "true" solution, including the cross-overhear z=4 (z is shown on the "x-axis").
