Context
I would like to solve a PDE on a boundary which is parametrized as a BSpline. I am trying to solve the force-free
Grad-Shafranov equation
on a boundary whose shape I do not know in advance.
Specifically I need to solve for the toroidal flux of the magnetic field above an accretion disc.
The Grad-Shafranov equation reads (in cylindrical coordinates)
R D[P[R, z], {R, 2}] + R D[P[R, z], {z, 2}] - D[P[R, z], R] == - R/2;
and I am seeking solution satisfying P==0
on a spline, see below.
This question is related to the physical context of that question
,
where we try in to explain astrophysical jets like this:
Eventually I would like to optimize the problem while changing the shape of the spline.
First attempt
I define my region via a BSpline:
ff0 = BSplineFunction[pts = {{1, 0}, {1.2, 2}, {0, 2}}]
So the upper envelope of the jet looks like this:
pl0 = ParametricPlot[ ff0[t] // Release, {t, 0, 1},
Frame -> False, Axes -> False, PlotPoints -> 15, ImageSize -> Small]
and the region like that:
pl = ParametricPlot[r ff0[t] // Release, {t, 0, 1}, {r, 0.01, 1},
Frame -> False, Axes -> False, PlotPoints -> 15, ImageSize -> Small]
I can then discretize both the boundary and the region:
Ω = DiscretizeGraphics[pl]
δΩ = DiscretizeGraphics[pl0, MaxCellMeasure -> 0.1]
and then solve for the PDE
eqn0 = R D[P[R, z], {R, 2}] + R D[P[R, z], {z, 2}] - D[P[R, z], R] == - R/2;
P0 = NDSolveValue[{eqn0,
DirichletCondition[P[R, z] == 0, R == 0],
DirichletCondition[P[R, z] == 0, {R, z} ∈ δΩ],
DirichletCondition[P[R, z] == E R^2 Log[1/R^2], z == 0]},
P, {R, z} ∈ Ω, Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 1/10000},
"IntegrationOrder" -> 3}}]
If I then try and plot the resulting PDE solution, P0
,
ContourPlot[P0[R, z], {R, z} ∈ Ω,
PlotLegends -> Automatic, PlotPoints -> 30,
ColorFunction -> "LightTemperatureMap", ImageSize -> Small,
PlotRange -> All,
FrameLabel -> {R, z},
AspectRatio -> 1]
Even though it seems happy, it satisfies very poorly the boundary on the spine:
Plot[ P0 @@ ff0[t], {t, 0, 1}, ImageSize -> Small]
This should be zero…
Second attempt
Following J. M., I have attempted using explicit splines and ParametricRegion
as follows:
pts = {{1, 0}, {1.8, 3}, {0, 2}};
{xu, yu} = Transpose[pts];
n = 2;m = Length[pts];
knots = {ConstantArray[0, n + 1], Range[m - n - 1]/(m - n),
ConstantArray[1, n + 1]} // Flatten;
fx[t_] = xu.Table[ BSplineBasis[{n, knots}, i - 1, t], {i, Length[pts]}];
fy[t_] = yu.Table[ BSplineBasis[{n, knots}, i - 1, t], {i, Length[pts]}];
Indeed
ParametricPlot[{fx[t], fy[t]}, {t, 0, 1}, Axes -> None, Frame -> True,
Epilog -> {Directive[AbsolutePointSize[5], Red], Point[pts]}]
seems to return the same spine; now I can define my region and triangulate it:
pr = ParametricRegion[{{r fx[t], r fy[t]}, 1 <= t <= 1 && 0 <= r <= 1}, {t, r}];
Ω = DiscretizeRegion[pr, MaxCellMeasure -> 0.001]
RegionPlot[Ω]
and similarly its boundary:
dpr = ParametricRegion[{{ fx[t], fy[t]}, 0 <= t <= 1}, t];
δΩ = DiscretizeRegion[dpr, MaxCellMeasure -> 0.001];
But applying the same PDE on these regions/boundary with these newly regions yields the same inaccuracies as before (boundary condition not satisfied properly on δΩ).
The problem might be with the second discretize region: indeed
Show[δΩ, Axes -> True]
presents some defect in the triangulation.
Note in particular the two points at the origin and at coordinate (0.9,-0.2)
.
Questions
Any suggestion on why it fails to satisfy the boundary?
Any suggestion on how to avoid going through
DiscretizeGraphics
?Any suggestion on how to specify
DirichletCondition
onBSplineFunction
?
I feel I am not using the most straightforward method here but…
Thanks!
BSplineFunction[]
and forming the correspondingParametricRegion[]
? $\endgroup$