From cf90d61abe22bc02b63bcc44d7a4ea7b170c3110 Mon Sep 17 00:00:00 2001 From: zhi Date: Thu, 21 Nov 2024 14:02:33 -0500 Subject: [PATCH 1/9] add dloga for grids --- pyro/mesh/patch.py | 15 ++++++++++++++- pyro/mesh/tests/test_patch.py | 21 +++++++++++++++++++++ 2 files changed, 35 insertions(+), 1 deletion(-) diff --git a/pyro/mesh/patch.py b/pyro/mesh/patch.py index 721bfc9d1..42e011496 100644 --- a/pyro/mesh/patch.py +++ b/pyro/mesh/patch.py @@ -220,6 +220,11 @@ def __init__(self, nx, ny, *, ng=1, self.Ay = self.Lx + # Spatial derivative of log(Area) + + self.dlogAx = np.zeros_like(self.Ax) + self.dlogAy = np.zeros_like(self.Ay) + # Volume of the cell. self.V = ArrayIndexer(np.full((self.qx, self.qy), self.dx * self.dy), @@ -235,7 +240,9 @@ def __str__(self): class SphericalPolar(Grid2d): """ This class defines a spherical polar grid. - This is technically a 2D geometry but assumes azimuthal symmetry. + This is technically a 3D geometry but assumes azimuthal symmetry, + and zero velocity in phi-direction. + Hence 2D. Define: r = x @@ -279,6 +286,12 @@ def __init__(self, nx, ny, *, ng=1, self.Ay = np.abs(np.pi * np.sin(self.yl2d) * (self.xr2d**2 - self.xl2d**2)) + # dlogAx = 1 / r^2 d( r^2 ) / dr = 2 / r + self.dlogAx = 2.0 / self.x2d + + # dlogAy = 1 / (r sin(theta)) d( sin(theta) )/dtheta = cot(theta) / r + self.dlogAy = 1.0 / (np.tan(self.y2d) * self.x2d) + # Returns an array of the volume of each cell. # dV = dL_r * dL_theta * dL_phi # = (dr) * (r * dtheta) * (r * sin(theta) * dphi) diff --git a/pyro/mesh/tests/test_patch.py b/pyro/mesh/tests/test_patch.py index 3385b788d..b3ffb3ac3 100644 --- a/pyro/mesh/tests/test_patch.py +++ b/pyro/mesh/tests/test_patch.py @@ -92,6 +92,12 @@ def test_Ax(self): def test_Ay(self): assert np.all(self.g.Ay.v() == 0.25) + def test_dlogAx(self): + assert np.all(self.g.dlogAx.v() == 0.0) + + def test_dlogAy(self): + assert np.all(self.g.dlogAy.v() == 0.0) + def test_V(self): assert np.all(self.g.V.v() == 0.1 * 0.25) @@ -135,6 +141,21 @@ def test_Ay(self): (self.g.xr2d.v()**2 - self.g.xl2d.v()**2)) assert_array_equal(self.g.Ay.jp(1), area_y_r) + def test_dlogAx(self): + i = 4 + r = self.g.xmin + (i + 0.5 - self.g.ng) * self.g.dx + dlogAx = 2.0 / r + assert (self.g.dlogAx[i,:] == dlogAx).all() + + def test_dlogAy(self): + i = 4 + j = 6 + r = self.g.xmin + (i + 0.5 - self.g.ng) * self.g.dx + tan = np.tan(self.g.ymin + (j + 0.5 - self.g.ng) * self.g.dy) + dlogAy = 1.0 / (r * tan) + + assert self.g.dlogAy[i,j] == dlogAy + def test_V(self): volume = np.abs(-2.0 * np.pi / 3.0 * (np.cos(self.g.yr2d.v()) - np.cos(self.g.yl2d.v())) * From 7e20ddd6f1f61cfe53083fb4fa48caee46fd2483 Mon Sep 17 00:00:00 2001 From: zhi Date: Thu, 21 Nov 2024 14:58:16 -0500 Subject: [PATCH 2/9] add primitive geometric source terms --- pyro/compressible/interface.py | 24 +++++++++++++++++++++++- pyro/compressible/unsplit_fluxes.py | 6 ++++-- 2 files changed, 27 insertions(+), 3 deletions(-) diff --git a/pyro/compressible/interface.py b/pyro/compressible/interface.py index c6bca9717..4fd0a2bc6 100644 --- a/pyro/compressible/interface.py +++ b/pyro/compressible/interface.py @@ -3,7 +3,7 @@ @njit(cache=True) -def states(idir, ng, dx, dt, +def states(idir, ng, dx, dloga, dt, irho, iu, iv, ip, ix, nspec, gamma, qv, dqv): r""" @@ -212,6 +212,28 @@ def states(idir, ng, dx, dt, q_l[i, j + 1, m] = q_l[i, j + 1, m] + sum_l q_r[i, j, m] = q_r[i, j, m] + sum_r + + # Geometric Source term from converting conserved-variable to primitive + # It's only there for non Cartesian coord. + + if idir == 1: + rho_source = -0.5 * dt * dloga[i,j] * q[iu] + + q_l[i + 1, j, irho] += rho_source + q_r[i, j, irho] += rho_source + + q_l[i + 1, j, ip] += rho_source * cs * cs + q_r[i, j, ip] += rho_source * cs * cs + + else: + rho_source = -0.5 * dt * dloga[i,j] * q[iv] + + q_l[i, j + 1, irho] += rho_source + q_r[i, j, irho] += rho_source + + q_l[i, j + 1, ip] += rho_source * cs * cs + q_r[i, j, ip] += rho_source * cs * cs + return q_l, q_r diff --git a/pyro/compressible/unsplit_fluxes.py b/pyro/compressible/unsplit_fluxes.py index 64c656a9b..443d35197 100644 --- a/pyro/compressible/unsplit_fluxes.py +++ b/pyro/compressible/unsplit_fluxes.py @@ -206,11 +206,13 @@ def interface_states(my_data, rp, ivars, tc, dt): tm_states = tc.timer("interfaceStates") tm_states.begin() - V_l, V_r = ifc.states(1, myg.ng, myg.Lx, dt, + _V_l, _V_r = ifc.states(1, myg.ng, myg.Lx, myg.dlogAx, dt, ivars.irho, ivars.iu, ivars.iv, ivars.ip, ivars.ix, ivars.naux, gamma, q, ldx) + V_l = ai.ArrayIndexer(d=_V_l, grid=myg) + V_r = ai.ArrayIndexer(d=_V_r, grid=myg) tm_states.end() @@ -225,7 +227,7 @@ def interface_states(my_data, rp, ivars, tc, dt): # left and right primitive variable states tm_states.begin() - _V_l, _V_r = ifc.states(2, myg.ng, myg.Ly, dt, + _V_l, _V_r = ifc.states(2, myg.ng, myg.Ly, myg.dlogAy, dt, ivars.irho, ivars.iu, ivars.iv, ivars.ip, ivars.ix, ivars.naux, gamma, From 3ea3428f90cab7618bd71b7ccdc4bd6202d90dc9 Mon Sep 17 00:00:00 2001 From: zhi Date: Thu, 21 Nov 2024 17:04:21 -0500 Subject: [PATCH 3/9] update comment --- pyro/mesh/patch.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyro/mesh/patch.py b/pyro/mesh/patch.py index 42e011496..a79a5ecde 100644 --- a/pyro/mesh/patch.py +++ b/pyro/mesh/patch.py @@ -220,7 +220,7 @@ def __init__(self, nx, ny, *, ng=1, self.Ay = self.Lx - # Spatial derivative of log(Area) + # Spatial derivative of log(area). It's zero for Cartesian self.dlogAx = np.zeros_like(self.Ax) self.dlogAy = np.zeros_like(self.Ay) From a4201208e4252ff6760eba5f333edb31b4dabe9f Mon Sep 17 00:00:00 2001 From: zhi Date: Thu, 21 Nov 2024 17:08:29 -0500 Subject: [PATCH 4/9] fix geom term --- pyro/compressible/interface.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyro/compressible/interface.py b/pyro/compressible/interface.py index 4fd0a2bc6..de7dd9bf3 100644 --- a/pyro/compressible/interface.py +++ b/pyro/compressible/interface.py @@ -217,7 +217,7 @@ def states(idir, ng, dx, dloga, dt, # It's only there for non Cartesian coord. if idir == 1: - rho_source = -0.5 * dt * dloga[i,j] * q[iu] + rho_source = -0.5 * dt * dloga[i,j] * q[irho] * q[iu] q_l[i + 1, j, irho] += rho_source q_r[i, j, irho] += rho_source @@ -226,7 +226,7 @@ def states(idir, ng, dx, dloga, dt, q_r[i, j, ip] += rho_source * cs * cs else: - rho_source = -0.5 * dt * dloga[i,j] * q[iv] + rho_source = -0.5 * dt * dloga[i,j] * q[irho] * q[iv] q_l[i, j + 1, irho] += rho_source q_r[i, j, irho] += rho_source From 6fad483ceb75bd14a450d9d52b025bc396497341 Mon Sep 17 00:00:00 2001 From: zhi Date: Thu, 21 Nov 2024 17:20:13 -0500 Subject: [PATCH 5/9] Fix comment --- pyro/mesh/patch.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyro/mesh/patch.py b/pyro/mesh/patch.py index 42e011496..0d7aa4b6e 100644 --- a/pyro/mesh/patch.py +++ b/pyro/mesh/patch.py @@ -220,7 +220,7 @@ def __init__(self, nx, ny, *, ng=1, self.Ay = self.Lx - # Spatial derivative of log(Area) + # Spatial derivative of log(area). It's zero self.dlogAx = np.zeros_like(self.Ax) self.dlogAy = np.zeros_like(self.Ay) From 26334fa1fc70f0a3083c94285296c2b3032b375e Mon Sep 17 00:00:00 2001 From: zhi Date: Thu, 21 Nov 2024 17:22:39 -0500 Subject: [PATCH 6/9] comment --- pyro/mesh/patch.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyro/mesh/patch.py b/pyro/mesh/patch.py index 0d7aa4b6e..a79a5ecde 100644 --- a/pyro/mesh/patch.py +++ b/pyro/mesh/patch.py @@ -220,7 +220,7 @@ def __init__(self, nx, ny, *, ng=1, self.Ay = self.Lx - # Spatial derivative of log(area). It's zero + # Spatial derivative of log(area). It's zero for Cartesian self.dlogAx = np.zeros_like(self.Ax) self.dlogAy = np.zeros_like(self.Ay) From c237dc75e605a5c577bac07387b0b737c550bb93 Mon Sep 17 00:00:00 2001 From: zhi Date: Thu, 21 Nov 2024 17:24:09 -0500 Subject: [PATCH 7/9] fix flake8 --- pyro/mesh/tests/test_patch.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyro/mesh/tests/test_patch.py b/pyro/mesh/tests/test_patch.py index b3ffb3ac3..6fa5986af 100644 --- a/pyro/mesh/tests/test_patch.py +++ b/pyro/mesh/tests/test_patch.py @@ -145,7 +145,7 @@ def test_dlogAx(self): i = 4 r = self.g.xmin + (i + 0.5 - self.g.ng) * self.g.dx dlogAx = 2.0 / r - assert (self.g.dlogAx[i,:] == dlogAx).all() + assert (self.g.dlogAx[i, :] == dlogAx).all() def test_dlogAy(self): i = 4 @@ -154,7 +154,7 @@ def test_dlogAy(self): tan = np.tan(self.g.ymin + (j + 0.5 - self.g.ng) * self.g.dy) dlogAy = 1.0 / (r * tan) - assert self.g.dlogAy[i,j] == dlogAy + assert self.g.dlogAy[i, j] == dlogAy def test_V(self): volume = np.abs(-2.0 * np.pi / 3.0 * From 18271c435e90f3ec00fd76343c6aef27638d71d6 Mon Sep 17 00:00:00 2001 From: zhi Date: Thu, 21 Nov 2024 17:38:17 -0500 Subject: [PATCH 8/9] initialize dloga as arrayidnexer --- pyro/mesh/patch.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/pyro/mesh/patch.py b/pyro/mesh/patch.py index a79a5ecde..7c930b6b7 100644 --- a/pyro/mesh/patch.py +++ b/pyro/mesh/patch.py @@ -222,8 +222,10 @@ def __init__(self, nx, ny, *, ng=1, # Spatial derivative of log(area). It's zero for Cartesian - self.dlogAx = np.zeros_like(self.Ax) - self.dlogAy = np.zeros_like(self.Ay) + self.dlogAx = ArrayIndexer(np.zeros_like(self.Ax), + grid=self) + self.dlogAy = ArrayIndexer(np.zeros_like(self.Ay), + grid=self) # Volume of the cell. From 6e8038954cd5c3c68ff4cd228d04198741caa277 Mon Sep 17 00:00:00 2001 From: zhi Date: Thu, 21 Nov 2024 18:00:24 -0500 Subject: [PATCH 9/9] flake8 --- pyro/compressible/interface.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/pyro/compressible/interface.py b/pyro/compressible/interface.py index de7dd9bf3..df553661e 100644 --- a/pyro/compressible/interface.py +++ b/pyro/compressible/interface.py @@ -212,12 +212,11 @@ def states(idir, ng, dx, dloga, dt, q_l[i, j + 1, m] = q_l[i, j + 1, m] + sum_l q_r[i, j, m] = q_r[i, j, m] + sum_r - # Geometric Source term from converting conserved-variable to primitive # It's only there for non Cartesian coord. if idir == 1: - rho_source = -0.5 * dt * dloga[i,j] * q[irho] * q[iu] + rho_source = -0.5 * dt * dloga[i, j] * q[irho] * q[iu] q_l[i + 1, j, irho] += rho_source q_r[i, j, irho] += rho_source @@ -226,7 +225,7 @@ def states(idir, ng, dx, dloga, dt, q_r[i, j, ip] += rho_source * cs * cs else: - rho_source = -0.5 * dt * dloga[i,j] * q[irho] * q[iv] + rho_source = -0.5 * dt * dloga[i, j] * q[irho] * q[iv] q_l[i, j + 1, irho] += rho_source q_r[i, j, irho] += rho_source