add an example showing how to do convergence testing (#228)

python-hydro · Sep 1, 2024 · 0edb083 · 0edb083
1 parent 771daa7
commit 0edb083
Show file tree

Hide file tree

Showing 5 changed files with 242 additions and 10 deletions.
diff --git a/docs/source/compressible-convergence.ipynb b/docs/source/compressible-convergence.ipynb
@@ -0,0 +1,224 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "7ee63722-fe1b-41eb-99c2-ff6d437b16bd",
+   "metadata": {},
+   "source": [
+    "# Convergence of the compressible solvers"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e5d60167-70bb-44d4-a139-e0dd5646a426",
+   "metadata": {},
+   "source": [
+    "We'll look at convergence of the 2nd order `compressible` and 4th order\n",
+    "`compressible_fv4`  solvers using the `acoustic_pulse` problem and doing simple\n",
+    "Richardson convergence testing."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 19,
+   "id": "0c19f42b-16f1-48a8-ba19-e07f5addabd1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from pyro.pyro_sim import Pyro"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "79b31069-25a3-4926-be7a-9a23b7a6fe3a",
+   "metadata": {},
+   "source": [
+    "We want to keep $\\Delta t / \\Delta x$ constant as we test convergence so we will use a fixed timestep, following:\n",
+    "\n",
+    "$$\\Delta t = 3\\times 10^{-3} \\frac{64}{N}$$\n",
+    "\n",
+    "where $N$ is the number of zones in a dimension."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 20,
+   "id": "90900ff2-27b5-4642-a1de-006a9a30d975",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def timestep(N):\n",
+    "    return 3.e-3 * 64.0 / N"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7f0e962e-3728-4e3b-a4a9-9b3cd78e888c",
+   "metadata": {},
+   "source": [
+    "## `compressible`"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fa8b7497-d938-4ca2-8db6-74978baa81b9",
+   "metadata": {},
+   "source": [
+    "We'll run the problem at several different resolutions and store the `Pyro` simulation objects in a list."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cc7c0964-e0cf-43f4-8ca8-3ea6ed11c9fd",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\u001b[1mpyro ...\u001b[0m\n",
+      "\u001b[1minitializing the acoustic pulse problem...\u001b[0m\n",
+      "\u001b[1mpyro ...\u001b[0m\n",
+      "\u001b[1minitializing the acoustic pulse problem...\u001b[0m\n",
+      "\u001b[1mpyro ...\u001b[0m\n",
+      "\u001b[1minitializing the acoustic pulse problem...\u001b[0m\n",
+      "\u001b[1mpyro ...\u001b[0m\n",
+      "\u001b[1minitializing the acoustic pulse problem...\u001b[0m\n"
+     ]
+    }
+   ],
+   "source": [
+    "sims = []\n",
+    "\n",
+    "for N in [32, 64, 128, 256]:\n",
+    "    dt = timestep(N)\n",
+    "    params = {\"driver.fix_dt\": dt, \"mesh.nx\": N, \"mesh.ny\": N, \"driver.verbose\": 0}\n",
+    "    p = Pyro(\"compressible\")\n",
+    "    p.initialize_problem(problem_name=\"acoustic_pulse\", inputs_dict=params)\n",
+    "    p.run_sim()\n",
+    "    sims.append(p)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a624a88d-efd3-404e-9664-6346e191d00c",
+   "metadata": {},
+   "source": [
+    "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",
+    "for a single variable."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9705ab17-81c6-4b8a-becd-6a9af75371e1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from itertools import pairwise\n",
+    "var = \"density\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "97d051b5-563a-40ea-a838-9b4f7832380f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for coarse, fine in pairwise(sims):\n",
+    "    cvar = coarse.get_var(var)\n",
+    "    fvar = fine.sim.cc_data.restrict(var)\n",
+    "    e = cvar - fvar\n",
+    "    print(f\"{fine.get_grid().nx:3} -> {coarse.get_grid().nx:3} : {e.norm()}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "8d455059-3c7c-4597-a966-05f850eed570",
+   "metadata": {},
+   "source": [
+    "We see that the error is dropping by a factor of ~4 each time, indicating 2nd order convergence."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2caa12bc-d273-4bdc-af66-681ee30dde17",
+   "metadata": {},
+   "source": [
+    "## `compressible_fv4`"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "94e1b170-18ab-414c-a9f4-186f267c2d10",
+   "metadata": {},
+   "source": [
+    "Now we'll do the same for the 4th order solver.  We need to change the Riemann solver\n",
+    "to "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dd7a64cb-992e-4e0f-96f7-c8c03c0ca3eb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sims = []\n",
+    "\n",
+    "for N in [32, 64, 128, 256]:\n",
+    "    dt = timestep(N)\n",
+    "    params = {\"driver.fix_dt\": dt, \"mesh.nx\": N, \"mesh.ny\": N, \"driver.verbose\": 0}\n",
+    "    p = Pyro(\"compressible_fv4\")\n",
+    "    p.initialize_problem(problem_name=\"acoustic_pulse\", inputs_dict=params)\n",
+    "    p.run_sim()\n",
+    "    sims.append(p)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f03120c8-bc1d-4f0d-b79f-e498c64076a3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for coarse, fine in pairwise(sims):\n",
+    "    cvar = coarse.get_var(var)\n",
+    "    fvar = fine.sim.cc_data.restrict(var)\n",
+    "    e = cvar - fvar\n",
+    "    print(f\"{fine.get_grid().nx:3} -> {coarse.get_grid().nx:3} : {e.norm()}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c2c2fba6-cc5f-4ed3-ba0c-fb2bdfd509f3",
+   "metadata": {},
+   "source": [
+    "Now we see that the convergence is close to 4th order, with the error decreasing close to a factor of 16."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/docs/source/index.rst b/docs/source/index.rst
@@ -47,6 +47,12 @@ pyro: a python hydro code
    swe_basics
    particles_basics
 
+.. toctree::
+   :maxdepth: 1
+   :caption: Examples
+
+   compressible_convergence.ipynb
+
 .. toctree::
    :maxdepth: 1
    :caption: Utilities

diff --git a/pyro/compressible/problems/acoustic_pulse.py b/pyro/compressible/problems/acoustic_pulse.py
@@ -2,7 +2,7 @@
 
 import numpy as np
 
-from pyro.mesh import fv
+from pyro.mesh import fv, patch
 from pyro.util import msg
 
 DEFAULT_INPUTS = "inputs.acoustic_pulse"
@@ -15,7 +15,7 @@ def init_data(myd, rp):
     msg.bold("initializing the acoustic pulse problem...")
 
     # make sure that we are passed a valid patch object
-    if not isinstance(myd, fv.FV2d):
+    if not isinstance(myd, (patch.CellCenterData2d, fv.FV2d)):
         print("ERROR: patch invalid in acoustic_pulse.py")
         print(myd.__class__)
         sys.exit()

diff --git a/pyro/compressible_fv4/_defaults b/pyro/compressible_fv4/_defaults
@@ -21,5 +21,7 @@ temporal_method = RK4     ; integration method (see mesh/integration.py)
 
 grav = 0.0                ; gravitational acceleration (in y-direction)
 
+riemann = CGF
+
 
 
diff --git a/pyro/pyro_sim.py b/pyro/pyro_sim.py
@@ -75,13 +75,6 @@ def __init__(self, solver_name, from_commandline=False):
         self.rp.load_params(self.pyro_home + "_defaults")
         self.rp.load_params(self.pyro_home + self.solver_name + "/_defaults")
 
-        # manually override the dovis default
-        # for Jupyter, we want runtime vis disabled by default
-        if self.from_commandline:
-            self.rp.set_param("vis.dovis", 1)
-        else:
-            self.rp.set_param("vis.dovis", 0)
-
         self.tc = profile.TimerCollection()
 
         self.is_initialized = False
@@ -125,9 +118,16 @@ def initialize_problem(self, problem_name, inputs_file=None, inputs_dict=None,
 
             self.rp.load_params(inputs_file, no_new=1)
 
+        # manually override the dovis default
+        # for Jupyter, we want runtime vis disabled by default
+        if self.from_commandline:
+            self.rp.set_param("vis.dovis", 1)
+        else:
+            self.rp.set_param("vis.dovis", 0)
+
         if inputs_dict is not None:
             for k, v in inputs_dict.items():
-                self.rp.params[k] = v
+                self.rp.set_param(k, v)
 
         # and any commandline overrides
         if other_commands is not None:
Original file line number	Diff line number	Diff line change
Expand Up		@@ -21,5 +21,7 @@ temporal_method = RK4 ; integration method (see mesh/integration.py)

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

		riemann = CGF