Generowanie siatki 1D dla rozwiązania PDE
Próbuję rozwiązać układ dwóch PDE, które są zależne od czasu i odległości (H [x, t] i P [x, t]). Rozwiązuję problem metodą linii, ale chcę samodzielnie wygenerować siatkę i wprowadzić ją w NDsolve. Siatka, którą chcę wygenerować, jest następująca
Potrzebuję siatki takiej jak ta, ponieważ wartości jednej z funkcji (P [x, t]) zmieniają się w czasie tylko bardzo blisko x = 0, podczas gdy H [x, t] zmienia się w całym regionie 0 <x < xmax. Poniżej znajduje się przykład kodu, którego używam
(* 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;
Czy mogę uzyskać pomoc w tworzeniu siatki i wprowadzaniu jej w NDSolve? z góry dziękuję !
Odpowiedzi
Oto alternatywne podejście z użyciem siatki stopniowanej.
Zdefiniuj funkcje pomocnicze dla siatki stopniowanej
Oto kilka funkcji, których używałem do tworzenia siatek anizotropowych 1d do 3D. Nie wszystkie funkcje są używane.
(*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]
Utwórz stopniowany poziomy segment siatki
Poniższe działania utworzą poziomy obszar siatki składający się ze 100 elementów, w którym początkowa szerokość elementu będzie równa 1/10000 długości domeny.
(*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"]]
Widać wykładniczy wzrost rozmiaru elementu wraz ze wzrostem liczby elementów.
Przekształć region siatki w siatkę elementu i rozwiąż PDE
Generalnie konwertuję na MeshRegion„ElementMesh”, aby w razie potrzeby zastosować znaczniki elementów i punktów.
(*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
Dodatek: przykłady siatki anizotropowej
Jak wspomniałem w poniższym komentarzu, poniższa lista punktowana pokazuje kilka przykładów, w których użyłem anizotropowego poczwórnego siatki, aby uchwycić ostre interfejsy, które w przeciwnym razie byłyby bardzo kosztowne obliczeniowo. Kod działa, ale nie jest optymalny, a niektóre funkcje zostały z czasem zmodyfikowane. Używaj na własne ryzyko
- 2D-stacjonarne
- Mathematica kontra MATLAB: dlaczego otrzymuję różne wyniki dla PDE z niestałym warunkiem brzegowym?
- Poprawa konwergencji rozwiązań typu mesh i NDSolve
- Przejściowe 2D
- Kontrolowanie dynamicznego rozmiaru kroku czasowego w NDSolveValue
- Jak modelować dyfuzję przez membranę?
- Transport masowy MES przy użyciu poczwórnej siatki
- NDSolve z układem równań z nieznanymi funkcjami zdefiniowanymi w różnych dziedzinach
- Siatkowanie 3D
- Utwórz stopniowaną siatkę
- 3D-stacjonarne
- Jak ulepszyć rozwiązanie FEM za pomocą NDSolve?
- Potencjał wektora MES 3D
Jeśli masz dostęp do innych narzędzi, takich jak COMSOL, które mają funkcjonalność warstwy granicznej, możesz importować siatki za pomocą funkcji zasobów FEMAddOns . Nie będzie działać dla siatek 3D, które wymagają dodatkowych typów elementów, takich jak pryzmaty i piramidy, które nie są obecnie obsługiwane w MES Mathematica .
A co z tym?
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 ">"]}] *)
Wizualizujmy to:
mesh2["Wireframe"["MeshElementIDStyle" -> Red]]
Czerwone cyfry oznaczają elementy siatki. Miejsce, w którym zachodzą na siebie, jest w rzeczywistości tym, w którym siatka jest 10 razy gęstsza (patrz zdjęcie poniżej):
Baw się dobrze!