PDE 솔루션을위한 1D 메시 생성

Dec 02 2020

나는 시간과 거리 (H [x, t]와 P [x, t])에 의존하는 두 개의 PDE 시스템을 풀려고합니다. 선 방법을 사용하여 문제를 해결하고 있지만 메시를 직접 생성하여 NDsolve에 도입하고 싶습니다. 생성하고자하는 메시는 다음과 같습니다.

함수 (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"]]

요소 수가 증가함에 따라 요소 크기가 기하 급수적으로 증가하는 것을 볼 수 있습니다.

메쉬 영역을 요소 메쉬로 변환하고 PDE 풀기

일반적으로 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

부록 : 이방성 메싱 예제

아래 주석에서 언급했듯이 아래 글 머리 기호 목록은 이방성 쿼드 메싱을 사용하여 계산 비용이 많이 드는 날카로운 인터페이스를 캡처 한 몇 가지 예를 보여줍니다. 코드는 작동하지만 최적이 아니며 일부 기능은 시간이 지남에 따라 수정되었습니다. 자신의 책임하에 사용

2D 고정
- Mathematica vs. MATLAB : 경계 조건이 일정하지 않은 PDE에 대해 다른 결과를 얻는 이유는 무엇입니까?
- 메시 및 NDSolve 솔루션 수렴 개선
2D 과도
- NDSolveValue에서 동적 시간 단계 크기 제어
- 멤브레인을 통한 확산을 모델링하는 방법은 무엇입니까?
- Quad Mesh를 사용한 Mass Transport FEM
- 다른 도메인에 정의 된 알 수없는 함수가있는 방정식 시스템이있는 NDSolve
3D- 메싱
- 그레이드 메쉬 생성
3D 고정
- NDSolve로 FEM 솔루션을 개선하는 방법은 무엇입니까?
- 3D FEM 벡터 잠재력

경계 레이어 기능이있는 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 배 더 조밀 한 곳입니다 (아래의 확대 된 이미지 참조).