Skip to main content

differential equations - How to specify PDE Boundary condition on a B-spline?


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:


Mathematica graphics


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]

Mathematica graphics


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]

Mathematica graphics


I can then discretize both the boundary and the region:



Ω = DiscretizeGraphics[pl]

Mathematica graphics


δΩ = DiscretizeGraphics[pl0, MaxCellMeasure -> 0.1]

Mathematica graphics


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]

Mathematica graphics


Even though it seems happy, it satisfies very poorly the boundary on the spine:


Plot[ P0 @@ ff0[t], {t, 0, 1}, ImageSize -> Small]

Mathematica graphics


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]}]

Mathematica graphics


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[Ω]


Mathematica graphics


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]


Mathematica graphics


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 on BSplineFunction?



I feel I am not using the most straightforward method here but…



Thanks!



Answer



The best way (as pointed out by J. M.) is to convert splines into implicit functions. The real issue you are having is that you'd need a second order mesh to get a decent solution. Note that DiscretizeGraphics and DiscretizeRegion create first order meshes. So you'd need to use ToElementMesh. We also would like to have a finer boundary resolution, thus use "MaxBoundaryCellMeasure". Another thing to think about is the way the boundary condition is specified on the spline. A better way to specify is to say "all boundary elements where R and z are not 0 instead of the code to rest for region member ship on the boundary with Element.


This then gives:


Needs["NDSolve`FEM`"]
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]}];
pr = ParametricRegion[{{r fx[t], r fy[t]}, -1 <= t <= 1 && 
     0 <= r <= 1}, {t, r}];
mesh = ToElementMesh[pr, "MaxBoundaryCellMeasure" -> 0.01];
mesh["Wireframe"]

enter image description here



Note that the mesh order is 2.


mesh["MeshOrder"]
2

From there we go to:


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 != 0 && z != 0], 

    DirichletCondition[
     P[R, z] == E R^2 Log[1/(R + $MachineEpsilon)^2], z == 0]
    }, P, {R, z} \[Element] mesh];

Note the $MachineEpsilon to avoid division by zero.


ContourPlot[P0[R, z], {R, z} \[Element] mesh, 
 PlotLegends -> Automatic, PlotPoints -> 30, 
 ColorFunction -> "LightTemperatureMap", PlotRange -> All, 
 FrameLabel -> {R, z}, AspectRatio -> 1]


enter image description here


And then this is about 2 order of magnitude better:


ff0 = BSplineFunction[pts];
Plot[P0 @@ ff0[t], {t, 0, 1}]

enter image description here


Which I hope is reasonable.


Note, that the boundary conditions are set to zero in the interpolating function:


bmesh = ToBoundaryMesh[mesh];
bc = bmesh["Coordinates"];

nodes = DeleteCases[bc, {x_ /; x < 10^-3, y_} | {x_, y_ /; y < 10^-3}];
MinMax[P0 @@@ nodes]
{-1.3877787807814457`*^-17, 2.7755575615628914`*^-17}

So what you see above is a an interpolation "limitiation" (it's "only" second order accurate). What I am not sure about is why it does not deteriorate further if the boundary is refined. Nevertheless, I think, it's OK to take the derivative of the interpolating function since doing that (currently V10.2) does not evaluates points that are not on the mesh.


Comments

Post a Comment

Popular posts from this blog

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...
Post a Comment
Read more

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...
Post a Comment
Read more

What is and isn't a valid variable specification for Manipulate?

I have an expression whose terms have arguments (representing subscripts), like this: myExpr = A[0] + V[1,T] I would like to put it inside a Manipulate to see its value as I move around the parameters. (The goal is eventually to plot it wrt one of the variables inside.) However, Mathematica complains when I set V[1,T] as a manipulated variable: Manipulate[Evaluate[myExpr], {A[0], 0, 1}, {V[1, T], 0, 1}] (*Manipulate::vsform: Manipulate argument {V[1,T],0,1} does not have the correct form for a variable specification. >> *) As a workaround, if I get rid of the symbol T inside the argument, it works fine: Manipulate[ Evaluate[myExpr /. T -> 15], {A[0], 0, 1}, {V[1, 15], 0, 1}] Why this behavior? Can anyone point me to the documentation that says what counts as a valid variable? And is there a way to get Manpiulate to accept an expression with a symbolic argument as a variable? Investigations I've done so far: I tried using variableQ from this answer , but it says V[1...
Post a Comment
Read more