diff --git a/.coveragerc b/.coveragerc index 52a1daba9..4b16222f8 100755 --- a/.coveragerc +++ b/.coveragerc @@ -1,8 +1,9 @@ [run] omit = examples/* test.py - incompressible/setup.py - compressible/setup.py - lm_atm/setup.py - advection_fv4/setup.py - plot.py \ No newline at end of file + plot.py + +[report] +omit = examples/* + test.py + plot.py diff --git a/.travis.yml b/.travis.yml index e7dd7010d..ffe2c9d0b 100755 --- a/.travis.yml +++ b/.travis.yml @@ -3,13 +3,6 @@ language: python python: - "3.6" -matrix: - include: - - env: - - MK_ARGS="" - - env: - - MK_ARGS=fortran - before_install: - export PATH=$(echo $PATH | tr ':' "\n" | sed '/\/opt\/python/d' | tr "\n" ":" | sed "s|::|:|g") @@ -26,7 +19,6 @@ install: before_script: - export PYTHONPATH=$PYTHONPATH:$(pwd) - export PYRO_HOME=$(pwd) - - ./mk.sh $MK_ARGS script: - flake8 . diff --git a/README.md b/README.md index 4005bae0d..8d5d26a5d 100644 --- a/README.md +++ b/README.md @@ -42,13 +42,13 @@ http://python-hydro.github.io/pyro2/ switch to python 3.x - There are a few steps to take to get things running. You need to - make sure you have `numpy`, `f2py`, `matplotlib`, and `h5py` + make sure you have `numpy`, `numba`, `matplotlib`, and `h5py` installed. On a Fedora system, this can be accomplished by doing: - `dnf install python3-numpy python3-numpy-f2py python3-matplotlib python3-matplotlib-tk python3-h5py` + `dnf install python3-numpy python3-numba python3-matplotlib python3-matplotlib-tk python3-h5py` (note, for older Fedora releases, replace `dnf` with `yum`. For - python 2.x, leave off the `2` in the package names.) + python 2.x, leave off the `3` in the package names.) - You also need to make sure gfortran is present on you system. On a Fedora system, it can be installed as: @@ -79,10 +79,6 @@ http://python-hydro.github.io/pyro2/ * Define the environment variable `PYRO_HOME` to point to the `pyro2/` directory (only needed for regression testing) - * Build the Fortran source. In `pyro2/` type - - `./mk.sh` - * Run a quick test of the advection solver: `./pyro.py advection smooth inputs.smooth` diff --git a/advection/tests/smooth_0040.h5 b/advection/tests/smooth_0040.h5 index ff14dd334..4e29e6e08 100644 Binary files a/advection/tests/smooth_0040.h5 and b/advection/tests/smooth_0040.h5 differ diff --git a/advection_fv4/fluxes.py b/advection_fv4/fluxes.py index a1b6af6b5..a86333ccb 100644 --- a/advection_fv4/fluxes.py +++ b/advection_fv4/fluxes.py @@ -1,4 +1,4 @@ -import advection_fv4.interface as interface_f +import advection_fv4.interface as interface import mesh.array_indexer as ai @@ -77,13 +77,13 @@ def fluxes(my_data, rp, dt): 1./12.*(a.jp(-2, buf=1) + a.jp(1, buf=1)) else: - a_l, a_r = interface_f.states(a, myg.qx, myg.qy, myg.ng, 1) + a_l, a_r = interface.states(a, myg.ng, 1) if u > 0: a_x = ai.ArrayIndexer(d=a_l, grid=myg) else: a_x = ai.ArrayIndexer(d=a_r, grid=myg) - a_l, a_r = interface_f.states(a, myg.qx, myg.qy, myg.ng, 2) + a_l, a_r = interface.states(a, myg.ng, 2) if v > 0: a_y = ai.ArrayIndexer(d=a_l, grid=myg) else: diff --git a/advection_fv4/interface.py b/advection_fv4/interface.py index 5f911a4b3..1cfcbd5d3 100644 --- a/advection_fv4/interface.py +++ b/advection_fv4/interface.py @@ -3,7 +3,7 @@ @njit(cache=True) -def states(a, qx, qy, ng, idir): +def states(a, ng, idir): r""" Predict the cell-centered state to the edges in one-dimension using the reconstructed, limited slopes. We use a fourth-order Godunov method. @@ -18,8 +18,6 @@ def states(a, qx, qy, ng, idir): ---------- a : ndarray The cell-centered state. - qx, qy : int - The dimensions of `a`. ng : int The number of ghost cells idir : int @@ -31,6 +29,8 @@ def states(a, qx, qy, ng, idir): The state predicted to the left and right edges. """ + qx, qy = a.shape + al = np.zeros((qx, qy)) ar = np.zeros((qx, qy)) @@ -250,8 +250,6 @@ def states_nolimit(a, qx, qy, ng, idir): ---------- a : ndarray The cell-centered state. - qx, qy : int - The dimensions of `a`. ng : int The number of ghost cells idir : int diff --git a/advection_fv4/interface_f.f90 b/advection_fv4/interface_f.f90 deleted file mode 100644 index d8efb8258..000000000 --- a/advection_fv4/interface_f.f90 +++ /dev/null @@ -1,317 +0,0 @@ -subroutine states(a, qx, qy, ng, idir, & - al, ar) - - implicit none - - integer, intent(in) :: idir - integer, intent(in) :: qx, qy, ng - - double precision, intent(inout) :: a(0:qx-1, 0:qy-1) - - double precision, intent(out) :: al(0:qx-1, 0:qy-1) - double precision, intent(out) :: ar(0:qx-1, 0:qy-1) - -!f2py depend(qx, qy) :: a -!f2py depend(qx, qy) :: al, ar -!f2py intent(in) :: a -!f2py intent(out) :: al, ar - - double precision :: a_int(0:qx-1, 0:qy-1) - double precision :: dafm(0:qx-1, 0:qy-1) - double precision :: dafp(0:qx-1, 0:qy-1) - double precision :: d2af(0:qx-1, 0:qy-1) - double precision :: d2ac(0:qx-1, 0:qy-1) - double precision :: d3a(0:qx-1, 0:qy-1) - - integer :: ilo, ihi, jlo, jhi - integer :: nx, ny - - double precision, parameter :: C2 = 1.25d0 - double precision, parameter :: C3 = 0.1d0 - - integer :: i, j - double precision :: rho, s - - double precision :: d2a_lim, d3a_min, d3a_max - - nx = qx - 2*ng; ny = qy - 2*ng - ilo = ng; ihi = ng+nx-1; jlo = ng; jhi = ng+ny-1 - - ! our convention here is that: - ! al(i,j) will be al_{i-1/2,j), - ! al(i+1,j) will be al_{i+1/2,j) - - ! we need interface values on all faces of the domain - if (idir == 1) then - - do j = jlo-1, jhi+1 - do i = ilo-2, ihi+3 - - ! interpolate to the edges - a_int(i,j) = (7.0d0/12.0d0)*(a(i-1,j) + a(i,j)) - (1.0d0/12.0d0)*(a(i-2,j) + a(i+1,j)) - - al(i,j) = a_int(i,j) - ar(i,j) = a_int(i,j) - - enddo - enddo - - do j = jlo-1, jhi+1 - do i = ilo-2, ihi+3 - ! these live on cell-centers - dafm(i,j) = a(i,j) - a_int(i,j) - dafp(i,j) = a_int(i+1,j) - a(i,j) - - ! these live on cell-centers - d2af(i,j) = 6.0d0*(a_int(i,j) - 2.0d0*a(i,j) + a_int(i+1,j)) - enddo - enddo - - do j = jlo-1, jhi+1 - do i = ilo-3, ihi+3 - d2ac(i,j) = a(i-1,j) - 2.0d0*a(i,j) + a(i+1,j) - enddo - enddo - - do j = jlo-1, jhi+1 - do i = ilo-2, ihi+3 - ! this lives on the interface - d3a(i,j) = d2ac(i,j) - d2ac(i-1,j) - enddo - enddo - - ! this is a look over cell centers, affecting - ! i-1/2,R and i+1/2,L - do j = jlo-1, jhi+1 - do i = ilo-1, ihi+1 - - ! limit? MC Eq. 24 and 25 - if (dafm(i,j) * dafp(i,j) <= 0.0d0 .or. & - (a(i,j) - a(i-2,j))*(a(i+2,j) - a(i,j)) <= 0.0d0) then - - ! we are at an extrema - - s = sign(1.0d0, d2ac(i,j)) - if (s == sign(1.0d0, d2ac(i-1,j)) .and. s == sign(1.0d0, d2ac(i+1,j)) .and. & - s == sign(1.0d0, d2af(i,j))) then - ! MC Eq. 26 - d2a_lim = s*min(abs(d2af(i,j)), C2*abs(d2ac(i-1,j)), & - C2*abs(d2ac(i,j)), C2*abs(d2ac(i+1,j))) - else - d2a_lim = 0.0d0 - endif - - if (abs(d2af(i,j)) <= 1.e-12*max(abs(a(i-2,j)), abs(a(i-1,j)), & - abs(a(i,j)), abs(a(i+1,j)), abs(a(i+2,j)))) then - rho = 0.0d0 - else - ! MC Eq. 27 - rho = d2a_lim/d2af(i,j) - endif - - if (rho < 1.0d0 - 1.d-12) then - ! we may need to limit -- these quantities are at cell-centers - d3a_min = min(d3a(i-1,j), d3a(i,j), d3a(i+1,j), d3a(i+2,j)) - d3a_max = max(d3a(i-1,j), d3a(i,j), d3a(i+1,j), d3a(i+2,j)) - - if (C3*max(abs(d3a_min), abs(d3a_max)) <= (d3a_max - d3a_min)) then - ! limit - if (dafm(i,j)*dafp(i,j) < 0.0d0) then - ! Eqs. 29, 30 - ar(i,j) = a(i,j) - rho*dafm(i,j) ! note: typo in Eq 29 - al(i+1,j) = a(i,j) + rho*dafp(i,j) - else if (abs(dafm(i,j)) >= 2.0d0*abs(dafp(i,j))) then - ! Eq. 31 - ar(i,j) = a(i,j) - 2.0d0*(1.0d0 - rho)*dafp(i,j) - rho*dafm(i,j) - else if (abs(dafp(i,j)) >= 2.0d0*abs(dafm(i,j))) then - ! Eq. 32 - al(i+1,j) = a(i,j) + 2.0d0*(1.0d0 - rho)*dafm(i,j) + rho*dafp(i,j) - endif - - endif - endif - - else - ! if Eqs. 24 or 25 didn't hold we still may need to limit - if (abs(dafm(i,j)) >= 2.0d0*abs(dafp(i,j))) then - ar(i,j) = a(i,j) - 2.0d0*dafp(i,j) - endif - if (abs(dafp(i,j)) >= 2.0d0*abs(dafm(i,j))) then - al(i+1,j) = a(i,j) + 2.0d0*dafm(i,j) - endif - endif - - enddo - enddo - - else if (idir == 2) then - - do j = jlo-2, jhi+3 - do i = ilo-1, ihi+1 - - ! interpolate to the edges - a_int(i,j) = (7.0d0/12.0d0)*(a(i,j-1) + a(i,j)) - (1.0d0/12.0d0)*(a(i,j-2) + a(i,j+1)) - - al(i,j) = a_int(i,j) - ar(i,j) = a_int(i,j) - - enddo - enddo - - do j = jlo-2, jhi+3 - do i = ilo-1, ihi+1 - ! these live on cell-centers - dafm(i,j) = a(i,j) - a_int(i,j) - dafp(i,j) = a_int(i,j+1) - a(i,j) - - ! these live on cell-centers - d2af(i,j) = 6.0d0*(a_int(i,j) - 2.0d0*a(i,j) + a_int(i,j+1)) - enddo - enddo - - do j = jlo-3, jhi+3 - do i = ilo-1, ihi+1 - d2ac(i,j) = a(i,j-1) - 2.0d0*a(i,j) + a(i,j+1) - enddo - enddo - - do j = jlo-2, jhi+3 - do i = ilo-1, ihi+1 - ! this lives on the interface - d3a(i,j) = d2ac(i,j) - d2ac(i,j-1) - enddo - enddo - - ! this is a look over cell centers, affecting - ! j-1/2,R and j+1/2,L - do j = jlo-1, jhi+1 - do i = ilo-1, ihi+1 - - ! limit? MC Eq. 24 and 25 - if (dafm(i,j) * dafp(i,j) <= 0.0d0 .or. & - (a(i,j) - a(i,j-2))*(a(i,j+2) - a(i,j)) <= 0.0d0) then - - ! we are at an extrema - - s = sign(1.0d0, d2ac(i,j)) - if (s == sign(1.0d0, d2ac(i,j-1)) .and. s == sign(1.0d0, d2ac(i,j+1)) .and. & - s == sign(1.0d0, d2af(i,j))) then - ! MC Eq. 26 - d2a_lim = s*min(abs(d2af(i,j)), C2*abs(d2ac(i,j-1)), & - C2*abs(d2ac(i,j)), C2*abs(d2ac(i,j+1))) - else - d2a_lim = 0.0d0 - endif - - if (abs(d2af(i,j)) <= 1.e-12*max(abs(a(i,j-2)), abs(a(i,j-1)), & - abs(a(i,j)), abs(a(i,j+1)), abs(a(i,j+2)))) then - rho = 0.0d0 - else - ! MC Eq. 27 - rho = d2a_lim/d2af(i,j) - endif - - if (rho < 1.0d0 - 1.d-12) then - ! we may need to limit -- these quantities are at cell-centers - d3a_min = min(d3a(i,j-1), d3a(i,j), d3a(i,j+1), d3a(i,j+2)) - d3a_max = max(d3a(i,j-1), d3a(i,j), d3a(i,j+1), d3a(i,j+2)) - - if (C3*max(abs(d3a_min), abs(d3a_max)) <= (d3a_max - d3a_min)) then - ! limit - if (dafm(i,j)*dafp(i,j) < 0.0d0) then - ! Eqs. 29, 30 - ar(i,j) = a(i,j) - rho*dafm(i,j) ! note: typo in Eq 29 - al(i,j+1) = a(i,j) + rho*dafp(i,j) - else if (abs(dafm(i,j)) >= 2.0d0*abs(dafp(i,j))) then - ! Eq. 31 - ar(i,j) = a(i,j) - 2.0d0*(1.0d0 - rho)*dafp(i,j) - rho*dafm(i,j) - else if (abs(dafp(i,j)) >= 2.0*abs(dafm(i,j))) then - ! Eq. 32 - al(i,j+1) = a(i,j) + 2.0d0*(1.0d0 - rho)*dafm(i,j) + rho*dafp(i,j) - endif - - endif - endif - - else - ! if Eqs. 24 or 25 didn't hold we still may need to limit - if (abs(dafm(i,j)) >= 2.0d0*abs(dafp(i,j))) then - ar(i,j) = a(i,j) - 2.0d0*dafp(i,j) - endif - if (abs(dafp(i,j)) >= 2.0d0*abs(dafm(i,j))) then - al(i,j+1) = a(i,j) + 2.0d0*dafm(i,j) - endif - endif - - enddo - enddo - - endif - -end subroutine states - - -subroutine states_nolimit(a, qx, qy, ng, idir, & - al, ar) - - implicit none - - integer, intent(in) :: idir - integer, intent(in) :: qx, qy, ng - - double precision, intent(inout) :: a(0:qx-1, 0:qy-1) - - double precision, intent(out) :: al(0:qx-1, 0:qy-1) - double precision, intent(out) :: ar(0:qx-1, 0:qy-1) - -!f2py depend(qx, qy) :: a -!f2py depend(qx, qy) :: al, ar -!f2py intent(in) :: a -!f2py intent(out) :: al, ar - - double precision :: a_int(0:qx-1, 0:qy-1) - - integer :: ilo, ihi, jlo, jhi - integer :: nx, ny - - - integer :: i, j - - nx = qx - 2*ng; ny = qy - 2*ng - ilo = ng; ihi = ng+nx-1; jlo = ng; jhi = ng+ny-1 - - ! our convention here is that: - ! al(i,j) will be al_{i-1/2,j), - ! al(i+1,j) will be al_{i+1/2,j) - - ! we need interface values on all faces of the domain - if (idir == 1) then - - do j = jlo-1, jhi+1 - do i = ilo-2, ihi+3 - - ! interpolate to the edges - a_int(i,j) = (7.0d0/12.0d0)*(a(i-1,j) + a(i,j)) - (1.0d0/12.0d0)*(a(i-2,j) + a(i+1,j)) - - al(i,j) = a_int(i,j) - ar(i,j) = a_int(i,j) - - enddo - enddo - - else if (idir == 2) then - - do j = jlo-2, jhi+3 - do i = ilo-1, ihi+1 - - ! interpolate to the edges - a_int(i,j) = (7.0d0/12.0d0)*(a(i,j-1) + a(i,j)) - (1.0d0/12.0d0)*(a(i,j-2) + a(i,j+1)) - - al(i,j) = a_int(i,j) - ar(i,j) = a_int(i,j) - - enddo - enddo - endif - -end subroutine states_nolimit diff --git a/advection_nonuniform/tests/slotted_0248.h5 b/advection_nonuniform/tests/slotted_0248.h5 index a34a8f015..b3c6613d1 100644 Binary files a/advection_nonuniform/tests/slotted_0248.h5 and b/advection_nonuniform/tests/slotted_0248.h5 differ diff --git a/advection_rk/tests/smooth_0081.h5 b/advection_rk/tests/smooth_0081.h5 index b41f0376c..6957986c0 100644 Binary files a/advection_rk/tests/smooth_0081.h5 and b/advection_rk/tests/smooth_0081.h5 differ diff --git a/advection_weno/__init__.py b/advection_weno/__init__.py index a06d1c51b..be169b672 100644 --- a/advection_weno/__init__.py +++ b/advection_weno/__init__.py @@ -9,21 +9,21 @@ The general flow of the solver when invoked through pyro.py is: - create grid +- create grid - initial conditions +- initial conditions - main loop +- main loop - fill ghost cells + * fill ghost cells - compute dt + * compute dt - compute fluxes + * compute fluxes - conservative update + * conservative update - output + * output """ diff --git a/advection_weno/fluxes.py b/advection_weno/fluxes.py index 2c6c57208..b89b78e9a 100644 --- a/advection_weno/fluxes.py +++ b/advection_weno/fluxes.py @@ -41,12 +41,13 @@ def fvs(q, order, u, alpha): def fluxes(my_data, rp, dt): - """ + r""" Construct the fluxes through the interfaces for the linear advection - equation: + equation + + .. math:: - a + u a + v a = 0 - t x y + a_t + u a_x + v a_y = 0 We use a high-order flux split WENO method to construct the interface fluxes. No Riemann problems are solved. The Lax-Friedrichs flux split diff --git a/compressible/interface.py b/compressible/interface.py index 12147d70c..87833ca2e 100644 --- a/compressible/interface.py +++ b/compressible/interface.py @@ -3,8 +3,8 @@ @njit(cache=True) -def states(idir, qx, qy, ng, dx, dt, - irho, iu, iv, ip, ix, nvar, nspec, +def states(idir, ng, dx, dt, + irho, iu, iv, ip, ix, nspec, gamma, qv, dqv): r""" predict the cell-centered state to the edges in one-dimension @@ -65,8 +65,6 @@ def states(idir, qx, qy, ng, dx, dt, ---------- idir : int Are we predicting to the edges in the x-direction (1) or y-direction (2)? - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells dx : float @@ -76,8 +74,6 @@ def states(idir, qx, qy, ng, dx, dt, irho, iu, iv, ip, ix : int Indices of the density, x-velocity, y-velocity, pressure and species in the state vector - nvar : int - The total number of variables in the state vector nspec : int The number of species gamma : float @@ -93,8 +89,10 @@ def states(idir, qx, qy, ng, dx, dt, State vector predicted to the left and right edges """ - q_l = np.zeros((qx, qy, nvar)) - q_r = np.zeros((qx, qy, nvar)) + qx, qy, nvar = qv.shape + + q_l = np.zeros_like(qv) + q_r = np.zeros_like(qv) nx = qx - 2 * ng ny = qy - 2 * ng @@ -217,8 +215,8 @@ def states(idir, qx, qy, ng, dx, dt, @njit(cache=True) -def riemann_cgf(idir, qx, qy, ng, - nvar, idens, ixmom, iymom, iener, irhoX, nspec, +def riemann_cgf(idir, ng, + idens, ixmom, iymom, iener, irhoX, nspec, lower_solid, upper_solid, gamma, U_l, U_r): r""" @@ -253,12 +251,8 @@ def riemann_cgf(idir, qx, qy, ng, ---------- idir : int Are we predicting to the edges in the x-direction (1) or y-direction (2)? - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells - nvar : int - The number of variables in the state vector nspec : int The number of species idens, ixmom, iymom, iener, irhoX : int @@ -277,6 +271,8 @@ def riemann_cgf(idir, qx, qy, ng, Conserved flux """ + qx, qy, nvar = U_l.shape + F = np.zeros((qx, qy, nvar)) smallc = 1.e-10 @@ -523,8 +519,8 @@ def riemann_cgf(idir, qx, qy, ng, @njit(cache=True) -def riemann_prim(idir, qx, qy, ng, - nvar, irho, iu, iv, ip, iX, nspec, +def riemann_prim(idir, ng, + irho, iu, iv, ip, iX, nspec, lower_solid, upper_solid, gamma, q_l, q_r): r""" @@ -563,12 +559,8 @@ def riemann_prim(idir, qx, qy, ng, ---------- idir : int Are we predicting to the edges in the x-direction (1) or y-direction (2)? - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells - nvar : int - The number of variables in the state vector nspec : int The number of species irho, iu, iv, ip, iX : int @@ -586,6 +578,8 @@ def riemann_prim(idir, qx, qy, ng, Primitive flux """ + qx, qy, nvar = q_l.shape + q_int = np.zeros((qx, qy, nvar)) smallc = 1.e-10 @@ -808,8 +802,8 @@ def riemann_prim(idir, qx, qy, ng, @njit(cache=True) -def riemann_hllc(idir, qx, qy, ng, - nvar, idens, ixmom, iymom, iener, irhoX, nspec, +def riemann_hllc(idir, ng, + idens, ixmom, iymom, iener, irhoX, nspec, lower_solid, upper_solid, gamma, U_l, U_r): r""" @@ -821,12 +815,8 @@ def riemann_hllc(idir, qx, qy, ng, ---------- idir : int Are we predicting to the edges in the x-direction (1) or y-direction (2)? - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells - nvar : int - The number of variables in the state vector nspec : int The number of species idens, ixmom, iymom, iener, irhoX : int @@ -845,6 +835,8 @@ def riemann_hllc(idir, qx, qy, ng, Conserved flux """ + qx, qy, nvar = U_l.shape + F = np.zeros((qx, qy, nvar)) smallc = 1.e-10 @@ -1001,7 +993,7 @@ def riemann_hllc(idir, qx, qy, ng, # R region U_state[:] = U_r[i, j, :] - F[i, j, :] = consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nvar, nspec, + F[i, j, :] = consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nspec, U_state) elif (S_r > 0.0 and S_c <= 0): @@ -1026,7 +1018,7 @@ def riemann_hllc(idir, qx, qy, ng, U_r[i, j, irhoX:irhoX + nspec] / rho_r # find the flux on the right interface - F[i, j, :] = consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nvar, nspec, + F[i, j, :] = consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nspec, U_r[i, j, :]) # correct the flux @@ -1054,7 +1046,7 @@ def riemann_hllc(idir, qx, qy, ng, U_l[i, j, irhoX:irhoX + nspec] / rho_l # find the flux on the left interface - F[i, j, :] = consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nvar, nspec, + F[i, j, :] = consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nspec, U_l[i, j, :]) # correct the flux @@ -1064,7 +1056,7 @@ def riemann_hllc(idir, qx, qy, ng, # L region U_state[:] = U_l[i, j, :] - F[i, j, :] = consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nvar, nspec, + F[i, j, :] = consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nspec, U_state) # we should deal with solid boundaries somehow here @@ -1073,7 +1065,7 @@ def riemann_hllc(idir, qx, qy, ng, @njit(cache=True) -def consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nvar, nspec, U_state): +def consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nspec, U_state): r""" Calculate the conservative flux. @@ -1086,8 +1078,6 @@ def consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nvar, nspec, U_stat idens, ixmom, iymom, iener, irhoX : int The indices of the density, x-momentum, y-momentum, internal energy density and species partial densities in the conserved state vector. - nvar : int - The number of variables in the state vector nspec : int The number of species U_state : ndarray @@ -1099,7 +1089,7 @@ def consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nvar, nspec, U_stat Conserved flux """ - F = np.zeros(nvar) + F = np.zeros_like(U_state) u = U_state[ixmom] / U_state[idens] v = U_state[iymom] / U_state[idens] @@ -1128,7 +1118,7 @@ def consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nvar, nspec, U_stat @njit(cache=True) -def artificial_viscosity(qx, qy, ng, dx, dy, +def artificial_viscosity(ng, dx, dy, cvisc, u, v): r""" Compute the artifical viscosity. Here, we compute edge-centered @@ -1159,8 +1149,6 @@ def artificial_viscosity(qx, qy, ng, dx, dy, Parameters ---------- - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells dx, dy : float @@ -1176,6 +1164,8 @@ def artificial_viscosity(qx, qy, ng, dx, dy, Artificial viscosity in the x- and y-directions """ + qx, qy = u.shape + avisco_x = np.zeros((qx, qy)) avisco_y = np.zeros((qx, qy)) diff --git a/compressible/interface_f.f90 b/compressible/interface_f.f90 deleted file mode 100644 index 5dabf97ca..000000000 --- a/compressible/interface_f.f90 +++ /dev/null @@ -1,1240 +0,0 @@ -subroutine states(idir, qx, qy, ng, dx, dt, & - irho, iu, iv, ip, ix, nvar, nspec, & - gamma, & - qv, dqv, & - q_l, q_r) - - implicit none - - integer, intent(in) :: idir - integer, intent(in) :: qx, qy, ng - double precision, intent(in) :: dx, dt - integer, intent(in) :: irho, iu, iv, ip, ix, nvar, nspec - double precision, intent(in) :: gamma - - ! 0-based indexing to match python - double precision, intent(inout) :: qv(0:qx-1, 0:qy-1, 0:nvar-1) - double precision, intent(inout) :: dqv(0:qx-1, 0:qy-1, 0:nvar-1) - - double precision, intent( out) :: q_l(0:qx-1, 0:qy-1, 0:nvar-1) - double precision, intent( out) :: q_r(0:qx-1, 0:qy-1, 0:nvar-1) - - !f2py depend(qx, qy, nvar) :: qv, dqv - !f2py depend(qx, qy, nvar) :: q_l, q_r - !f2py intent(in) :: qv, dqv - !f2py intent(out) :: q_l, q_r - - ! predict the cell-centered state to the edges in one-dimension - ! using the reconstructed, limited slopes. - ! - ! We follow the convection here that V_l[i] is the left state at the - ! i-1/2 interface and V_l[i+1] is the left state at the i+1/2 - ! interface. - ! - ! - ! We need the left and right eigenvectors and the eigenvalues for - ! the system projected along the x-direction - ! - ! Taking our state vector as Q = (rho, u, v, p)^T, the eigenvalues - ! are u - c, u, u + c. - ! - ! We look at the equations of hydrodynamics in a split fashion -- - ! i.e., we only consider one dimension at a time. - ! - ! Considering advection in the x-direction, the Jacobian matrix for - ! the primitive variable formulation of the Euler equations - ! projected in the x-direction is: - ! - ! / u r 0 0 \ - ! | 0 u 0 1/r | - ! A = | 0 0 u 0 | - ! \ 0 rc^2 0 u / - ! - ! The right eigenvectors are - ! - ! / 1 \ / 1 \ / 0 \ / 1 \ - ! |-c/r | | 0 | | 0 | | c/r | - ! r1 = | 0 | r2 = | 0 | r3 = | 1 | r4 = | 0 | - ! \ c^2 / \ 0 / \ 0 / \ c^2 / - ! - ! In particular, we see from r3 that the transverse velocity (v in - ! this case) is simply advected at a speed u in the x-direction. - ! - ! The left eigenvectors are - ! - ! l1 = ( 0, -r/(2c), 0, 1/(2c^2) ) - ! l2 = ( 1, 0, 0, -1/c^2 ) - ! l3 = ( 0, 0, 1, 0 ) - ! l4 = ( 0, r/(2c), 0, 1/(2c^2) ) - ! - ! The fluxes are going to be defined on the left edge of the - ! computational zones - ! - ! | | | | - ! | | | | - ! -+------+------+------+------+------+------+-- - ! | i-1 | i | i+1 | - ! ^ ^ ^ - ! q_l,i q_r,i q_l,i+1 - ! - ! q_r,i and q_l,i+1 are computed using the information in zone i,j. - - integer :: ilo, ihi, jlo, jhi - integer :: nx, ny - integer :: i, j, n, m - - double precision :: dq(0:nvar-1), q(0:nvar-1) - double precision :: lvec(0:nvar-1,0:nvar-1), rvec(0:nvar-1,0:nvar-1) - double precision :: eval(0:nvar-1) - double precision :: betal(0:nvar-1), betar(0:nvar-1) - - double precision :: dtdx, dtdx4 - double precision :: cs - - double precision :: sum, sum_l, sum_r, factor - - integer :: ns - - nx = qx - 2*ng; ny = qy - 2*ng - ilo = ng; ihi = ng+nx-1; jlo = ng; jhi = ng+ny-1 - - ns = nvar - nspec - - dtdx = dt/dx - dtdx4 = 0.25d0*dtdx - - ! this is the loop over zones. For zone i, we see q_l[i+1] and q_r[i] - do j = jlo-2, jhi+2 - do i = ilo-2, ihi+2 - - dq(:) = dqv(i,j,:) - q(:) = qv(i,j,:) - - cs = sqrt(gamma*q(ip)/q(irho)) - - lvec(:,:) = 0.0d0 - rvec(:,:) = 0.0d0 - eval(:) = 0.0d0 - - ! compute the eigenvalues and eigenvectors - if (idir == 1) then - eval(:) = [q(iu) - cs, q(iu), q(iu), q(iu) + cs] - - lvec(0,0:ns-1) = [ 0.0d0, -0.5d0*q(irho)/cs, 0.0d0, 0.5d0/(cs*cs) ] - lvec(1,0:ns-1) = [ 1.0d0, 0.0d0, 0.0d0, -1.0d0/(cs*cs) ] - lvec(2,0:ns-1) = [ 0.0d0, 0.0d0, 1.0d0, 0.0d0 ] - lvec(3,0:ns-1) = [ 0.0d0, 0.5d0*q(irho)/cs, 0.0d0, 0.5d0/(cs*cs) ] - - rvec(0,0:ns-1) = [1.0d0, -cs/q(irho), 0.0d0, cs*cs ] - rvec(1,0:ns-1) = [1.0d0, 0.0d0, 0.0d0, 0.0d0 ] - rvec(2,0:ns-1) = [0.0d0, 0.0d0, 1.0d0, 0.0d0 ] - rvec(3,0:ns-1) = [1.0d0, cs/q(irho), 0.0d0, cs*cs ] - - ! now the species -- they only have a 1 in their corresponding slot - eval(ns:) = q(iu) - do n = ix, ix-1+nspec - lvec(n,n) = 1.0 - rvec(n,n) = 1.0 - enddo - else - eval = [q(iv) - cs, q(iv), q(iv), q(iv) + cs] - - lvec(0,0:ns-1) = [ 0.0d0, 0.0d0, -0.5d0*q(irho)/cs, 0.5d0/(cs*cs) ] - lvec(1,0:ns-1) = [ 1.0d0, 0.0d0, 0.0d0, -1.0d0/(cs*cs) ] - lvec(2,0:ns-1) = [ 0.0d0, 1.0d0, 0.0d0, 0.0d0 ] - lvec(3,0:ns-1) = [ 0.0d0, 0.0d0, 0.5d0*q(irho)/cs, 0.5d0/(cs*cs) ] - - rvec(0,0:ns-1) = [1.0d0, 0.0d0, -cs/q(irho), cs*cs ] - rvec(1,0:ns-1) = [1.0d0, 0.0d0, 0.0d0, 0.0d0 ] - rvec(2,0:ns-1) = [0.0d0, 1.0d0, 0.0d0, 0.0d0 ] - rvec(3,0:ns-1) = [1.0d0, 0.0d0, cs/q(irho), cs*cs ] - - ! now the species -- they only have a 1 in their corresponding slot - eval(ns:) = q(iv) - do n = ix, ix-1+nspec - lvec(n,n) = 1.0 - rvec(n,n) = 1.0 - enddo - - endif - - ! define the reference states - if (idir == 1) then - ! this is one the right face of the current zone, - ! so the fastest moving eigenvalue is eval[3] = u + c - factor = 0.5d0*(1.0d0 - dtdx*max(eval(3), 0.0d0)) - q_l(i+1,j,:) = q(:) + factor*dq(:) - - ! left face of the current zone, so the fastest moving - ! eigenvalue is eval[3] = u - c - factor = 0.5d0*(1.0d0 + dtdx*min(eval(0), 0.0d0)) - q_r(i, j,:) = q(:) - factor*dq(:) - - else - - factor = 0.5d0*(1.0d0 - dtdx*max(eval(3), 0.0d0)) - q_l(i,j+1,:) = q(:) + factor*dq(:) - - factor = 0.5d0*(1.0d0 + dtdx*min(eval(0), 0.0d0)) - q_r(i,j, :) = q(:) - factor*dq(:) - - endif - - ! compute the Vhat functions - do m = 0, nvar-1 - sum = dot_product(lvec(m,:),dq(:)) - - betal(m) = dtdx4*(eval(3) - eval(m))*(sign(1.0d0,eval(m)) + 1.0d0)*sum - betar(m) = dtdx4*(eval(0) - eval(m))*(1.0d0 - sign(1.0d0,eval(m)))*sum - enddo - - ! construct the states - do m = 0, nvar-1 - sum_l = dot_product(betal(:),rvec(:,m)) - sum_r = dot_product(betar(:),rvec(:,m)) - - if (idir == 1) then - q_l(i+1,j,m) = q_l(i+1,j,m) + sum_l - q_r(i, j,m) = q_r(i, j,m) + sum_r - else - 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 - endif - - enddo - - enddo - enddo - -end subroutine states - - -subroutine riemann_cgf(idir, qx, qy, ng, & - nvar, idens, ixmom, iymom, iener, irhoX, nspec, & - lower_solid, upper_solid, & - gamma, U_l, U_r, F) - - implicit none - - integer, intent(in) :: idir - integer, intent(in) :: qx, qy, ng - integer, intent(in) :: nvar, idens, ixmom, iymom, iener, irhoX, nspec - integer, intent(in) :: lower_solid, upper_solid - double precision, intent(in) :: gamma - - ! 0-based indexing to match python - double precision, intent(inout) :: U_l(0:qx-1,0:qy-1,0:nvar-1) - double precision, intent(inout) :: U_r(0:qx-1,0:qy-1,0:nvar-1) - double precision, intent( out) :: F(0:qx-1,0:qy-1,0:nvar-1) - -!f2py depend(qx, qy, nvar) :: U_l, U_r, F -!f2py intent(in) :: U_l, U_r -!f2py intent(out) :: F - - ! Solve riemann shock tube problem for a general equation of - ! state using the method of Colella, Glaz, and Ferguson. See - ! Almgren et al. 2010 (the CASTRO paper) for details. - ! - ! The Riemann problem for the Euler's equation produces 4 regions, - ! separated by the three characteristics (u - cs, u, u + cs): - ! - ! - ! u - cs t u u + cs - ! \ ^ . / - ! \ *L | . *R / - ! \ | . / - ! \ | . / - ! L \ | . / R - ! \ | . / - ! \ |. / - ! \|./ - ! ----------+----------------> x - ! - ! We care about the solution on the axis. The basic idea is to use - ! estimates of the wave speeds to figure out which region we are in, - ! and then use jump conditions to evaluate the state there. - ! - ! Only density jumps across the u characteristic. All primitive - ! variables jump across the other two. Special attention is needed - ! if a rarefaction spans the axis. - - integer :: ilo, ihi, jlo, jhi - integer :: nx, ny - integer :: i, j - - double precision, parameter :: smallc = 1.e-10 - double precision, parameter :: smallrho = 1.e-10 - double precision, parameter :: smallp = 1.e-10 - - double precision :: rho_l, un_l, ut_l, p_l, rhoe_l - double precision :: rho_r, un_r, ut_r, p_r, rhoe_r - double precision :: xn(nspec) - double precision :: rhostar_l, rhostar_r, rhoestar_l, rhoestar_r - double precision :: ustar, pstar, cstar_l, cstar_r - double precision :: lambda_l, lambdastar_l, lambda_r, lambdastar_r - double precision :: W_l, W_r, c_l, c_r, sigma - double precision :: alpha - - double precision :: rho_state, un_state, ut_state, p_state, rhoe_state - - - nx = qx - 2*ng; ny = qy - 2*ng - ilo = ng; ihi = ng+nx-1; jlo = ng; jhi = ng+ny-1 - - do j = jlo-1, jhi+1 - do i = ilo-1, ihi+1 - - ! primitive variable states - rho_l = U_l(i,j,idens) - - ! un = normal velocity; ut = transverse velocity - if (idir == 1) then - un_l = U_l(i,j,ixmom)/rho_l - ut_l = U_l(i,j,iymom)/rho_l - else - un_l = U_l(i,j,iymom)/rho_l - ut_l = U_l(i,j,ixmom)/rho_l - endif - - rhoe_l = U_l(i,j,iener) - 0.5*rho_l*(un_l**2 + ut_l**2) - - p_l = rhoe_l*(gamma - 1.0d0) - p_l = max(p_l, smallp) - - rho_r = U_r(i,j,idens) - - if (idir == 1) then - un_r = U_r(i,j,ixmom)/rho_r - ut_r = U_r(i,j,iymom)/rho_r - else - un_r = U_r(i,j,iymom)/rho_r - ut_r = U_r(i,j,ixmom)/rho_r - endif - - rhoe_r = U_r(i,j,iener) - 0.5*rho_r*(un_r**2 + ut_r**2) - - p_r = rhoe_r*(gamma - 1.0d0) - p_r = max(p_r, smallp) - - - ! define the Lagrangian sound speed - W_l = max(smallrho*smallc, sqrt(gamma*p_l*rho_l)) - W_r = max(smallrho*smallc, sqrt(gamma*p_r*rho_r)) - - ! and the regular sound speeds - c_l = max(smallc, sqrt(gamma*p_l/rho_l)) - c_r = max(smallc, sqrt(gamma*p_r/rho_r)) - - ! define the star states - pstar = (W_l*p_r + W_r*p_l + W_l*W_r*(un_l - un_r))/(W_l + W_r) - pstar = max(pstar, smallp) - ustar = (W_l*un_l + W_r*un_r + (p_l - p_r))/(W_l + W_r) - - ! now compute the remaining state to the left and right - ! of the contact (in the star region) - rhostar_l = rho_l + (pstar - p_l)/c_l**2 - rhostar_r = rho_r + (pstar - p_r)/c_r**2 - - rhoestar_l = rhoe_l + & - (pstar - p_l)*(rhoe_l/rho_l + p_l/rho_l)/c_l**2 - rhoestar_r = rhoe_r + & - (pstar - p_r)*(rhoe_r/rho_r + p_r/rho_r)/c_r**2 - - cstar_l = max(smallc,sqrt(gamma*pstar/rhostar_l)) - cstar_r = max(smallc,sqrt(gamma*pstar/rhostar_r)) - - ! figure out which state we are in, based on the location of - ! the waves - if (ustar > 0.0d0) then - - ! contact is moving to the right, we need to understand - ! the L and *L states - - ! Note: transverse velocity only jumps across contact - ut_state = ut_l - - ! define eigenvalues - lambda_l = un_l - c_l - lambdastar_l = ustar - cstar_l - - if (pstar > p_l) then - ! the wave is a shock -- find the shock speed - sigma = (lambda_l + lambdastar_l)/2.0d0 - - if (sigma > 0.0d0) then - ! shock is moving to the right -- solution is L state - rho_state = rho_l - un_state = un_l - p_state = p_l - rhoe_state = rhoe_l - - else - ! solution is *L state - rho_state = rhostar_l - un_state = ustar - p_state = pstar - rhoe_state = rhoestar_l - endif - - else - ! the wave is a rarefaction - if (lambda_l < 0.0d0 .and. lambdastar_l < 0.0d0) then - ! rarefaction fan is moving to the left -- solution is - ! *L state - rho_state = rhostar_l - un_state = ustar - p_state = pstar - rhoe_state = rhoestar_l - - else if (lambda_l > 0.0d0 .and. lambdastar_l > 0.0d0) then - ! rarefaction fan is moving to the right -- solution is - ! L state - rho_state = rho_l - un_state = un_l - p_state = p_l - rhoe_state = rhoe_l - - else - ! rarefaction spans x/t = 0 -- interpolate - alpha = lambda_l/(lambda_l - lambdastar_l) - - rho_state = alpha*rhostar_l + (1.0d0 - alpha)*rho_l - un_state = alpha*ustar + (1.0d0 - alpha)*un_l - p_state = alpha*pstar + (1.0d0 - alpha)*p_l - rhoe_state = alpha*rhoestar_l + (1.0d0 - alpha)*rhoe_l - endif - - endif - - else if (ustar < 0) then - - ! contact moving left, we need to understand the R and *R - ! states - - ! Note: transverse velocity only jumps across contact - ut_state = ut_r - - ! define eigenvalues - lambda_r = un_r + c_r - lambdastar_r = ustar + cstar_r - - if (pstar > p_r) then - ! the wave if a shock -- find the shock speed - sigma = (lambda_r + lambdastar_r)/2.0d0 - - if (sigma > 0.0d0) then - ! shock is moving to the right -- solution is *R state - rho_state = rhostar_r - un_state = ustar - p_state = pstar - rhoe_state = rhoestar_r - - else - ! solution is R state - rho_state = rho_r - un_state = un_r - p_state = p_r - rhoe_state = rhoe_r - endif - - else - ! the wave is a rarefaction - if (lambda_r < 0.0d0 .and. lambdastar_r < 0.0d0) then - ! rarefaction fan is moving to the left -- solution is - ! R state - rho_state = rho_r - un_state = un_r - p_state = p_r - rhoe_state = rhoe_r - - else if (lambda_r > 0.0d0 .and. lambdastar_r > 0.0d0) then - ! rarefaction fan is moving to the right -- solution is - ! *R state - rho_state = rhostar_r - un_state = ustar - p_state = pstar - rhoe_state = rhoestar_r - - else - ! rarefaction spans x/t = 0 -- interpolate - alpha = lambda_r/(lambda_r - lambdastar_r) - - rho_state = alpha*rhostar_r + (1.0d0 - alpha)*rho_r - un_state = alpha*ustar + (1.0d0 - alpha)*un_r - p_state = alpha*pstar + (1.0d0 - alpha)*p_r - rhoe_state = alpha*rhoestar_r + (1.0d0 - alpha)*rhoe_r - - endif - - endif - - else ! ustar == 0 - - rho_state = 0.5*(rhostar_l + rhostar_r) - un_state = ustar - ut_state = 0.5*(ut_l + ut_r) - p_state = pstar - rhoe_state = 0.5*(rhoestar_l + rhoestar_r) - - endif - - ! species now - if (nspec > 0) then - if (ustar > 0.0) then - xn(:) = U_l(i,j,irhoX:irhoX-1+nspec)/U_l(i,j,idens) - - else if (ustar < 0.0) then - xn(:) = U_r(i,j,irhoX:irhoX-1+nspec)/U_r(i,j,idens) - else - xn(:) = 0.5d0*(U_l(i,j,irhoX:irhoX-1+nspec)/U_l(i,j,idens) + & - U_r(i,j,irhoX:irhoX-1+nspec)/U_r(i,j,idens)) - endif - endif - - ! are we on a solid boundary? - if (idir == 1) then - if (i == ilo .and. lower_solid == 1) then - un_state = 0.0 - endif - - if (i == ihi+1 .and. upper_solid == 1) then - un_state = 0.0 - endif - - else if (idir == 2) then - if (j == jlo .and. lower_solid == 1) then - un_state = 0.0 - endif - - if (j == jhi+1 .and. upper_solid == 1) then - un_state = 0.0 - endif - - endif - - ! compute the fluxes - F(i,j,idens) = rho_state*un_state - - if (idir == 1) then - F(i,j,ixmom) = rho_state*un_state**2 + p_state - F(i,j,iymom) = rho_state*ut_state*un_state - else - F(i,j,ixmom) = rho_state*ut_state*un_state - F(i,j,iymom) = rho_state*un_state**2 + p_state - endif - - F(i,j,iener) = rhoe_state*un_state + & - 0.5*rho_state*(un_state**2 + ut_state**2)*un_state + & - p_state*un_state - - if (nspec > 0) then - F(i,j,irhoX:irhoX-1+nspec) = xn(:)*rho_state*un_state - endif - - enddo - enddo - -end subroutine riemann_cgf - - -subroutine riemann_prim(idir, qx, qy, ng, & - nvar, irho, iu, iv, ip, iX, nspec, & - lower_solid, upper_solid, & - gamma, q_l, q_r, q_int) - - ! this is like riemann_cgf, except that it works on a primitive - ! variable input state and returns the primitive variable interface - ! state - - implicit none - - integer, intent(in) :: idir - integer, intent(in) :: qx, qy, ng - integer, intent(in) :: nvar, irho, iu, iv, ip, iX, nspec - integer, intent(in) :: lower_solid, upper_solid - double precision, intent(in) :: gamma - - ! 0-based indexing to match python - double precision, intent(inout) :: q_l(0:qx-1,0:qy-1,0:nvar-1) - double precision, intent(inout) :: q_r(0:qx-1,0:qy-1,0:nvar-1) - double precision, intent( out) :: q_int(0:qx-1,0:qy-1,0:nvar-1) - -!f2py depend(qx, qy, nvar) :: q_l, q_r, q_int -!f2py intent(in) :: q_l, q_r -!f2py intent(out) :: q_int - - ! Solve riemann shock tube problem for a general equation of - ! state using the method of Colella, Glaz, and Ferguson. See - ! Almgren et al. 2010 (the CASTRO paper) for details. - ! - ! The Riemann problem for the Euler's equation produces 4 regions, - ! separated by the three characteristics (u - cs, u, u + cs): - ! - ! - ! u - cs t u u + cs - ! \ ^ . / - ! \ *L | . *R / - ! \ | . / - ! \ | . / - ! L \ | . / R - ! \ | . / - ! \ |. / - ! \|./ - ! ----------+----------------> x - ! - ! We care about the solution on the axis. The basic idea is to use - ! estimates of the wave speeds to figure out which region we are in, - ! and then use jump conditions to evaluate the state there. - ! - ! Only density jumps across the u characteristic. All primitive - ! variables jump across the other two. Special attention is needed - ! if a rarefaction spans the axis. - - integer :: ilo, ihi, jlo, jhi - integer :: nx, ny - integer :: i, j - - double precision, parameter :: smallc = 1.e-10 - double precision, parameter :: smallrho = 1.e-10 - double precision, parameter :: smallp = 1.e-10 - - double precision :: rho_l, un_l, ut_l, p_l - double precision :: rho_r, un_r, ut_r, p_r - double precision :: xn(nspec) - double precision :: rhostar_l, rhostar_r - double precision :: ustar, pstar, cstar_l, cstar_r - double precision :: lambda_l, lambdastar_l, lambda_r, lambdastar_r - double precision :: W_l, W_r, c_l, c_r, sigma - double precision :: alpha - - double precision :: rho_state, un_state, ut_state, p_state - - nx = qx - 2*ng; ny = qy - 2*ng - ilo = ng; ihi = ng+nx-1; jlo = ng; jhi = ng+ny-1 - - do j = jlo-1, jhi+1 - do i = ilo-1, ihi+1 - - ! primitive variable states - rho_l = q_l(i,j,irho) - - ! un = normal velocity; ut = transverse velocity - if (idir == 1) then - un_l = q_l(i,j,iu) - ut_l = q_l(i,j,iv) - else - un_l = q_l(i,j,iv) - ut_l = q_l(i,j,iu) - endif - - p_l = q_l(i,j,ip) - p_l = max(p_l, smallp) - - rho_r = q_r(i,j,irho) - - if (idir == 1) then - un_r = q_r(i,j,iu) - ut_r = q_r(i,j,iv) - else - un_r = q_r(i,j,iv) - ut_r = q_r(i,j,iu) - endif - - p_r = q_r(i,j,ip) - p_r = max(p_r, smallp) - - - ! define the Lagrangian sound speed - W_l = max(smallrho*smallc, sqrt(gamma*p_l*rho_l)) - W_r = max(smallrho*smallc, sqrt(gamma*p_r*rho_r)) - - ! and the regular sound speeds - c_l = max(smallc, sqrt(gamma*p_l/rho_l)) - c_r = max(smallc, sqrt(gamma*p_r/rho_r)) - - ! define the star states - pstar = (W_l*p_r + W_r*p_l + W_l*W_r*(un_l - un_r))/(W_l + W_r) - pstar = max(pstar, smallp) - ustar = (W_l*un_l + W_r*un_r + (p_l - p_r))/(W_l + W_r) - - ! now compute the remaining state to the left and right - ! of the contact (in the star region) - rhostar_l = rho_l + (pstar - p_l)/c_l**2 - rhostar_r = rho_r + (pstar - p_r)/c_r**2 - - cstar_l = max(smallc,sqrt(gamma*pstar/rhostar_l)) - cstar_r = max(smallc,sqrt(gamma*pstar/rhostar_r)) - - ! figure out which state we are in, based on the location of - ! the waves - if (ustar > 0.0d0) then - - ! contact is moving to the right, we need to understand - ! the L and *L states - - ! Note: transverse velocity only jumps across contact - ut_state = ut_l - - ! define eigenvalues - lambda_l = un_l - c_l - lambdastar_l = ustar - cstar_l - - if (pstar > p_l) then - ! the wave is a shock -- find the shock speed - sigma = (lambda_l + lambdastar_l)/2.0d0 - - if (sigma > 0.0d0) then - ! shock is moving to the right -- solution is L state - rho_state = rho_l - un_state = un_l - p_state = p_l - - else - ! solution is *L state - rho_state = rhostar_l - un_state = ustar - p_state = pstar - endif - - else - ! the wave is a rarefaction - if (lambda_l < 0.0d0 .and. lambdastar_l < 0.0d0) then - ! rarefaction fan is moving to the left -- solution is - ! *L state - rho_state = rhostar_l - un_state = ustar - p_state = pstar - - else if (lambda_l > 0.0d0 .and. lambdastar_l > 0.0d0) then - ! rarefaction fan is moving to the right -- solution is - ! L state - rho_state = rho_l - un_state = un_l - p_state = p_l - - else - ! rarefaction spans x/t = 0 -- interpolate - alpha = lambda_l/(lambda_l - lambdastar_l) - - rho_state = alpha*rhostar_l + (1.0d0 - alpha)*rho_l - un_state = alpha*ustar + (1.0d0 - alpha)*un_l - p_state = alpha*pstar + (1.0d0 - alpha)*p_l - endif - - endif - - else if (ustar < 0) then - - ! contact moving left, we need to understand the R and *R - ! states - - ! Note: transverse velocity only jumps across contact - ut_state = ut_r - - ! define eigenvalues - lambda_r = un_r + c_r - lambdastar_r = ustar + cstar_r - - if (pstar > p_r) then - ! the wave if a shock -- find the shock speed - sigma = (lambda_r + lambdastar_r)/2.0d0 - - if (sigma > 0.0d0) then - ! shock is moving to the right -- solution is *R state - rho_state = rhostar_r - un_state = ustar - p_state = pstar - - else - ! solution is R state - rho_state = rho_r - un_state = un_r - p_state = p_r - endif - - else - ! the wave is a rarefaction - if (lambda_r < 0.0d0 .and. lambdastar_r < 0.0d0) then - ! rarefaction fan is moving to the left -- solution is - ! R state - rho_state = rho_r - un_state = un_r - p_state = p_r - - else if (lambda_r > 0.0d0 .and. lambdastar_r > 0.0d0) then - ! rarefaction fan is moving to the right -- solution is - ! *R state - rho_state = rhostar_r - un_state = ustar - p_state = pstar - - else - ! rarefaction spans x/t = 0 -- interpolate - alpha = lambda_r/(lambda_r - lambdastar_r) - - rho_state = alpha*rhostar_r + (1.0d0 - alpha)*rho_r - un_state = alpha*ustar + (1.0d0 - alpha)*un_r - p_state = alpha*pstar + (1.0d0 - alpha)*p_r - - endif - - endif - - else ! ustar == 0 - - rho_state = 0.5*(rhostar_l + rhostar_r) - un_state = ustar - ut_state = 0.5*(ut_l + ut_r) - p_state = pstar - - endif - - ! species now - if (nspec > 0) then - if (ustar > 0.0) then - xn(:) = q_l(i,j,iX:iX-1+nspec) - - else if (ustar < 0.0) then - xn(:) = q_r(i,j,iX:iX-1+nspec) - else - xn(:) = 0.5d0*(q_l(i,j,iX:iX-1+nspec) + q_r(i,j,iX:iX-1+nspec)) - endif - endif - - ! are we on a solid boundary? - if (idir == 1) then - if (i == ilo .and. lower_solid == 1) then - un_state = 0.0 - endif - - if (i == ihi+1 .and. upper_solid == 1) then - un_state = 0.0 - endif - - else if (idir == 2) then - if (j == jlo .and. lower_solid == 1) then - un_state = 0.0 - endif - - if (j == jhi+1 .and. upper_solid == 1) then - un_state = 0.0 - endif - - endif - - q_int(i,j,irho) = rho_state - if (idir == 1) then - q_int(i,j,iu) = un_state - q_int(i,j,iv) = ut_state - else - q_int(i,j,iu) = ut_state - q_int(i,j,iv) = un_state - endif - q_int(i,j,ip) = p_state - - if (nspec > 0) then - q_int(i,j,iX:iX-1+nspec) = xn(:) - endif - - enddo - enddo - -end subroutine riemann_prim - - -subroutine riemann_HLLC(idir, qx, qy, ng, & - nvar, idens, ixmom, iymom, iener, irhoX, nspec, & - lower_solid, upper_solid, & - gamma, U_l, U_r, F) - - - implicit none - - integer, intent(in) :: idir - integer, intent(in) :: qx, qy, ng - integer, intent(in) :: nvar, idens, ixmom, iymom, iener, irhoX, nspec - integer, intent(in) :: lower_solid, upper_solid - double precision, intent(in) :: gamma - - ! 0-based indexing to match python - double precision, intent(inout) :: U_l(0:qx-1,0:qy-1,0:nvar-1) - double precision, intent(inout) :: U_r(0:qx-1,0:qy-1,0:nvar-1) - double precision, intent( out) :: F(0:qx-1,0:qy-1,0:nvar-1) - -!f2py depend(qx, qy, nvar) :: U_l, U_r -!f2py intent(in) :: U_l, U_r -!f2py intent(out) :: F - - ! this is the HLLC Riemann solver. The implementation follows - ! directly out of Toro's book. Note: this does not handle the - ! transonic rarefaction. - - integer :: ilo, ihi, jlo, jhi - integer :: nx, ny - integer :: i, j - - double precision, parameter :: smallc = 1.e-10 - double precision, parameter :: smallrho = 1.e-10 - double precision, parameter :: smallp = 1.e-10 - - double precision :: rho_l, un_l, ut_l, rhoe_l, p_l - double precision :: rho_r, un_r, ut_r, rhoe_r, p_r - double precision :: xn(nspec) - - double precision :: rhostar_l, rhostar_r, rho_avg - double precision :: ustar, pstar - double precision :: Q, p_min, p_max, p_lr, p_guess - double precision :: factor, factor2 - double precision :: g_l, g_r, A_l, B_l, A_r, B_r, z - double precision :: S_l, S_r, S_c - double precision :: c_l, c_r, c_avg - - double precision :: U_state(0:nvar-1) - double precision :: HLLCfactor - - - - nx = qx - 2*ng; ny = qy - 2*ng - ilo = ng; ihi = ng+nx-1; jlo = ng; jhi = ng+ny-1 - - do j = jlo-1, jhi+1 - do i = ilo-1, ihi+1 - - ! primitive variable states - rho_l = U_l(i,j,idens) - - ! un = normal velocity; ut = transverse velocity - if (idir == 1) then - un_l = U_l(i,j,ixmom)/rho_l - ut_l = U_l(i,j,iymom)/rho_l - else - un_l = U_l(i,j,iymom)/rho_l - ut_l = U_l(i,j,ixmom)/rho_l - endif - - rhoe_l = U_l(i,j,iener) - 0.5*rho_l*(un_l**2 + ut_l**2) - - p_l = rhoe_l*(gamma - 1.0d0) - p_l = max(p_l, smallp) - - rho_r = U_r(i,j,idens) - - if (idir == 1) then - un_r = U_r(i,j,ixmom)/rho_r - ut_r = U_r(i,j,iymom)/rho_r - else - un_r = U_r(i,j,iymom)/rho_r - ut_r = U_r(i,j,ixmom)/rho_r - endif - - rhoe_r = U_r(i,j,iener) - 0.5*rho_r*(un_r**2 + ut_r**2) - - p_r = rhoe_r*(gamma - 1.0d0) - p_r = max(p_r, smallp) - - - ! compute the sound speeds - c_l = max(smallc, sqrt(gamma*p_l/rho_l)) - c_r = max(smallc, sqrt(gamma*p_r/rho_r)) - - ! Estimate the star quantities -- use one of three methods to - ! do this -- the primitive variable Riemann solver, the two - ! shock approximation, or the two rarefaction approximation. - ! Pick the method based on the pressure states at the - ! interface. - - p_max = max(p_l, p_r) - p_min = min(p_l, p_r) - - Q = p_max/p_min - - rho_avg = 0.5*(rho_l + rho_r) - c_avg = 0.5*(c_l + c_r) - - ! primitive variable Riemann solver (Toro, 9.3) - factor = rho_avg*c_avg - factor2 = rho_avg/c_avg - - pstar = 0.5*(p_l + p_r) + 0.5*(un_l - un_r)*factor - ustar = 0.5*(un_l + un_r) + 0.5*(p_l - p_r)/factor - - rhostar_l = rho_l + (un_l - ustar)*factor2 - rhostar_r = rho_r + (ustar - un_r)*factor2 - - if (Q > 2 .and. (pstar < p_min .or. pstar > p_max)) then - - ! use a more accurate Riemann solver for the estimate here - - if (pstar < p_min) then - - ! 2-rarefaction Riemann solver - z = (gamma - 1.0d0)/(2.0d0*gamma) - p_lr = (p_l/p_r)**z - - ustar = (p_lr*un_l/c_l + un_r/c_r + & - 2.0d0*(p_lr - 1.0d0)/(gamma - 1.0d0)) / & - (p_lr/c_l + 1.0d0/c_r) - - pstar = 0.5d0*(p_l*(1.0d0 + (gamma - 1.0d0)*(un_l - ustar)/ & - (2.0d0*c_l) )**(1.0d0/z) + & - p_r*(1.0d0 + (gamma - 1.0d0)*(ustar - un_r)/ & - (2.0d0*c_r) )**(1.0d0/z) ) - - rhostar_l = rho_l*(pstar/p_l)**(1.0d0/gamma) - rhostar_r = rho_r*(pstar/p_r)**(1.0d0/gamma) - - else - - ! 2-shock Riemann solver - A_r = 2.0/((gamma + 1.0d0)*rho_r) - B_r = p_r*(gamma - 1.0d0)/(gamma + 1.0d0) - - A_l = 2.0/((gamma + 1.0d0)*rho_l) - B_l = p_l*(gamma - 1.0d0)/(gamma + 1.0d0) - - ! guess of the pressure - p_guess = max(0.0d0, pstar) - - g_l = sqrt(A_l / (p_guess + B_l)) - g_r = sqrt(A_r / (p_guess + B_r)) - - pstar = (g_l*p_l + g_r*p_r - (un_r - un_l))/(g_l + g_r) - - ustar = 0.5*(un_l + un_r) + & - 0.5*( (pstar - p_r)*g_r - (pstar - p_l)*g_l) - - rhostar_l = rho_l*(pstar/p_l + (gamma-1.0d0)/(gamma+1.0d0))/ & - ( (gamma-1.0d0)/(gamma+1.0d0)*(pstar/p_l) + 1.0d0) - - rhostar_r = rho_r*(pstar/p_r + (gamma-1.0d0)/(gamma+1.0d0))/ & - ( (gamma-1.0d0)/(gamma+1.0d0)*(pstar/p_r) + 1.0d0) - - endif - endif - - ! estimate the nonlinear wave speeds - - if (pstar <= p_l) then - ! rarefaction - S_l = un_l - c_l - else - ! shock - S_l = un_l - c_l*sqrt(1.0d0 + ((gamma+1.0d0)/(2.0d0*gamma))* & - (pstar/p_l - 1.0d0)) - endif - - if (pstar <= p_r) then - ! rarefaction - S_r = un_r + c_r - else - ! shock - S_r = un_r + c_r*sqrt(1.0d0 + ((gamma+1.0d0)/(2.0d0/gamma))* & - (pstar/p_r - 1.0d0)) - endif - - ! We could just take S_c = u_star as the estimate for the - ! contact speed, but we can actually do this more accurately - ! by using the Rankine-Hugonoit jump conditions across each - ! of the waves (see Toro 10.58, Batten et al. SIAM - ! J. Sci. and Stat. Comp., 18:1553 (1997) - S_c = (p_r - p_l + rho_l*un_l*(S_l - un_l) - rho_r*un_r*(S_r - un_r))/ & - (rho_l*(S_l - un_l) - rho_r*(S_r - un_r)) - - - ! figure out which region we are in and compute the state and - ! the interface fluxes using the HLLC Riemann solver - if (S_r <= 0.0d0) then - ! R region - U_state(:) = U_r(i,j,:) - - call consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nvar, nspec, & - U_state, F(i,j,:)) - - else if (S_r > 0.0d0 .and. S_c <= 0) then - ! R* region - HLLCfactor = rho_r*(S_r - un_r)/(S_r - S_c) - - U_state(idens) = HLLCfactor - - if (idir == 1) then - U_state(ixmom) = HLLCfactor*S_c - U_state(iymom) = HLLCfactor*ut_r - else - U_state(ixmom) = HLLCfactor*ut_r - U_state(iymom) = HLLCfactor*S_c - endif - - U_state(iener) = HLLCfactor*(U_r(i,j,iener)/rho_r + & - (S_c - un_r)*(S_c + p_r/(rho_r*(S_r - un_r)))) - - ! species - if (nspec > 0) then - U_state(irhoX:irhoX-1+nspec) = HLLCfactor*U_r(i,j,irhoX:irhoX-1+nspec)/rho_r - endif - - ! find the flux on the right interface - call consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nvar, nspec, & - U_r(i,j,:), F(i,j,:)) - - ! correct the flux - F(i,j,:) = F(i,j,:) + S_r*(U_state(:) - U_r(i,j,:)) - - else if (S_c > 0.0d0 .and. S_l < 0.0) then - ! L* region - HLLCfactor = rho_l*(S_l - un_l)/(S_l - S_c) - - U_state(idens) = HLLCfactor - - if (idir == 1) then - U_state(ixmom) = HLLCfactor*S_c - U_state(iymom) = HLLCfactor*ut_l - else - U_state(ixmom) = HLLCfactor*ut_l - U_state(iymom) = HLLCfactor*S_c - endif - - U_state(iener) = HLLCfactor*(U_l(i,j,iener)/rho_l + & - (S_c - un_l)*(S_c + p_l/(rho_l*(S_l - un_l)))) - - ! species - if (nspec > 0) then - U_state(irhoX:irhoX-1+nspec) = HLLCfactor*U_l(i,j,irhoX:irhoX-1+nspec)/rho_l - endif - - ! find the flux on the left interface - call consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nvar, nspec, & - U_l(i,j,:), F(i,j,:)) - - ! correct the flux - F(i,j,:) = F(i,j,:) + S_l*(U_state(:) - U_l(i,j,:)) - - else - ! L region - U_state(:) = U_l(i,j,:) - - call consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nvar, nspec, & - U_state, F(i,j,:)) - - endif - - ! we should deal with solid boundaries somehow here - - enddo - enddo -end subroutine riemann_HLLC - -subroutine consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nvar, nspec, U_state, F) - - integer, intent(in) :: idir - double precision, intent(in) :: gamma - integer, intent(in) :: idens, ixmom, iymom, iener, irhoX, nvar, nspec - double precision, intent(in) :: U_state(0:nvar-1) - double precision, intent(out) :: F(0:nvar-1) - - double precision :: p, u, v - - u = U_state(ixmom)/U_state(idens) - v = U_state(iymom)/U_state(idens) - - p = (U_state(iener) - 0.5d0*U_state(idens)*(u*u + v*v))*(gamma - 1.0d0) - - if (idir == 1) then - F(idens) = U_state(idens)*u - F(ixmom) = U_state(ixmom)*u + p - F(iymom) = U_state(iymom)*u - F(iener) = (U_state(iener) + p)*u - if (nspec > 0) then - F(irhoX:irhoX-1+nspec) = U_state(irhoX:irhoX-1+nspec)*u - endif - else - F(idens) = U_state(idens)*v - F(ixmom) = U_state(ixmom)*v - F(iymom) = U_state(iymom)*v + p - F(iener) = (U_state(iener) + p)*v - if (nspec > 0) then - F(irhoX:irhoX-1+nspec) = U_state(irhoX:irhoX-1+nspec)*v - endif - endif - -end subroutine consFlux - - -subroutine artificial_viscosity(qx, qy, ng, dx, dy, & - cvisc, u, v, avisco_x, avisco_y) - - implicit none - integer, intent(in) :: qx, qy, ng - double precision, intent(in) :: dx, dy - double precision, intent(in) :: cvisc - - ! 0-based indexing to match python - double precision, intent(in) :: u(0:qx-1, 0:qy-1) - double precision, intent(in) :: v(0:qx-1, 0:qy-1) - double precision, intent(out) :: avisco_x(0:qx-1, 0:qy-1) - double precision, intent(out) :: avisco_y(0:qx-1, 0:qy-1) - -!f2py depend(qx, qy) :: u, v -!f2py depend(qx, qy) :: avisco_x, avisco_y -!f2py intent(in) :: u, v -!f2py intent(out) :: avisco_x, avisco_y - - ! compute the artifical viscosity. Here, we compute edge-centered - ! approximations to the divergence of the velocity. This follows - ! directly Colella & Woodward (1984) Eq. 4.5 - ! - ! data locations: - ! - ! j+3/2--+---------+---------+---------+ - ! | | | | - ! j+1 + | | | - ! | | | | - ! j+1/2--+---------+---------+---------+ - ! | | | | - ! j + X | | - ! | | | | - ! j-1/2--+---------+----Y----+---------+ - ! | | | | - ! j-1 + | | | - ! | | | | - ! j-3/2--+---------+---------+---------+ - ! | | | | | | | - ! i-1 i i+1 - ! i-3/2 i-1/2 i+1/2 i+3/2 - ! - ! X is the location of avisco_x(i,j) - ! Y is the location of avisco_y(i,j) - - integer :: ilo, ihi, jlo, jhi - integer :: nx, ny - - integer :: i, j - - double precision :: divU_x, divU_y - - nx = qx - 2*ng; ny = qy - 2*ng - ilo = ng; ihi = ng+nx-1; jlo = ng; jhi = ng+ny-1 - - do j = jlo-1, jhi+1 - do i = ilo-1, ihi+1 - - ! start by computing the divergence on the x-interface. The - ! x-difference is simply the difference of the cell-centered - ! x-velocities on either side of the x-interface. For the - ! y-difference, first average the four cells to the node on - ! each end of the edge, and then difference these to find the - ! edge centered y difference. - divU_x = (u(i,j) - u(i-1,j))/dx + & - 0.25d0*(v(i,j+1) + v(i-1,j+1) - v(i,j-1) - v(i-1,j-1))/dy - - avisco_x(i,j) = cvisc*max(-divU_x*dx, 0.0d0) - - ! now the y-interface value - divU_y = 0.25d0*(u(i+1,j) + u(i+1,j-1) - u(i-1,j) - u(i-1,j-1))/dx + & - (v(i,j) - v(i,j-1))/dy - - avisco_y(i,j) = cvisc*max(-divU_y*dy, 0.0d0) - - enddo - enddo - -end subroutine artificial_viscosity diff --git a/compressible/tests/rt_0945.h5 b/compressible/tests/rt_0945.h5 index c5c071de9..bde7bd9f4 100644 Binary files a/compressible/tests/rt_0945.h5 and b/compressible/tests/rt_0945.h5 differ diff --git a/compressible/tests/sod_x_0076.h5 b/compressible/tests/sod_x_0076.h5 index 84bbecb2c..af3814e93 100644 Binary files a/compressible/tests/sod_x_0076.h5 and b/compressible/tests/sod_x_0076.h5 differ diff --git a/compressible/unsplit_fluxes.py b/compressible/unsplit_fluxes.py index 1d74fb00b..21dd983ff 100644 --- a/compressible/unsplit_fluxes.py +++ b/compressible/unsplit_fluxes.py @@ -218,9 +218,9 @@ def unsplit_fluxes(my_data, my_aux, rp, ivars, solid, tc, dt): tm_states = tc.timer("interfaceStates") tm_states.begin() - V_l, V_r = ifc.states(1, myg.qx, myg.qy, myg.ng, myg.dx, dt, + V_l, V_r = ifc.states(1, myg.ng, myg.dx, dt, ivars.irho, ivars.iu, ivars.iv, ivars.ip, ivars.ix, - ivars.nvar, ivars.naux, + ivars.naux, gamma, q, ldx) @@ -237,9 +237,9 @@ def unsplit_fluxes(my_data, my_aux, rp, ivars, solid, tc, dt): # left and right primitive variable states tm_states.begin() - _V_l, _V_r = ifc.states(2, myg.qx, myg.qy, myg.ng, myg.dy, dt, + _V_l, _V_r = ifc.states(2, myg.ng, myg.dy, dt, ivars.irho, ivars.iu, ivars.iv, ivars.ip, ivars.ix, - ivars.nvar, ivars.naux, + ivars.naux, gamma, q, ldy) V_l = ai.ArrayIndexer(d=_V_l, grid=myg) @@ -295,13 +295,13 @@ def unsplit_fluxes(my_data, my_aux, rp, ivars, solid, tc, dt): else: msg.fail("ERROR: Riemann solver undefined") - _fx = riemannFunc(1, myg.qx, myg.qy, myg.ng, - ivars.nvar, ivars.idens, ivars.ixmom, ivars.iymom, ivars.iener, ivars.irhox, ivars.naux, + _fx = riemannFunc(1, myg.ng, + ivars.idens, ivars.ixmom, ivars.iymom, ivars.iener, ivars.irhox, ivars.naux, solid.xl, solid.xr, gamma, U_xl, U_xr) - _fy = riemannFunc(2, myg.qx, myg.qy, myg.ng, - ivars.nvar, ivars.idens, ivars.ixmom, ivars.iymom, ivars.iener, ivars.irhox, ivars.naux, + _fy = riemannFunc(2, myg.ng, + ivars.idens, ivars.ixmom, ivars.iymom, ivars.iener, ivars.irhox, ivars.naux, solid.yl, solid.yr, gamma, U_yl, U_yr) @@ -394,13 +394,13 @@ def unsplit_fluxes(my_data, my_aux, rp, ivars, solid, tc, dt): tm_riem.begin() - _fx = riemannFunc(1, myg.qx, myg.qy, myg.ng, - ivars.nvar, ivars.idens, ivars.ixmom, ivars.iymom, ivars.iener, ivars.irhox, ivars.naux, + _fx = riemannFunc(1, myg.ng, + ivars.idens, ivars.ixmom, ivars.iymom, ivars.iener, ivars.irhox, ivars.naux, solid.xl, solid.xr, gamma, U_xl, U_xr) - _fy = riemannFunc(2, myg.qx, myg.qy, myg.ng, - ivars.nvar, ivars.idens, ivars.ixmom, ivars.iymom, ivars.iener, ivars.irhox, ivars.naux, + _fy = riemannFunc(2, myg.ng, + ivars.idens, ivars.ixmom, ivars.iymom, ivars.iener, ivars.irhox, ivars.naux, solid.yl, solid.yr, gamma, U_yl, U_yr) @@ -414,8 +414,7 @@ def unsplit_fluxes(my_data, my_aux, rp, ivars, solid, tc, dt): # ========================================================================= cvisc = rp.get_param("compressible.cvisc") - _ax, _ay = ifc.artificial_viscosity( - myg.qx, myg.qy, myg.ng, myg.dx, myg.dy, + _ax, _ay = ifc.artificial_viscosity(myg.ng, myg.dx, myg.dy, cvisc, q.v(n=ivars.iu, buf=myg.ng), q.v(n=ivars.iv, buf=myg.ng)) avisco_x = ai.ArrayIndexer(d=_ax, grid=myg) diff --git a/compressible_fv4/fluxes.py b/compressible_fv4/fluxes.py index f4578315b..d36bfb2a5 100644 --- a/compressible_fv4/fluxes.py +++ b/compressible_fv4/fluxes.py @@ -2,7 +2,7 @@ import numpy as np -import advection_fv4.interface as interface_f +import advection_fv4.interface as interface import compressible as comp import compressible.interface as cf import mesh.reconstruction as reconstruction @@ -102,7 +102,7 @@ def fluxes(myd, rp, ivars, solid, tc): else: for n in range(ivars.nq): - q_l[:, :, n], q_r[:, :, n] = interface_f.states(q_avg[:, :, n], myg.qx, myg.qy, myg.ng, idir) + q_l[:, :, n], q_r[:, :, n] = interface.states(q_avg[:, :, n], myg.ng, idir) # apply flattening for n in range(ivars.nq): @@ -117,8 +117,8 @@ def fluxes(myd, rp, ivars, solid, tc): q_r.v(n=n, buf=2)[:, :] = xi.v(buf=2)*q_r.v(n=n, buf=2) + \ (1.0 - xi.v(buf=2))*q_avg.v(n=n, buf=2) - _q = cf.riemann_prim(idir, myg.qx, myg.qy, myg.ng, - ivars.nq, ivars.irho, ivars.iu, ivars.iv, ivars.ip, ivars.ix, ivars.naux, + _q = cf.riemann_prim(idir, myg.ng, + ivars.irho, ivars.iu, ivars.iv, ivars.ip, ivars.ix, ivars.naux, 0, 0, gamma, q_l, q_r) diff --git a/compressible_rk/fluxes.py b/compressible_rk/fluxes.py index 027e5848b..54427b7d1 100644 --- a/compressible_rk/fluxes.py +++ b/compressible_rk/fluxes.py @@ -19,7 +19,7 @@ """ -import compressible.interface as interface_f +import compressible.interface as interface import compressible as comp import mesh.reconstruction as reconstruction import mesh.array_indexer as ai @@ -168,19 +168,19 @@ def fluxes(my_data, rp, ivars, solid, tc): riemann = rp.get_param("compressible.riemann") if riemann == "HLLC": - riemannFunc = interface_f.riemann_hllc + riemannFunc = interface.riemann_hllc elif riemann == "CGF": - riemannFunc = interface_f.riemann_cgf + riemannFunc = interface.riemann_cgf else: msg.fail("ERROR: Riemann solver undefined") - _fx = riemannFunc(1, myg.qx, myg.qy, myg.ng, - ivars.nvar, ivars.idens, ivars.ixmom, ivars.iymom, ivars.iener, ivars.irhox, ivars.naux, + _fx = riemannFunc(1, myg.ng, + ivars.idens, ivars.ixmom, ivars.iymom, ivars.iener, ivars.irhox, ivars.naux, solid.xl, solid.xr, gamma, U_xl, U_xr) - _fy = riemannFunc(2, myg.qx, myg.qy, myg.ng, - ivars.nvar, ivars.idens, ivars.ixmom, ivars.iymom, ivars.iener, ivars.irhox, ivars.naux, + _fy = riemannFunc(2, myg.ng, + ivars.idens, ivars.ixmom, ivars.iymom, ivars.iener, ivars.irhox, ivars.naux, solid.yl, solid.yr, gamma, U_yl, U_yr) @@ -194,8 +194,7 @@ def fluxes(my_data, rp, ivars, solid, tc): # ========================================================================= cvisc = rp.get_param("compressible.cvisc") - _ax, _ay = interface_f.artificial_viscosity( - myg.qx, myg.qy, myg.ng, myg.dx, myg.dy, + _ax, _ay = interface.artificial_viscosity(myg.ng, myg.dx, myg.dy, cvisc, q.v(n=ivars.iu, buf=myg.ng), q.v(n=ivars.iv, buf=myg.ng)) avisco_x = ai.ArrayIndexer(d=_ax, grid=myg) diff --git a/compressible_rk/tests/rt_1835.h5 b/compressible_rk/tests/rt_1835.h5 index a67a5dd63..d35dc13fd 100644 Binary files a/compressible_rk/tests/rt_1835.h5 and b/compressible_rk/tests/rt_1835.h5 differ diff --git a/diffusion/tests/gaussian_0164.h5 b/diffusion/tests/gaussian_0164.h5 index 9120518ba..dd27eac19 100644 Binary files a/diffusion/tests/gaussian_0164.h5 and b/diffusion/tests/gaussian_0164.h5 differ diff --git a/docs/source/ack.rst b/docs/source/ack.rst index 75917ce5d..041aa17ba 100644 --- a/docs/source/ack.rst +++ b/docs/source/ack.rst @@ -47,5 +47,6 @@ using the Numeric array package and handwritten C extensions for the compute-intensive kernels. It was ported to numarray when that replaced Numeric, and continued to use C extensions. This version "pyro2" was resurrected beginning in 2012 and rewritten for numpy -using f2py, and brought up to date. +using f2py, and brought up to date. Most recently we've dropped +f2py and are using numba for the compute-intensive kernels. diff --git a/docs/source/design.rst b/docs/source/design.rst index d6556a89a..83504a309 100644 --- a/docs/source/design.rst +++ b/docs/source/design.rst @@ -1,12 +1,11 @@ Design ideas ============ -pyro is written primarily in python (by default, we expect python 3), -with a few low-level routines written in Fortran for performance. The +pyro is written entirely in python (by default, we expect python 3), +with a few low-level routines compiled *just-in-time* by `numba` for performance. The ``numpy`` package is used for representing arrays throughout the python code and the ``matplotlib`` library is used for -visualization. We use ``f2py`` (part of NumPy) to interface with some -Fortran code. Finally, ``pytest`` is used for unit testing of some +visualization. Finally, ``pytest`` is used for unit testing of some components. All solvers are written for a 2-d grid. This gives a good balance @@ -128,31 +127,10 @@ The overall structure is: modes. -Fortran -------- - -Fortran is used to speed up some critical portions of the code, and in -many cases, provides more clarity than trying to write optimized -python code using array operations in numpy. The Fortran code -seemlessly integrates into python using f2py. - -Wherever Fortran is used, we enforce the following design rule: the -Fortran functions must be completely self-contained, with all -information coming through the interface. No external dependencies -are allowed. Each pyro module will have (at most) a single Fortran -file and can be compiled into a library via a single f2py command line -invocation. - -A single script, ``mk.sh``, in the top-level directory run with the ``fortran`` -argument will compile all the Fortran source. - - Numba ----- -Instead of using Fortran to speed up performance-critical parts of the code, -``pyro`` can also use python-only versions of these functions compiled with -``numba``. Numba is a *just-in-time compiler* for python. When a call is first +``numba`` is used to speed up some critical portions of the code. Numba is a *just-in-time compiler* for python. When a call is first made to a function decorated with Numba's ``@njit`` decorator, it is compiled to machine code 'just-in-time' for it to be executed. Once compiled, it can then run at (near-to) native machine code speed. @@ -163,9 +141,9 @@ time you run the same bit of code, Numba will use the saved version rather than compiling the code again, saving some compilation time at the start of the simulation. -``pyro`` will use the ``numba`` functions (rather than the Fortran modules) by -default. +.. note:: + Because we have chosen to cache the compiled code, Numba will save it in the ``__pycache__`` directories. If you change the code, a new version will be compiled and saved, but the old version will not be deleted. Over time, you may end up with many unneeded files saved in the ``__pycache__`` directories. To clean up these files, you can run ``./mk.sh clean`` in the main ``pyro2`` directory. Main driver ----------- diff --git a/docs/source/installation.rst b/docs/source/installation.rst index 65a3d4e64..9b6e4f8d4 100644 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -16,14 +16,6 @@ The following python packages are required: * ``numba`` * ``pytest`` (for unit tests) -By default, ``pyro`` will run using code written entirely in python (with the help -of ``numba`` to speed up some functions). However, some of the functions also have -Fortran versions which run several times faster than their python counterparts. -This can be useful for running more intensive calculations (e.g. with higher -resolution or a large number of timesteps). To compile the Fortran functions, -you will need the python package ``f2py`` (part of NumPy) and have ``gFortran`` -installed on your computer for ``f2py`` to be able to compile the code. - The following steps are needed before running pyro: * add ``pyro/`` to your ``PYTHONPATH`` environment variable. For @@ -33,20 +25,6 @@ The following steps are needed before running pyro: export PYTHONPATH="/path/to/pyro/:${PYTHONPATH}" -* build the modules by running the ``mk.sh`` script. To use the python-only version, it - should be sufficient to just do: - - .. code-block:: none - - ./mk.sh - - - Alternatively, to build the faster Fortran modules do: - - .. code-block:: none - - ./mk.sh fortran - * define the environment variable ``PYRO_HOME`` to point to the ``pyro2/`` directory @@ -54,11 +32,6 @@ The following steps are needed before running pyro: export PYRO_HOME="/path/to/pyro/" -.. note:: - - If you have built the Fortran modules and want to go back to using the python-only - modules, you will need to first run ``./mk.sh clean``. - Quick test ---------- diff --git a/incompressible/incomp_interface.py b/incompressible/incomp_interface.py index c9361e1d1..8f4164d77 100644 --- a/incompressible/incomp_interface.py +++ b/incompressible/incomp_interface.py @@ -3,7 +3,7 @@ @njit(cache=True) -def mac_vels(qx, qy, ng, dx, dy, dt, +def mac_vels(ng, dx, dy, dt, u, v, ldelta_ux, ldelta_vx, ldelta_uy, ldelta_vy, @@ -13,8 +13,6 @@ def mac_vels(qx, qy, ng, dx, dy, dt, Parameters ---------- - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells dx, dy : float @@ -38,7 +36,7 @@ def mac_vels(qx, qy, ng, dx, dy, dt, # get the full u and v left and right states (including transverse # terms) on both the x- and y-interfaces - u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = get_interface_states(qx, qy, ng, dx, dy, dt, + u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = get_interface_states(ng, dx, dy, dt, u, v, ldelta_ux, ldelta_vx, ldelta_uy, ldelta_vy, @@ -47,15 +45,15 @@ def mac_vels(qx, qy, ng, dx, dy, dt, # Riemann problem -- this follows Burger's equation. We don't use # any input velocity for the upwinding. Also, we only care about # the normal states here (u on x and v on y) - u_MAC = riemann_and_upwind(qx, qy, ng, u_xl, u_xr) - v_MAC = riemann_and_upwind(qx, qy, ng, v_yl, v_yr) + u_MAC = riemann_and_upwind(ng, u_xl, u_xr) + v_MAC = riemann_and_upwind(ng, v_yl, v_yr) return u_MAC, v_MAC # xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx @njit(cache=True) -def states(qx, qy, ng, dx, dy, dt, +def states(ng, dx, dy, dt, u, v, ldelta_ux, ldelta_vx, ldelta_uy, ldelta_vy, @@ -68,8 +66,6 @@ def states(qx, qy, ng, dx, dy, dt, Parameters ---------- - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells dx, dy : float @@ -95,7 +91,7 @@ def states(qx, qy, ng, dx, dy, dt, # get the full u and v left and right states (including transverse # terms) on both the x- and y-interfaces - u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = get_interface_states(qx, qy, ng, dx, dy, dt, + u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = get_interface_states(ng, dx, dy, dt, u, v, ldelta_ux, ldelta_vx, ldelta_uy, ldelta_vy, @@ -103,17 +99,17 @@ def states(qx, qy, ng, dx, dy, dt, # upwind using the MAC velocity to determine which state exists on # the interface - u_xint = upwind(qx, qy, ng, u_xl, u_xr, u_MAC) - v_xint = upwind(qx, qy, ng, v_xl, v_xr, u_MAC) - u_yint = upwind(qx, qy, ng, u_yl, u_yr, v_MAC) - v_yint = upwind(qx, qy, ng, v_yl, v_yr, v_MAC) + u_xint = upwind(ng, u_xl, u_xr, u_MAC) + v_xint = upwind(ng, v_xl, v_xr, u_MAC) + u_yint = upwind(ng, u_yl, u_yr, v_MAC) + v_yint = upwind(ng, v_yl, v_yr, v_MAC) return u_xint, v_xint, u_yint, v_yint # xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx @njit(cache=True) -def get_interface_states(qx, qy, ng, dx, dy, dt, +def get_interface_states(ng, dx, dy, dt, u, v, ldelta_ux, ldelta_vx, ldelta_uy, ldelta_vy, @@ -124,8 +120,6 @@ def get_interface_states(qx, qy, ng, dx, dy, dt, Parameters ---------- - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells dx, dy : float @@ -148,6 +142,8 @@ def get_interface_states(qx, qy, ng, dx, dy, dt, y-interfaces """ + qx, qy = u.shape + u_xl = np.zeros((qx, qy)) u_xr = np.zeros((qx, qy)) u_yl = np.zeros((qx, qy)) @@ -200,19 +196,19 @@ def get_interface_states(qx, qy, ng, dx, dy, dt, # now get the normal advective velocities on the interfaces by solving # the Riemann problem. - uhat_adv = riemann(qx, qy, ng, u_xl, u_xr) - vhat_adv = riemann(qx, qy, ng, v_yl, v_yr) + uhat_adv = riemann(ng, u_xl, u_xr) + vhat_adv = riemann(ng, v_yl, v_yr) # now that we have the advective velocities, upwind the left and right # states using the appropriate advective velocity. # on the x-interfaces, we upwind based on uhat_adv - u_xint = upwind(qx, qy, ng, u_xl, u_xr, uhat_adv) - v_xint = upwind(qx, qy, ng, v_xl, v_xr, uhat_adv) + u_xint = upwind(ng, u_xl, u_xr, uhat_adv) + v_xint = upwind(ng, v_xl, v_xr, uhat_adv) # on the y-interfaces, we upwind based on vhat_adv - u_yint = upwind(qx, qy, ng, u_yl, u_yr, vhat_adv) - v_yint = upwind(qx, qy, ng, v_yl, v_yr, vhat_adv) + u_yint = upwind(ng, u_yl, u_yr, vhat_adv) + v_yint = upwind(ng, v_yl, v_yr, vhat_adv) # at this point, these states are the `hat' states -- they only # considered the normal to the interface portion of the predictor. @@ -261,15 +257,13 @@ def get_interface_states(qx, qy, ng, dx, dy, dt, # xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx @njit(cache=True) -def upwind(qx, qy, ng, q_l, q_r, s): +def upwind(ng, q_l, q_r, s): r""" upwind the left and right states based on the specified input velocity, s. The resulting interface state is q_int Parameters ---------- - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells q_l, q_r : ndarray @@ -283,6 +277,8 @@ def upwind(qx, qy, ng, q_l, q_r, s): Upwinded state """ + qx, qy = s.shape + nx = qx - 2 * ng ny = qy - 2 * ng ilo = ng @@ -290,7 +286,7 @@ def upwind(qx, qy, ng, q_l, q_r, s): jlo = ng jhi = ng + ny - q_int = np.zeros((qx, qy)) + q_int = np.zeros_like(s) for i in range(ilo - 1, ihi + 1): for j in range(jlo - 1, jhi + 1): @@ -307,7 +303,7 @@ def upwind(qx, qy, ng, q_l, q_r, s): # xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx @njit(cache=True) -def riemann(qx, qy, ng, q_l, q_r): +def riemann(ng, q_l, q_r): """ Solve the Burger's Riemann problem given the input left and right states and return the state on the interface. @@ -316,8 +312,6 @@ def riemann(qx, qy, ng, q_l, q_r): Parameters ---------- - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells q_l, q_r : ndarray @@ -329,6 +323,8 @@ def riemann(qx, qy, ng, q_l, q_r): Interface state """ + qx, qy = q_l.shape + nx = qx - 2 * ng ny = qy - 2 * ng ilo = ng @@ -336,7 +332,7 @@ def riemann(qx, qy, ng, q_l, q_r): jlo = ng jhi = ng + ny - s = np.zeros((qx, qy)) + s = np.zeros_like(q_l) for i in range(ilo - 1, ihi + 1): for j in range(jlo - 1, jhi + 1): @@ -353,7 +349,7 @@ def riemann(qx, qy, ng, q_l, q_r): # xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx @njit(cache=True) -def riemann_and_upwind(qx, qy, ng, q_l, q_r): +def riemann_and_upwind(ng, q_l, q_r): r""" First solve the Riemann problem given q_l and q_r to give the velocity on the interface and: use this velocity to upwind to @@ -364,8 +360,6 @@ def riemann_and_upwind(qx, qy, ng, q_l, q_r): Parameters ---------- - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells q_l, q_r : ndarray @@ -377,5 +371,5 @@ def riemann_and_upwind(qx, qy, ng, q_l, q_r): Upwinded state """ - s = riemann(qx, qy, ng, q_l, q_r) - return upwind(qx, qy, ng, q_l, q_r, s) + s = riemann(ng, q_l, q_r) + return upwind(ng, q_l, q_r, s) diff --git a/incompressible/incomp_interface_f.f90 b/incompressible/incomp_interface_f.f90 deleted file mode 100644 index 0f0567e0c..000000000 --- a/incompressible/incomp_interface_f.f90 +++ /dev/null @@ -1,372 +0,0 @@ -subroutine mac_vels(qx, qy, ng, dx, dy, dt, & - u, v, & - ldelta_ux, ldelta_vx, & - ldelta_uy, ldelta_vy, & - gradp_x, gradp_y, & - u_MAC, v_MAC) - - implicit none - - integer, intent(in) :: qx, qy, ng - double precision, intent(in) :: dx, dy, dt - - ! 0-based indexing to match python - double precision, intent(inout) :: u(0:qx-1, 0:qy-1) - double precision, intent(inout) :: v(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: ldelta_ux(0:qx-1, 0:qy-1) - double precision, intent(inout) :: ldelta_vx(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: ldelta_uy(0:qx-1, 0:qy-1) - double precision, intent(inout) :: ldelta_vy(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: gradp_x(0:qx-1, 0:qy-1) - double precision, intent(inout) :: gradp_y(0:qx-1, 0:qy-1) - - double precision, intent( out) :: u_MAC(0:qx-1, 0:qy-1) - double precision, intent( out) :: v_MAC(0:qx-1, 0:qy-1) - -!f2py depend(qx, qy) :: u, v -!f2py depend(qx, qy) :: ldelta_ux, ldelta_vx, ldelta_uy, ldelta_vy -!f2py depend(qx, qy) :: gradp_x, gradp_y -!f2py depend(qx, qy) :: u_MAC, v_MAC -!f2py intent(in) :: u, v, gradp_x, gradp_y -!f2py intent(in) :: ldelta_ux, ldelta_vx, ldelta_uy, ldelta_vy -!f2py intent(out) :: u_MAC, v_MAC - - - double precision :: u_xl(0:qx-1, 0:qy-1), u_xr(0:qx-1, 0:qy-1) - double precision :: u_yl(0:qx-1, 0:qy-1), u_yr(0:qx-1, 0:qy-1) - double precision :: v_xl(0:qx-1, 0:qy-1), v_xr(0:qx-1, 0:qy-1) - double precision :: v_yl(0:qx-1, 0:qy-1), v_yr(0:qx-1, 0:qy-1) - - - ! get the full u and v left and right states (including transverse - ! terms) on both the x- and y-interfaces - call get_interface_states(qx, qy, ng, dx, dy, dt, & - u, v, & - ldelta_ux, ldelta_vx, & - ldelta_uy, ldelta_vy, & - gradp_x, gradp_y, & - u_xl, u_xr, u_yl, u_yr, & - v_xl, v_xr, v_yl, v_yr) - - - ! Riemann problem -- this follows Burger's equation. We don't use - ! any input velocity for the upwinding. Also, we only care about - ! the normal states here (u on x and v on y) - call riemann_and_upwind(qx, qy, ng, u_xl, u_xr, u_MAC) - call riemann_and_upwind(qx, qy, ng, v_yl, v_yr, v_MAC) - -end subroutine mac_vels - - -!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx -subroutine states(qx, qy, ng, dx, dy, dt, & - u, v, & - ldelta_ux, ldelta_vx, & - ldelta_uy, ldelta_vy, & - gradp_x, gradp_y, & - u_MAC, v_MAC, & - u_xint, v_xint, u_yint, v_yint) - - ! this is similar to mac_vels, but it predicts the interface states - ! of both u and v on both interfaces, using the MAC velocities to - ! do the upwinding. - - implicit none - - integer, intent(in) :: qx, qy, ng - double precision, intent(in) :: dx, dy, dt - - ! 0-based indexing to match python - double precision, intent(inout) :: u(0:qx-1, 0:qy-1) - double precision, intent(inout) :: v(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: ldelta_ux(0:qx-1, 0:qy-1) - double precision, intent(inout) :: ldelta_vx(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: ldelta_uy(0:qx-1, 0:qy-1) - double precision, intent(inout) :: ldelta_vy(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: gradp_x(0:qx-1, 0:qy-1) - double precision, intent(inout) :: gradp_y(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: u_MAC(0:qx-1, 0:qy-1) - double precision, intent(inout) :: v_MAC(0:qx-1, 0:qy-1) - - double precision, intent(out) :: u_xint(0:qx-1, 0:qy-1), u_yint(0:qx-1, 0:qy-1) - double precision, intent(out) :: v_xint(0:qx-1, 0:qy-1), v_yint(0:qx-1, 0:qy-1) - -!f2py depend(qx, qy) :: u, v -!f2py depend(qx, qy) :: ldelta_ux, ldelta_vx, ldelta_uy, ldelta_vy -!f2py depend(qx, qy) :: gradp_x, gradp_y -!f2py depend(qx, qy) :: u_MAC, v_MAC -!f2py intent(in) :: u, v, gradp_x, gradp_y -!f2py intent(in) :: ldelta_ux, ldelta_vx, ldelta_uy, ldelta_vy -!f2py intent(in) :: u_MAC, v_MAC -!f2py intent(out) :: u_xint, v_xint, u_yint, v_yint - - double precision :: u_xl(0:qx-1, 0:qy-1), u_xr(0:qx-1, 0:qy-1) - double precision :: u_yl(0:qx-1, 0:qy-1), u_yr(0:qx-1, 0:qy-1) - double precision :: v_xl(0:qx-1, 0:qy-1), v_xr(0:qx-1, 0:qy-1) - double precision :: v_yl(0:qx-1, 0:qy-1), v_yr(0:qx-1, 0:qy-1) - - - ! get the full u and v left and right states (including transverse - ! terms) on both the x- and y-interfaces - call get_interface_states(qx, qy, ng, dx, dy, dt, & - u, v, & - ldelta_ux, ldelta_vx, & - ldelta_uy, ldelta_vy, & - gradp_x, gradp_y, & - u_xl, u_xr, u_yl, u_yr, & - v_xl, v_xr, v_yl, v_yr) - - - ! upwind using the MAC velocity to determine which state exists on - ! the interface - call upwind(qx, qy, ng, u_xl, u_xr, u_MAC, u_xint) - call upwind(qx, qy, ng, v_xl, v_xr, u_MAC, v_xint) - call upwind(qx, qy, ng, u_yl, u_yr, v_MAC, u_yint) - call upwind(qx, qy, ng, v_yl, v_yr, v_MAC, v_yint) - -end subroutine states - - -!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx -subroutine get_interface_states(qx, qy, ng, dx, dy, dt, & - u, v, & - ldelta_ux, ldelta_vx, & - ldelta_uy, ldelta_vy, & - gradp_x, gradp_y, & - u_xl, u_xr, u_yl, u_yr, & - v_xl, v_xr, v_yl, v_yr) - - ! Compute the unsplit predictions of u and v on both the x- and - ! y-interfaces. This includes the transverse terms. - - implicit none - - integer, intent(in) :: qx, qy, ng - double precision, intent(in) :: dx, dy, dt - - ! 0-based indexing to match python - double precision, intent(inout) :: u(0:qx-1, 0:qy-1) - double precision, intent(inout) :: v(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: ldelta_ux(0:qx-1, 0:qy-1) - double precision, intent(inout) :: ldelta_vx(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: ldelta_uy(0:qx-1, 0:qy-1) - double precision, intent(inout) :: ldelta_vy(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: gradp_x(0:qx-1, 0:qy-1) - double precision, intent(inout) :: gradp_y(0:qx-1, 0:qy-1) - - double precision, intent(out) :: u_xl(0:qx-1, 0:qy-1), u_xr(0:qx-1, 0:qy-1) - double precision, intent(out) :: u_yl(0:qx-1, 0:qy-1), u_yr(0:qx-1, 0:qy-1) - double precision, intent(out) :: v_xl(0:qx-1, 0:qy-1), v_xr(0:qx-1, 0:qy-1) - double precision, intent(out) :: v_yl(0:qx-1, 0:qy-1), v_yr(0:qx-1, 0:qy-1) - - integer :: ilo, ihi, jlo, jhi - integer :: nx, ny - integer :: i, j - - double precision :: uhat_adv(0:qx-1, 0:qy-1), vhat_adv(0:qx-1, 0:qy-1) - - double precision :: u_xint(0:qx-1, 0:qy-1), u_yint(0:qx-1, 0:qy-1) - double precision :: v_xint(0:qx-1, 0:qy-1), v_yint(0:qx-1, 0:qy-1) - - double precision :: dtdx, dtdy - double precision :: ubar, vbar, uv_x, vu_y, uu_x, vv_y - - nx = qx - 2*ng; ny = qy - 2*ng - ilo = ng; ihi = ng+nx-1; jlo = ng; jhi = ng+ny-1 - - - ! first predict u and v to both interfaces, considering only the normal - ! part of the predictor. These are the 'hat' states. - - - dtdx = dt/dx - dtdy = dt/dy - - do j = jlo-2, jhi+2 - do i = ilo-2, ihi+2 - - ! u on x-edges - u_xl(i+1,j) = u(i,j) + 0.5d0*(1.0d0 - dtdx*u(i,j))*ldelta_ux(i,j) - u_xr(i ,j) = u(i,j) - 0.5d0*(1.0d0 + dtdx*u(i,j))*ldelta_ux(i,j) - - ! v on x-edges - v_xl(i+1,j) = v(i,j) + 0.5d0*(1.0d0 - dtdx*u(i,j))*ldelta_vx(i,j) - v_xr(i ,j) = v(i,j) - 0.5d0*(1.0d0 + dtdx*u(i,j))*ldelta_vx(i,j) - - ! u on y-edges - u_yl(i,j+1) = u(i,j) + 0.5d0*(1.0d0 - dtdy*v(i,j))*ldelta_uy(i,j) - u_yr(i,j ) = u(i,j) - 0.5d0*(1.0d0 + dtdy*v(i,j))*ldelta_uy(i,j) - - ! v on y-edges - v_yl(i,j+1) = v(i,j) + 0.5d0*(1.0d0 - dtdy*v(i,j))*ldelta_vy(i,j) - v_yr(i,j ) = v(i,j) - 0.5d0*(1.0d0 + dtdy*v(i,j))*ldelta_vy(i,j) - - enddo - enddo - - - ! now get the normal advective velocities on the interfaces by solving - ! the Riemann problem. - call riemann(qx, qy, ng, u_xl, u_xr, uhat_adv) - call riemann(qx, qy, ng, v_yl, v_yr, vhat_adv) - - - ! now that we have the advective velocities, upwind the left and right - ! states using the appropriate advective velocity. - - ! on the x-interfaces, we upwind based on uhat_adv - call upwind(qx, qy, ng, u_xl, u_xr, uhat_adv, u_xint) - call upwind(qx, qy, ng, v_xl, v_xr, uhat_adv, v_xint) - - ! on the y-interfaces, we upwind based on vhat_adv - call upwind(qx, qy, ng, u_yl, u_yr, vhat_adv, u_yint) - call upwind(qx, qy, ng, v_yl, v_yr, vhat_adv, v_yint) - - ! at this point, these states are the `hat' states -- they only - ! considered the normal to the interface portion of the predictor. - - - ! add the transverse flux differences to the preliminary interface states - do j = jlo-2, jhi+2 - do i = ilo-2, ihi+2 - - ubar = 0.5d0*(uhat_adv(i,j) + uhat_adv(i+1,j)) - vbar = 0.5d0*(vhat_adv(i,j) + vhat_adv(i,j+1)) - - ! v du/dy is the transerse term for the u states on x-interfaces - vu_y = vbar*(u_yint(i,j+1) - u_yint(i,j)) - - u_xl(i+1,j) = u_xl(i+1,j) - 0.5*dtdy*vu_y - 0.5*dt*gradp_x(i,j) - u_xr(i ,j) = u_xr(i ,j) - 0.5*dtdy*vu_y - 0.5*dt*gradp_x(i,j) - - ! v dv/dy is the transverse term for the v states on x-interfaces - vv_y = vbar*(v_yint(i,j+1) - v_yint(i,j)) - - v_xl(i+1,j) = v_xl(i+1,j) - 0.5*dtdy*vv_y - 0.5*dt*gradp_y(i,j) - v_xr(i ,j) = v_xr(i ,j) - 0.5*dtdy*vv_y - 0.5*dt*gradp_y(i,j) - - ! u dv/dx is the transverse term for the v states on y-interfaces - uv_x = ubar*(v_xint(i+1,j) - v_xint(i,j)) - - v_yl(i,j+1) = v_yl(i,j+1) - 0.5*dtdx*uv_x - 0.5*dt*gradp_y(i,j) - v_yr(i,j ) = v_yr(i,j ) - 0.5*dtdx*uv_x - 0.5*dt*gradp_y(i,j) - - ! u du/dx is the transverse term for the u states on y-interfaces - uu_x = ubar*(u_xint(i+1,j) - u_xint(i,j)) - - u_yl(i,j+1) = u_yl(i,j+1) - 0.5*dtdx*uu_x - 0.5*dt*gradp_x(i,j) - u_yr(i,j ) = u_yr(i,j ) - 0.5*dtdx*uu_x - 0.5*dt*gradp_x(i,j) - - enddo - enddo - -end subroutine get_interface_states - - -!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx -subroutine upwind(qx, qy, ng, q_l, q_r, s, q_int) - - ! upwind the left and right states based on the specified input - ! velocity, s. The resulting interface state is q_int - - implicit none - - integer :: qx, qy, ng - double precision :: q_l(0:qx-1, 0:qy-1), q_r(0:qx-1, 0:qy-1) - double precision :: s(0:qx-1, 0:qy-1) - double precision :: q_int(0:qx-1, 0:qy-1) - - integer :: ilo, ihi, jlo, jhi - integer :: nx, ny - integer :: i, j - - nx = qx - 2*ng; ny = qy - 2*ng - ilo = ng; ihi = ng+nx-1; jlo = ng; jhi = ng+ny-1 - - do j = jlo-1, jhi+1 - do i = ilo-1, ihi+1 - - if (s(i,j) > 0.0d0) then - q_int(i,j) = q_l(i,j) - else if (s(i,j) == 0.0d0) then - q_int(i,j) = 0.5d0*(q_l(i,j) + q_r(i,j)) - else - q_int(i,j) = q_r(i,j) - endif - - enddo - enddo - -end subroutine upwind - - -!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx -subroutine riemann(qx, qy, ng, q_l, q_r, s) - - ! Solve the Burger's Riemann problem given the input left and right - ! states and return the state on the interface. - ! - ! This uses the expressions from Almgren, Bell, and Szymczak 1996. - - implicit none - - integer :: qx, qy, ng - double precision :: q_l(0:qx-1, 0:qy-1), q_r(0:qx-1, 0:qy-1) - double precision :: s(0:qx-1, 0:qy-1) - - integer :: ilo, ihi, jlo, jhi - integer :: nx, ny - integer :: i, j - - nx = qx - 2*ng; ny = qy - 2*ng - ilo = ng; ihi = ng+nx-1; jlo = ng; jhi = ng+ny-1 - - do j = jlo-1, jhi+1 - do i = ilo-1, ihi+1 - - if (q_l(i,j) > 0.0d0 .and. q_l(i,j) + q_r(i,j) > 0.0d0) then - s(i,j) = q_l(i,j) - else if (q_l(i,j) <= 0.0d0 .and. q_r(i,j) >= 0.0d0) then - s(i,j) = 0.0d0 - else - s(i,j) = q_r(i,j) - endif - - enddo - enddo - -end subroutine riemann - - -!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx -subroutine riemann_and_upwind(qx, qy, ng, q_l, q_r, q_int) - - ! First solve the Riemann problem given q_l and q_r to give the - ! velocity on the interface and then use this velocity to upwind to - ! determine the state (q_l, q_r, or a mix) on the interface). - ! - ! This differs from upwind, above, in that we don't take in a - ! velocity to upwind with). - - implicit none - - integer :: qx, qy, ng - double precision :: q_l(0:qx-1, 0:qy-1), q_r(0:qx-1, 0:qy-1) - double precision :: q_int(0:qx-1, 0:qy-1) - - double precision :: s(0:qx-1, 0:qy-1) - - call riemann(qx, qy, ng, q_l, q_r, s) - call upwind(qx, qy, ng, q_l, q_r, s, q_int) - -end subroutine riemann_and_upwind diff --git a/incompressible/problems/converge.py b/incompressible/problems/converge.py index ce7602a4f..d6ed8196a 100644 --- a/incompressible/problems/converge.py +++ b/incompressible/problems/converge.py @@ -1,15 +1,22 @@ -""" +r""" Initialize a smooth incompressible convergence test. Here, the velocities are initialized as -u(x,y) = 1 - 2 cos(2 pi x) sin(2 pi y) -v(x,y) = 1 + 2 sin(2 pi x) cos(2 pi y) +.. math:: + + u(x,y) = 1 - 2 \cos(2 \pi x) \sin(2 \pi y) + + v(x,y) = 1 + 2 \sin(2 \pi x) \cos(2 \pi y) and the exact solution at some later time t is then -u(x,y,t) = 1 - 2 cos(2 pi (x - t)) sin(2 pi (y - t)) -v(x,y,t) = 1 + 2 sin(2 pi (x - t)) cos(2 pi (y - t)) -p(x,y,t) = -cos(4 pi (x - t)) - cos(4 pi (y - t)) +.. math:: + + u(x,y,t) = 1 - 2 \cos(2 \pi (x - t)) \sin(2 \pi (y - t)) + + v(x,y,t) = 1 + 2 \sin(2 \pi (x - t)) \cos(2 \pi (y - t)) + + p(x,y,t) = -\cos(4 \pi (x - t)) - \cos(4 \pi (y - t)) The numerical solution can be compared to the exact solution to measure the convergence rate of the algorithm. diff --git a/incompressible/simulation.py b/incompressible/simulation.py index 7a22107e4..18fced6fc 100644 --- a/incompressible/simulation.py +++ b/incompressible/simulation.py @@ -5,7 +5,7 @@ import numpy as np import matplotlib.pyplot as plt -import incompressible.incomp_interface as incomp_interface_f +import incompressible.incomp_interface as incomp_interface import mesh.reconstruction as reconstruction import mesh.patch as patch import mesh.array_indexer as ai @@ -222,8 +222,7 @@ def evolve(self): if self.verbose > 0: print(" making MAC velocities") - _um, _vm = incomp_interface_f.mac_vels(myg.qx, myg.qy, myg.ng, - myg.dx, myg.dy, self.dt, + _um, _vm = incomp_interface.mac_vels(myg.ng, myg.dx, myg.dy, self.dt, u, v, ldelta_ux, ldelta_vx, ldelta_uy, ldelta_vy, @@ -288,8 +287,7 @@ def evolve(self): print(" making u, v edge states") _ux, _vx, _uy, _vy = \ - incomp_interface_f.states(myg.qx, myg.qy, myg.ng, - myg.dx, myg.dy, self.dt, + incomp_interface.states(myg.ng, myg.dx, myg.dy, self.dt, u, v, ldelta_ux, ldelta_vx, ldelta_uy, ldelta_vy, diff --git a/lm_atm/LM_atm_interface.py b/lm_atm/LM_atm_interface.py index 165ff4df0..51e12e7cd 100644 --- a/lm_atm/LM_atm_interface.py +++ b/lm_atm/LM_atm_interface.py @@ -3,14 +3,12 @@ @njit(cache=True) -def is_symmetric_pair(qx, qy, ng, nodal, sl, sr): +def is_symmetric_pair(ng, nodal, sl, sr): r""" Are sl and sr symmetric about an axis parallel with the y-axis in the center of domain the x-direction? Parameters ---------- - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells nodal: bool @@ -24,6 +22,8 @@ def is_symmetric_pair(qx, qy, ng, nodal, sl, sr): Are they symmetric? (1 = yes, 0 = no) """ + qx, qy = sl.shape + nx = qx - 2 * ng ny = qy - 2 * ng ilo = ng @@ -67,14 +67,12 @@ def is_symmetric_pair(qx, qy, ng, nodal, sl, sr): @njit(cache=True) -def is_symmetric(qx, qy, ng, nodal, s): +def is_symmetric(ng, nodal, s): r""" Is the left half of s the mirror image of the right half? Parameters ---------- - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells nodal: bool @@ -88,18 +86,16 @@ def is_symmetric(qx, qy, ng, nodal, s): Is it symmetric? (1 = yes, 0 = no) """ - return is_symmetric_pair(qx, qy, ng, nodal, s, s) + return is_symmetric_pair(ng, nodal, s, s) @njit(cache=True) -def is_asymmetric_pair(qx, qy, ng, nodal, sl, sr): +def is_asymmetric_pair(ng, nodal, sl, sr): r""" Are sl and sr asymmetric about an axis parallel with the y-axis in the center of domain the x-direction? Parameters ---------- - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells nodal: bool @@ -113,6 +109,8 @@ def is_asymmetric_pair(qx, qy, ng, nodal, sl, sr): Are they asymmetric? (1 = yes, 0 = no) """ + qx, qy = sl.shape + nx = qx - 2 * ng ny = qy - 2 * ng ilo = ng @@ -157,14 +155,12 @@ def is_asymmetric_pair(qx, qy, ng, nodal, sl, sr): @njit(cache=True) -def is_asymmetric(qx, qy, ng, nodal, s): +def is_asymmetric(ng, nodal, s): """ Is the left half of s asymmetric to the right half? Parameters ---------- - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells nodal: bool @@ -178,11 +174,11 @@ def is_asymmetric(qx, qy, ng, nodal, s): Is it asymmetric? (1 = yes, 0 = no) """ - return is_asymmetric_pair(qx, qy, ng, nodal, s, s) + return is_asymmetric_pair(ng, nodal, s, s) @njit(cache=True) -def mac_vels(qx, qy, ng, dx, dy, dt, +def mac_vels(ng, dx, dy, dt, u, v, ldelta_ux, ldelta_vx, ldelta_uy, ldelta_vy, @@ -193,8 +189,6 @@ def mac_vels(qx, qy, ng, dx, dy, dt, Parameters ---------- - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells dx, dy : float @@ -246,7 +240,7 @@ def mac_vels(qx, qy, ng, dx, dy, dt, # get the full u and v left and right states (including transverse # terms) on both the x- and y-interfaces - u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = get_interface_states(qx, qy, ng, dx, dy, dt, + u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = get_interface_states(ng, dx, dy, dt, u, v, ldelta_ux, ldelta_vx, ldelta_uy, ldelta_vy, @@ -261,8 +255,8 @@ def mac_vels(qx, qy, ng, dx, dy, dt, # Riemann problem -- this follows Burger's equation. We don't use # any input velocity for the upwinding. Also, we only care about # the normal states here (u on x and v on y) - u_MAC = riemann_and_upwind(qx, qy, ng, u_xl, u_xr) - v_MAC = riemann_and_upwind(qx, qy, ng, v_yl, v_yr) + u_MAC = riemann_and_upwind(ng, u_xl, u_xr) + v_MAC = riemann_and_upwind(ng, v_yl, v_yr) # print *, 'checking U_MAC' # if (not is_asymmetric(qx, qy, ng, .true., u_MAC) == 1): @@ -274,7 +268,7 @@ def mac_vels(qx, qy, ng, dx, dy, dt, # xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx @njit(cache=True) -def states(qx, qy, ng, dx, dy, dt, +def states(ng, dx, dy, dt, u, v, ldelta_ux, ldelta_vx, ldelta_uy, ldelta_vy, @@ -288,8 +282,6 @@ def states(qx, qy, ng, dx, dy, dt, Parameters ---------- - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells dx, dy : float @@ -317,7 +309,7 @@ def states(qx, qy, ng, dx, dy, dt, # get the full u and v left and right states (including transverse # terms) on both the x- and y-interfaces - u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = get_interface_states(qx, qy, ng, dx, dy, dt, + u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = get_interface_states(ng, dx, dy, dt, u, v, ldelta_ux, ldelta_vx, ldelta_uy, ldelta_vy, @@ -326,17 +318,17 @@ def states(qx, qy, ng, dx, dy, dt, # upwind using the MAC velocity to determine which state exists on # the interface - u_xint = upwind(qx, qy, ng, u_xl, u_xr, u_MAC) - v_xint = upwind(qx, qy, ng, v_xl, v_xr, u_MAC) - u_yint = upwind(qx, qy, ng, u_yl, u_yr, v_MAC) - v_yint = upwind(qx, qy, ng, v_yl, v_yr, v_MAC) + u_xint = upwind(ng, u_xl, u_xr, u_MAC) + v_xint = upwind(ng, v_xl, v_xr, u_MAC) + u_yint = upwind(ng, u_yl, u_yr, v_MAC) + v_yint = upwind(ng, v_yl, v_yr, v_MAC) return u_xint, v_xint, u_yint, v_yint # xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx @njit(cache=True) -def rho_states(qx, qy, ng, dx, dy, dt, +def rho_states(ng, dx, dy, dt, rho, u_MAC, v_MAC, ldelta_rx, ldelta_ry): r""" @@ -345,8 +337,6 @@ def rho_states(qx, qy, ng, dx, dy, dt, Parameters ---------- - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells dx, dy : float @@ -364,6 +354,8 @@ def rho_states(qx, qy, ng, dx, dy, dt, rho predicted to the interfaces """ + qx, qy = rho.shape + rho_xl = np.zeros((qx, qy)) rho_xr = np.zeros((qx, qy)) rho_yl = np.zeros((qx, qy)) @@ -395,8 +387,8 @@ def rho_states(qx, qy, ng, dx, dy, dt, (1.0 + dtdy * v_MAC[i, j]) * ldelta_ry[i, j] # we upwind based on the MAC velocities - rho_xint = upwind(qx, qy, ng, rho_xl, rho_xr, u_MAC) - rho_yint = upwind(qx, qy, ng, rho_yl, rho_yr, v_MAC) + rho_xint = upwind(ng, rho_xl, rho_xr, u_MAC) + rho_yint = upwind(ng, rho_yl, rho_yr, v_MAC) # now add the transverse term and the non-advective part of the normal # divergence @@ -425,15 +417,15 @@ def rho_states(qx, qy, ng, dx, dy, dt, rho_yr[i, j] = rho_yr[i, j] - 0.5 * dt * (rhou_x + rho[i, j] * v_y) # finally upwind the full states - rho_xint = upwind(qx, qy, ng, rho_xl, rho_xr, u_MAC) - rho_yint = upwind(qx, qy, ng, rho_yl, rho_yr, v_MAC) + rho_xint = upwind(ng, rho_xl, rho_xr, u_MAC) + rho_yint = upwind(ng, rho_yl, rho_yr, v_MAC) return rho_xint, rho_yint # xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx @njit(cache=True) -def get_interface_states(qx, qy, ng, dx, dy, dt, +def get_interface_states(ng, dx, dy, dt, u, v, ldelta_ux, ldelta_vx, ldelta_uy, ldelta_vy, @@ -448,8 +440,6 @@ def get_interface_states(qx, qy, ng, dx, dy, dt, Parameters ---------- - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells dx, dy : float @@ -474,6 +464,8 @@ def get_interface_states(qx, qy, ng, dx, dy, dt, y-interfaces """ + qx, qy = u.shape + u_xl = np.zeros((qx, qy)) u_xr = np.zeros((qx, qy)) u_yl = np.zeros((qx, qy)) @@ -531,19 +523,19 @@ def get_interface_states(qx, qy, ng, dx, dy, dt, # now get the normal advective velocities on the interfaces by solving # the Riemann problem. - uhat_adv = riemann(qx, qy, ng, u_xl, u_xr) - vhat_adv = riemann(qx, qy, ng, v_yl, v_yr) + uhat_adv = riemann(ng, u_xl, u_xr) + vhat_adv = riemann(ng, v_yl, v_yr) # now that we have the advective velocities, upwind the left and right # states using the appropriate advective velocity. # on the x-interfaces, we upwind based on uhat_adv - u_xint = upwind(qx, qy, ng, u_xl, u_xr, uhat_adv) - v_xint = upwind(qx, qy, ng, v_xl, v_xr, uhat_adv) + u_xint = upwind(ng, u_xl, u_xr, uhat_adv) + v_xint = upwind(ng, v_xl, v_xr, uhat_adv) # on the y-interfaces, we upwind based on vhat_adv - u_yint = upwind(qx, qy, ng, u_yl, u_yr, vhat_adv) - v_yint = upwind(qx, qy, ng, v_yl, v_yr, vhat_adv) + u_yint = upwind(ng, u_yl, u_yr, vhat_adv) + v_yint = upwind(ng, v_yl, v_yr, vhat_adv) # at this point, these states are the `hat' states -- they only # considered the normal to the interface portion of the predictor. @@ -592,15 +584,13 @@ def get_interface_states(qx, qy, ng, dx, dy, dt, # xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx @njit(cache=True) -def upwind(qx, qy, ng, q_l, q_r, s): +def upwind(ng, q_l, q_r, s): r""" upwind the left and right states based on the specified input velocity, s. The resulting interface state is q_int Parameters ---------- - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells q_l, q_r : ndarray @@ -614,7 +604,9 @@ def upwind(qx, qy, ng, q_l, q_r, s): Upwinded state """ - q_int = np.zeros((qx, qy)) + qx, qy = s.shape + + q_int = np.zeros_like(s) nx = qx - 2 * ng ny = qy - 2 * ng @@ -638,7 +630,7 @@ def upwind(qx, qy, ng, q_l, q_r, s): # xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx @njit(cache=True) -def riemann(qx, qy, ng, q_l, q_r): +def riemann(ng, q_l, q_r): """ Solve the Burger's Riemann problem given the input left and right states and return the state on the interface. @@ -647,8 +639,6 @@ def riemann(qx, qy, ng, q_l, q_r): Parameters ---------- - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells q_l, q_r : ndarray @@ -660,6 +650,8 @@ def riemann(qx, qy, ng, q_l, q_r): Interface state """ + qx, qy = q_l.shape + nx = qx - 2 * ng ny = qy - 2 * ng ilo = ng @@ -684,7 +676,7 @@ def riemann(qx, qy, ng, q_l, q_r): # xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx @njit(cache=True) -def riemann_and_upwind(qx, qy, ng, q_l, q_r): +def riemann_and_upwind(ng, q_l, q_r): r""" First solve the Riemann problem given q_l and q_r to give the velocity on the interface and: use this velocity to upwind to @@ -695,8 +687,6 @@ def riemann_and_upwind(qx, qy, ng, q_l, q_r): Parameters ---------- - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells q_l, q_r : ndarray @@ -708,5 +698,5 @@ def riemann_and_upwind(qx, qy, ng, q_l, q_r): Upwinded state """ - s = riemann(qx, qy, ng, q_l, q_r) - return upwind(qx, qy, ng, q_l, q_r, s) + s = riemann(ng, q_l, q_r) + return upwind(ng, q_l, q_r, s) diff --git a/lm_atm/LM_atm_interface_f.f90 b/lm_atm/LM_atm_interface_f.f90 deleted file mode 100644 index dd2500d06..000000000 --- a/lm_atm/LM_atm_interface_f.f90 +++ /dev/null @@ -1,670 +0,0 @@ -function is_symmetric_pair(qx, qy, ng, nodal, sl, sr) result (sym) - - implicit none - - integer :: sym - integer :: qx, qy, ng - double precision :: sl(0:qx-1, 0:qy-1) - double precision :: sr(0:qx-1, 0:qy-1) - logical :: nodal - - integer :: i, j - integer :: nx, ny, ilo, ihi, jlo, jhi - integer :: il, ir - - nx = qx - 2*ng; ny = qy - 2*ng - ilo = ng; ihi = ng+nx-1; jlo = ng; jhi = ng+ny-1 - - sym = 1 - - if (.not. nodal) then - loop_outer: do j = jlo, jhi - do i = 1, nx/2 - il = ilo + i - 1 - ir = ihi - i + 1 - - if (.not. sl(il,j) == sr(ir,j)) then - sym = 0 - exit loop_outer - endif - enddo - enddo loop_outer - - else - - loop_nodal_outer: do j = jlo, jhi - do i = 1, nx/2 - il = ilo + i - 1 - ir = ihi - i + 2 - - if (.not. sl(il,j) == sr(ir,j)) then - sym = 0 - exit loop_nodal_outer - endif - enddo - enddo loop_nodal_outer - endif - -end function is_symmetric_pair - -function is_symmetric(qx, qy, ng, nodal, s) result (sym) - - implicit none - - integer :: sym - integer :: qx, qy, ng - double precision :: s(0:qx-1, 0:qy-1) - logical :: nodal - - integer :: is_symmetric_pair - - sym = is_symmetric_pair(qx, qy, ng, nodal, s, s) - -end function is_symmetric - - -function is_asymmetric_pair(qx, qy, ng, nodal, sl, sr) result (asym) - - implicit none - - integer :: asym - integer :: qx, qy, ng - double precision :: sl(0:qx-1, 0:qy-1) - double precision :: sr(0:qx-1, 0:qy-1) - logical :: nodal - - integer :: i, j - integer :: nx, ny, ilo, ihi, jlo, jhi - integer :: il, ir - - nx = qx - 2*ng; ny = qy - 2*ng - ilo = ng; ihi = ng+nx-1; jlo = ng; jhi = ng+ny-1 - - asym = 1 - - if (.not. nodal) then - loop_aouter: do j = jlo, jhi - do i = 1, nx/2 - il = ilo + i - 1 - ir = ihi - i + 1 - - print *, il, ir, sl(il,j), -sr(ir,j) - if (.not. sl(il,j) == -sr(ir,j)) then - asym = 0 - exit loop_aouter - endif - enddo - enddo loop_aouter - - else - - loop_nodal_aouter: do j = jlo, jhi - do i = 1, nx/2 - il = ilo + i - 1 - ir = ihi - i + 2 - - !print *, il, ir, sl(il,j), -sr(ir,j) - if (.not. sl(il,j) == -sr(ir,j)) then - asym = 0 - exit loop_nodal_aouter - endif - enddo - enddo loop_nodal_aouter - endif -end function is_asymmetric_pair - -function is_asymmetric(qx, qy, ng, nodal, s) result (asym) - - implicit none - - integer :: asym - integer :: qx, qy, ng - double precision :: s(0:qx-1, 0:qy-1) - logical :: nodal - - integer :: is_asymmetric_pair - - asym = is_asymmetric_pair(qx, qy, ng, nodal, s, s) - -end function is_asymmetric - -subroutine mac_vels(qx, qy, ng, dx, dy, dt, & - u, v, & - ldelta_ux, ldelta_vx, & - ldelta_uy, ldelta_vy, & - gradp_x, gradp_y, & - source, & - u_MAC, v_MAC) - - implicit none - - integer, intent(in) :: qx, qy, ng - double precision, intent(in) :: dx, dy, dt - - ! 0-based indexing to match python - double precision, intent(inout) :: u(0:qx-1, 0:qy-1) - double precision, intent(inout) :: v(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: ldelta_ux(0:qx-1, 0:qy-1) - double precision, intent(inout) :: ldelta_vx(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: ldelta_uy(0:qx-1, 0:qy-1) - double precision, intent(inout) :: ldelta_vy(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: gradp_x(0:qx-1, 0:qy-1) - double precision, intent(inout) :: gradp_y(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: source(0:qx-1, 0:qy-1) - - double precision, intent( out) :: u_MAC(0:qx-1, 0:qy-1) - double precision, intent( out) :: v_MAC(0:qx-1, 0:qy-1) - - integer :: is_symmetric, is_asymmetric, is_asymmetric_pair - -!f2py depend(qx, qy) :: u, v -!f2py depend(qx, qy) :: ldelta_ux, ldelta_vx, ldelta_uy, ldelta_vy -!f2py depend(qx, qy) :: gradp_x, gradp_y -!f2py depend(qx, qy) :: source -!f2py depend(qx, qy) :: u_MAC, v_MAC -!f2py intent(in) :: u, v, gradp_x, gradp_y -!f2py intent(in) :: ldelta_ux, ldelta_vx, ldelta_uy, ldelta_vy -!f2py intent(in) :: source -!f2py intent(out) :: u_MAC, v_MAC - - - double precision :: u_xl(0:qx-1, 0:qy-1), u_xr(0:qx-1, 0:qy-1) - double precision :: u_yl(0:qx-1, 0:qy-1), u_yr(0:qx-1, 0:qy-1) - double precision :: v_xl(0:qx-1, 0:qy-1), v_xr(0:qx-1, 0:qy-1) - double precision :: v_yl(0:qx-1, 0:qy-1), v_yr(0:qx-1, 0:qy-1) - - integer :: i, j - integer :: nx, ny, ilo, ihi, jlo, jhi - nx = qx - 2*ng; ny = qy - 2*ng - ilo = ng; ihi = ng+nx-1; jlo = ng; jhi = ng+ny-1 - - - ! assertions - ! print *, "checking ldelta_ux" - ! if (.not. is_asymmetric(qx, qy, ng, .false., ldelta_ux) == 1) then - ! stop 'ldelta_ux not asymmetric' - ! endif - - ! print *, "checking ldelta_uy" - ! if (.not. is_symmetric(qx, qy, ng, .false., ldelta_uy) == 1) then - ! stop 'ldelta_uy not symmetric' - ! endif - - ! print *, "checking ldelta_vx" - ! if (.not. is_symmetric(qx, qy, ng, .false., ldelta_vx) == 1) then - ! stop 'ldelta_vx not symmetric' - ! endif - - ! print *, "checking ldelta_vy" - ! if (.not. is_symmetric(qx, qy, ng, .false., ldelta_vy) == 1) then - ! stop 'ldelta_vy not symmetric' - ! endif - - ! print *, "checking gradp_x" - ! if (.not. is_asymmetric(qx, qy, ng, .false., gradp_x) == 1) then - ! stop 'gradp_x not asymmetric' - ! endif - - - ! get the full u and v left and right states (including transverse - ! terms) on both the x- and y-interfaces - call get_interface_states(qx, qy, ng, dx, dy, dt, & - u, v, & - ldelta_ux, ldelta_vx, & - ldelta_uy, ldelta_vy, & - gradp_x, gradp_y, & - source, & - u_xl, u_xr, u_yl, u_yr, & - v_xl, v_xr, v_yl, v_yr) - - ! print *, 'checking u_xl' - ! if (.not. is_asymmetric_pair(qx, qy, ng, .true., u_xl, u_xr) == 1) then - ! stop 'u_xl/r not asymmetric' - ! endif - - - ! Riemann problem -- this follows Burger's equation. We don't use - ! any input velocity for the upwinding. Also, we only care about - ! the normal states here (u on x and v on y) - call riemann_and_upwind(qx, qy, ng, u_xl, u_xr, u_MAC) - call riemann_and_upwind(qx, qy, ng, v_yl, v_yr, v_MAC) - - ! print *, 'checking U_MAC' - ! if (.not. is_asymmetric(qx, qy, ng, .true., u_MAC) == 1) then - ! stop 'u_MAC not asymmetric' - ! endif - -end subroutine mac_vels - - -!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx -subroutine states(qx, qy, ng, dx, dy, dt, & - u, v, & - ldelta_ux, ldelta_vx, & - ldelta_uy, ldelta_vy, & - gradp_x, gradp_y, & - source, & - u_MAC, v_MAC, & - u_xint, v_xint, u_yint, v_yint) - - ! this is similar to mac_vels, but it predicts the interface states - ! of both u and v on both interfaces, using the MAC velocities to - ! do the upwinding. - - implicit none - - integer, intent(in) :: qx, qy, ng - double precision, intent(in) :: dx, dy, dt - - ! 0-based indexing to match python - double precision, intent(inout) :: u(0:qx-1, 0:qy-1) - double precision, intent(inout) :: v(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: ldelta_ux(0:qx-1, 0:qy-1) - double precision, intent(inout) :: ldelta_vx(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: ldelta_uy(0:qx-1, 0:qy-1) - double precision, intent(inout) :: ldelta_vy(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: gradp_x(0:qx-1, 0:qy-1) - double precision, intent(inout) :: gradp_y(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: source(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: u_MAC(0:qx-1, 0:qy-1) - double precision, intent(inout) :: v_MAC(0:qx-1, 0:qy-1) - - double precision, intent(out) :: u_xint(0:qx-1, 0:qy-1), u_yint(0:qx-1, 0:qy-1) - double precision, intent(out) :: v_xint(0:qx-1, 0:qy-1), v_yint(0:qx-1, 0:qy-1) - -!f2py depend(qx, qy) :: u, v -!f2py depend(qx, qy) :: ldelta_ux, ldelta_vx, ldelta_uy, ldelta_vy -!f2py depend(qx, qy) :: gradp_x, gradp_y -!f2py depend(qx, qy) :: source -!f2py depend(qx, qy) :: u_MAC, v_MAC -!f2py intent(in) :: u, v, gradp_x, gradp_y -!f2py intent(in) :: ldelta_ux, ldelta_vx, ldelta_uy, ldelta_vy -!f2py intent(in) :: source -!f2py intent(in) :: u_MAC, v_MAC -!f2py intent(out) :: u_xint, v_xint, u_yint, v_yint - - double precision :: u_xl(0:qx-1, 0:qy-1), u_xr(0:qx-1, 0:qy-1) - double precision :: u_yl(0:qx-1, 0:qy-1), u_yr(0:qx-1, 0:qy-1) - double precision :: v_xl(0:qx-1, 0:qy-1), v_xr(0:qx-1, 0:qy-1) - double precision :: v_yl(0:qx-1, 0:qy-1), v_yr(0:qx-1, 0:qy-1) - - - ! get the full u and v left and right states (including transverse - ! terms) on both the x- and y-interfaces - call get_interface_states(qx, qy, ng, dx, dy, dt, & - u, v, & - ldelta_ux, ldelta_vx, & - ldelta_uy, ldelta_vy, & - gradp_x, gradp_y, & - source, & - u_xl, u_xr, u_yl, u_yr, & - v_xl, v_xr, v_yl, v_yr) - - - ! upwind using the MAC velocity to determine which state exists on - ! the interface - call upwind(qx, qy, ng, u_xl, u_xr, u_MAC, u_xint) - call upwind(qx, qy, ng, v_xl, v_xr, u_MAC, v_xint) - call upwind(qx, qy, ng, u_yl, u_yr, v_MAC, u_yint) - call upwind(qx, qy, ng, v_yl, v_yr, v_MAC, v_yint) - -end subroutine states - - -!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx -subroutine rho_states(qx, qy, ng, dx, dy, dt, & - rho, u_MAC, v_MAC, & - ldelta_rx, ldelta_ry, & - rho_xint, rho_yint) - - ! this predicts rho to the interfaces. We use the MAC velocities to do - ! the upwinding - - implicit none - - integer, intent(in) :: qx, qy, ng - double precision, intent(in) :: dx, dy, dt - - ! 0-based indexing to match python - double precision, intent(inout) :: rho(0:qx-1, 0:qy-1) - double precision, intent(inout) :: u_MAC(0:qx-1, 0:qy-1) - double precision, intent(inout) :: v_MAC(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: ldelta_rx(0:qx-1, 0:qy-1) - double precision, intent(inout) :: ldelta_ry(0:qx-1, 0:qy-1) - - double precision, intent(out) :: rho_xint(0:qx-1, 0:qy-1), rho_yint(0:qx-1, 0:qy-1) - -!f2py depend(qx, qy) :: rho, u_MAC, v_MAC -!f2py depend(qx, qy) :: ldelta_rx, ldelta_ry -!f2py intent(in) :: rho, u_MAC, v_MAC -!f2py intent(in) :: ldelta_rx, ldelta_ry -!f2py intent(out) :: rho_xint, rho_yint - - double precision :: rho_xl(0:qx-1, 0:qy-1), rho_xr(0:qx-1, 0:qy-1) - double precision :: rho_yl(0:qx-1, 0:qy-1), rho_yr(0:qx-1, 0:qy-1) - - integer :: ilo, ihi, jlo, jhi - integer :: nx, ny - integer :: i, j - - double precision :: dtdx, dtdy - double precision :: u_x, v_y, rhov_y, rhou_x - - nx = qx - 2*ng; ny = qy - 2*ng - ilo = ng; ihi = ng+nx-1; jlo = ng; jhi = ng+ny-1 - - dtdx = dt/dx - dtdy = dt/dy - - do j = jlo-2, jhi+2 - do i = ilo-2, ihi+2 - - ! u on x-edges - rho_xl(i+1,j) = rho(i,j) + 0.5d0*(1.0d0 - dtdx*u_MAC(i+1,j))*ldelta_rx(i,j) - rho_xr(i ,j) = rho(i,j) - 0.5d0*(1.0d0 + dtdx*u_MAC(i,j))*ldelta_rx(i,j) - - ! u on y-edges - rho_yl(i,j+1) = rho(i,j) + 0.5d0*(1.0d0 - dtdy*v_MAC(i,j+1))*ldelta_ry(i,j) - rho_yr(i,j ) = rho(i,j) - 0.5d0*(1.0d0 + dtdy*v_MAC(i,j))*ldelta_ry(i,j) - - enddo - enddo - - - ! we upwind based on the MAC velocities - call upwind(qx, qy, ng, rho_xl, rho_xr, u_MAC, rho_xint) - call upwind(qx, qy, ng, rho_yl, rho_yr, v_MAC, rho_yint) - - - ! now add the transverse term and the non-advective part of the normal - ! divergence - do j = jlo-2, jhi+2 - do i = ilo-2, ihi+2 - - u_x = (u_MAC(i+1,j) - u_MAC(i,j))/dx - v_y = (v_MAC(i,j+1) - v_MAC(i,j))/dy - - ! (rho v)_y is the transverse term for the x-interfaces - ! rho u_x is the non-advective piece for the x-interfaces - rhov_y = (rho_yint(i,j+1)*v_MAC(i,j+1) - rho_yint(i,j)*v_MAC(i,j))/dy - - rho_xl(i+1,j) = rho_xl(i+1,j) - 0.5*dt*(rhov_y + rho(i,j)*u_x) - rho_xr(i ,j) = rho_xr(i ,j) - 0.5*dt*(rhov_y + rho(i,j)*u_x) - - ! (rho u)_x is the transverse term for the y-interfaces - ! rho v_y is the non-advective piece for the y-interfaces - rhou_x = (rho_xint(i+1,j)*u_MAC(i+1,j) - rho_xint(i,j)*u_MAC(i,j))/dx - - rho_yl(i,j+1) = rho_yl(i,j+1) - 0.5*dt*(rhou_x + rho(i,j)*v_y) - rho_yr(i,j ) = rho_yr(i,j ) - 0.5*dt*(rhou_x + rho(i,j)*v_y) - - enddo - enddo - - ! finally upwind the full states - call upwind(qx, qy, ng, rho_xl, rho_xr, u_MAC, rho_xint) - call upwind(qx, qy, ng, rho_yl, rho_yr, v_MAC, rho_yint) - -end subroutine rho_states - - -!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx -subroutine get_interface_states(qx, qy, ng, dx, dy, dt, & - u, v, & - ldelta_ux, ldelta_vx, & - ldelta_uy, ldelta_vy, & - gradp_x, gradp_y, & - source, & - u_xl, u_xr, u_yl, u_yr, & - v_xl, v_xr, v_yl, v_yr) - - ! Compute the unsplit predictions of u and v on both the x- and - ! y-interfaces. This includes the transverse terms. - - ! note that the gradp_x, gradp_y should have any coefficients - ! already included (e.g. beta_0/rho) - - implicit none - - integer, intent(in) :: qx, qy, ng - double precision, intent(in) :: dx, dy, dt - - ! 0-based indexing to match python - double precision, intent(inout) :: u(0:qx-1, 0:qy-1) - double precision, intent(inout) :: v(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: ldelta_ux(0:qx-1, 0:qy-1) - double precision, intent(inout) :: ldelta_vx(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: ldelta_uy(0:qx-1, 0:qy-1) - double precision, intent(inout) :: ldelta_vy(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: gradp_x(0:qx-1, 0:qy-1) - double precision, intent(inout) :: gradp_y(0:qx-1, 0:qy-1) - - double precision, intent(inout) :: source(0:qx-1, 0:qy-1) - - double precision, intent(out) :: u_xl(0:qx-1, 0:qy-1), u_xr(0:qx-1, 0:qy-1) - double precision, intent(out) :: u_yl(0:qx-1, 0:qy-1), u_yr(0:qx-1, 0:qy-1) - double precision, intent(out) :: v_xl(0:qx-1, 0:qy-1), v_xr(0:qx-1, 0:qy-1) - double precision, intent(out) :: v_yl(0:qx-1, 0:qy-1), v_yr(0:qx-1, 0:qy-1) - - integer :: ilo, ihi, jlo, jhi - integer :: nx, ny - integer :: i, j - - double precision :: uhat_adv(0:qx-1, 0:qy-1), vhat_adv(0:qx-1, 0:qy-1) - - double precision :: u_xint(0:qx-1, 0:qy-1), u_yint(0:qx-1, 0:qy-1) - double precision :: v_xint(0:qx-1, 0:qy-1), v_yint(0:qx-1, 0:qy-1) - - double precision :: dtdx, dtdy - double precision :: ubar, vbar, uv_x, vu_y, uu_x, vv_y - - integer :: is_asymmetric_pair - - nx = qx - 2*ng; ny = qy - 2*ng - ilo = ng; ihi = ng+nx-1; jlo = ng; jhi = ng+ny-1 - - - ! first predict u and v to both interfaces, considering only the normal - ! part of the predictor. These are the 'hat' states. - - - dtdx = dt/dx - dtdy = dt/dy - - do j = jlo-2, jhi+2 - do i = ilo-2, ihi+2 - - ! u on x-edges - u_xl(i+1,j) = u(i,j) + 0.5d0*(1.0d0 - dtdx*u(i,j))*ldelta_ux(i,j) - u_xr(i ,j) = u(i,j) - 0.5d0*(1.0d0 + dtdx*u(i,j))*ldelta_ux(i,j) - - ! v on x-edges - v_xl(i+1,j) = v(i,j) + 0.5d0*(1.0d0 - dtdx*u(i,j))*ldelta_vx(i,j) - v_xr(i ,j) = v(i,j) - 0.5d0*(1.0d0 + dtdx*u(i,j))*ldelta_vx(i,j) - - ! u on y-edges - u_yl(i,j+1) = u(i,j) + 0.5d0*(1.0d0 - dtdy*v(i,j))*ldelta_uy(i,j) - u_yr(i,j ) = u(i,j) - 0.5d0*(1.0d0 + dtdy*v(i,j))*ldelta_uy(i,j) - - ! v on y-edges - v_yl(i,j+1) = v(i,j) + 0.5d0*(1.0d0 - dtdy*v(i,j))*ldelta_vy(i,j) - v_yr(i,j ) = v(i,j) - 0.5d0*(1.0d0 + dtdy*v(i,j))*ldelta_vy(i,j) - - enddo - enddo - - ! print *, 'checking u_xl in states' - ! if (.not. is_asymmetric_pair(qx, qy, ng, .true., u_xl, u_xr) == 1) then - ! stop 'u_xl/r not asymmetric' - ! endif - - - ! now get the normal advective velocities on the interfaces by solving - ! the Riemann problem. - call riemann(qx, qy, ng, u_xl, u_xr, uhat_adv) - call riemann(qx, qy, ng, v_yl, v_yr, vhat_adv) - - - - ! now that we have the advective velocities, upwind the left and right - ! states using the appropriate advective velocity. - - ! on the x-interfaces, we upwind based on uhat_adv - call upwind(qx, qy, ng, u_xl, u_xr, uhat_adv, u_xint) - call upwind(qx, qy, ng, v_xl, v_xr, uhat_adv, v_xint) - - ! on the y-interfaces, we upwind based on vhat_adv - call upwind(qx, qy, ng, u_yl, u_yr, vhat_adv, u_yint) - call upwind(qx, qy, ng, v_yl, v_yr, vhat_adv, v_yint) - - ! at this point, these states are the `hat' states -- they only - ! considered the normal to the interface portion of the predictor. - - - ! add the transverse flux differences to the preliminary interface states - do j = jlo-1, jhi+1 - do i = ilo-1, ihi+1 - - ubar = 0.5d0*(uhat_adv(i,j) + uhat_adv(i+1,j)) - vbar = 0.5d0*(vhat_adv(i,j) + vhat_adv(i,j+1)) - - ! v du/dy is the transerse term for the u states on x-interfaces - vu_y = vbar*(u_yint(i,j+1) - u_yint(i,j)) - - u_xl(i+1,j) = u_xl(i+1,j) - 0.5*dtdy*vu_y - 0.5*dt*gradp_x(i,j) - u_xr(i ,j) = u_xr(i ,j) - 0.5*dtdy*vu_y - 0.5*dt*gradp_x(i,j) - - ! v dv/dy is the transverse term for the v states on x-interfaces - vv_y = vbar*(v_yint(i,j+1) - v_yint(i,j)) - - v_xl(i+1,j) = v_xl(i+1,j) - 0.5*dtdy*vv_y - 0.5*dt*gradp_y(i,j) + 0.5d0*dt*source(i,j) - v_xr(i ,j) = v_xr(i ,j) - 0.5*dtdy*vv_y - 0.5*dt*gradp_y(i,j) + 0.5d0*dt*source(i,j) - - ! u dv/dx is the transverse term for the v states on y-interfaces - uv_x = ubar*(v_xint(i+1,j) - v_xint(i,j)) - - v_yl(i,j+1) = v_yl(i,j+1) - 0.5*dtdx*uv_x - 0.5*dt*gradp_y(i,j) + 0.5d0*dt*source(i,j) - v_yr(i,j ) = v_yr(i,j ) - 0.5*dtdx*uv_x - 0.5*dt*gradp_y(i,j) + 0.5d0*dt*source(i,j) - - ! u du/dx is the transverse term for the u states on y-interfaces - uu_x = ubar*(u_xint(i+1,j) - u_xint(i,j)) - - u_yl(i,j+1) = u_yl(i,j+1) - 0.5*dtdx*uu_x - 0.5*dt*gradp_x(i,j) - u_yr(i,j ) = u_yr(i,j ) - 0.5*dtdx*uu_x - 0.5*dt*gradp_x(i,j) - - enddo - enddo - -end subroutine get_interface_states - - -!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx -subroutine upwind(qx, qy, ng, q_l, q_r, s, q_int) - - ! upwind the left and right states based on the specified input - ! velocity, s. The resulting interface state is q_int - - implicit none - - integer :: qx, qy, ng - double precision :: q_l(0:qx-1, 0:qy-1), q_r(0:qx-1, 0:qy-1) - double precision :: s(0:qx-1, 0:qy-1) - double precision :: q_int(0:qx-1, 0:qy-1) - - integer :: ilo, ihi, jlo, jhi - integer :: nx, ny - integer :: i, j - - nx = qx - 2*ng; ny = qy - 2*ng - ilo = ng; ihi = ng+nx-1; jlo = ng; jhi = ng+ny-1 - - do j = jlo-1, jhi+2 - do i = ilo-1, ihi+2 - - if (s(i,j) > 0.0d0) then - q_int(i,j) = q_l(i,j) - else if (s(i,j) == 0.0d0) then - q_int(i,j) = 0.5d0*(q_l(i,j) + q_r(i,j)) - else - q_int(i,j) = q_r(i,j) - endif - - enddo - enddo - -end subroutine upwind - - -!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx -subroutine riemann(qx, qy, ng, q_l, q_r, s) - - ! Solve the Burger's Riemann problem given the input left and right - ! states and return the state on the interface. - ! - ! This uses the expressions from Almgren, Bell, and Szymczak 1996. - - implicit none - - integer :: qx, qy, ng - double precision :: q_l(0:qx-1, 0:qy-1), q_r(0:qx-1, 0:qy-1) - double precision :: s(0:qx-1, 0:qy-1) - - integer :: ilo, ihi, jlo, jhi - integer :: nx, ny - integer :: i, j - - nx = qx - 2*ng; ny = qy - 2*ng - ilo = ng; ihi = ng+nx-1; jlo = ng; jhi = ng+ny-1 - - do j = jlo-1, jhi+2 - do i = ilo-1, ihi+2 - - if (q_l(i,j) > 0.0d0 .and. q_l(i,j) + q_r(i,j) > 0.0d0) then - s(i,j) = q_l(i,j) - else if (q_l(i,j) <= 0.0d0 .and. q_r(i,j) >= 0.0d0) then - s(i,j) = 0.0d0 - else - s(i,j) = q_r(i,j) - endif - - enddo - enddo - -end subroutine riemann - - -!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx -subroutine riemann_and_upwind(qx, qy, ng, q_l, q_r, q_int) - - ! First solve the Riemann problem given q_l and q_r to give the - ! velocity on the interface and then use this velocity to upwind to - ! determine the state (q_l, q_r, or a mix) on the interface). - ! - ! This differs from upwind, above, in that we don't take in a - ! velocity to upwind with). - - implicit none - - integer :: qx, qy, ng - double precision :: q_l(0:qx-1, 0:qy-1), q_r(0:qx-1, 0:qy-1) - double precision :: q_int(0:qx-1, 0:qy-1) - - double precision :: s(0:qx-1, 0:qy-1) - - call riemann(qx, qy, ng, q_l, q_r, s) - call upwind(qx, qy, ng, q_l, q_r, s, q_int) - -end subroutine riemann_and_upwind diff --git a/lm_atm/simulation.py b/lm_atm/simulation.py index d0f1917f4..5e18ecae4 100644 --- a/lm_atm/simulation.py +++ b/lm_atm/simulation.py @@ -5,7 +5,7 @@ import numpy as np import matplotlib.pyplot as plt -import lm_atm.LM_atm_interface as lm_interface_f +import lm_atm.LM_atm_interface as lm_interface import mesh.reconstruction as reconstruction import mesh.boundary as bnd import mesh.patch as patch @@ -362,8 +362,7 @@ def evolve(self): source.v()[:, :] = rhoprime.v()*g/rho.v() self.aux_data.fill_BC("source_y") - _um, _vm = lm_interface_f.mac_vels(myg.qx, myg.qy, myg.ng, - myg.dx, myg.dy, self.dt, + _um, _vm = lm_interface.mac_vels(myg.ng, myg.dx, myg.dy, self.dt, u, v, ldelta_ux, ldelta_vx, ldelta_uy, ldelta_vy, @@ -447,8 +446,7 @@ def evolve(self): # --------------------------------------------------------------------- # predict rho to the edges and do its conservative update # --------------------------------------------------------------------- - _rx, _ry = lm_interface_f.rho_states(myg.qx, myg.qy, myg.ng, - myg.dx, myg.dy, self.dt, + _rx, _ry = lm_interface.rho_states(myg.ng, myg.dx, myg.dy, self.dt, rho, u_MAC, v_MAC, ldelta_rx, ldelta_ry) @@ -483,8 +481,7 @@ def evolve(self): self.aux_data.fill_BC("coeff") _ux, _vx, _uy, _vy = \ - lm_interface_f.states(myg.qx, myg.qy, myg.ng, - myg.dx, myg.dy, self.dt, + lm_interface.states(myg.ng, myg.dx, myg.dy, self.dt, u, v, ldelta_ux, ldelta_vx, ldelta_uy, ldelta_vy, diff --git a/mk.sh b/mk.sh index d19e9c460..40a56ac04 100755 --- a/mk.sh +++ b/mk.sh @@ -14,46 +14,7 @@ if [ "$1" == "clean" ]; then - rm -rf advection_fv4/*.so - rm -rf compressible/*.so - rm -rf incompressible/*.so - rm -rf lm_atm/*.so - rm -rf mesh/*.so - rm -rf swe/*.so - find . -name "*.pyc" -exec rm -f {} \; find . -type d -name "__pycache__" -exec rm -rf {} \; - find . -type d -name "build" -exec rm -rf {} \; - - # move numba interface files back - regex="([a-z0-9A-Z_/]*/)_([a-z0-9A-Z_]*interface.py)" - for f in $(find . -name "*interface.py") - do - if [[ $f =~ $regex ]] - then - mv ".${BASH_REMATCH[1]}_${BASH_REMATCH[2]}" ".${BASH_REMATCH[1]}${BASH_REMATCH[2]}" - fi - done - -elif [ "$1" == "fortran" ]; then - - if [ "$1" == "debug" ]; then - FFLAGS="-fbounds-check -fbacktrace -Wuninitialized -Wunused -ffpe-trap=invalid -finit-real=snan" - else - FFLAGS="-C" - fi - - # move numba interface files out of the way - regex="([a-z0-9A-Z_/]*/)([a-z0-9A-Z_]*interface.py)" - for f in $(find . -name "*interface.py") - do - if [[ $f =~ $regex ]] - then - mv ".${BASH_REMATCH[1]}${BASH_REMATCH[2]}" ".${BASH_REMATCH[1]}_${BASH_REMATCH[2]}" - fi - done - - ${PYTHON} setup_fortran.py config_fc --f90flags "${FFLAGS}" build_ext - else ${PYTHON} setup.py build diff --git a/multigrid/MG.py b/multigrid/MG.py index 09ebd57fb..42b1df0e3 100644 --- a/multigrid/MG.py +++ b/multigrid/MG.py @@ -1,16 +1,18 @@ -""" +r""" The multigrid module provides a framework for solving elliptic problems. A multigrid object is just a list of grids, from the finest mesh down (by factors of two) to a single interior zone (each grid has the same number of guardcells). The main multigrid class is setup to solve a constant-coefficient -Helmholtz equation:: +Helmholtz equation + +.. math:: - (alpha - beta L) phi = f + (\alpha - \beta L) \phi = f -where L is the Laplacian and alpha and beta are constants. If alpha = -0 and beta = -1, then this is the Poisson equation. +where :math:`L` is the Laplacian and :math:`\alpha` and :math:`\beta` are constants. If :math:`\alpha = +0` and :math:`\beta = -1`, then this is the Poisson equation. We support Dirichlet or Neumann BCs, or a periodic domain. @@ -19,7 +21,7 @@ a = multigrid.CellCenterMG2d(nx, ny, verbose=1, alpha=alpha, beta=beta) this creates the multigrid object a, with a finest grid of nx by ny -zones and the default boundary condition types. alpha and beta are +zones and the default boundary condition types. :math:`\alpha` and :math:`\beta` are the coefficients of the Helmholtz equation. Setting verbose = 1 causing debugging information to be output, so you can see the residual errors in each of the V-cycles. @@ -502,9 +504,9 @@ def init_zeros(self): v[:, :] = 0.0 def init_RHS(self, data): - """ + r""" Initialize the right hand side, f, of the Helmholtz equation - (alpha - beta L) phi = f + :math:`(\alpha - \beta L) \phi = f` Parameters ---------- diff --git a/multigrid/__init__.py b/multigrid/__init__.py index 8e7feab78..0b5787b76 100644 --- a/multigrid/__init__.py +++ b/multigrid/__init__.py @@ -1,17 +1,23 @@ -"""This is the pyro multigrid solver. THere are several versions. +r"""This is the pyro multigrid solver. There are several versions. MG implements a second-order discretization of a constant-coefficient -Helmholtz equation:: +Helmholtz equation: - (alpha - beta L) phi = f +.. math:: -variable_coeff_MG implements a variable-coefficient Poisson equation:: + (\alpha - \beta L) \phi = f - div { eta grad phi } = f +variable_coeff_MG implements a variable-coefficient Poisson equation -general_MG implements a more general elliptic equation:: +.. math:: - alpha phi + div { beta grad phi } + gamma . grad phi = f + \nabla \cdot { \eta \nabla \phi } = f + +general_MG implements a more general elliptic equation + +.. math:: + + \alpha \phi + \nabla \cdot { \beta \nabla \phi } + \gamma \cdot \nabla \phi = f All use pure V-cycles to solve elliptic problems diff --git a/multigrid/general_MG.py b/multigrid/general_MG.py index 233aa5b8f..d415842ae 100644 --- a/multigrid/general_MG.py +++ b/multigrid/general_MG.py @@ -1,8 +1,10 @@ -""" +r""" This multigrid solver is build from multigrid/MG.py -and implements a more general solver for an equation of the form:: +and implements a more general solver for an equation of the form + +.. math:: - alpha phi + div { beta grad phi } + gamma . grad phi = f + \alpha \phi + \nabla\cdot { \beta \nabla \phi } + \gamma \cdot \nabla \phi = f where alpha, beta, and gamma are defined on the same grid as phi. These should all come in as cell-centered quantities. The solver @@ -24,18 +26,18 @@ class GeneralMG2d(MG.CellCenterMG2d): - """ + r""" this is a multigrid solver that supports our general elliptic equation. - we need to accept a coefficient CellCenterData2d object with - fields defined for alpha, beta, gamma_x, and gamma_y on the + we need to accept a coefficient ``CellCenterData2d`` object with + fields defined for ``alpha``, ``beta``, ``gamma_x``, and ``gamma_y`` on the fine level. We then restrict this data through the MG hierarchy (and average beta to the edges). - we need a new compute_residual() and smooth() function, that + we need a ``new compute_residual()`` and ``smooth()`` function, that understands these coeffs. """ diff --git a/multigrid/tests/mg_general_poisson_inhomogeneous.h5 b/multigrid/tests/mg_general_poisson_inhomogeneous.h5 index a794dd5bd..340cfaabe 100644 Binary files a/multigrid/tests/mg_general_poisson_inhomogeneous.h5 and b/multigrid/tests/mg_general_poisson_inhomogeneous.h5 differ diff --git a/multigrid/tests/mg_poisson_dirichlet.h5 b/multigrid/tests/mg_poisson_dirichlet.h5 index bdd99fad3..78faae305 100644 Binary files a/multigrid/tests/mg_poisson_dirichlet.h5 and b/multigrid/tests/mg_poisson_dirichlet.h5 differ diff --git a/multigrid/tests/mg_vc_poisson_dirichlet.h5 b/multigrid/tests/mg_vc_poisson_dirichlet.h5 index 0f226b771..0d6a63d51 100644 Binary files a/multigrid/tests/mg_vc_poisson_dirichlet.h5 and b/multigrid/tests/mg_vc_poisson_dirichlet.h5 differ diff --git a/multigrid/tests/mg_vc_poisson_periodic.h5 b/multigrid/tests/mg_vc_poisson_periodic.h5 index b87337bc6..3ea6a8b83 100644 Binary files a/multigrid/tests/mg_vc_poisson_periodic.h5 and b/multigrid/tests/mg_vc_poisson_periodic.h5 differ diff --git a/multigrid/variable_coeff_MG.py b/multigrid/variable_coeff_MG.py index 090054ac1..5e5c85983 100644 --- a/multigrid/variable_coeff_MG.py +++ b/multigrid/variable_coeff_MG.py @@ -1,10 +1,12 @@ -""" +r""" This multigrid solver is build from multigrid/MG.py and implements -a variable coefficient solver for an equation of the form:: +a variable coefficient solver for an equation of the form + +.. math:: - div { eta grad phi } = f + \nabla\cdot { \eta \nabla \phi } = f -where eta is defined on the same grid as phi. +where :math:`\eta` is defined on the same grid as :math:`\phi`. A cell-centered discretization is used throughout. """ @@ -21,14 +23,14 @@ class VarCoeffCCMG2d(MG.CellCenterMG2d): - """ + r""" this is a multigrid solver that supports variable coefficients we need to accept a coefficient array, coeffs, defined at each level. We can do this at the fine level and restrict it down the MG grids once. - we need a new compute_residual() and smooth() function, that + we need a new ``compute_residual()`` and ``smooth()`` function, that understands coeffs. """ diff --git a/paper/paper.md b/paper/paper.md index 4c3f0aea6..b54be8bb1 100644 --- a/paper/paper.md +++ b/paper/paper.md @@ -44,7 +44,12 @@ and a number of new contributors. pyro's functionality can now be accessed directly through a `Pyro()` class, in addition to the original commandline script interface. This new interface in particular allows for easy use within Jupyter notebooks. We also now use HDF5 -for output instead of python's `pickle()` function. +for output instead of python's `pickle()` function. Previously, we used Fortran +to speed up some performance-critical portions of the code. These routines +could be called by the main python code by first compiling them using `f2py`. +In the new version, we have replaced these Fortran routines by python functions +that are compiled at runtime by `numba`. Consequently, pyro is now written +entirely in python. The original goal of pyro was to learn hydrodynamics methods through example, and it still serves this goal. At Stony Brook, pyro is used diff --git a/particles/particles.py b/particles/particles.py index b399ce018..5683c1c88 100644 --- a/particles/particles.py +++ b/particles/particles.py @@ -210,13 +210,13 @@ def array_generate_particles(self, pos_array, init_array=None): self.particles[(ix, iy)] = Particle(x, y) def update_particles(self, dt, u=None, v=None): - """ + r""" Update the particles on the grid. This is based off the - AdvectWithUcc function in AMReX, which used the midpoint + ``AdvectWithUcc`` function in AMReX, which used the midpoint method to advance particles using the cell-centered velocity. We will explicitly pass in u and v if they cannot be accessed from the - sim_data using get_var("velocity"). + ``sim_data`` using ``get_var("velocity")``. Parameters ---------- diff --git a/setup.py b/setup.py index 41270a9bd..fdc6b0ae9 100644 --- a/setup.py +++ b/setup.py @@ -4,6 +4,6 @@ # # Note the setup.cfg directs the build to be done in-place. -from numpy.distutils.core import setup +from setuptools import setup setup(name='pyro') diff --git a/setup_fortran.py b/setup_fortran.py deleted file mode 100644 index 984c9bc36..000000000 --- a/setup_fortran.py +++ /dev/null @@ -1,15 +0,0 @@ -# this works for python 2 or 3 directly. To build, do: -# -# python setup.py build_ext -# -# Note the setup.cfg directs the build to be done in-place. - -from numpy.distutils.core import setup, Extension - -f_modules = ["compressible.interface", "advection_fv4.interface", - "lm_atm.LM_atm_interface", "incompressible.incomp_interface", "swe.interface"] - -ext_modules = [Extension(module, [module.replace('.', '/') + '_f.f90']) - for module in f_modules] - -setup(ext_modules=ext_modules) diff --git a/solver-test.ipynb b/solver-test.ipynb index 9fc36f014..c5d2c651e 100644 --- a/solver-test.ipynb +++ b/solver-test.ipynb @@ -88,14 +88,14 @@ "output_type": "stream", "text": [ "\u001b[33moutputting...\u001b[0m\n", - "main: 0.39704036712646484\n", - " vis: 0.1215507984161377\n" + "main: 0.3755648136138916\n", + " vis: 0.11605525016784668\n" ] }, { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 3, @@ -226,14 +226,14 @@ "\u001b[33mparameter particles.particle_generator never used\u001b[0m\n", "\u001b[33mparameter advection.u never used\u001b[0m\n", "\u001b[33mparameter advection.v never used\u001b[0m\n", - "main: 0.31522274017333984\n", - " vis: 0.10975933074951172\n" + "main: 0.3078031539916992\n", + " vis: 0.10725545883178711\n" ] }, { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 6, @@ -358,15 +358,15 @@ "\u001b[33moutputting...\u001b[0m\n", "\u001b[33mparameter particles.n_particles never used\u001b[0m\n", "\u001b[33mparameter particles.particle_generator never used\u001b[0m\n", - "main: 0.34589266777038574\n", - " evolve: 0.05175304412841797\n", - " vis: 0.12535333633422852\n" + "main: 1.5392191410064697\n", + " evolve: 1.2681171894073486\n", + " vis: 0.10479354858398438\n" ] }, { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 9, @@ -491,15 +491,15 @@ "\u001b[33moutputting...\u001b[0m\n", "\u001b[33mparameter particles.n_particles never used\u001b[0m\n", "\u001b[33mparameter particles.particle_generator never used\u001b[0m\n", - "main: 0.33099818229675293\n", - " evolve: 0.0028314590454101562\n", - " vis: 0.12505078315734863\n" + "main: 0.3059535026550293\n", + " evolve: 0.003072977066040039\n", + " vis: 0.11471986770629883\n" ] }, { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 12, @@ -568,7 +568,7 @@ "solver = \"compressible\"\n", "problem_name = \"rt\"\n", "param_file = \"inputs.rt\"\n", - "other_commands = [\"driver.max_steps=1\", \"mesh.nx=8\", \"mesh.ny=24\", \"driver.verbose=0\"]" + "other_commands = [\"driver.max_steps=1\", \"mesh.nx=8\", \"mesh.ny=24\", \"driver.verbose=0\", \"compressible.riemann=CGF\"]" ] }, { @@ -612,7 +612,7 @@ }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAADjCAYAAABEvCkpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xm4HHWd7/H355ycAEFCwqIGEgUVNxhhFMEZrysuwKOiXhyCW2DwcscRUZ/RER23Yca5bs84eHEZRmLADRFFc68ooI7bKEjgsgUEwyIciUYkLBGznfO5f1Qd7HT6dHeqqrt+3f19PU8/6a6qrv7151T6W8uvqmSbEEIIoYixuhsQQghhcEURCSGEUFgUkRBCCIVFEQkhhFBYFJEQQgiFRREJIYRQWBSREEIIhUURCT0hyZIe1zTsA5K+UFebQgjViyISQgihsKEvIpJ2k3SWpPWS1kl6W91tCqHXJN0u6V2SbsiX/c9J2rnudo0KSftI+pqk30m6TdKpdbepV4a+iADfAG4BHgksBT4m6ZH1NimEvngN8GLgscDjgffU25zRIGkM+D/ANcC+wBHAWyW9uNaG9chQFxFJLwGw/WHbm2x/H/g12X+oEIbdmbbvtH0P8EHg+LobNCKeDuxt+3Tbm23fCvwH2Urs0JlTdwN67GXAN2de5GsIuwO/ra1Fo2MKmGgaNgFsqaEto+rOhue/AvapqyEj5tHAPpLubRg2Dvy4pvb01LAXkcOBHza8fj5wt+2bamrPKLkD2A+4sWHY/sDNtbRmNC1peP4o4K66GjJi7gRus31A3Q3ph6HdnSVpAjgAOFbSzpIOBD4FvLPelo2MrwDvkbRY0pikFwAvBS6ouV2j5E15/nsA7yb7m4Te+zlwv6R3StpF0rikgyQ9ve6G9cLQFhHgScDtwPVku6++AXzQdvyI9cfpwE+BnwDrgY8Ar7F9fa2tGi1fAi4Bbs0f/1xvc0aD7SmyFaZDgNuAu4HPku1KHzoa1ptSSXot8Arb/73utoTQb5JuB95g+7t1tyUMt2HeEjmYbffHhxBCqFhfi0h+bOLnkq6RtFrSP/bw454C/KKH86+EpCWS/lPSjXkmb2kxjSR9QtIaSddKemrDuGWSfpk/lvW39SGEUdfX3VmSBOxqe0N+4PsnwFtsX9a3RiRG0iJgke2rJO0GXAm83PYNDdMcDbwZOJqsx9kZtg/PD5iuAg4FnL/3abbX9/t7hBBGU1+3RJzZkL+cyB/DeVCmS7bX2r4qf/4A2S64fZsmOwY4N8/vMmBBXnxeDFxq+568cFwKHNnH5ocQRlzfzxORNE62xvw44JO2L+93G6oi6WmdprF95Q7Mbz/gz4HmTPZl2xPHJvNhsw0fCVXnH7oX2dcntez7XkTy7m+HSFoAXCjpoOZun5JOBk4G2HWenvbEx83tdzNndeW1m+62vTfA7vPHVi16xHjb6SVtABpPbjzL9lktpnsY8DXgrbbvbx7dYtZuM7ywxuzHGX/aPOaXmV2lHmD9Q9lD1/lv856URfb1GpT8U8u+tjPWbd8r6Qdku1+ubxp3FnAWwKEH7+yfX7xk+xnUZHzRml/NPH/sfhNcfvHittNPLLrlJtuHtpsmPz70NeCLtr/eYpJJtj37eDHZ2ceTwHObhv+gbYM6aMx+vvbw4TqizOwq9V1f8KvG113m/6u2EyQksq/XoORfRfaSlgDnkl2Ydpps5faMhvFvBz5Kdg2wu9vNq69FRNLewJa8gOwCvAD4cD/bUCVjNrncpaDyzgZnAzfa/tdZJlsJnCLpPLID6/fZXivpYuBfJC3Mp3sR8K5SDRogVeQfions61NR9luBv2vs0CPpUts35AXmhWSXLuqo31sii4Bz8uMiY8D5tv9vn9tQGWO2eLrsbJ4JvA64TtLV+bB3k13rCNufAS4i65m1BngQODEfd4+kfwKuyN93en7F1pFQUf6hgMi+PlVkb3stsDZ//oCkmQ49NwAfB/6ehovXttPXImL7WrIDx0NhGtjoqVLzsP0TWh/baJzGwJtmGbccWF6qEQOqivxDMZF9farOvrFDj6SXAb+2fU22k6SzYb+Kb08Zs2W0eyjXapjz15xxxhfuWXcz/uR3274c5uwBkNBEOh162Pynp11mv5ekVQ2vO3boIdvF9Q9ku8W7FkWkBBu2DPH/o9RF/vWJ7OvTZfZ372iHHkl/Rna7hpmtkMXAVZIOs/2b2eYTRaQEIzZ6mC8/lrbIvz6RfX2qyL5Vhx7b1wEPb5jmduDQpHpnDRsDW+I/Um0i//pE9vWpKPuWHXpsX7SjM4oiUsI0YjPtT/oJvTPU+c+ZA3su7DxdvzQdExnq7AFJaG7z3Z1r1HBMpIrsu+zQs18384oiUoIRG6cjwrpE/vWJ7OuTWvbptGQAGdgyxGtjqYv86xPZ1ye17KOIlDBtsdEJbfKOmMi/PpF9fVLLPopICUZscURYl8i/PpF9fVLLPp2WDCAjNjudzcpRM8z5e3yMrXvuWnczZjXM2QMwJrTzTnW34k82/OlpatlHESlhGrFxOp3NylET+dcnsq9PatlHESkhO3M0IqxL5F+fyL4+qWWfTksGULZvMp3NylET+dcnsq9PatlHESkhu/xAOpuVo2aY8/ccsXlhQhcAbDLM2QP5BRjT/H6pZR9FpITU1ghGTeRfn8i+PqllH0WkhOwaNun8MUdN5F+fqrKX9Bbgf5BdguM/bP9b6ZkOudSW+7iCWgnTFpumJ9o+Qu9Ukb+kIyXdJGmNpNNajH+2pKskbZV0bMPwQyT9TNJqSddKOq7ir5e0irI/iKyAHAYcDLxE0gE9bvrAS+13J7ZESjBiy3Q6awSjpmz++W2aP0l2P+lJ4ApJK23f0DDZHcAJwNub3v4g8Hrbv5S0D9k9qi+2fW/hBg2Qipb9JwGX2X4QQNIPgVcAHyk742GW2u9OFJESjNhU8gCXpOXAS4B1tg9qMf4dwGvyl3PI/uPtnd9f/XbgAWAK2NrpJjTDpoL8DwPW2L4VQNJ5wDFk95nOPsO+PR+3zU2tbd/c8PwuSeuAvYFKisj0uNi4IJ0fimZdZt/p7nrXAx+UtCfwR+BooHH6+kiwU5odG6r43alSFJESbNgyXXqP4ArgTODc1p/hjwIfBZD0UuBttu9pmOR5nW4aM6wqyH9f4M6G15PA4Ts6E0mHAXOBW8o0ZpB0mX3bu+vZvlHSh4FLyc7JvobsFq2hjYp+dyoTRaSEKnpJ2P6RpP26nPx44MulPnCIdJl/u7XhVvdT2KGbvkpaBHweWGZ7utP0w6KqHkK2zya7wx6S/oWskIc2onfWEDFic5+u6y9pHnAkcMo2TYBLJBn496ZdBUOvy/zbrQ1PAksaXi8G7ur28yXNB74FvMf2Zd2+bxhUtexLerjtdZIeBbwS+IvSMx1yVWQvaQnZ3o9HAtNkK1dnSNoD+AqwH3A78Fe217ebVxSRErLNylJrwjvipcB/Ne3Kema+P/7hwKWSfmH7RwXmPZC6zL+dK4ADJO0P/BpYCry6mzdKmgtcCJxr+6tlGtGKx2HT7unssmhWQfYzvpYfE9kCvKnTD1bfSHgizZ/HirLfCvyd7ask7UbWMeRSsk4k37P9oby34mnAO9vNKM2UBoQRWzvf67jtfuEdsJSmXVm278r/XSfpQrIDxaNTRLrLf/b321slnQJcDIwDy22vlnQ6sMr2SklPJysWC4GXSvpH2wcCfwU8G9hT0gn5LE+wffX2nzR8ymb/0HzsZ1XQnJFSRfa21wJr8+cPSLqR7BjhMcBz88nOAX5AFJHescXmqd5HKGl34DnAaxuG7QqM5QvArsCLgNN73piEVJG/7YuAi5qGva/h+RVku7ma3/cF4AulPnyA9WvZD9vrMvuu94Dkx2T/HLgceEReYLC9Nt/L0VYsBSUYSq8RSPoyWeXfS9Ik8H5gAsD2Z/LJXgFcYvsPDW99BHChJMj+jl+y/Z1SjRkwVeQfions69Nl9l3tAZH0MOBrwFtt35//nuyQKCIlTCM2T5XunXV8F9OsIOsK3DjsVrKzfEdWFfmHYiL7+lSVvaQJsgLyRdtfzwf/VtKifCtkEbCu03z6WkRm6xHQzzZUyrE2Vqshzt/jsGV+3a1oY4izB7KTDeemc0LfNirIXtkmx9nAjbb/tWHUSmAZ8KH83292mle/t0Ra9ghouszEwDCwNaGTfkZN5F+fyL4+FWX/TOB1wHWSZjqDvJuseJwv6SSyS/68qtOM+lpE2vQIGNAiktY1bEZN5F+fyL4+VWRv+ye0PtkW4IgdmVdtx0SaegQMJBumYm2sNpF/fSL7+qSWfS1FpLlHQIvxJwMnAzxq35SP/SupP2YVGrPfmXk1t6aT4cq/Mfs5Cxayef4OXYGlz4Yre2ha9id2xxOpbmmllX3ff6Fn6RGwjbw/81kAhx68c7L/k1K7EFoVGrOfrz2SzR6GL//G7HdevCSy77PG/Heft0+y+aeWfb97Z83WI2AgObE1glET+dcnsq9Patn3e0ukZY+A/KzhgWPD1FQ6f8xRE/nXJ7KvT2rZ97t3VrseAQNpanqovs7AifzrE9nXJ6XsUz5qnTwjphParBw1w5y/x2Hr/HRvTzLM2QNYwnPT/HlMLfs0UxoUhumE1ghGTuRfn8i+PollH0WkpJT+mKMo8q9PZF+flLKPIlKCDU5os3LURP71iezrk1r2UURKEdNT6awRVE3j44zP3/2h19N/+GPH94ztukvHaQrPZ7t73g1x/uOG+VvqbkUb1WQv6W3AG8guCXUdcKLtjaVnXJZgeiKdH+ptpbXcp5rSYDB4Wm0foYci//pUkL2kfYFTgUNtH0R2d8mlPW754EtsuY8tkbIcP1S1ivzrU032c4BdJG0B5gF3VTHToZfQch9FpAwDCW1WjpzIvz7dZd/2Fq22fy3pY2SXHP8j2d07L6m8rcMmseU+dmeV5On2j9BbZfOXdKSkmyStkXRai/E7SfpKPv7y/OrTSJqQdI6k6yTdKOldVX+31HWR/d22D214bHOPb0kLgWOA/YF9gF0lvbbvX2QApfS7E1siZRhUco1A0nLgJcC6fL9w8/jnkt1d7LZ80Ndtn56POxI4g2xf8mdtf6hUY5occOADXHTxDx96feD//tsWX2Dbl6tP+VTH+R54Zov5NF3ubvWbt5/P+KLt31Mmf0njwCeBFwKTwBWSVjbdJO0kYL3tx0laCnwYOI7sZj072f4zSfOAGyR92fbthRvUYGxsmofN79wBoTYVLPvAC4DbbP8OQNLXgb8EvlB2xqWNiem5iV7Ft5rsZ/3tkfRm4BSymwh+y/bft5tPbImUIpju8OhsBXBkh2l+bPuQ/DFTQGZ+AI8CngwcL+nJJb7MACqd/2HAGtu32t4MnEe2ZtzoGOCc/PkFwBH5hURNtuY8B9gF2Axsd1uD4VXJsn8H8AxJ8/JMjwBu7Gmzh0Il2UOL3x5JzyNb5p9i+0DgY51mUriISDol3xwdbdMdHh3Y/hFwT4FP7uYHcPh1zn8vSasaHic3vHtf4M6G15P5MFpNY3srcB+wJ1lB+QPZnTrvAD5mu8jfcXCVX/YvJ8vxKrLuvWPkl2IPHZTMHmb97Xkj8CHbm/Jp1nWaT5ndWY8k2/y/ClgOXGw72Wvw90RFm5Vd+AtJ15D1XHm77dW0/gE8vB+NSUZ3+d9t+9BZxrV6c/MyPNs0hwFTZPvyFwI/lvRd27d2atBQqGjZt/1+4P3lGzRCusu+baeGNh4PPEvSB4GNZL83V7R7Q+EiYvs9kt4LvAg4EThT0vnA2bZvKTrfZjdfO48X73NIVbOrwJptXqlz2Sz6x5xxFfBo2xskHQ18AziA7n4Ae6/IJ1bYyi7yb2cSWNLwejHbdzGdmWYy33W1O9na26uB79jeAqyT9F/AoUAlRWR8bJoFu9R/zl07JbNP2sOWbOBZZ1xWdzMe8p9P2fZ1F9m3W3lqZw7ZStEzgKcD50t6TLsNhFLHRPIZ/yZ/bM0//AJJHykz34Fhutk32baHSsePsO+3vSF/fhEwIWkvuvsBHG7d5d/OFcABkvaXNJfsRLeVTdOsBJblz48Fvp8v93cAz1dmV7L/dL+o6qslr3z2oajeZj9J1nnHtn9OtnNsr3ZvKLwlIulUsv9cdwOfBd5he4ukMeCXQNsj+sNCUz2ev/RI4Le2LekwssL/e+Be8h9A4NdkP4Cv7m1r0lMmf9tbJZ0CXEzWw2257dWSTgdW2V5JdifOz0taQ7YFMnNG9SeBzwHXk20Vfs72tcVbM3h6veyH2fUw+28Azwd+IOnxwFyy3/hZlTkmshfwStu/ahxoe1rSS0rMd2DIoJJ9siV9GXgu2W6vSbL9wxMAtj9Dtvb7RklbyU7IWpqvCbf8ASzXmsFSRf751t1FTcPe1/B8I1l33ub3bWg1fFRUkX0opqrsZ/ntWQ4sl3Q9WY/DZZ2OdZc5JvK+NuNGppte2YOLto/vMP5M4MxZxm33Azhq+tSxIbQQ2denok4Ns/327NAJn3GyYRlDvjb2y9W7cfSBz3vo9ZINq7afaGzbhfnozzxv+2maLHnwyu0HTm+7stN6Ptt2ahjm/MdlHjZ3U93NmN0QZw9Z/ruNJ9qxIbHso4iUldAfcyRF/vWJ7OuTUPZRREpKaY1gFEX+9Yns65NS9lFEykhss3LkRP71iezrk1j2UURKEGn9MavmqSmm1m93O8G2pjZVsx9/av3mjtMMc/7jY9MsmJvuBRiHOXuAMaaZN5bmManUso8iUoajr3ytIv/6RPb1SSz7KCIlpbRGMIoi//pE9vVJKfsoImUktkYwciL/+kT29Uks+ygiJaW0RjCKIv/6RPb1SSn7vhaRTnfxGzgmqf7aI2eI8x9X2gfWhzl7gLHETzZMKft+39lwBZ3v4jcwBIxNt3+E3on86xPZ1ye17PtaRErcxS9dFdxhLJQQ+denZPaSniDp6obH/ZLe2sMWD4+Elvs4JlJGYif9jJzIvz7VXEH5JuAQAEnjZLc0uLB024ZdYst9kkUkvw/2yQA7M6/m1rQ3llAviSoMUvYwXPk3Zj9/0S7sMfGHmlvUXsXZHwHc0nxriX5qzH+vfSbYNdGTDSGt5b7fx0S6YvusmTsBTrBT3c2Z3cwBrkQ2K6swMNnD0OXfmP28hUOR/V6SVjU8Tm4zx6XAl3va5g4a8999jyTXrzOJLfcJJ5W+1C4/MGoi//p0mX1X9/nOb038MuBd5Vs2/Kpa7lv1lpX0UeClZDekugU40fa97ebT1y2R/E5aPwOeIGlS0kn9/PzKGcam3PYReijyr0+12R8FXGX7tz1q7XCpLvsVbN9b9lLgINtPAW6mi8Le795Zx9teZHvC9mLbZ/fz83tB0+0fobfK5i/pSEk3SVoj6bQW43eS9JV8/OWS9msa/yhJGyS9varvNCgqXPaPp+ZdWYOmiuxb9Za1fYntrfnLy4DFneYTu7PKqODyA51OwJT0GuCd+csNwBttX5OPux14AJgCtnaz62ColMw/7xH0SeCFwCRwhaSVtm9omOwkYL3tx0laCnwYOK5h/MeBbxdvRWtzNM0ecxI+sF7RpTckzSPL/3+Wn1t1xphmt7FET/bs32VP/hr4SqeJooiUkO2bLL3LZAXZPdTPnWX8bcBzbK+XdBRwFnB4w/jn2b67bCMGUQX5HwassX0rgKTzgGOAxiJyDPCB/PkFwJmSZNuSXg7cCiT8a98bFS372H4Q2LP0jEZIl9nvJanxftZn2T6r68+Q/gHYCnyx07RRRMqopq/8j5p3kTSN/2nDy642L0dG+fz3Be5seD3JtgV6m2lsb5V0H7CnpD+SbSG+EBi5XVmpnaswUrrLvqtODa1IWka2d+QI2x2rVRSRksa2dp6mQiex7a4TA5dIMvDvO7KmMSy6yL/dGplaTN/8n2a2af4R+LjtDVKrSYZfn5f90KBX2Us6kmzl6Dn5VmJHUUTKcO83K2dIeh5ZEflvDYOfafsuSQ8HLpX0i/xg2WjoLv92a2STwJKG14uBu2aZZlLSHGB3soORhwPHSvoIsACYlrTR9pk7+C1aGmeKPeZsqGJWvdFd9gNrTGZXdb67Zi0qyj7vLftcst+oSeD9ZL2xdiL7PQG4zPbftJtPFJESquwr3/ZzpKcAnwWOsv37meG278r/XSfpQrJ9/CNTRCroL38FcICk/ckuubEUeHXTNCuBZWRd048Fvp9v4j/roXZIHwA2VFVABkGco1OfqrK3fXyLwTvcYzaKSBk26vG5CJIeBXwdeJ3tmxuG7wqM2X4gf/4i4PSeNiY1JfPPj3GcAlwMjAPLba+WdDqwyvZKsv9Un5e0hmwLZGkFLR98fVj2wywSyz6KSBmm9B9zlk3KCQDbnwHeR9Z75VP55uVMV95HABfmw+YAX7L9nVKNGTQV5G/7IuCipmHva3i+EXhVh3l8oFQjBlEF2YeCEss+ikhJFfyItdqkbBz/BuANLYbfChxc6sOHQEr/mUZNZF+flLKPIlJGdHOs1xDnP6Ep9p5zf93NmN0QZw8wjtltLOUD63U34k+iiJQg0lojGDWRf30i+/qkln0UkTJsNJXQKsGoifzrE9nXJ7Hso4iUNMx95QdB5F+fyL4+KWUfRaSMxHpJjJwhzn9c0+w5lvAluYY4e4AxzLw+XeVwhyWWfRSRklLarBxFkX99Ivv6pJR9FJESlNhJP6Mm8q9PZF+f1LKPIlKGga3prBGMnMi/PpF9fRLLPopISZpO5485iiL/+kT29Ukp+ygiZdiQ0B9z5Axx/hNMs/d4onfWg8qyl7SA7OKiB5GtY/+17Z+VnnFJ4xK7jSV6if/ElvsoIiUpoc3KURT516ei7M8AvmP7WElzgXlVzHTYpbTcRxEpw4aEekmMnMi/PhVkL2k+8GzghGyW3gwkeq2RhCS23EcRKSuhzcqRFPnXp3P2nW7I9hjgd8DnJB0MXAm8xXbCJ8gkoppdiW8ju7irgeuAE/OrVu+QKCJl2LA17hFamyHOf1xij7Gxupsxu+6y73RDtjnAU4E3275c0hnAacB7K2plYWOIeZqouxmtVbDcS9oXOBV4su0/Sjqf7F45K3Z0XlFEyjBJbVaOnMi/PtVkPwlM2r48f30BWREJ7VS33M8BdpG0hexYVPOtobueSSgsrV4Soyfyr0/57G3/RtKdkp5g+ybgCOCGSpo31CrJ/teSPgbcAfwRuMT2JUXmFUWkDDO0u1MGQuRfn+qyfzPwxbxn1q3AiVXMdKh1l33b41GSFgLHAPsD9wJflfRa21/Y0eZEESnDxlOJXqRtFET+9akoe9tXA+2Om4Rm3WXf6XjUC4DbbP8OQNLXgb8E0i8iko4k6xs+DnzW9of63YbK2LCl9AGu5cBLgHW2D2oxXmR5HQ08CJxg+6p83DLgPfmk/2z7nFKNGTTV5N92eZS0E3Au8DTg98Bxtm/Px70LOAmYAk61fXGpxjSYwxgLxxM+ZaKC7FM2hpg3NrfuZrRWTfZ3AM+QNI9sd9YRwKr2b2mtr0VE0jjwSeCFZAfVrpC00vbA7getYG1sBXAm2Q9VK0cBB+SPw4FPA4dL2gN4P9lanIEr8yzXl23QICmTf5fL40nAetuPk7QU+DBwnKQnk/VmORDYB/iupMfbHplNo9gKrE/Z7PPecBcAVwFbgf8HnNX+Xa31uw/hYcAa27fmJxadR7ZfbjDNnPTT7tFxFv4RcE+bSY4BznXmMmCBpEXAi4FLbd+TF45LgSMr+FaDo3z+3SyPxwAzW3gXAEfkW4fHAOfZ3mT7NmBNPr/RUMGyHwqqKHvb77f9RNsH2X6d7U1FmtPv3Vn7Anc2vJ4kW7seSFPewn1bftvrj2mV2b5tho+MCvLvZnl8aBrbWyXdB+yZD7+s6b0jk3+flv3QQmrZ97uItLqi2XYXxpd0MnBy/nLTd33B9T1t1Y55wsyTP/DAvZdPXbq20/QdztrtZLbMuspyRw1K9lBJ/t1k2Lf8m7MfX7Rm0LN/dA/bU7lByT+17PtdRCaBJQ2vF9PiBJf8P/lZAJJWdehl0FeNP0i2F/bhI2fLbBJ4btPwH5T9sEHJHirJv5vlcWaaSUlzgN3Jdj92tSzviBHLPjmDkn9q2ff7mMgVwAGS9s/7hS8FVva5DYNmJfB6ZZ4B3Gd7LXAx8CJJC/M+3y/Kh4XudbM8rgSW5c+PBb5v2/nwpZJ2krQ/WceHn/ep3SEko69bIvk+5VPIfuzGgeW2V/ezDamR9GWyLYq9JE2S9biaALD9GeAisu69a8i6+J6Yj7tH0j+R/RACnG673QH60GS25VHS6cAq2yuBs4HPS1pDtgWyNH/v6vx6QzeQ9W550yj1zAphhrKVqnRJOnkHjyH0VGrt6aXUvmtq7eml1L5rau3ptdS+b2rtaZR8EQkhhJCuhK81HUIIIXXJFhFJR0q6SdIaSbVfHlrScknrJKXU7a9nUso/sq+1LZF9ve1JPv8ki0jD5SiOAp4MHJ9fZqJOKxiRM8ITzH8FkX1dVhDZ12kFieefZBEhwcujdHF5kmGSVP6RfWTfJ0llD4ORf6pFZOQv6VGzyL8+kX19IvsCUi0iPbmkR+ha5F+fyL4+kX0BqRaRyi8pEXZI5F+fyL4+kX0BqRaRuDxKvSL/+kT29YnsC0iyiNjeCsxcjuJG4Py6L4+SX57kZ2RXhZ2UdFKd7eml1PKP7CP7fkgtexiM/OOM9RBCCIUluSUSQghhMEQRCSGEUFgUkRBCCIVFEQkhhFBYFJEQQgiFRREJIYRQWBSREEIIhQ11EZH0dEnXStpZ0q6SVks6qO52jYLIvl6Rf31GLfuhP9lQ0j8DOwO7AJO2/1fNTRoZkX29Iv/6jFL2o1BE5pJdE2cj8Je2p2pu0siI7OsV+ddnlLIf6t1ZuT2AhwG7ka0ZhP6J7OsV+ddnZLIfhS2RlWR3KNsfWGT7lJqbNDIi+3pF/vUZpezn1N2AXpL0emCr7S/l90/+qaTn2/5+3W0bdpF9vSL/+oxa9kO/JRJCCKF3RuGYSAghhB6JIhJCCKGwKCIhhBAKiyISQgihsCgiIYQQCosiEkIIobAoIiGEEAqLIhJCCKGsIQjHAAAEJklEQVSwKCIhhBAKiyISQgihsCgiIYQQCosiEkIIobAoIiGEEAqLIhJCCKGwKCIhhBAKiyISQgihsCgiIYQQCosiEkIIobAoIiGEEAqLIhJCCKGwKCIhhBAKiyISQgihsCgiIYQQCosiEkIIobAoIiGEEAqLIhJCCKGwKCIhhBAKiyISQgihsCgiIYQQCosiEkIIobAoIiGEEAqLIhJCCKGwKCIhhBAKiyISQgihsCgiIYQQCosiEkIIobAoIiGEEAqLIhJCCKGwKCIhhBAKiyISQgihsCgiIYQQCosiEkIIobAoIiGEEAqLIhJCCKGwKCIhhBAKiyISQgihsCgiIYQQCosiEkIIobAoIiGEEAprW0QkLZD0t1V/qDKfkLRG0rWSnjrLdE+TdF0+3SckKR++h6RLJf0y/3dhu/lKOkTSzyStzocf19SWD0q6WdKNkk7Nhz8xf88mSW9vkcsFkn6Rv+cvqs4ohBAGQactkQVA5UUEOAo4IH+cDHx6luk+nY+fmfbIfPhpwPdsHwB8L3/dbr4PAq+3fWA+j3+TtCAfdwKwBHii7ScB5+XD7wFOBT7Wol1nAN+x/UTgYODGbr94CCEMk05F5EPAYyVdLemjFX7uMcC5zlwGLJC0qHGC/PV82z+zbeBc4OUN7z8nf35O0/Dt5mv7Ztu/BLB9F7AO2Dt/zxuB021P5+PXzfxr+wpgS1O75gPPBs7Op9ts+94KMgkhhIEzp8P404CDbB/SaqSkHwO7tRj1dtvfbTPffYE7G15P5sPWNk0z2WIagEfYXgtge62kh3c7X0mHAXOBW/JBjwWOk/QK4HfAqTMFZxaPyaf7nKSDgSuBt9j+Q5v3hBDCUOpURNqy/ayCb1Wr2RWYZofmm2/dfB5YNrPlAewEbLR9qKRXAsuBdt9rDvBU4M22L5d0BlmxfW+HtoUQwtAp1TtL0o/zXV3Njxd0eOsk2XGIGYuBu1pMs3iWaX47s/sr/3ddp/nmu6G+Bbwn39XV+Dlfy59fCDyli7ZP2r48f30BWVEJIYSR06mIPEDr3VVAtiVi+5AWj3a7sgBWAq/Pe0Y9A7hvZvdUw7zXAg9IekbeK+v1wDcb3r8sf76safh285U0l6xAnGv7q01t+Qbw/Pz5c4Cb2zXc9m+AOyU9IR90BHBDh+8bQghDSdkx6zYTSF8iWzv/tu13VPKhWVE4k6yn1IPAibZX5eOunjkGI+lQYAWwC/Btsl1IlrQncD7wKOAO4FW275ltvpJeC3wOWN3QjBNsX5330vpiPq8NwN/YvkbSI4FVwHxgOh/3ZNv3SzoE+CzZsZVb889ZX0U2IYQwSDoWkRBCCGE2ccZ6CCGEwqKIhBBCKCyKSAghhMKiiIQQQigsikgIIYTCooiEEEIoLIpICCGEwqKIhBBCKOz/A5ARba6YGR6kAAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAADjCAYAAABEvCkpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XuYJHV97/H3Z2ZngUWWXS7qwq6CihohQhTBxKOCKAKPinpIWLwtBA8nRkR9QiImRg2JiReeGHzwEiLrgjdEFLPnhAgYY9REkIXDbUFwuQgjKythuay4t5nP+aNqsLe3p7u3qrrr193f1/P0M91V1dW/+UxNf+vyqyrZJoQQQihirO4GhBBCGFxRREIIIRQWRSSEEEJhUURCCCEUFkUkhBBCYVFEQgghFBZFJIQQQmFRREJPSLKkZzUN+7CkL9XVphBC9aKIhBBCKGzoi4ik3SSdL2m9pHWS3lt3m0LoNUn3SHq/pFvzZf8Lknauu12jQtI+kr4h6ZeS7pZ0Rt1t6pWhLyLAt4A7gacCS4FzJD213iaF0BdvBl4NPBN4NvCBepszGiSNAf8HuBHYFzgKeI+kV9fasB4Z6iIi6TUAtj9me5Pt7wI/J/uHCmHYnWf7PtsPAR8BTqq7QSPiRcDets+2vdn2XcA/ka3EDp05dTegx14H/PPMi3wNYXfggdpaNDqmgImmYRPAlhraMqrua3j+M2CfuhoyYp4O7CPp4YZh48APampPTw17ETkc+I+G168AHrR9e03tGSX3AvsBtzUM2x+4o5bWjKYlDc+fBtxfV0NGzH3A3bYPqLsh/TC0u7MkTQAHACdI2lnSgcBngPfV27KR8TXgA5IWSxqT9ErgtcClNbdrlLwzz38P4M/J/iah934MPCrpfZJ2kTQu6SBJL6q7Yb0wtEUE+C3gHuAWst1X3wI+Yju+xPrjbOC/gB8C64GPA2+2fUutrRotXwGuBO7KH39Tb3NGg+0pshWmQ4C7gQeBz5PtSh86GtabUkl6C/AG2/+z7raE0G+S7gHebvs7dbclDLdh3hI5mG33x4cQQqhYX4tIfmzix5JulLRa0l/18OOeD/ykh/OvhKQlkv5d0m15Ju9uMY0kfUrSGkk3SXpBw7hlkn6aP5b1t/UhhFHX191ZkgTsantDfuD7h8C7bV/dt0YkRtIiYJHt6yXtBlwHvN72rQ3THAe8CziOrMfZubYPzw+YrgIOBZy/94W21/f79wghjKa+bok4syF/OZE/hvOgTJdsr7V9ff78MbJdcPs2TXY8cFGe39XAgrz4vBq4yvZDeeG4Cjimj80PIYy4vp8nImmcbI35WcCnbV/T7zZURdILO01j+7odmN9+wO8AzZnsy7Ynjk3mw2YbPhKqzj90L7KvT2rZ972I5N3fDpG0ALhM0kHN3T4lnQacBrDrPL3wuc+a2+9mzuq6mzY9aHtvgN3nj61a9JTxttNL2gA0ntx4vu3zW0z3JOAbwHtsP9o8usWs3WZ4YY3ZjzP+wnnMLzO7Sj3G+ieyh67z3+Y9KYvs6zUo+aeWfW1nrNt+WNL3yHa/3NI07nzgfIBDD97ZP75iyfYzqMn4ojU/m3n+zP0muOaKxW2nn1h05+22D203TX586BvAl21/s8Ukk2x79vFisrOPJ4EjmoZ/r22DOmjMfr728OE6qszsKvUdX/qzxtdd5v+zthMkJLKv16DkX0X2kpYAF5FdmHaabOX23IbxZwKfILsG2IPt5tXXIiJpb2BLXkB2AV4JfKyfbaiSMZtc7lJQeWeDC4DbbP/9LJOtBE6XdDHZgfVHbK+VdAXwt5IW5tMdDby/VIMGSBX5h2Ii+/pUlP1W4E8aO/RIusr2rXmBeRXZpYs66veWyCLgwvy4yBhwie3/2+c2VMaYLZ4uO5uXAG8FbpZ0Qz7sz8mudYTtzwGXk/XMWgM8DpySj3tI0l8D1+bvOzu/YutIqCj/UEBkX58qsre9FlibP39M0kyHnluBTwJ/RsPFa9vpaxGxfRPZgeOhMA1s9FSpedj+Ia2PbTROY+Cds4xbDiwv1YgBVUX+oZjIvj5VZ9/YoUfS64Cf274x20nS2bBfxbenjNky2j2UazXM+WvOOOML96y7Gb/xy21fDnP2AEhoIp0OPWz+zdMus99L0qqG1x079JDt4voLst3iXYsiUoINW4b4/yh1kX99Ivv6dJn9gzvaoUfSb5PdrmFmK2QxcL2kw2z/Yrb5RBEpwYiNHubLj6Ut8q9PZF+fKrJv1aHH9s3AkxumuQc4NKneWcPGwJb4R6pN5F+fyL4+FWXfskOP7ct3dEZRREqYRmym/Uk/oXeGOv85c2DPhZ2n65emYyJDnT0gCc1tvrtzjRqOiVSRfZcdevbrZl5RREowYuN0RFiXyL8+kX19Uss+nZYMIANbhnhtLHWRf30i+/qkln0UkRKmLTY6oU3eERP51yeyr09q2UcRKcGILY4I6xL51yeyr09q2afTkgFkxGans1k5aoY5f4+PsXXPXetuxqyGOXsAxoR23qnuVvzGht88TS37KCIlTCM2TqezWTlqIv/6RPb1SS37KCIlZGeORoR1ifzrE9nXJ7Xs02nJAMr2TaazWTlqIv/6RPb1SS37KCIlZJcfSGezctQMc/6eIzYvTOgCgE2GOXsgvwBjmr9fatlHESkhtTWCURP51yeyr09q2UcRKSG7hk06f8xRE/nXp6rsJb0b+F9kl+D4J9v/UHqmQy615T6uoFbCtMWm6Ym2j9A7VeQv6RhJt0taI+msFuNfJul6SVslndAw/BBJP5K0WtJNkk6s+NdLWkXZH0RWQA4DDgZeI+mAHjd94KX2vRNbIiUYsWU6nTWCUVM2//w2zZ8mu5/0JHCtpJW2b22Y7F7gZODMprc/DrzN9k8l7UN2j+orbD9cuEEDpKJl/7eAq20/DiDpP4A3AB8vO+Nhltr3ThSREozYVPIAl6TlwGuAdbYPajH+T4E35y/nkP3j7Z3fX/0e4DFgCtja6SY0w6aC/A8D1ti+C0DSxcDxZPeZzj7Dvicft81NrW3f0fD8fknrgL2BSorI9LjYuCCdL4pmXWbf6e56twAfkbQn8GvgOKBx+vpIsFOaHRuq+N6pUhSREmzYMl16j+AK4Dzgotaf4U8AnwCQ9FrgvbYfapjkyE43jRlWFeS/L3Bfw+tJ4PAdnYmkw4C5wJ1lGjNIusy+7d31bN8m6WPAVWTnZN9IdovW0EZF3zuViSJSQhW9JGx/X9J+XU5+EvDVUh84RLrMv93acKv7KezQTV8lLQK+CCyzPd1p+mFRVQ8h2xeQ3WEPSX9LVshDG9E7a4gYsblP1/WXNA84Bjh9mybAlZIM/GPTroKh12X+7daGJ4ElDa8XA/d3+/mS5gP/AnzA9tXdvm8YVLXsS3qy7XWSnga8Efjd0jMdclVkL2kJ2d6PpwLTZCtX50raA/gasB9wD/AHtte3m1cUkRKyzcpSa8I74rXAfzbtynpJvj/+ycBVkn5i+/sF5j2Qusy/nWuBAyTtD/wcWAq8qZs3SpoLXAZcZPvrZRrRisdh0+7p7LJoVkH2M76RHxPZAryz0xdW30h4Is2vx4qy3wr8ie3rJe1G1jHkKrJOJP9m+6N5b8WzgPe1m1GaKQ0II7Z2vtdx2/3CO2ApTbuybN+f/1wn6TKyA8WjU0S6y3/299tbJZ0OXAGMA8ttr5Z0NrDK9kpJLyIrFguB10r6K9sHAn8AvAzYU9LJ+SxPtn3D9p80fMpm/8R87JdW0JyRUkX2ttcCa/Pnj0m6jewY4fHAEflkFwLfI4pI79hi81TvI5S0O/By4C0Nw3YFxvIFYFfgaODsnjcmIVXkb/ty4PKmYR9seH4t2W6u5vd9CfhSqQ8fYP1a9sP2usy+6z0g+THZ3wGuAZ6SFxhsr833crQVS0EJhtJrBJK+Slb595I0CXwImACw/bl8sjcAV9r+VcNbnwJcJgmyv+NXbH+7VGMGTBX5h2Ii+/p0mX1Xe0AkPQn4BvAe24/m3yc7JIpICdOIzVOle2ed1MU0K8i6AjcOu4vsLN+RVUX+oZjIvj5VZS9pgqyAfNn2N/PBD0halG+FLALWdZpPX4vIbD0C+tmGSjnWxmo1xPl7HLbMr7sVbQxx9kB2suHcdE7o20YF2Svb5LgAuM323zeMWgksAz6a//znTvPq95ZIyx4BTZeZGBgGtiZ00s+oifzrE9nXp6LsXwK8FbhZ0kxnkD8nKx6XSDqV7JI/v99pRn0tIm16BAxoEUnrGjajJvKvT2Rfnyqyt/1DWp9sC3DUjsyrtmMiTT0CBpINU7E2VpvIvz6RfX1Sy76WItLcI6DF+NOA0wCetm/Kx/6V1B+zCo3Z78y8mlvTyXDl35j9nAUL2Tx/h67A0mfDlT00LfsTu+OJVLe00sq+79/Qs/QI2Eben/l8gEMP3jnZ/6TULoRWhcbs52uPZLOH4cu/MfudFy+J7PusMf/d5+2TbP6pZd/v3lmz9QgYSE5sjWDURP71iezrk1r2/d4SadkjID9reODYMDWVzh9z1ET+9Yns65Na9v3undWuR8BAmpoeql9n4ET+9Yns65NS9ikftU6eEdMJbVaOmmHO3+OwdX66tycZ5uwBLOG5aX49ppZ9mikNCsN0QmsEIyfyr09kX5/Eso8iUlJKf8xRFPnXJ7KvT0rZRxEpwQYntFk5aiL/+kT29Ukt+ygipYjpqXTWCKqm8XHG5+/+xOvpX/2643vGdt2l4zSF57PdPe+GOP9xw/wtdbeijWqyl/Re4O1kl4S6GTjF9sbSMy5LMD2Rzhf1ttJa7lNNaTAYPK22j9BDkX99Kshe0r7AGcChtg8iu7vk0h63fPAlttzHlkhZji+qWkX+9akm+znALpK2APOA+6uY6dBLaLmPIlKGgYQ2K0dO5F+f7rJve4tW2z+XdA7ZJcd/TXb3zisrb+uwSWy5j91ZJXm6/SP0Vtn8JR0j6XZJaySd1WL8TpK+lo+/Jr/6NJImJF0o6WZJt0l6f9W/W+q6yP5B24c2PLa5x7ekhcDxwP7APsCukt7S919kAKX0vRNbImUYVHKNQNJy4DXAuny/cPP4I8juLnZ3Puibts/Oxx0DnEu2L/nztj9aqjFNDjjwMS6/4j+eeH3geX/c8T2rT/9Mx2mKzmd8UdOAkvlLGgc+DbwKmASulbSy6SZppwLrbT9L0lLgY8CJZDfr2cn2b0uaB9wq6au27yncoAZjY9M8aX7nDgi1qWDZB14J3G37lwCSvgn8HvClsjMubUxMz030Kr7VZD/rd4+kdwGnk91E8F9s/1m7+cSWSCmC6Q6PzlYAx3SY5ge2D8kfMwVk5gvwWOB5wEmSnlfilxlApfM/DFhj+y7bm4GLydaMGx0PXJg/vxQ4Kr+QqMnWnOcAuwCbge1uazC8Kln27wVeLGlenulRwG09bfZQqCR7aPHdI+lIsmX++bYPBM7pNJPCRUTS6fnm6Gib7vDowPb3gYcKfHI3X4DDr3P+e0la1fA4reHd+wL3NbyezIfRahrbW4FHgD3JCsqvyO7UeS9wju0if8fBVX7Zv4Ysx+vJuveOkV+KPXRQMnuY9bvnHcBHbW/Kp1nXaT5ldmc9lWzz/3pgOXCF7WSvwd8TFW1WduF3Jd1I1nPlTNuraf0FeHg/GpOM7vJ/0Pahs4xr9ebmZXi2aQ4Dpsj25S8EfiDpO7bv6tSgoVDRsm/7Q8CHyjdohHSXfdtODW08G3ippI8AG8m+b65t94bCRcT2ByT9JXA0cApwnqRLgAts31l0vs3uuGker97nkKpmV4E127xS57JZ9I8543rg6bY3SDoO+BZwAN19AVarqrlX2Mou8m9nEljS8Hox23cxnZlmMt91tTvZ2tubgG/b3gKsk/SfwKFAJUVkfGyaBbvUf85dOyWzT9qTlmzgpedeXXcznvDvz9/2dRfZt1t5amcO2UrRi4EXAZdIeka7DYRSx0TyGf8if2zNP/xSSR8vM9+BYbrZN9m2h0rHj7Aftb0hf345MCFpL7r7Ahxu3eXfzrXAAZL2lzSX7ES3lU3TrASW5c9PAL6bL/f3Aq9QZleyf7qfVPWrJa989qGo3mY/SdZ5x7Z/TLZzbK92byi8JSLpDLJ/rgeBzwN/anuLpDHgp0DbI/rDQlM9nr/0VOAB25Z0GFnh/2/gYfIvQODnZF+Ab+pta9JTJn/bWyWdDlxB1sNtue3Vks4GVtleSXYnzi9KWkO2BTJzRvWngS8At5BtFX7B9k3FWzN4er3sh9n1MPtvAa8Avifp2cBcsu/4WZU5JrIX8EbbP2scaHta0mtKzHdgyKCSfbIlfRU4gmy31yTZ/uEJANufI1v7fYekrWQnZC3N14RbfgGWa81gqSL/fOvu8qZhH2x4vpGsO2/z+za0Gj4qqsg+FFNV9rN89ywHlku6hazH4bJOx7rLHBP5YJtxI9NNr+zBRdsndRh/HnDeLOO2+wIcNX3q2BBaiOzrU1Gnhtm+e3bohM842bCMIV8b++nq3TjuwCOfeL1kw6rtJxrbdmE+7nNHbj9NkyWPX7f9wOltV3Zaz2fbTg3DnP+4zJPmbqq7GbMb4uwhy3+38UQ7NiSWfRSRshL6Y46kyL8+kX19Eso+ikhJKa0RjKLIvz6RfX1Syj6KSBmJbVaOnMi/PpF9fRLLPopICSKtP2bVPDXF1PrtbifY1tSmavbjT63f3HGaYc5/fGyaBXPTvQDjMGcPMMY088bSPCaVWvZRRMpw9JWvVeRfn8i+PollH0WkpJTWCEZR5F+fyL4+KWUfRaSMxNYIRk7kX5/Ivj6JZR9FpKSU1ghGUeRfn8i+Pill39ci0ukufgPHJNVfe+QMcf7jSvvA+jBnDzCW+MmGKWXf7zsbrqDzXfwGhoCx6faP0DuRf30i+/qkln1fi0iJu/ilq4I7jIUSIv/6lMxe0nMk3dDweFTSe3rY4uGR0HIfx0TKSOykn5ET+denmiso3w4cAiBpnOyWBpeVbtuwS2y5T7KI5PfBPg1gZ+bV3Jr2xhLqJVGFQcoehiv/xuznL9qFPSZ+VXOL2qs4+6OAO5tvLdFPjfnvtc8EuyZ6siGktdz3+5hIV2yfP3MnwAl2qrs5s5s5wJXIZmUVBiZ7GLr8G7Oft3Aost9L0qqGx2lt5rgU+GpP29xBY/6775Hk+nUmseU+4aTSl9rlB0ZN5F+fLrPv6j7f+a2JXwe8v3zLhl9Vy32r3rKSPgG8luyGVHcCp9h+uN18+rolkt9J60fAcyRNSjq1n59fOcPYlNs+Qg9F/vWpNvtjgettP9Cj1g6X6rJfwfa9Za8CDrL9fOAOuijs/e6ddZLtRbYnbC+2fUE/P78XNN3+EXqrbP6SjpF0u6Q1ks5qMX4nSV/Lx18jab+m8U+TtEHSmVX9ToOiwmX/JGrelTVoqsi+VW9Z21fa3pq/vBpY3Gk+sTurjAouP9DpBExJbwbel7/cALzD9o35uHuAx4ApYGs3uw6GSsn88x5BnwZeBUwC10paafvWhslOBdbbfpakpcDHgBMbxn8S+NfirWhtjqbZY07CB9YruvSGpHlk+f/v8nOrzhjT7DaW6Mme/bvsyR8CX+s0URSRErJ9k6V3mawgu4f6RbOMvxt4ue31ko4FzgcObxh/pO0HyzZiEFWQ/2HAGtt3AUi6GDgeaCwixwMfzp9fCpwnSbYt6fXAXUDC3/a9UdGyj+3HgT1Lz2iEdJn9XpIa72d9vu3zu/4M6S+ArcCXO00bRaSMavrKf795F0nT+P9qeNnV5uXIKJ//vsB9Da8n2bZAbzON7a2SHgH2lPRrsi3EVwEjtysrtXMVRkp32XfVqaEVScvI9o4cZbtjtYoiUtLY1s7TVOhUtt11YuBKSQb+cUfWNIZFF/m3WyNTi+mb/2lmm+avgE/a3iC1mmT49XnZDw16lb2kY8hWjl6ebyV2FEWkDPd+s3KGpCPJisj/aBj8Etv3S3oycJWkn+QHy0ZDd/m3WyObBJY0vF4M3D/LNJOS5gC7kx2MPBw4QdLHgQXAtKSNts/bwd+ipXGm2GPOhipm1RvdZT+wxmR2Vee7a9aiouzz3rJHkH1HTQIfIuuNtRPZ9wnA1bb/qN18ooiUUGVf+bafIz0f+DxwrO3/nhlu+/785zpJl5Ht4x+ZIlJBf/lrgQMk7U92yY2lwJuaplkJLCPrmn4C8N18E/+lT7RD+jCwoaoCMgjiHJ36VJW97ZNaDN7hHrNRRMqwUY/PRZD0NOCbwFtt39EwfFdgzPZj+fOjgbN72pjUlMw/P8ZxOnAFMA4st71a0tnAKtsryf6pvihpDdkWyNIKWj74+rDsh1kkln0UkTJM6T/mLJuUEwC2Pwd8kKz3ymfyzcuZrrxPAS7Lh80BvmL726UaM2gqyN/25cDlTcM+2PB8I/D7Hebx4VKNGEQVZB8KSiz7KCIlVfAl1mqTsnH824G3txh+F3BwqQ8fAin9M42ayL4+KWUfRaSM6OZYryHOf0JT7D3n0bqbMbshzh5gHLPbWMoH1utuxG9EESlBpLVGMGoi//pE9vVJLfsoImXYaCqhVYJRE/nXJ7KvT2LZRxEpaZj7yg+CyL8+kX19Uso+ikgZifWSGDlDnP+4ptlzLOFLcg1x9gBjmHl9usrhDkss+ygiJaW0WTmKIv/6RPb1SSn7KCIlKLGTfkZN5F+fyL4+qWUfRaQMA1vTWSMYOZF/fSL7+iSWfRSRkjSdzh9zFEX+9Yns65NS9lFEyrAhoT/myBni/CeYZu/xRO+sB5VlL2kB2cVFDyJbx/5D2z8qPeOSxiV2G0v0Ev+JLfdRREpSQpuVoyjyr09F2Z8LfNv2CZLmAvOqmOmwS2m5jyJShg0J9ZIYOZF/fSrIXtJ84GXAydksvRlI9FojCUlsuY8iUlZCm5UjKfKvT+fsO92Q7RnAL4EvSDoYuA54t+2ET5BJRDW7Et9LdnFXAzcDp+RXrd4hUUTKsGFr3CO0NkOc/7jEHmNjdTdjdt1l3+mGbHOAFwDvsn2NpHOBs4C/rKiVhY0h5mmi7ma0VsFyL2lf4AzgebZ/LekSsnvlrNjReUURKcMktVk5ciL/+lST/SQwafua/PWlZEUktFPdcj8H2EXSFrJjUc23hu56JqGwtHpJjJ7Ivz7ls7f9C0n3SXqO7duBo4BbK2neUKsk+59LOge4F/g1cKXtK4vMK4pIGWZod6cMhMi/PtVl/y7gy3nPrLuAU6qY6VDrLvu2x6MkLQSOB/YHHga+Lukttr+0o82JIlKGjacSvUjbKIj861NR9rZvANodNwnNusu+0/GoVwJ32/4lgKRvAr8HpF9EJB1D1jd8HPi87Y/2uw2VsWFL6QNcy4HXAOtsH9RivMjyOg54HDjZ9vX5uGXAB/JJ/8b2haUaM2iqyb/t8ihpJ+Ai4IXAfwMn2r4nH/d+4FRgCjjD9hWlGtNgDmMsHE/4lIkKsk/ZGGLe2Ny6m9FaNdnfC7xY0jyy3VlHAavav6W1vhYRSePAp4FXkR1Uu1bSStsDux+0grWxFcB5ZF9UrRwLHJA/Dgc+CxwuaQ/gQ2RrcQauy7NcX7ZBg6RM/l0uj6cC620/S9JS4GPAiZKeR9ab5UBgH+A7kp5te2Q2jWIrsD5ls897w10KXA9sBf4fcH77d7XW7z6EhwFrbN+Vn1h0Mdl+ucE0c9JPu0fHWfj7wENtJjkeuMiZq4EFkhYBrwausv1QXjiuAo6p4LcaHOXz72Z5PB6Y2cK7FDgq3zo8HrjY9ibbdwNr8vmNhgqW/VBQRdnb/pDt59o+yPZbbW8q0px+787aF7iv4fUk2dr1QJryFh7Z8kCvP6ZVZvu2GT4yKsi/m+XxiWlsb5X0CLBnPvzqpveOTP59WvZDC6ll3+8i0uqKZttdGF/SacBp+ctN3/Glt/S0VTvmOTNPfsVjD18zddXaTtN3OGu3k9ky6yrLHTUo2UMl+XeTYd/yb85+fNGaQc/+6T1sT+UGJf/Usu93EZkEljS8XkyLE1zyf/LzASSt6tDLoK8av5BsL+zDR86W2SRwRNPw75X9sEHJHirJv5vlcWaaSUlzgN3Jdj92tSzviBHLPjmDkn9q2ff7mMi1wAGS9s/7hS8FVva5DYNmJfA2ZV4MPGJ7LXAFcLSkhXmf76PzYaF73SyPK4Fl+fMTgO/adj58qaSdJO1P1vHhx31qdwjJ6OuWSL5P+XSyL7txYLnt1f1sQ2okfZVsi2IvSZNkPa4mAGx/DricrHvvGrIuvqfk4x6S9NdkX4QAZ9tud4A+NJlteZR0NrDK9krgAuCLktaQbYEszd+7Or/e0K1kvVveOUo9s0KYoWylKl2STtvBYwg9lVp7eim13zW19vRSar9rau3ptdR+39Ta0yj5IhJCCCFdCV9rOoQQQuqSLSKSjpF0u6Q1kmq/PLSk5ZLWSUqp21/PpJR/ZF9rWyL7etuTfP5JFpGGy1EcCzwPOCm/zESdVjAiZ4QnmP8KIvu6rCCyr9MKEs8/ySJCgpdH6eLyJMMkqfwj+8i+T5LKHgYj/1SLyMhf0qNmkX99Ivv6RPYFpFpEenJJj9C1yL8+kX19IvsCUi0ilV9SIuyQyL8+kX19IvsCUi0icXmUekX+9Yns6xPZF5BkEbG9FZi5HMVtwCV1Xx4lvzzJj8iuCjsp6dQ629NLqeUf2Uf2/ZBa9jAY+ccZ6yGEEApLckskhBDCYIgiEkIIobAoIiGEEAqLIhJCCKGwKCIhhBAKiyISQgihsCgiIYQQChvqIiLpRZJukrSzpF0lrZZ0UN3tGgWRfb0i//qMWvZDf7KhpL8BdgZ2ASZt/13NTRoZkX29Iv/6jFL2o1BE5pJdE2cj8Hu2p2pu0siI7OsV+ddnlLIf6t1ZuT2AJwG7ka0ZhP6J7OsV+ddnZLIfhS2RlWR3KNsfWGT79JqbNDIi+3pF/vUZpezn1N2AXpL0NmCr7a/k90/+L0mvsP3duts27CL7ekX+9Rm17Id+SySEEELvjMIxkRBCCD0SRSSEEEJhUURCCCEUFkUkhBBCYVFEQgghFBZFJIQQQmFRREIIIRQWRSSEEEJhUURCCCEUFkUkhBA5xa2dAAAEG0lEQVRCYVFEQgghFBZFJIQQQmFRREIIIRQWRSSEEEJhUURCCCEUFkUkhBBCYVFEQgghFBZFJIQQQmFRREIIIRQWRSSEEEJhUURCCCEUFkUkhBBCYVFEQgghFBZFJIQQQmFRREIIIRQWRSSEEEJhUURCCCEUFkUkhBBCYVFEQgghFBZFJIQQQmFRREIIIRQWRSSEEEJhUURCCCEUFkUkhBBCYVFEQgghFBZFJIQQQmFRREIIIRQWRSSEEEJhUURCCCEUFkUkhBBCYVFEQgghFBZFJIQQQmFRREIIIRQWRSSEEEJhUURCCCEUFkUkhBBCYVFEQgghFBZFJIQQQmFRREIIIRTWtohIWiDpj6v+UGU+JWmNpJskvWCW6V4o6eZ8uk9JUj58D0lXSfpp/nNhu/lKOkTSjyStzoef2NSWj0i6Q9Jtks7Ihz83f88mSWe2yOVSST/J3/O7VWcUQgiDoNOWyAKg8iICHAsckD9OAz47y3SfzcfPTHtMPvws4N9sHwD8W/663XwfB95m+8B8Hv8gaUE+7mRgCfBc278FXJwPfwg4AzinRbvOBb5t+7nAwcBt3f7iIYQwTDoVkY8Cz5R0g6RPVPi5xwMXOXM1sEDSosYJ8tfzbf/ItoGLgNc3vP/C/PmFTcO3m6/tO2z/FMD2/cA6YO/8Pe8AzrY9nY9fN/PT9rXAlqZ2zQdeBlyQT7fZ9sMVZBJCCANnTofxZwEH2T6k1UhJPwB2azHqTNvfaTPffYH7Gl5P5sPWNk0z2WIagKfYXgtge62kJ3c7X0mHAXOBO/NBzwROlPQG4JfAGTMFZxbPyKf7gqSDgeuAd9v+VZv3hBDCUOpURNqy/dKCb1Wr2RWYZofmm2/dfBFYNrPlAewEbLR9qKQ3AsuBdr/XHOAFwLtsXyPpXLJi+5cd2hZCCEOnVO8sST/Id3U1P17Z4a2TZMchZiwG7m8xzeJZpnlgZvdX/nNdp/nmu6H+BfhAvqur8XO+kT+/DHh+F22ftH1N/vpSsqISQggjp1MReYzWu6uAbEvE9iEtHu12ZQGsBN6W94x6MfDIzO6phnmvBR6T9OK8V9bbgH9ueP+y/PmypuHbzVfSXLICcZHtrze15VvAK/LnLwfuaNdw278A7pP0nHzQUcCtHX7fEEIYSsqOWbeZQPoK2dr5v9r+00o+NCsK55H1lHocOMX2qnzcDTPHYCQdCqwAdgH+lWwXkiXtCVwCPA24F/h92w/NNl9JbwG+AKxuaMbJtm/Ie2l9OZ/XBuCPbN8o6anAKmA+MJ2Pe57tRyUdAnye7NjKXfnnrK8imxBCGCQdi0gIIYQwmzhjPYQQQmFRREIIIRQWRSSEEEJhUURCCCEUFkUkhBBCYVFEQgghFBZFJIQQQmFRREIIIRT2/wFztGyuif/0zwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] @@ -625,7 +625,7 @@ { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 15, @@ -757,7 +757,7 @@ { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 18, @@ -871,7 +871,7 @@ { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 21, @@ -1014,7 +1014,7 @@ { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 24, @@ -1146,7 +1146,7 @@ { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 27, @@ -1272,7 +1272,7 @@ { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 30, @@ -1386,7 +1386,7 @@ { "data": { "text/plain": [ - "" + "" ] }, "execution_count": 33, diff --git a/swe/interface.py b/swe/interface.py index c18430cbd..08f51ae15 100644 --- a/swe/interface.py +++ b/swe/interface.py @@ -3,8 +3,8 @@ @njit(cache=True) -def states(idir, qx, qy, ng, dx, dt, - ih, iu, iv, ix, nvar, nspec, +def states(idir, ng, dx, dt, + ih, iu, iv, ix, nspec, g, qv, dqv): r""" @@ -19,7 +19,7 @@ def states(idir, qx, qy, ng, dx, dt, We need the left and right eigenvectors and the eigenvalues for the system projected along the x-direction - Taking our state vector as :math:`Q = (rho, u, v, p)^T`, the eigenvalues + Taking our state vector as :math:`Q = (\rho, u, v, p)^T`, the eigenvalues are :math:`u - c`, :math:`u`, :math:`u + c`. We look at the equations of hydrodynamics in a split fashion -- @@ -61,8 +61,6 @@ def states(idir, qx, qy, ng, dx, dt, ---------- idir : int Are we predicting to the edges in the x-direction (1) or y-direction (2)? - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells dx : float @@ -72,8 +70,6 @@ def states(idir, qx, qy, ng, dx, dt, ih, iu, iv, ix : int Indices of the height, x-velocity, y-velocity and species in the state vector - nvar : int - The total number of variables in the state vector nspec : int The number of species g : float @@ -89,8 +85,11 @@ def states(idir, qx, qy, ng, dx, dt, State vector predicted to the left and right edges """ - q_l = np.zeros((qx, qy, nvar)) - q_r = np.zeros((qx, qy, nvar)) + qx, qy, nvar = qv.shape + + q_l = np.zeros_like(qv) + q_r = np.zeros_like(qv) + lvec = np.zeros((nvar, nvar)) rvec = np.zeros((nvar, nvar)) e_val = np.zeros(nvar) @@ -209,8 +208,8 @@ def states(idir, qx, qy, ng, dx, dt, @njit(cache=True) -def riemann_roe(idir, qx, qy, ng, - nvar, ih, ixmom, iymom, ihX, nspec, +def riemann_roe(idir, ng, + ih, ixmom, iymom, ihX, nspec, lower_solid, upper_solid, g, U_l, U_r): r""" @@ -221,12 +220,8 @@ def riemann_roe(idir, qx, qy, ng, ---------- idir : int Are we predicting to the edges in the x-direction (1) or y-direction (2)? - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells - nvar : int - The number of variables in the state vector ih, ixmom, iymom, ihX : int The indices of the height, x-momentum, y-momentum and height*species fractions in the conserved state vector. nspec : int @@ -243,6 +238,9 @@ def riemann_roe(idir, qx, qy, ng, out : ndarray Conserved flux """ + + qx, qy, nvar = U_l.shape + F = np.zeros((qx, qy, nvar)) smallc = 1.e-10 @@ -328,9 +326,9 @@ def riemann_roe(idir, qx, qy, ng, K_roe[n, :] = 0.0 K_roe[n, n] = 1.0 - F[i, j, :] = consFlux(idir, g, ih, ixmom, iymom, ihX, nvar, nspec, + F[i, j, :] = consFlux(idir, g, ih, ixmom, iymom, ihX, nspec, U_l[i, j, :]) - F_r = consFlux(idir, g, ih, ixmom, iymom, ihX, nvar, nspec, + F_r = consFlux(idir, g, ih, ixmom, iymom, ihX, nspec, U_r[i, j, :]) F[i, j, :] = 0.5 * (F[i, j, :] + F_r) @@ -358,8 +356,8 @@ def riemann_roe(idir, qx, qy, ng, @njit(cache=True) -def riemann_hllc(idir, qx, qy, ng, - nvar, ih, ixmom, iymom, ihX, nspec, +def riemann_hllc(idir, ng, + ih, ixmom, iymom, ihX, nspec, lower_solid, upper_solid, g, U_l, U_r): r""" @@ -371,12 +369,8 @@ def riemann_hllc(idir, qx, qy, ng, ---------- idir : int Are we predicting to the edges in the x-direction (1) or y-direction (2)? - qx, qy : int - The dimensions of the grid. ng : int The number of ghost cells - nvar : int - The number of variables in the state vector ih, ixmom, iymom, ihX : int The indices of the height, x-momentum, y-momentum and height*species fractions in the conserved state vector. nspec : int @@ -394,6 +388,8 @@ def riemann_hllc(idir, qx, qy, ng, Conserved flux """ + qx, qy, nvar = U_l.shape + F = np.zeros((qx, qy, nvar)) smallc = 1.e-10 @@ -469,7 +465,7 @@ def riemann_hllc(idir, qx, qy, ng, # R region U_state[:] = U_r[i, j, :] - F[i, j, :] = consFlux(idir, g, ih, ixmom, iymom, ihX, nvar, nspec, + F[i, j, :] = consFlux(idir, g, ih, ixmom, iymom, ihX, nspec, U_state) elif (S_r > 0.0 and S_c <= 0): @@ -491,7 +487,7 @@ def riemann_hllc(idir, qx, qy, ng, U_r[i, j, ihX:ihX + nspec] / h_r # find the flux on the right interface - F[i, j, :] = consFlux(idir, g, ih, ixmom, iymom, ihX, nvar, nspec, + F[i, j, :] = consFlux(idir, g, ih, ixmom, iymom, ihX, nspec, U_r[i, j, :]) # correct the flux @@ -516,7 +512,7 @@ def riemann_hllc(idir, qx, qy, ng, U_l[i, j, ihX:ihX + nspec] / h_l # find the flux on the left interface - F[i, j, :] = consFlux(idir, g, ih, ixmom, iymom, ihX, nvar, nspec, + F[i, j, :] = consFlux(idir, g, ih, ixmom, iymom, ihX, nspec, U_l[i, j, :]) # correct the flux @@ -526,13 +522,13 @@ def riemann_hllc(idir, qx, qy, ng, # L region U_state[:] = U_l[i, j, :] - F[i, j, :] = consFlux(idir, g, ih, ixmom, iymom, ihX, nvar, nspec, + F[i, j, :] = consFlux(idir, g, ih, ixmom, iymom, ihX, nspec, U_state) return F @njit(cache=True) -def consFlux(idir, g, ih, ixmom, iymom, ihX, nvar, nspec, U_state): +def consFlux(idir, g, ih, ixmom, iymom, ihX, nspec, U_state): r""" Calculate the conserved flux for the shallow water equations. In the x-direction, this is given by:: @@ -549,8 +545,6 @@ def consFlux(idir, g, ih, ixmom, iymom, ihX, nvar, nspec, U_state): Graviational acceleration ih, ixmom, iymom, ihX : int The indices of the height, x-momentum, y-momentum, height*species fraction in the conserved state vector. - nvar : int - The number of variables in the state vector nspec : int The number of species U_state : ndarray @@ -562,7 +556,7 @@ def consFlux(idir, g, ih, ixmom, iymom, ihX, nvar, nspec, U_state): Conserved flux """ - F = np.zeros(nvar) + F = np.zeros_like(U_state) u = U_state[ixmom] / U_state[ih] v = U_state[iymom] / U_state[ih] diff --git a/swe/interface_f.f90 b/swe/interface_f.f90 deleted file mode 100644 index e2c5c0d15..000000000 --- a/swe/interface_f.f90 +++ /dev/null @@ -1,594 +0,0 @@ -subroutine states(idir, qx, qy, ng, dx, dt, & - ih, iu, iv, ix, nvar, nspec, & - g, & - qv, dqv, & - q_l, q_r) - - implicit none - - integer, intent(in) :: idir - integer, intent(in) :: qx, qy, ng - double precision, intent(in) :: dx, dt - integer, intent(in) :: ih, iu, iv, ix, nvar, nspec - double precision, intent(in) :: g - - ! 0-based indexing to match python - double precision, intent(inout) :: qv(0:qx-1, 0:qy-1, 0:nvar-1) - double precision, intent(inout) :: dqv(0:qx-1, 0:qy-1, 0:nvar-1) - - double precision, intent( out) :: q_l(0:qx-1, 0:qy-1, 0:nvar-1) - double precision, intent( out) :: q_r(0:qx-1, 0:qy-1, 0:nvar-1) - -!f2py depend(qx, qy, nvar) :: qv, dqv -!f2py depend(qx, qy, nvar) :: q_l, q_r -!f2py intent(in) :: qv, dqv -!f2py intent(out) :: q_l, q_r - - ! predict the cell-centered state to the edges in one-dimension - ! using the reconstructed, limited slopes. - ! - ! We follow the convection here that V_l[i] is the left state at the - ! i-1/2 interface and V_l[i+1] is the left state at the i+1/2 - ! interface. - ! - ! - ! We need the left and right eigenvectors and the eigenvalues for - ! the system projected along the x-direction - ! - ! Taking our state vector as Q = (rho, u, v, p)^T, the eigenvalues - ! are u - c, u, u + c. - ! - ! We look at the equations of hydrodynamics in a split fashion -- - ! i.e., we only consider one dimension at a time. - ! - ! Considering advection in the x-direction, the Jacobian matrix for - ! the primitive variable formulation of the Euler equations - ! projected in the x-direction is: - ! - ! / u 0 0 \ - ! | g u 0 | - ! A = \ 0 0 u / - ! - ! The right eigenvectors are - ! - ! / h \ / 0 \ / h \ - ! r1 = | -c | r2 = | 0 | r3 = | c | - ! \ 0 / \ 1 / \ 0 / - ! - ! The left eigenvectors are - ! - ! l1 = ( 1/(2h), -h/(2hc), 0 ) - ! l2 = ( 0, 0, 1 ) - ! l3 = ( -1/(2h), -h/(2hc), 0 ) - ! - ! The fluxes are going to be defined on the left edge of the - ! computational zones - ! - ! | | | | - ! | | | | - ! -+------+------+------+------+------+------+-- - ! | i-1 | i | i+1 | - ! ^ ^ ^ - ! q_l,i q_r,i q_l,i+1 - ! - ! q_r,i and q_l,i+1 are computed using the information in zone i,j. - - integer :: ilo, ihi, jlo, jhi - integer :: nx, ny - integer :: i, j, n, m - - double precision :: dq(0:nvar-1), q(0:nvar-1) - double precision :: lvec(0:nvar-1,0:nvar-1), rvec(0:nvar-1,0:nvar-1) - double precision :: eval(0:nvar-1) - double precision :: betal(0:nvar-1), betar(0:nvar-1) - - double precision :: dtdx, dtdx3 - double precision :: cs - - double precision :: sum, sum_l, sum_r, factor - - integer :: ns - - nx = qx - 2*ng; ny = qy - 2*ng - ilo = ng; ihi = ng+nx-1; jlo = ng; jhi = ng+ny-1 - - ns = nvar - nspec - - dtdx = dt/dx - dtdx3 = 0.33333d0*dtdx - - ! this is the loop over zones. For zone i, we see q_l[i+1] and q_r[i] - do j = jlo-2, jhi+2 - do i = ilo-2, ihi+2 - - dq(:) = dqv(i,j,:) - q(:) = qv(i,j,:) - - cs = sqrt(g*q(ih)) - - lvec(:,:) = 0.0d0 - rvec(:,:) = 0.0d0 - eval(:) = 0.0d0 - - ! compute the eigenvalues and eigenvectors - if (idir == 1) then - eval(0:ns-1) = [q(iu) - cs, q(iu), q(iu) + cs] - - lvec(0,0:ns-1) = [ cs, -q(ih), 0.0d0 ] - lvec(1,0:ns-1) = [ 0.0d0, 0.0d0, 1.0d0 ] - lvec(2,0:ns-1) = [ cs, q(ih), 0.0d0 ] - - rvec(0,0:ns-1) = [ q(ih), -cs, 0.0d0 ] - rvec(1,0:ns-1) = [ 0.0d0, 0.0d0, 1.0d0 ] - rvec(2,0:ns-1) = [ q(ih), cs, 0.0d0 ] - - ! now the species -- they only have a 1 in their corresponding slot - eval(ns:) = q(iu) - do n = ix, ix-1+nspec - lvec(n,n) = 1.0 - rvec(n,n) = 1.0 - enddo - - ! multiply by scaling factors - lvec(0,:) = lvec(0,:) * 0.50d0 / (cs * q(ih)) - lvec(2,:) = -lvec(2,:) * 0.50d0 / (cs * q(ih)) - else - eval(0:ns-1) = [q(iv) - cs, q(iv), q(iv) + cs] - - lvec(0,0:ns-1) = [ cs, 0.0d0, -q(ih) ] - lvec(1,0:ns-1) = [ 0.0d0, 1.0d0, 0.0d0 ] - lvec(2,0:ns-1) = [ cs, 0.0d0, q(ih) ] - - rvec(0,0:ns-1) = [ q(ih), 0.0d0, -cs ] - rvec(1,0:ns-1) = [ 0.0d0, 1.0d0, 0.0d0 ] - rvec(2,0:ns-1) = [ q(ih), 0.0d0, cs ] - - ! now the species -- they only have a 1 in their corresponding slot - eval(ns:) = q(iv) - do n = ix, ix-1+nspec - lvec(n,n) = 1.0 - rvec(n,n) = 1.0 - enddo - - ! multiply by scaling factors - lvec(0,:) = lvec(0,:) * 0.50d0 / (cs * q(ih)) - lvec(2,:) = -lvec(2,:) * 0.50d0 / (cs * q(ih)) - - endif - - ! define the reference states - if (idir == 1) then - ! this is one the right face of the current zone, - ! so the fastest moving eigenvalue is eval[2] = u + c - factor = 0.5d0*(1.0d0 - dtdx*max(eval(2), 0.0d0)) - q_l(i+1,j,:) = q(:) + factor*dq(:) - - ! left face of the current zone, so the fastest moving - ! eigenvalue is eval[3] = u - c - factor = 0.5d0*(1.0d0 + dtdx*min(eval(0), 0.0d0)) - q_r(i, j,:) = q(:) - factor*dq(:) - - else - - factor = 0.5d0*(1.0d0 - dtdx*max(eval(2), 0.0d0)) - q_l(i,j+1,:) = q(:) + factor*dq(:) - - factor = 0.5d0*(1.0d0 + dtdx*min(eval(0), 0.0d0)) - q_r(i,j, :) = q(:) - factor*dq(:) - - endif - - ! compute the Vhat functions - do m = 0, nvar-1 - sum = dot_product(lvec(m,:),dq(:)) - - betal(m) = dtdx3*(eval(2) - eval(m))*(sign(1.0d0,eval(m)) + 1.0d0)*sum - betar(m) = dtdx3*(eval(0) - eval(m))*(1.0d0 - sign(1.0d0,eval(m)))*sum - enddo - - ! construct the states - do m = 0, nvar-1 - sum_l = dot_product(betal(:),rvec(:,m)) - sum_r = dot_product(betar(:),rvec(:,m)) - - if (idir == 1) then - q_l(i+1,j,m) = q_l(i+1,j,m) + sum_l - q_r(i, j,m) = q_r(i, j,m) + sum_r - else - 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 - endif - - enddo - - enddo - enddo - -end subroutine states - - -subroutine riemann_Roe(idir, qx, qy, ng, & - nvar, ih, ixmom, iymom, ihX, nspec, & - lower_solid, upper_solid, & - g, U_l, U_r, F) - - implicit none - - integer, intent(in) :: idir - integer, intent(in) :: qx, qy, ng - integer, intent(in) :: nvar, ih, ixmom, iymom, ihX, nspec - integer, intent(in) :: lower_solid, upper_solid - double precision, intent(in) :: g - - ! 0-based indexing to match python - double precision, intent(inout) :: U_l(0:qx-1,0:qy-1,0:nvar-1) - double precision, intent(inout) :: U_r(0:qx-1,0:qy-1,0:nvar-1) - double precision, intent( out) :: F(0:qx-1,0:qy-1,0:nvar-1) - -!f2py depend(qx, qy, nvar) :: U_l, U_r -!f2py intent(in) :: U_l, U_r -!f2py intent(out) :: F - - ! This is the Roe Riemann solver with entropy fix. The implementation - ! follows Toro's SWE book and the clawpack 2d SWE Roe solver. - - integer :: ilo, ihi, jlo, jhi - integer :: nx, ny - integer :: i, j, n, m - integer :: ns - - double precision, parameter :: smallc = 1.e-10 - double precision, parameter :: tol = 0.1d-1 ! entropy fix parameter - ! Note that I've basically assumed that cfl = 0.1 here to get away with - ! not passing dx/dt or cfl to this function. If this isn't the case, will need - ! to pass one of these to the function or else things will go wrong. - - double precision :: h_l, un_l, ut_l, c_l - double precision :: h_r, un_r, ut_r, c_r - double precision :: h_star, u_star, c_star - double precision :: xn(nspec) - - double precision :: U_roe(0:nvar-1), c_roe, un_roe - double precision :: lambda_roe(0:nvar-1), K_roe(0:nvar-1, 0:nvar-1) - double precision :: alpha_roe(0:nvar-1), delta(0:nvar-1), F_r(0:nvar-1) - - nx = qx - 2*ng; ny = qy - 2*ng - ilo = ng; ihi = ng+nx-1; jlo = ng; jhi = ng+ny-1 - ns = nvar - nspec - - do j = jlo-1, jhi+1 - do i = ilo-1, ihi+1 - - ! primitive variable states - h_l = U_l(i,j,ih) - - ! un = normal velocity; ut = transverse velocity - if (idir == 1) then - un_l = U_l(i,j,ixmom)/h_l - ut_l = U_l(i,j,iymom)/h_l - else - un_l = U_l(i,j,iymom)/h_l - ut_l = U_l(i,j,ixmom)/h_l - endif - - h_r = U_r(i,j,ih) - - if (idir == 1) then - un_r = U_r(i,j,ixmom)/h_r - ut_r = U_r(i,j,iymom)/h_r - else - un_r = U_r(i,j,iymom)/h_r - ut_r = U_r(i,j,ixmom)/h_r - endif - - ! compute the sound speeds - c_l = max(smallc, sqrt(g*h_l)) - c_r = max(smallc, sqrt(g*h_r)) - - ! Calculate the Roe averages - U_roe = (U_l(i,j,:)/sqrt(h_l) + U_r(i,j,:)/sqrt(h_r)) / & - (sqrt(h_l) + sqrt(h_r)) - - U_roe(ih) = sqrt(h_l * h_r) - c_roe = sqrt(0.5d0 * (c_l**2 + c_r**2)) - - delta(:) = U_r(i,j,:)/h_r - U_l(i,j,:)/h_l - delta(ih) = h_r - h_l - - ! evalues and right evectors - if (idir == 1) then - un_roe = U_roe(ixmom) - else - un_roe = U_roe(iymom) - endif - - K_roe(:,:) = 0.0d0 - - lambda_roe(0:2) = [un_roe - c_roe, un_roe, un_roe + c_roe] - if (idir == 1) then - alpha_roe(0:2) = [0.5d0*(delta(ih) - U_roe(ih)/c_roe*delta(ixmom)), & - U_roe(ih) * delta(iymom), & - 0.5d0*(delta(ih) + U_roe(ih)/c_roe*delta(ixmom))] - - K_roe(0, 0:2) = [1.0d0, un_roe - c_roe, U_roe(iymom)] - K_roe(1, 0:2) = [0.0d0, 0.0d0, 1.0d0] - K_roe(2, 0:2) = [1.0d0, un_roe + c_roe, U_roe(iymom)] - else - alpha_roe(0:2) = [0.5d0*(delta(ih) - U_roe(ih)/c_roe*delta(iymom)), & - U_roe(ih) * delta(ixmom), & - 0.5d0*(delta(ih) + U_roe(ih)/c_roe*delta(iymom))] - - K_roe(0, 0:2) = [1.0d0, U_roe(ixmom), un_roe - c_roe] - K_roe(1, 0:2) = [0.0d0, 1.0d0, 0.0d0] - K_roe(2, 0:2) = [1.0d0, U_roe(ixmom), un_roe + c_roe] - endif - - lambda_roe(ns:) = un_roe - alpha_roe(ns:) = U_roe(ih) * delta(ns:) - do n = ns, nvar-1 - K_roe(n,:) = 0.0d0 - K_roe(n,n) = 1.0d0 - enddo - - call consFlux(idir, g, ih, ixmom, iymom, ihX, nvar, nspec, & - U_l(i,j,:), F(i,j,:)) - call consFlux(idir, g, ih, ixmom, iymom, ihX, nvar, nspec, & - U_r(i,j,:), F_r) - - F(i,j,:) = 0.5d0 * (F(i,j,:) + F_r) - - h_star = 1.0d0 / g * (0.5d0 * (c_l + c_r) + 0.25d0 * (un_l - un_r))**2 - u_star = 0.5d0 * (un_l + un_r) + c_l - c_r - - c_star = sqrt(g * h_star) - - ! modified evalues for entropy fix - if (abs(lambda_roe(0)) < tol) then - lambda_roe(0) = lambda_roe(0) * (u_star - c_star - lambda_roe(0)) / & - (u_star - c_star - (un_l - c_l)) - endif - if (abs(lambda_roe(2)) < tol) then - lambda_roe(2) = lambda_roe(2) * (u_star + c_star - lambda_roe(2)) / & - (u_star + c_star - (un_r + c_r)) - endif - - do n = 0, nvar-1 - do m = 0, nvar-1 - F(i,j,n) = F(i,j,n) - & - 0.5d0 * alpha_roe(m) * abs(lambda_roe(m)) * K_roe(m,n) - enddo - enddo - - enddo - enddo -end subroutine riemann_Roe - - -subroutine riemann_HLLC(idir, qx, qy, ng, & - nvar, ih, ixmom, iymom, ihX, nspec, & - lower_solid, upper_solid, & - g, U_l, U_r, F) - - implicit none - - integer, intent(in) :: idir - integer, intent(in) :: qx, qy, ng - integer, intent(in) :: nvar, ih, ixmom, iymom, ihX, nspec - integer, intent(in) :: lower_solid, upper_solid - double precision, intent(in) :: g - - ! 0-based indexing to match python - double precision, intent(inout) :: U_l(0:qx-1,0:qy-1,0:nvar-1) - double precision, intent(inout) :: U_r(0:qx-1,0:qy-1,0:nvar-1) - double precision, intent( out) :: F(0:qx-1,0:qy-1,0:nvar-1) - -!f2py depend(qx, qy, nvar) :: U_l, U_r -!f2py intent(in) :: U_l, U_r -!f2py intent(out) :: F - - ! this is the HLLC Riemann solver. The implementation follows - ! directly out of Toro's book. Note: this does not handle the - ! transonic rarefaction. - - integer :: ilo, ihi, jlo, jhi - integer :: nx, ny - integer :: i, j - - double precision, parameter :: smallc = 1.e-10 - double precision, parameter :: smallrho = 1.e-10 - double precision, parameter :: smallp = 1.e-10 - - double precision :: h_l, un_l, ut_l - double precision :: h_r, un_r, ut_r - double precision :: xn(nspec) - - double precision :: hstar_l, hstar_r, h_avg - double precision :: hstar, ustar, u_avg - double precision :: S_l, S_r, S_c - double precision :: c_l, c_r, c_avg - - double precision :: U_state(0:nvar-1) - double precision :: HLLCfactor - - nx = qx - 2*ng; ny = qy - 2*ng - ilo = ng; ihi = ng+nx-1; jlo = ng; jhi = ng+ny-1 - - do j = jlo-1, jhi+1 - do i = ilo-1, ihi+1 - - ! primitive variable states - h_l = U_l(i,j,ih) - - ! un = normal velocity; ut = transverse velocity - if (idir == 1) then - un_l = U_l(i,j,ixmom)/h_l - ut_l = U_l(i,j,iymom)/h_l - else - un_l = U_l(i,j,iymom)/h_l - ut_l = U_l(i,j,ixmom)/h_l - endif - - h_r = U_r(i,j,ih) - - if (idir == 1) then - un_r = U_r(i,j,ixmom)/h_r - ut_r = U_r(i,j,iymom)/h_r - else - un_r = U_r(i,j,iymom)/h_r - ut_r = U_r(i,j,ixmom)/h_r - endif - - ! compute the sound speeds - c_l = max(smallc, sqrt(g*h_l)) - c_r = max(smallc, sqrt(g*h_r)) - - ! Estimate the star quantities -- use one of three methods to - ! do this -- the primitive variable Riemann solver, the two - ! shock approximation, or the two rarefaction approximation. - ! Pick the method based on the pressure states at the - ! interface. - - h_avg = 0.5*(h_l + h_r) - c_avg = 0.5*(c_l + c_r) - u_avg = 0.5*(un_l + un_r) - - hstar = h_avg - 0.25d0 * (un_r - un_l) * h_avg / c_avg - ustar = u_avg - (h_r - h_l) * c_avg / h_avg - - ! estimate the nonlinear wave speeds - - if (hstar <= h_l) then - ! rarefaction - S_l = un_l - c_l - else - ! shock - S_l = un_l - c_l * sqrt(0.5d0 * (hstar+h_l) * hstar) / h_l - endif - - if (hstar <= h_r) then - ! rarefaction - S_r = un_r + c_r - else - ! shock - S_r = un_r + c_r*sqrt(0.5d0 * (hstar+h_r) * hstar) / h_r - endif - - S_c = (S_l*h_r*(un_r-S_r) - S_r*h_l*(un_l-S_l)) / & - (h_r*(un_r-S_r) - h_l*(un_l-S_l)) - - ! figure out which region we are in and compute the state and - ! the interface fluxes using the HLLC Riemann solver - if (S_r <= 0.0d0) then - ! R region - U_state(:) = U_r(i,j,:) - - call consFlux(idir, g, ih, ixmom, iymom, ihX, nvar, nspec, & - U_state, F(i,j,:)) - - else if (S_r > 0.0d0 .and. S_c <= 0) then - ! R* region - HLLCfactor = h_r*(S_r - un_r)/(S_r - S_c) - - U_state(ih) = HLLCfactor - - if (idir == 1) then - U_state(ixmom) = HLLCfactor*S_c - U_state(iymom) = HLLCfactor*ut_r - else - U_state(ixmom) = HLLCfactor*ut_r - U_state(iymom) = HLLCfactor*S_c - endif - - ! species - if (nspec > 0) then - U_state(ihX:ihX-1+nspec) = HLLCfactor*U_r(i,j,ihX:ihX-1+nspec)/h_r - endif - - ! find the flux on the right interface - call consFlux(idir, g, ih, ixmom, iymom, ihX, nvar, nspec, & - U_r(i,j,:), F(i,j,:)) - - ! correct the flux - F(i,j,:) = F(i,j,:) + S_r*(U_state(:) - U_r(i,j,:)) - - else if (S_c > 0.0d0 .and. S_l < 0.0) then - ! L* region - HLLCfactor = h_l*(S_l - un_l)/(S_l - S_c) - - U_state(ih) = HLLCfactor - - if (idir == 1) then - U_state(ixmom) = HLLCfactor*S_c - U_state(iymom) = HLLCfactor*ut_l - else - U_state(ixmom) = HLLCfactor*ut_l - U_state(iymom) = HLLCfactor*S_c - endif - - ! species - if (nspec > 0) then - U_state(ihX:ihX-1+nspec) = HLLCfactor*U_l(i,j,ihX:ihX-1+nspec)/h_l - endif - - ! find the flux on the left interface - call consFlux(idir, g, ih, ixmom, iymom, ihX, nvar, nspec, & - U_l(i,j,:), F(i,j,:)) - - ! correct the flux - F(i,j,:) = F(i,j,:) + S_l*(U_state(:) - U_l(i,j,:)) - - else - ! L region - U_state(:) = U_l(i,j,:) - - call consFlux(idir, g, ih, ixmom, iymom, ihX, nvar, nspec, & - U_state, F(i,j,:)) - - endif - - ! we should deal with solid boundaries somehow here - - enddo - enddo -end subroutine riemann_HLLC - - -subroutine consFlux(idir, g, ih, ixmom, iymom, ihX, nvar, nspec, U_state, F) - - implicit none - - integer, intent(in) :: idir - double precision, intent(in) :: g - integer, intent(in) :: ih, ixmom, iymom, ihX, nvar, nspec - double precision, intent(in) :: U_state(0:nvar-1) - double precision, intent(out) :: F(0:nvar-1) - -! Calculate the conserved flux for the shallow water equations. In the -! x-direction, this is given by -! -! / hu \ -! F = | hu^2 + gh^2/2 | -! \ huv / - - double precision :: u, v - - u = U_state(ixmom)/U_state(ih) - v = U_state(iymom)/U_state(ih) - - if (idir == 1) then - F(ih) = U_state(ih)*u - F(ixmom) = U_state(ixmom)*u + 0.5d0 * g * U_state(ih)**2 - F(iymom) = U_state(iymom)*u - if (nspec > 0) then - F(ihX:ihX-1+nspec) = U_state(ihX:ihX-1+nspec)*u - endif - else - F(ih) = U_state(ih)*v - F(ixmom) = U_state(ixmom)*v - F(iymom) = U_state(iymom)*v + 0.5d0 * g * U_state(ih)**2 - if (nspec > 0) then - F(ihX:ihX-1+nspec) = U_state(ihX:ihX-1+nspec)*v - endif - endif - -end subroutine consFlux diff --git a/swe/tests/dam_x_0081.h5 b/swe/tests/dam_x_0081.h5 index 2f86bdf14..d163692df 100644 Binary files a/swe/tests/dam_x_0081.h5 and b/swe/tests/dam_x_0081.h5 differ diff --git a/swe/unsplit_fluxes.py b/swe/unsplit_fluxes.py index b22fd07eb..5011b5201 100644 --- a/swe/unsplit_fluxes.py +++ b/swe/unsplit_fluxes.py @@ -212,9 +212,9 @@ def unsplit_fluxes(my_data, my_aux, rp, ivars, solid, tc, dt): tm_states = tc.timer("interfaceStates") tm_states.begin() - V_l, V_r = ifc.states(1, myg.qx, myg.qy, myg.ng, myg.dx, dt, + V_l, V_r = ifc.states(1, myg.ng, myg.dx, dt, ivars.ih, ivars.iu, ivars.iv, ivars.ix, - ivars.nvar, ivars.naux, + ivars.naux, g, q, ldx) @@ -231,9 +231,9 @@ def unsplit_fluxes(my_data, my_aux, rp, ivars, solid, tc, dt): # left and right primitive variable states tm_states.begin() - _V_l, _V_r = ifc.states(2, myg.qx, myg.qy, myg.ng, myg.dy, dt, + _V_l, _V_r = ifc.states(2, myg.ng, myg.dy, dt, ivars.ih, ivars.iu, ivars.iv, ivars.ix, - ivars.nvar, ivars.naux, + ivars.naux, g, q, ldy) V_l = ai.ArrayIndexer(d=_V_l, grid=myg) @@ -260,13 +260,13 @@ def unsplit_fluxes(my_data, my_aux, rp, ivars, solid, tc, dt): else: msg.fail("ERROR: Riemann solver undefined") - _fx = riemannFunc(1, myg.qx, myg.qy, myg.ng, - ivars.nvar, ivars.ih, ivars.ixmom, ivars.iymom, ivars.ihx, ivars.naux, + _fx = riemannFunc(1, myg.ng, + ivars.ih, ivars.ixmom, ivars.iymom, ivars.ihx, ivars.naux, solid.xl, solid.xr, g, U_xl, U_xr) - _fy = riemannFunc(2, myg.qx, myg.qy, myg.ng, - ivars.nvar, ivars.ih, ivars.ixmom, ivars.iymom, ivars.ihx, ivars.naux, + _fy = riemannFunc(2, myg.ng, + ivars.ih, ivars.ixmom, ivars.iymom, ivars.ihx, ivars.naux, solid.yl, solid.yr, g, U_yl, U_yr) @@ -359,13 +359,13 @@ def unsplit_fluxes(my_data, my_aux, rp, ivars, solid, tc, dt): tm_riem.begin() - _fx = riemannFunc(1, myg.qx, myg.qy, myg.ng, - ivars.nvar, ivars.ih, ivars.ixmom, ivars.iymom, ivars.ihx, ivars.naux, + _fx = riemannFunc(1, myg.ng, + ivars.ih, ivars.ixmom, ivars.iymom, ivars.ihx, ivars.naux, solid.xl, solid.xr, g, U_xl, U_xr) - _fy = riemannFunc(2, myg.qx, myg.qy, myg.ng, - ivars.nvar, ivars.ih, ivars.ixmom, ivars.iymom, ivars.ihx, ivars.naux, + _fy = riemannFunc(2, myg.ng, + ivars.ih, ivars.ixmom, ivars.iymom, ivars.ihx, ivars.naux, solid.yl, solid.yr, g, U_yl, U_yr)