Generación de malla 1D para solución PDE
Estoy tratando de resolver un sistema de dos PDE que dependen del tiempo y la distancia (H [x, t] y P [x, t]). Estoy resolviendo el problema usando el método de líneas pero quiero generar yo mismo la malla e introducirla en NDsolve. La malla que quiero generar es la siguiente
Necesito una malla como esta porque los valores de una de las funciones (P [x, t]) cambian con el tiempo solo muy cerca de x = 0, mientras que H [x, t] cambia en toda la región 0 <x < xmax. A continuación se muestra un ejemplo del código que estoy usando
(* 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;
¿Puedo obtener ayuda sobre cómo crear la malla e introducirla en NDSolve? gracias por adelantado !
Respuestas
Aquí hay un enfoque alternativo que utiliza una malla graduada.
Definir algunas funciones auxiliares para una malla graduada
Aquí hay algunas funciones que he usado para crear mallas anisotrópicas 1d a 3D. No se utilizan todas las funciones.
(*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]
Crear un segmento de malla horizontal graduado
Lo siguiente creará una región de malla horizontal de 100 elementos donde el ancho del elemento inicial es 1/10000 de la longitud del dominio.
(*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"]]
Se puede ver el crecimiento exponencial del tamaño del elemento a medida que aumenta el número del elemento.
Convierta la región de la malla en una malla de elementos y resuelva el PDE
Por lo general, convierto el MeshRegiona 'ElementMesh' para poder aplicar marcadores de elementos y puntos si es necesario.
(*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
Apéndice: ejemplos de mallado anisotrópico
Como mencioné en el comentario a continuación, la lista de viñetas a continuación muestra varios ejemplos en los que he utilizado cuadrícula anisotrópica para capturar interfaces nítidas que de otro modo serían muy costosas desde el punto de vista informático. El código es funcional, pero no óptimo y algunas de las funciones se han modificado con el tiempo. Úselo bajo su propio riesgo
- 2D-estacionario
- Mathematica vs.MATLAB: ¿por qué obtengo resultados diferentes para PDE con una condición de límite no constante?
- Mejora de la convergencia de soluciones de malla y NDSolve
- Transitorio 2D
- Controlar el tamaño del paso de tiempo dinámico en NDSolveValue
- ¿Cómo modelar la difusión a través de una membrana?
- Transporte masivo FEM usando Quad Mesh
- NDSolve con un sistema de ecuaciones con funciones desconocidas definidas en diferentes dominios
- Mallado 3D
- Crear malla graduada
- 3D-estacionario
- ¿Cómo mejorar la solución FEM con NDSolve?
- Potencial de vector 3D FEM
Si tiene acceso a otras herramientas, como COMSOL, que tienen funcionalidad de capa límite, puede importar mallas a través de la función de recursos FEMAddOns . No funcionará para mallas 3D que requieran tipos de elementos adicionales como prismas y pirámides que actualmente no son compatibles con FEM de Mathematica .
¿Qué pasa con esto?
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 ">"]}] *)
Visualicémoslo:
mesh2["Wireframe"["MeshElementIDStyle" -> Red]]
Las figuras rojas indican los elementos de la malla. El lugar donde se superponen es, de hecho, donde la malla es 10 veces más densa (vea la imagen ampliada a continuación):
¡Que te diviertas!