PDEソリューションの1Dメッシュ生成

Dec 02 2020

時間と距離に依存する2つのPDEのシステム(H [x、t]とP [x、t])を解こうとしています。線の方法を使用して問題を解決していますが、メッシュを自分で生成してNDsolveに導入したいと思います。生成したいメッシュは以下の通りです

関数の1つ(P [x、t])の値はx = 0に非常に近い時間でのみ変化するのに対し、H [x、t]は領域全体で変化するため、このようなメッシュが必要です0 <x < xmax。以下は私が使用しているコードの例です

(* Constants *)
f = 38.94; logL = -2;  
Ls = 10^logL; a = 0.5;            
C1 = 1*^-5; dH = 1*^-6;       
Ea = 0.1; 
tmax = 40; (* Time in seconds *)
xmax = 10 Sqrt[dH] Sqrt[tmax]; (* Maximum distance to simulate. cm *)


(* PDE system *)
 eqsH = {D[H[x, t], t] - dH D[H[x, t], x, x] ==  NeumannValue[Ls Exp[a f Ea ] P[x, t] - Ls Exp[-a f Ea ] H[x, t], 
    x == 0], H[x, 0] == 1};
eqsP = {D[P[x, t], t] == NeumannValue[-Ls Exp[a f Ea ] P[x, t] + Ls Exp[-a f Ea ] H[x, t], 
     x == 0], P[x, 0] == 1};

(*Solution of the differential equations*)
prec = 7;
msf = 0.001; 

sol = NDSolve[{eqsH, eqsP}, {H, P}, {x, 0, xmax}, {t, 0, tmax}, 
     AccuracyGoal -> prec, PrecisionGoal -> prec, 
     Method -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"FiniteElement"}}] // First // 
   Quiet;

メッシュを作成してNDSolveに導入する方法についてサポートをいただけますか?前もって感謝します !

回答

9 TimLaska Dec 02 2020 at 23:33

これは、段階的メッシュを使用する別のアプローチです。

傾斜メッシュのいくつかのヘルパー関数を定義します

これは、1Dから3Dの異方性メッシュを作成するために使用したいくつかの関数です。すべての機能が使用されるわけではありません。

(*Import required FEM package*)
Needs["NDSolve`FEM`"];
(* Define Some Helper Functions For Structured Meshes*)
pointsToMesh[data_] :=
  MeshRegion[Transpose[{data}], 
   Line@Table[{i, i + 1}, {i, Length[data] - 1}]];
unitMeshGrowth[n_, r_] := 
 Table[(r^(j/(-1 + n)) - 1.)/(r - 1.), {j, 0, n - 1}]
meshGrowth[x0_, xf_, n_, r_] := (xf - x0) unitMeshGrowth[n, r] + x0
firstElmHeight[x0_, xf_, n_, r_] := 
 Abs@First@Differences@meshGrowth[x0, xf, n, r]
lastElmHeight[x0_, xf_, n_, r_] := 
 Abs@Last@Differences@meshGrowth[x0, xf, n, r]
findGrowthRate[x0_, xf_, n_, fElm_] := 
 Quiet@Abs@
   FindRoot[firstElmHeight[x0, xf, n, r] - fElm, {r, 1.0001, 100000}, 
     Method -> "Brent"][[1, 2]]
meshGrowthByElm[x0_, xf_, n_, fElm_] := 
 N@Sort@Chop@meshGrowth[x0, xf, n, findGrowthRate[x0, xf, n, fElm]]
meshGrowthByElm0[len_, n_, fElm_] := meshGrowthByElm[0, len, n, fElm]
flipSegment[l_] := (#1 - #2) & @@ {First[#], #} &@Reverse[l];
leftSegmentGrowth[len_, n_, fElm_] := meshGrowthByElm0[len, n, fElm]
rightSegmentGrowth[len_, n_, fElm_] := Module[{seg},
  seg = leftSegmentGrowth[len, n, fElm];
  flipSegment[seg]
  ]
reflectRight[pts_] := With[{rt = ReflectionTransform[{1}, {Last@pts}]},
  Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
reflectLeft[pts_] := 
 With[{rt = ReflectionTransform[{-1}, {First@pts}]},
  Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
extendMesh[mesh_, newmesh_] := Union[mesh, Max@mesh + newmesh]

段階的な水平メッシュセグメントを作成する

以下は、初期要素幅がドメイン長の1/10000である100要素の水平メッシュ領域を作成します。

(*Create a graded horizontal mesh segment*)
(*Initial element width is 1/10000 the domain length*)
seg = leftSegmentGrowth[xmax, 100, xmax/10000];
Print["Horizontal segment"]
rh = pointsToMesh@seg
(*Convert mesh region to element mesh*)
(*Extract Coords from horizontal region*)
crd = MeshCoordinates[rh];
(*Create element mesh*)
mesh = ToElementMesh[crd];
Print["ListPlot of exponential growth of element size"]
ListPlot[Transpose@mesh["Coordinates"]]

要素数が増えるにつれて、要素サイズが指数関数的に増加することがわかります。

メッシュ領域を要素メッシュに変換し、偏微分方程式を解きます

私は通常、MeshRegionを「ElementMesh」に変換して、必要に応じて要素マーカーとポイントマーカーを適用できるようにします。

(*Solve PDE on graded mesh*)
{Hfun, Pfun} = 
  NDSolveValue[{eqsH, eqsP}, {H, P}, x ∈ mesh, {t, 0, tmax}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"FiniteElement"}}];
(*Animate Hfun solution*)
imgs = Plot[Hfun[x, #], x ∈ mesh, 
     PlotRange -> {0.9999999, 1.0018}] & /@ Subdivide[0, tmax, 120];
Print["Animation of Hfun solution"]
ListAnimate@imgs

付録:異方性メッシュの例

以下のコメントでほのめかしたように、以下の箇条書きリストは、異方性クワッドメッシュを使用して、計算コストが非常に高くなるシャープなインターフェイスをキャプチャしたいくつかの例を示しています。コードは機能していますが、最適ではなく、一部の機能は時間の経過とともに変更されています。自己責任

  1. 2D-静止
    • MathematicaとMATLAB:境界条件が一定でない偏微分方程式で異なる結果が得られるのはなぜですか?
    • メッシュとNDSolveソリューションの収束の改善
  2. 2D-トランジェント
    • NDSolveValueで動的な時間ステップサイズを制御する
    • 膜を介した拡散をモデル化する方法は?
    • クワッドメッシュを使用した物質移動FEM
    • 異なるドメインで定義された未知の関数を持つ方程式システムを使用したNDSolve
  3. 3D-メッシュ
    • 傾斜メッシュを作成する
  4. 3D-静止
    • NDSolveでFEMソリューションを改善する方法は?
    • 3DFEMベクトルポテンシャル

境界層機能を備えたCOMSOLなどの他のツールにアクセスできる場合は、FEMAddOnsリソース関数を介してメッシュをインポートできます。現在MathematicaのFEMでサポートされていないプリズムやピラミッドのような追加の要素タイプを必要とする3Dメッシュでは機能しません。

5 AlexeiBoulbitch Dec 02 2020 at 22:36

これはどうですか?

lst1 = Partition[
  Join[Table[0.01*i, {i, 0, 5}], Table[0.1*i, {i, 0, 15}]], 1];
lst2 = Table[{i, i + 1}, {i, 1, Length[lst1] - 1}];

<< NDSolve`FEM`

mesh2 = ToElementMesh["Coordinates" -> lst1, 
  "MeshElements" -> {LineElement[lst2]}]

(*  ElementMesh[{{0., 1.5}}, {LineElement["<" 21 ">"]}]  *)

それを視覚化しましょう:

mesh2["Wireframe"["MeshElementIDStyle" -> Red]]

赤い数字はメッシュ要素を示しています。それらが重なる場所は、実際にはメッシュが10倍密度が高い場所です(下の拡大画像を参照)。

楽しんで!