Génération de maillage 1D pour solution PDE
J'essaie de résoudre un système de deux PDE qui dépendent du temps et de la distance (H [x, t] et P [x, t]). Je résous le problème en utilisant la méthode des lignes mais je veux générer moi-même le maillage et l'introduire dans NDsolve. Le maillage que je souhaite générer est le suivant
J'ai besoin d'un maillage comme celui-ci car les valeurs de l'une des fonctions (P [x, t]) ne changent avec le temps que très près de x = 0, alors que H [x, t] change sur toute la région 0 <x < xmax. Voici un exemple du code que j'utilise
(* 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;
Puis-je obtenir de l'aide pour créer le maillage et l'introduire dans NDSolve? Merci d'avance !
Réponses
Voici une autre approche utilisant un maillage gradué.
Définir des fonctions d'aide pour un maillage gradué
Voici quelques fonctions que j'ai utilisées pour créer des maillages anisotropes 1d à 3D. Toutes les fonctions ne sont pas utilisées.
(*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]
Créer un segment de maillage horizontal gradué
Ce qui suit va créer une région de maillage horizontal de 100 éléments où la largeur initiale de l'élément est 1/10000 de la longueur du domaine.
(*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"]]
On peut voir la croissance exponentielle de la taille de l'élément à mesure que le nombre d'élément augmente.
Convertissez la région maillée en maillage d'élément et résolvez le PDE
Je convertis généralement le MeshRegionen un 'ElementMesh' pour que je puisse appliquer des marqueurs d'élément et de point si nécessaire.
(*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
Annexe: Exemples de maillage anisotrope
Comme je l'ai mentionné dans le commentaire ci-dessous, la liste à puces ci-dessous montre plusieurs exemples dans lesquels j'ai utilisé un maillage quadruplé anisotrope pour capturer des interfaces nettes qui seraient autrement très coûteuses en calcul. Le code est fonctionnel, mais pas optimal et certaines des fonctions ont été modifiées au fil du temps. À utiliser à vos risques et périls
- 2D-stationnaire
- Mathematica vs MATLAB: pourquoi est-ce que j'obtiens des résultats différents pour PDE avec une condition aux limites non constante?
- Amélioration de la convergence des solutions de maillage et NDSolve
- 2D-transitoire
- Contrôle de la taille des pas de temps dynamiques dans NDSolveValue
- Comment modéliser la diffusion à travers une membrane?
- Transport de masse MEF à l'aide de quadruple maillage
- NDSolve avec un système d'équation avec des fonctions inconnues définies sur différents domaines
- Maillage 3D
- Créer un maillage gradué
- 3D-stationnaire
- Comment améliorer la solution FEM avec NDSolve?
- Potentiel vectoriel FEM 3D
Si vous avez accès à d'autres outils, tels que COMSOL, qui ont une fonctionnalité de couche limite, vous pouvez importer des maillages via la fonction de ressource FEMAddOns . Cela ne fonctionnera pas pour les maillages 3D qui nécessitent des types d'éléments supplémentaires tels que les prismes et les pyramides qui ne sont actuellement pas pris en charge dans le FEM de Mathematica .
Et ça?
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 ">"]}] *)
Laissez-nous le visualiser:
mesh2["Wireframe"["MeshElementIDStyle" -> Red]]
Les chiffres rouges indiquent les éléments du maillage. L'endroit où ils se chevauchent est en fait celui où le maillage est 10 fois plus dense (voir l'image agrandie ci-dessous):
S'amuser!