Generalize ablastr::coarsen::average to any staggering

Source/Utils/check_interp_points_and_weights.py

@@ -25,48 +25,108 @@
 # Source/ablastr/coarsen/sample.(H/.cpp)
 # -------------------------------------------------------------------------------
 
-import sys
 
 import numpy as np
 
+# Coarsening: map fine grid to coarse grid
+
 
 # Fine grid limits (without ghost cells)
-def fine_grid_limits(sf):
+# The find grid limits are arbitrarily chosen here to act as an example
+def coarsening_fine_grid_limits(sc, sf, cr):
     if sf == 0:  # cell-centered
-        iimin = 0
-        iimax = 7
+        ii_min = 0
+        ii_max = 4 * cr - 1
     elif sf == 1:  # nodal
-        iimin = 0
-        iimax = 8
-    return [iimin, iimax]
+        ii_min = 0
+        ii_max = 4 * cr
+    return [ii_min, ii_max]
 
 
 # Coarse grid limits (without ghost cells)
-def coarse_grid_limits(sc, sf, iimin, iimax):
-    imin = int(iimin / cr)
-    imax = int(iimax / cr) - (1 - sc) * sf + (1 - sf) * sc
-    return [imin, imax]
+def coarsening_coarse_grid_limits(sc, sf, cr, ii_min, ii_max):
+    return coarsening_coarse_grid_limits_brute_force(sc, sf, cr, ii_min, ii_max)
+
+
+def coarsening_coarse_grid_limits_brute_force(sc, sf, cr, ii_min, ii_max):
+    # Find coarse grid limits given fine grid limits using brute force scan of a
+    # by checking ii values produced by coarsening_points_and_weights for a large range of i values
+    i_range_start = (ii_min // cr) - 100
+    i_range_end = (ii_max // cr) + 100
+
+    i_min = i_range_end
+    i_max = i_range_start
+
+    for i in range(i_range_start, i_range_end + 1):
+        num_ii_pts, ii_start, weights = coarsening_points_and_weights(i, sc, sf, cr)
+        ii_end = ii_start + num_ii_pts - 1
+        if ii_min <= ii_end:
+            i_min = min(i_min, i)
+        if ii_max >= ii_start:
+            i_max = max(i_max, i)
+    return [i_min, i_max]
 
 
 # Coarsening for MR: interpolation points and weights
 def coarsening_points_and_weights(i, sc, sf, cr):
-    if cr == 1:
-        numpts = 1
-        idxmin = i
-    elif cr >= 2:
-        numpts = cr * (1 - sf) * (1 - sc) + (2 * (cr - 1) + 1) * sf * sc
-        idxmin = i * cr * (1 - sf) * (1 - sc) + (i * cr - cr + 1) * sf * sc
-    weights = np.zeros(numpts)
-    for ir in range(numpts):
-        ii = idxmin + ir
-        weights[ir] = (1 / cr) * (1 - sf) * (1 - sc) + (
-            (abs(cr - abs(ii - i * cr))) / (cr * cr)
-        ) * sf * sc
-    return [numpts, idxmin, weights]
+    two_ii_start = -cr * sc + sf - 1
+    if two_ii_start % 2 == 0:
+        ii_start = i * cr + two_ii_start // 2
+        num_ii_pts = cr + 1
+        weights = np.zeros(num_ii_pts)
+        weights[0] = 1.0 / (2 * cr)
+        weights[num_ii_pts - 1] = 1.0 / (2 * cr)
+        for ir in range(1, num_ii_pts - 1):
+            weights[ir] = 1.0 / cr
+    else:
+        ii_start = i * cr + (two_ii_start + 1) // 2
+        num_ii_pts = cr
+        weights = np.zeros(num_ii_pts)
+        for ir in range(0, num_ii_pts):
+            weights[ir] = 1.0 / cr
+
+    return [num_ii_pts, ii_start, weights]
+
+
+# Refinement: map coarse grid to fine grid
+
+
+def refinement_coarse_grid_limits(sc, sf, cr):
+    i_min = 0
+    i_max = 3
+    return [i_min, i_max]
+
+
+def refinement_fine_grid_limits(sc, sf, cr, i_min, i_max):
+    ii_range_start = i_min * cr - 100 * cr
+    ii_range_end = i_max * cr + 100 * cr
+
+    # print("ii_range_start={} and ii_range_end={}".format(ii_range_start,ii_range_end))
+
+    ii_min = ii_range_end
+    ii_max = ii_range_start
+
+    # print(" Before ii_min={} and ii_max={}".format(ii_min,ii_max))
+
+    for ii in range(ii_range_start, ii_range_end + 1):
+        num_i_pts, i_start, weights = refinement_points_and_weights(
+            ii, sc, sf, cr, ii_min, ii_max
+        )
+        i_end = i_start + num_i_pts - 1
+        if i_min <= i_end:
+            ii_min = min(ii_min, ii)
+        if i_max >= i_start:
+            ii_max = max(ii_max, ii)
+
+    # print(" After ii_min={} and ii_max={}".format(ii_min,ii_max))
+
+    return [ii_min, ii_max]
 
 
 # Refinement for MR: interpolation points and weights
-def refinement_points_and_weights(ii, sc, sf, cr):
+def refinement_points_and_weights(ii, sc, sf, cr, iimin, iimax):
+    idxmin = 0
+    numpts = 0
     if cr == 1:
         numpts = 1
         idxmin = ii
@@ -108,10 +168,7 @@
 # -------------------------------------------------------------------------------
 
 # Input coarsening ratio
-cr = int(input("\n Select coarsening ratio (cr=1,2,4): cr="))
-if cr != 1 and cr != 2 and cr != 4:
-    print()
-    sys.exit("coarsening ratio cr={} is not valid".format(cr))
+cr = int(input("\n Select coarsening ratio (cr>=1): cr="))
 
 # Loop over possible staggering of coarse and fine grid (cell-centered or nodal)
 for sc in [0, 1]:
@@ -129,8 +186,8 @@
             print(" nodal          *")
         print(" **************************************************")
 
-        iimin, iimax = fine_grid_limits(sf)
-        imin, imax = coarse_grid_limits(sc, sf, iimin, iimax)
+        iimin, iimax = coarsening_fine_grid_limits(sc, sf, cr)
+        imin, imax = coarsening_coarse_grid_limits(sc, sf, cr, iimin, iimax)
 
         print(
             "\n Min and max index on coarse grid:  imin={}  imax={}".format(imin, imax)
@@ -139,77 +196,103 @@
             " Min and max index on fine   grid: iimin={} iimax={}".format(iimin, iimax)
         )
 
-        # Number of grid points
-        nc = imax - imin + 1
-        nf = iimax - iimin + 1
-
-        print("\n Number of points on coarse grid: nc={}".format(nc))
-        print(" Number of points on fine   grid: nf={}".format(nf))
-
-        if sf != sc:
-            print(
-                "\n WARNING: sc={} not equal to sf={}, not implemented for MR, continue ...".format(
-                    sc, sf
-                )
-            )
-            continue
-
         print("\n Coarsening for MR: check interpolation points and weights")
         print(" ---------------------------------------------------------")
 
         # Coarsening for MR: interpolation points and weights
-        for i in range(nc):  # index on coarse grid
+        for i in range(imin, imax + 1):  # index on coarse grid
             numpts, idxmin, weights = coarsening_points_and_weights(i, sc, sf, cr)
             print(
                 "\n Find value at i={} by interpolating over the following points and weights:".format(
                     i
                 )
             )
+            wtotal = 0.0
             for ir in range(numpts):  # interpolation points and weights
                 ii = idxmin + ir
+                wtotal += weights[ir]
                 print(" ({},{})".format(ii, weights[ir]), end="")
                 if not (ir == numpts - 1):
                     print(" ", end="")
             print()
+            if abs(wtotal - 1.0) > 1e-9:
+                print(
+                    "\n ERROR: total weight wtotal={} should be 1 for coarse index i={}".format(
+                        wtotal, i
+                    )
+                )
 
         # Coarsening for MR: check conservation properties
-        for ii in range(nf):  # index on fine grid
+        for ii in range(iimin, iimax + 1):  # index on fine grid
             ws = 0.0
-            for i in range(nc):  # index on coarse grid
-                numpts, idxmin, weights = coarsening_points_and_weights(i, sc, sf, cr)
-                for ir in range(numpts):  # interpolation points and weights
-                    jj = idxmin + ir
+            for i in range(imin, imax + 1):  # index on coarse grid
+                num_ii_pts, ii_start, weights = coarsening_points_and_weights(
+                    i, sc, sf, cr
+                )
+                for ir in range(num_ii_pts):  # interpolation points and weights
+                    jj = ii_start + ir
                     if jj == ii:  # interpolation point matches point on fine grid
                         ws += weights[ir]
-            if ws != 1.0 / cr:
-                print("\n ERROR: sum of weights ws={} should be 1/cr".format(ws))
+            if abs(ws - 1.0 / cr) > 1e-9:
+                print(
+                    "\n ERROR: sum of weights ws={} should be 1/cr={} for ii={}".format(
+                        ws, 1.0 / cr, ii
+                    )
+                )
 
         print("\n Refinement for MR: check interpolation points and weights")
         print(" ---------------------------------------------------------")
 
+        if sf != sc:
+            print(
+                "\n WARNING: sc={} not equal to sf={}, not implemented for Refinement for MR, continue ...".format(
+                    sc, sf
+                )
+            )
+            continue
+
+        imin, imax = refinement_coarse_grid_limits(sc, sf, cr)
+        iimin, iimax = refinement_fine_grid_limits(sc, sf, cr, imin, imax)
+
+        # Number of grid points
+        print(
+            "\n Min and max index on coarse grid:  imin={}  imax={}".format(imin, imax)
+        )
+        print(
+            " Min and max index on fine   grid: iimin={} iimax={}".format(iimin, iimax)
+        )
+
         # Refinement for MR: interpolation points and weights
-        for ii in range(nf):  # index on fine grid
-            numpts, idxmin, weights = refinement_points_and_weights(ii, sc, sf, cr)
+        for ii in range(iimin, iimax + 1):  # index on fine grid
+            num_i_pts, i_start, weights = refinement_points_and_weights(
+                ii, sc, sf, cr, iimin, iimax
+            )
             print(
                 "\n Find value at ii={} by interpolating over the following points and weights:".format(
                     ii
                 )
             )
-            for ir in range(numpts):  # interpolation points and weights
-                i = idxmin + ir
+            for ir in range(num_i_pts):  # interpolation points and weights
+                i = i_start + ir
                 print(" ({},{})".format(i, weights[ir]), end="")
-                if not (ir == numpts - 1):
+                if not (ir == num_i_pts - 1):
                     print(" ", end="")
             print()
 
         # Refinement for MR: check conservation properties
-        for i in range(nc):  # index on coarse grid
+        for i in range(imin, imax + 1):  # index on coarse grid
             ws = 0.0
-            for ii in range(nf):  # index on fine grid
-                numpts, idxmin, weights = refinement_points_and_weights(ii, sc, sf, cr)
-                for ir in range(numpts):  # interpolation points and weights
-                    jj = idxmin + ir
-                    if jj == i:  # interpolation point matches point on coarse grid
+            for ii in range(iimin, iimax + 1):  # index on fine grid
+                num_i_pts, idxmin, weights = refinement_points_and_weights(
+                    ii, sc, sf, cr, iimin, iimax
+                )
+                for ir in range(num_i_pts):  # interpolation points and weights
+                    j = idxmin + ir
+                    if j == i:  # interpolation point matches point on coarse grid
                         ws += weights[ir]
-            if ws != cr:
-                print("\n ERROR: sum of weights ws={} should be cr".format(ws))
+            if abs(ws - cr) > 1e-9:
+                print(
+                    "\n ERROR: sum of weights ws={} should be cr={} for i={}".format(
+                        ws, cr, i
+                    )
+                )