numba and fortran don't quite agree for incompressible and LM solvers (as you'd expect), so have made separate solver-test notebooks for them. This is not ideal, but can't find a way to pass in a relative error value to nbval

Author: harpolea

.travis.yml

@@ -5,8 +5,12 @@ python:
 
 matrix:
     include:
-    - env: MK_ARGS=""
-    - env: MK_ARGS=fortran
+    - env:
+        - MK_ARGS=""
+        - NB_IGNORE="solver-test-fortran.ipynb"
+    - env:
+        - MK_ARGS=fortran
+        - NB_IGNORE="solver-test.ipynb"
 
 before_install:
   - export PATH=$(echo $PATH | tr ':' "\n" | sed '/\/opt\/python/d' | tr "\n" ":" | sed "s|::|:|g")
@@ -28,7 +32,7 @@ before_script:
 
 script:
   - flake8 --ignore=W504,F403,E226,E126,E127,E128,E129,E241,E741 .
-  - pytest -v --cov=. --cov-config .coveragerc --nbval --ignore=docs --ignore=./multigrid/variable_coeff_elliptic.ipynb --ignore=examples/mesh --ignore=examples/multigrid
+  - pytest -v --cov=. --cov-config .coveragerc --nbval --ignore=docs --ignore=./multigrid/variable_coeff_elliptic.ipynb --ignore=examples/mesh --ignore=examples/multigrid --ignore=$NB_IGNORE
 #   - travis-sphinx build
 
 # after_success:

advection_fv4/interface.py

@@ -32,8 +32,8 @@ def states(a, qx, qy, ng, idir):
     # we need interface values on all faces of the domain
     if (idir == 1):
 
-        for j in range(jlo - 1, jhi + 2):
-            for i in range(ilo - 2, ihi + 4):
+        for j in range(jlo - 1, jhi + 1):
+            for i in range(ilo - 2, ihi + 3):
 
                 # interpolate to the edges
                 a_int[i, j] = (7.0 / 12.0) * (a[i - 1, j] + a[i, j]) - \
@@ -42,8 +42,8 @@ def states(a, qx, qy, ng, idir):
                 al[i, j] = a_int[i, j]
                 ar[i, j] = a_int[i, j]
 
-        for j in range(jlo - 1, jhi + 2):
-            for i in range(ilo - 2, ihi + 4):
+        for j in range(jlo - 1, jhi + 1):
+            for i in range(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]
@@ -52,19 +52,19 @@ def states(a, qx, qy, ng, idir):
                 d2af[i, j] = 6.0 * (a_int[i, j] - 2.0 *
                                     a[i, j] + a_int[i + 1, j])
 
-        for j in range(jlo - 1, jhi + 2):
-            for i in range(ilo - 3, ihi + 4):
+        for j in range(jlo - 1, jhi + 1):
+            for i in range(ilo - 3, ihi + 3):
                 d2ac[i, j] = a[i - 1, j] - 2.0 * a[i, j] + a[i + 1, j]
 
-        for j in range(jlo - 1, jhi + 2):
-            for i in range(ilo - 2, ihi + 4):
+        for j in range(jlo - 1, jhi + 1):
+            for i in range(ilo - 2, ihi + 3):
                 # this lives on the interface
                 d3a[i, j] = d2ac[i, j] - d2ac[i - 1, j]
 
         # this is a look over cell centers, affecting
         # i-1/2,R and i+1/2,L
-        for j in range(jlo - 1, jhi + 2):
-            for i in range(ilo - 1, ihi + 2):
+        for j in range(jlo - 1, jhi + 1):
+            for i in range(ilo - 1, ihi + 1):
 
                 # limit? MC Eq. 24 and 25
                 if (dafm[i, j] * dafp[i, j] <= 0.0 or
@@ -121,8 +121,8 @@ def states(a, qx, qy, ng, idir):
 
     elif (idir == 2):
 
-        for j in range(jlo - 2, jhi + 4):
-            for i in range(ilo - 1, ihi + 2):
+        for j in range(jlo - 2, jhi + 3):
+            for i in range(ilo - 1, ihi + 1):
 
                 # interpolate to the edges
                 a_int[i, j] = (7.0 / 12.0) * (a[i, j - 1] + a[i, j]) - \
@@ -131,8 +131,8 @@ def states(a, qx, qy, ng, idir):
                 al[i, j] = a_int[i, j]
                 ar[i, j] = a_int[i, j]
 
-        for j in range(jlo - 2, jhi + 4):
-            for i in range(ilo - 1, ihi + 2):
+        for j in range(jlo - 2, jhi + 3):
+            for i in range(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]
@@ -141,19 +141,19 @@ def states(a, qx, qy, ng, idir):
                 d2af[i, j] = 6.0 * (a_int[i, j] - 2.0 *
                                     a[i, j] + a_int[i, j + 1])
 
-        for j in range(jlo - 3, jhi + 4):
-            for i in range(ilo - 1, ihi + 2):
+        for j in range(jlo - 3, jhi + 3):
+            for i in range(ilo - 1, ihi + 1):
                 d2ac[i, j] = a[i, j - 1] - 2.0 * a[i, j] + a[i, j + 1]
 
-        for j in range(jlo - 2, jhi + 3):
-            for i in range(ilo - 1, ihi + 2):
+        for j in range(jlo - 2, jhi + 2):
+            for i in range(ilo - 1, ihi + 1):
                 # this lives on the interface
                 d3a[i, j] = d2ac[i, j] - d2ac[i, j - 1]
 
         # this is a look over cell centers, affecting
         # j-1/2,R and j+1/2,L
-        for j in range(jlo - 1, jhi + 2):
-            for i in range(ilo - 1, ihi + 2):
+        for j in range(jlo - 1, jhi + 1):
+            for i in range(ilo - 1, ihi + 1):
 
                 # limit? MC Eq. 24 and 25
                 if (dafm[i, j] * dafp[i, j] <= 0.0 or
@@ -232,8 +232,8 @@ def states_nolimit(a, qx, qy, ng, idir):
     # we need interface values on all faces of the domain
     if (idir == 1):
 
-        for j in range(jlo - 1, jhi + 2):
-            for i in range(ilo - 2, ihi + 4):
+        for j in range(jlo - 1, jhi + 1):
+            for i in range(ilo - 2, ihi + 3):
 
                 # interpolate to the edges
                 a_int[i, j] = (7.0 / 12.0) * (a[i - 1, j] + a[i, j]) - \
@@ -244,8 +244,8 @@ def states_nolimit(a, qx, qy, ng, idir):
 
     elif (idir == 2):
 
-        for j in range(jlo - 2, jhi + 4):
-            for i in range(ilo - 1, ihi + 2):
+        for j in range(jlo - 2, jhi + 3):
+            for i in range(ilo - 1, ihi + 1):
 
                 # interpolate to the edges
                 a_int[i, j] = (7.0 / 12.0) * (a[i, j - 1] + a[i, j]) - \

compressible/interface.py

@@ -87,8 +87,8 @@ def states(idir, qx, qy, ng, dx, dt,
     betar = np.zeros(nvar)
 
     # this is the loop over zones.  For zone i, we see q_l[i+1] and q_r[i]
-    for j in range(jlo - 2, jhi + 3):
-        for i in range(ilo - 2, ihi + 3):
+    for j in range(jlo - 2, jhi + 2):
+        for i in range(ilo - 2, ihi + 2):
 
             dq[:] = dqv[i, j, :]
             q[:] = qv[i, j, :]
@@ -166,7 +166,7 @@ def states(idir, qx, qy, ng, dx, dt,
 
             # compute the Vhat functions
             for m in range(nvar):
-                sum = np.dot(lvec[m, :], dq[:])
+                sum = np.dot(lvec[m, :], dq)
 
                 betal[m] = dtdx4 * (e_val[3] - e_val[m]) * \
                     (np.copysign(1.0, e_val[m]) + 1.0) * sum
@@ -175,8 +175,8 @@ def states(idir, qx, qy, ng, dx, dt,
 
             # construct the states
             for m in range(nvar):
-                sum_l = np.dot(betal[:], rvec[:, m])
-                sum_r = np.dot(betar[:], rvec[:, m])
+                sum_l = np.dot(betal, rvec[:, m])
+                sum_r = np.dot(betar, rvec[:, m])
 
                 if (idir == 1):
                     q_l[i + 1, j, m] = q_l[i + 1, j, m] + sum_l
@@ -237,8 +237,8 @@ def riemann_cgf(idir, qx, qy, ng,
     jlo = ng
     jhi = ng + ny
 
-    for j in range(jlo - 1, jhi + 2):
-        for i in range(ilo - 1, ihi + 2):
+    for j in range(jlo - 1, jhi + 1):
+        for i in range(ilo - 1, ihi + 1):
 
             # primitive variable states
             rho_l = U_l[i, j, idens]
@@ -522,8 +522,8 @@ def riemann_prim(idir, qx, qy, ng,
     jlo = ng
     jhi = ng + ny
 
-    for j in range(jlo - 1, jhi + 2):
-        for i in range(ilo - 1, ihi + 2):
+    for j in range(jlo - 1, jhi + 1):
+        for i in range(ilo - 1, ihi + 1):
 
             # primitive variable states
             rho_l = q_l[i, j, irho]
@@ -755,8 +755,8 @@ def riemann_hllc(idir, qx, qy, ng,
     jlo = ng
     jhi = ng + ny
 
-    for j in range(jlo - 1, jhi + 2):
-        for i in range(ilo - 1, ihi + 2):
+    for j in range(jlo - 1, jhi + 1):
+        for i in range(ilo - 1, ihi + 1):
 
             # primitive variable states
             rho_l = U_l[i, j, idens]
@@ -976,8 +976,7 @@ def consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nvar, nspec, U_stat
     u = U_state[ixmom] / U_state[idens]
     v = U_state[iymom] / U_state[idens]
 
-    p = (U_state[iener] - 0.5 * U_state[idens] *
-         (u * u + v * v)) * (gamma - 1.0)
+    p = (U_state[iener] - 0.5 * U_state[idens] * (u * u + v * v)) * (gamma - 1.0)
 
     if (idir == 1):
         F[idens] = U_state[idens] * u
@@ -1041,8 +1040,8 @@ def artificial_viscosity(qx, qy, ng, dx, dy,
     jlo = ng
     jhi = ng + ny
 
-    for j in range(jlo - 1, jhi + 2):
-        for i in range(ilo - 1, ihi + 2):
+    for j in range(jlo - 1, jhi + 1):
+        for i in range(ilo - 1, ihi + 1):
 
                 # start by computing the divergence on the x-interface.  The
                 # x-difference is simply the difference of the cell-centered

incompressible/incomp_interface.py

@@ -109,8 +109,8 @@ def get_interface_states(qx, qy, ng, dx, dy, dt,
     dtdx = dt / dx
     dtdy = dt / dy
 
-    for j in range(jlo - 2, jhi + 3):
-        for i in range(ilo - 2, ihi + 3):
+    for j in range(jlo - 2, jhi + 2):
+        for i in range(ilo - 2, ihi + 2):
 
             # u on x-edges
             u_xl[i + 1, j] = u[i, j] + 0.5 * \
@@ -156,8 +156,8 @@ def get_interface_states(qx, qy, ng, dx, dy, dt,
     # considered the normal to the interface portion of the predictor.
 
     # add the transverse flux differences to the preliminary interface states
-    for j in range(jlo - 2, jhi + 3):
-        for i in range(ilo - 2, ihi + 3):
+    for j in range(jlo - 2, jhi + 2):
+        for i in range(ilo - 2, ihi + 2):
 
             ubar = 0.5 * (uhat_adv[i, j] + uhat_adv[i + 1, j])
             vbar = 0.5 * (vhat_adv[i, j] + vhat_adv[i, j + 1])
@@ -212,8 +212,8 @@ def upwind(qx, qy, ng, q_l, q_r, s, q_int):
     jlo = ng
     jhi = ng + ny
 
-    for j in range(jlo - 1, jhi + 2):
-        for i in range(ilo - 1, ihi + 2):
+    for j in range(jlo - 1, jhi + 1):
+        for i in range(ilo - 1, ihi + 1):
 
             if (s[i, j] > 0.0):
                 q_int[i, j] = q_l[i, j]
@@ -240,8 +240,8 @@ def riemann(qx, qy, ng, q_l, q_r, s):
     jlo = ng
     jhi = ng + ny
 
-    for j in range(jlo - 1, jhi + 2):
-        for i in range(ilo - 1, ihi + 2):
+    for j in range(jlo - 1, jhi + 1):
+        for i in range(ilo - 1, ihi + 1):
 
             if (q_l[i, j] > 0.0 and q_l[i, j] + q_r[i, j] > 0.0):
                 s[i, j] = q_l[i, j]

lm_atm/LM_atm_interface.py

@@ -235,8 +235,8 @@ def rho_states(qx, qy, ng, dx, dy, dt,
     dtdx = dt / dx
     dtdy = dt / dy
 
-    for j in range(jlo - 2, jhi + 3):
-        for i in range(ilo - 2, ihi + 3):
+    for j in range(jlo - 2, jhi + 2):
+        for i in range(ilo - 2, ihi + 2):
 
             # u on x-edges
             rho_xl[i + 1, j] = rho[i, j] + 0.5 * \
@@ -256,8 +256,8 @@ def rho_states(qx, qy, ng, dx, dy, dt,
 
     # now add the transverse term and the non-advective part of the normal
     # divergence
-    for j in range(jlo - 2, jhi + 3):
-        for i in range(ilo - 2, ihi + 3):
+    for j in range(jlo - 2, jhi + 2):
+        for i in range(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
@@ -333,8 +333,8 @@ def get_interface_states(qx, qy, ng, dx, dy, dt,
     dtdx = dt / dx
     dtdy = dt / dy
 
-    for j in range(jlo - 2, jhi + 3):
-        for i in range(ilo - 2, ihi + 3):
+    for j in range(jlo - 2, jhi + 2):
+        for i in range(ilo - 2, ihi + 2):
 
             # u on x-edges
             u_xl[i + 1, j] = u[i, j] + 0.5 * \
@@ -385,8 +385,8 @@ def get_interface_states(qx, qy, ng, dx, dy, dt,
     # considered the normal to the interface portion of the predictor.
 
     # add the transverse flux differences to the preliminary interface states
-    for j in range(jlo - 1, jhi + 2):
-        for i in range(ilo - 1, ihi + 2):
+    for j in range(jlo - 1, jhi + 1):
+        for i in range(ilo - 1, ihi + 1):
 
             ubar = 0.5 * (uhat_adv[i, j] + uhat_adv[i + 1, j])
             vbar = 0.5 * (vhat_adv[i, j] + vhat_adv[i, j + 1])
@@ -440,8 +440,8 @@ def upwind(qx, qy, ng, q_l, q_r, s, q_int):
     jlo = ng
     jhi = ng + ny
 
-    for j in range(jlo - 1, jhi + 3):
-        for i in range(ilo - 1, ihi + 3):
+    for j in range(jlo - 1, jhi + 2):
+        for i in range(ilo - 1, ihi + 2):
 
             if (s[i, j] > 0.0):
                 q_int[i, j] = q_l[i, j]
@@ -467,8 +467,8 @@ def riemann(qx, qy, ng, q_l, q_r, s):
     jlo = ng
     jhi = ng + ny
 
-    for j in range(jlo - 1, jhi + 3):
-        for i in range(ilo - 1, ihi + 3):
+    for j in range(jlo - 1, jhi + 2):
+        for i in range(ilo - 1, ihi + 2):
 
             if (q_l[i, j] > 0.0 and q_l[i, j] + q_r[i, j] > 0.0):
                 s[i, j] = q_l[i, j]

setup_fortran.py

@@ -6,8 +6,10 @@
 
 from numpy.distutils.core import setup, Extension
 
-f_modules = ["compressible.interface", "advection_fv4.interface", "lm_atm.LM_atm_interface", "incompressible.incomp_interface", "swe.interface"]
+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]
+ext_modules = [Extension(module, [module.replace('.', '/') + '_f.f90'])
+               for module in f_modules]
 
 setup(ext_modules=ext_modules)