Generasi mesh 1D untuk solusi PDE
Saya mencoba menyelesaikan sistem dua PDE yang bergantung pada waktu dan jarak (H [x, t] dan P [x, t]). Saya memecahkan masalah menggunakan metode garis tetapi saya ingin membuat sendiri mesh dan memperkenalkannya di NDsolve. Jala yang ingin saya buat adalah sebagai berikut
Saya membutuhkan mesh seperti ini karena nilai untuk salah satu fungsi (P [x, t]) berubah dengan waktu hanya sangat dekat dengan x = 0, sedangkan H [x, t] berubah di seluruh wilayah 0 <x < xmax. Di bawah ini adalah contoh kode yang saya gunakan
(* 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;
Bisakah saya mendapatkan bantuan tentang cara membuat mesh dan memperkenalkannya di NDSolve? Terima kasih sebelumnya !
Jawaban
Berikut adalah pendekatan alternatif menggunakan jaring bertingkat.
Tentukan beberapa fungsi pembantu untuk mesh bergradasi
Berikut adalah beberapa fungsi yang saya gunakan untuk membuat mesh anisotropik 1d ke 3D. Tidak semua fungsi digunakan.
(*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]
Buat segmen jaring horizontal bergradasi
Berikut ini akan membuat daerah mesh horizontal dari 100 elemen dengan lebar elemen awal 1/10000 dari panjang domain.
(*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"]]
Seseorang dapat melihat pertumbuhan eksponensial dari ukuran elemen dengan bertambahnya jumlah elemen.
Ubah wilayah mesh menjadi mesh elemen dan selesaikan PDE
Saya biasanya mengonversi MeshRegionmenjadi 'ElementMesh'so sehingga saya dapat menerapkan elemen dan penanda titik jika diperlukan.
(*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
Lampiran: Contoh meshing anisotropik
Seperti yang saya singgung dalam komentar di bawah, daftar peluru di bawah ini menunjukkan beberapa contoh di mana saya telah menggunakan meshing quad anisotropik untuk menangkap antarmuka tajam yang jika tidak akan menjadi sangat mahal secara komputasi. Kode tersebut berfungsi, tetapi tidak optimal dan beberapa fungsi telah dimodifikasi seiring waktu. Gunakan dengan resiko Anda sendiri
- 2D-Stasioner
- Mathematica vs. MATLAB: Mengapa saya mendapatkan hasil yang berbeda untuk PDE dengan kondisi batas tidak konstan?
- Meningkatkan konvergensi solusi mesh dan NDSolve
- 2D-Transien
- Mengontrol ukuran langkah waktu dinamis di NDSolveValue
- Bagaimana memodelkan difusi melalui membran?
- FEM Transportasi Massal Menggunakan Quad Mesh
- Selesaikan dengan sistem persamaan dengan fungsi yang tidak diketahui yang ditentukan pada domain yang berbeda
- 3D-Meshing
- Buat jaring bertingkat
- 3D-Stasioner
- Bagaimana cara Meningkatkan Solusi FEM dengan NDSolve?
- Potensi Vektor 3D FEM
Jika Anda memiliki akses ke alat lain, seperti COMSOL, yang memiliki fungsionalitas lapisan batas, Anda dapat mengimpor mesh melalui fungsi sumber daya FEMAddOns . Ini tidak akan berfungsi untuk mesh 3D yang memerlukan jenis elemen tambahan seperti prisma dan piramida yang saat ini tidak didukung di FEM Mathematica .
Bagaimana dengan ini?
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 ">"]}] *)
Mari kita visualisasikan:
mesh2["Wireframe"["MeshElementIDStyle" -> Red]]
Gambar merah menunjukkan elemen jaring. Tempat di mana mereka tumpang tindih sebenarnya adalah tempat di mana mesh 10 kali lebih padat (lihat gambar yang diledakkan di bawah):
Selamat bersenang-senang!