Skip to content

Commit 0edb083

Browse files
zingalezingale
authored
add an example showing how to do convergence testing (#228)
1 parent 771daa7 commit 0edb083

File tree

5 files changed

+242
-10
lines changed

5 files changed

+242
-10
lines changed

‎docs/source/compressible-convergence.ipynb‎

Lines changed: 224 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,224 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"id": "7ee63722-fe1b-41eb-99c2-ff6d437b16bd",
6+
"metadata": {},
7+
"source": [
8+
"# Convergence of the compressible solvers"
9+
]
10+
},
11+
{
12+
"cell_type": "markdown",
13+
"id": "e5d60167-70bb-44d4-a139-e0dd5646a426",
14+
"metadata": {},
15+
"source": [
16+
"We'll look at convergence of the 2nd order `compressible` and 4th order\n",
17+
"`compressible_fv4` solvers using the `acoustic_pulse` problem and doing simple\n",
18+
"Richardson convergence testing."
19+
]
20+
},
21+
{
22+
"cell_type": "code",
23+
"execution_count": 19,
24+
"id": "0c19f42b-16f1-48a8-ba19-e07f5addabd1",
25+
"metadata": {},
26+
"outputs": [],
27+
"source": [
28+
"from pyro.pyro_sim import Pyro"
29+
]
30+
},
31+
{
32+
"cell_type": "markdown",
33+
"id": "79b31069-25a3-4926-be7a-9a23b7a6fe3a",
34+
"metadata": {},
35+
"source": [
36+
"We want to keep $\\Delta t / \\Delta x$ constant as we test convergence so we will use a fixed timestep, following:\n",
37+
"\n",
38+
"$$\\Delta t = 3\\times 10^{-3} \\frac{64}{N}$$\n",
39+
"\n",
40+
"where $N$ is the number of zones in a dimension."
41+
]
42+
},
43+
{
44+
"cell_type": "code",
45+
"execution_count": 20,
46+
"id": "90900ff2-27b5-4642-a1de-006a9a30d975",
47+
"metadata": {},
48+
"outputs": [],
49+
"source": [
50+
"def timestep(N):\n",
51+
" return 3.e-3 * 64.0 / N"
52+
]
53+
},
54+
{
55+
"cell_type": "markdown",
56+
"id": "7f0e962e-3728-4e3b-a4a9-9b3cd78e888c",
57+
"metadata": {},
58+
"source": [
59+
"## `compressible`"
60+
]
61+
},
62+
{
63+
"cell_type": "markdown",
64+
"id": "fa8b7497-d938-4ca2-8db6-74978baa81b9",
65+
"metadata": {},
66+
"source": [
67+
"We'll run the problem at several different resolutions and store the `Pyro` simulation objects in a list."
68+
]
69+
},
70+
{
71+
"cell_type": "code",
72+
"execution_count": null,
73+
"id": "cc7c0964-e0cf-43f4-8ca8-3ea6ed11c9fd",
74+
"metadata": {},
75+
"outputs": [
76+
{
77+
"name": "stdout",
78+
"output_type": "stream",
79+
"text": [
80+
"\u001b[1mpyro ...\u001b[0m\n",
81+
"\u001b[1minitializing the acoustic pulse problem...\u001b[0m\n",
82+
"\u001b[1mpyro ...\u001b[0m\n",
83+
"\u001b[1minitializing the acoustic pulse problem...\u001b[0m\n",
84+
"\u001b[1mpyro ...\u001b[0m\n",
85+
"\u001b[1minitializing the acoustic pulse problem...\u001b[0m\n",
86+
"\u001b[1mpyro ...\u001b[0m\n",
87+
"\u001b[1minitializing the acoustic pulse problem...\u001b[0m\n"
88+
]
89+
}
90+
],
91+
"source": [
92+
"sims = []\n",
93+
"\n",
94+
"for N in [32, 64, 128, 256]:\n",
95+
" dt = timestep(N)\n",
96+
" params = {\"driver.fix_dt\": dt, \"mesh.nx\": N, \"mesh.ny\": N, \"driver.verbose\": 0}\n",
97+
" p = Pyro(\"compressible\")\n",
98+
" p.initialize_problem(problem_name=\"acoustic_pulse\", inputs_dict=params)\n",
99+
" p.run_sim()\n",
100+
" sims.append(p)"
101+
]
102+
},
103+
{
104+
"cell_type": "markdown",
105+
"id": "a624a88d-efd3-404e-9664-6346e191d00c",
106+
"metadata": {},
107+
"source": [
108+
"Now we want to loop over each adjacent pair of simulations, coarsen the finer resolution simulation and compute the norm of the difference. We'll do this\n",
109+
"for a single variable."
110+
]
111+
},
112+
{
113+
"cell_type": "code",
114+
"execution_count": null,
115+
"id": "9705ab17-81c6-4b8a-becd-6a9af75371e1",
116+
"metadata": {},
117+
"outputs": [],
118+
"source": [
119+
"from itertools import pairwise\n",
120+
"var = \"density\""
121+
]
122+
},
123+
{
124+
"cell_type": "code",
125+
"execution_count": null,
126+
"id": "97d051b5-563a-40ea-a838-9b4f7832380f",
127+
"metadata": {},
128+
"outputs": [],
129+
"source": [
130+
"for coarse, fine in pairwise(sims):\n",
131+
" cvar = coarse.get_var(var)\n",
132+
" fvar = fine.sim.cc_data.restrict(var)\n",
133+
" e = cvar - fvar\n",
134+
" print(f\"{fine.get_grid().nx:3} -> {coarse.get_grid().nx:3} : {e.norm()}\")"
135+
]
136+
},
137+
{
138+
"cell_type": "markdown",
139+
"id": "8d455059-3c7c-4597-a966-05f850eed570",
140+
"metadata": {},
141+
"source": [
142+
"We see that the error is dropping by a factor of ~4 each time, indicating 2nd order convergence."
143+
]
144+
},
145+
{
146+
"cell_type": "markdown",
147+
"id": "2caa12bc-d273-4bdc-af66-681ee30dde17",
148+
"metadata": {},
149+
"source": [
150+
"## `compressible_fv4`"
151+
]
152+
},
153+
{
154+
"cell_type": "markdown",
155+
"id": "94e1b170-18ab-414c-a9f4-186f267c2d10",
156+
"metadata": {},
157+
"source": [
158+
"Now we'll do the same for the 4th order solver. We need to change the Riemann solver\n",
159+
"to "
160+
]
161+
},
162+
{
163+
"cell_type": "code",
164+
"execution_count": null,
165+
"id": "dd7a64cb-992e-4e0f-96f7-c8c03c0ca3eb",
166+
"metadata": {},
167+
"outputs": [],
168+
"source": [
169+
"sims = []\n",
170+
"\n",
171+
"for N in [32, 64, 128, 256]:\n",
172+
" dt = timestep(N)\n",
173+
" params = {\"driver.fix_dt\": dt, \"mesh.nx\": N, \"mesh.ny\": N, \"driver.verbose\": 0}\n",
174+
" p = Pyro(\"compressible_fv4\")\n",
175+
" p.initialize_problem(problem_name=\"acoustic_pulse\", inputs_dict=params)\n",
176+
" p.run_sim()\n",
177+
" sims.append(p)"
178+
]
179+
},
180+
{
181+
"cell_type": "code",
182+
"execution_count": null,
183+
"id": "f03120c8-bc1d-4f0d-b79f-e498c64076a3",
184+
"metadata": {},
185+
"outputs": [],
186+
"source": [
187+
"for coarse, fine in pairwise(sims):\n",
188+
" cvar = coarse.get_var(var)\n",
189+
" fvar = fine.sim.cc_data.restrict(var)\n",
190+
" e = cvar - fvar\n",
191+
" print(f\"{fine.get_grid().nx:3} -> {coarse.get_grid().nx:3} : {e.norm()}\")"
192+
]
193+
},
194+
{
195+
"cell_type": "markdown",
196+
"id": "c2c2fba6-cc5f-4ed3-ba0c-fb2bdfd509f3",
197+
"metadata": {},
198+
"source": [
199+
"Now we see that the convergence is close to 4th order, with the error decreasing close to a factor of 16."
200+
]
201+
}
202+
],
203+
"metadata": {
204+
"kernelspec": {
205+
"display_name": "Python 3 (ipykernel)",
206+
"language": "python",
207+
"name": "python3"
208+
},
209+
"language_info": {
210+
"codemirror_mode": {
211+
"name": "ipython",
212+
"version": 3
213+
},
214+
"file_extension": ".py",
215+
"mimetype": "text/x-python",
216+
"name": "python",
217+
"nbconvert_exporter": "python",
218+
"pygments_lexer": "ipython3",
219+
"version": "3.12.4"
220+
}
221+
},
222+
"nbformat": 4,
223+
"nbformat_minor": 5
224+
}

‎docs/source/index.rst‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -47,6 +47,12 @@ pyro: a python hydro code
4747
swe_basics
4848
particles_basics
4949

50+
.. toctree::
51+
:maxdepth: 1
52+
:caption: Examples
53+
54+
compressible_convergence.ipynb
55+
5056
.. toctree::
5157
:maxdepth: 1
5258
:caption: Utilities

‎pyro/compressible/problems/acoustic_pulse.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22

33
import numpy as np
44

5-
from pyro.mesh import fv
5+
from pyro.mesh import fv, patch
66
from pyro.util import msg
77

88
DEFAULT_INPUTS = "inputs.acoustic_pulse"
@@ -15,7 +15,7 @@ def init_data(myd, rp):
1515
msg.bold("initializing the acoustic pulse problem...")
1616

1717
# make sure that we are passed a valid patch object
18-
if not isinstance(myd, fv.FV2d):
18+
if not isinstance(myd, (patch.CellCenterData2d, fv.FV2d)):
1919
print("ERROR: patch invalid in acoustic_pulse.py")
2020
print(myd.__class__)
2121
sys.exit()

‎pyro/compressible_fv4/_defaults‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -21,5 +21,7 @@ temporal_method = RK4 ; integration method (see mesh/integration.py)
2121

2222
grav = 0.0 ; gravitational acceleration (in y-direction)
2323

24+
riemann = CGF
25+
2426

2527

‎pyro/pyro_sim.py‎

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -75,13 +75,6 @@ def __init__(self, solver_name, from_commandline=False):
7575
self.rp.load_params(self.pyro_home + "_defaults")
7676
self.rp.load_params(self.pyro_home + self.solver_name + "/_defaults")
7777

78-
# manually override the dovis default
79-
# for Jupyter, we want runtime vis disabled by default
80-
if self.from_commandline:
81-
self.rp.set_param("vis.dovis", 1)
82-
else:
83-
self.rp.set_param("vis.dovis", 0)
84-
8578
self.tc = profile.TimerCollection()
8679

8780
self.is_initialized = False
@@ -125,9 +118,16 @@ def initialize_problem(self, problem_name, inputs_file=None, inputs_dict=None,
125118

126119
self.rp.load_params(inputs_file, no_new=1)
127120

121+
# manually override the dovis default
122+
# for Jupyter, we want runtime vis disabled by default
123+
if self.from_commandline:
124+
self.rp.set_param("vis.dovis", 1)
125+
else:
126+
self.rp.set_param("vis.dovis", 0)
127+
128128
if inputs_dict is not None:
129129
for k, v in inputs_dict.items():
130-
self.rp.params[k] = v
130+
self.rp.set_param(k, v)
131131

132132
# and any commandline overrides
133133
if other_commands is not None:

0 commit comments

Comments
 (0)