diff --git a/Source/Utils/check_interp_points_and_weights.py b/Source/Utils/check_interp_points_and_weights.py index 2f5a6e13b96..1e6a62e911e 100644 --- a/Source/Utils/check_interp_points_and_weights.py +++ b/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 @@ def refinement_points_and_weights(ii, sc, sf, cr): # ------------------------------------------------------------------------------- # 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 @@ def refinement_points_and_weights(ii, sc, sf, cr): 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 @@ def refinement_points_and_weights(ii, sc, sf, cr): " 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 + ) + ) diff --git a/Source/Utils/check_interp_points_and_weights_dev.py b/Source/Utils/check_interp_points_and_weights_dev.py new file mode 100644 index 00000000000..2f5a6e13b96 --- /dev/null +++ b/Source/Utils/check_interp_points_and_weights_dev.py @@ -0,0 +1,215 @@ +# Copyright 2020 Edoardo Zoni +# +# This file is part of WarpX +# +# License: BSD-3-Clause-LBNL + +# ------------------------------------------------------------------------------- +# Compute interpolation points and weights for coarsening and refinement in IO +# and MR applications in 1D (extensions to 2D and 3D are trivial). Weights are +# computed in order to guarantee total charge conservation for both cell-centered +# data (equal weights) and nodal data (weights depend on distance between points +# on fine and coarse grids). +# +# Notation: +# - index i refers to points on coarse grid +# - index ii refers to points on fine grid +# - sc denotes the staggering of the coarse data +# - sf denotes the staggering of the fine data +# - cr denotes the coarsening ratio (must be cr=1,2,4 only) +# +# For MR applications only the cases sc=sf=0 and sc=sf=1 are considered. Terms +# multiplied by (1-sf)*(1-sc) are ON for cell-centered data and OFF for nodal data, +# while terms multiplied by sf*sc are ON for nodal data and OFF for cell-centered +# data. C++ implementation in Source/ablastr/coarsen/average.(H/.cpp) and +# Source/ablastr/coarsen/sample.(H/.cpp) +# ------------------------------------------------------------------------------- + +import sys + +import numpy as np + + +# Fine grid limits (without ghost cells) +def fine_grid_limits(sf): + if sf == 0: # cell-centered + iimin = 0 + iimax = 7 + elif sf == 1: # nodal + iimin = 0 + iimax = 8 + return [iimin, iimax] + + +# 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] + + +# 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] + + +# Refinement for MR: interpolation points and weights +def refinement_points_and_weights(ii, sc, sf, cr): + if cr == 1: + numpts = 1 + idxmin = ii + elif cr >= 2: + if ii % cr == 0: + numpts = (1 - sf) * (1 - sc) + sf * sc + elif ii % cr != 0: + numpts = (1 - sf) * (1 - sc) + 2 * sf * sc + idxmin = (ii // cr) * (1 - sf) * (1 - sc) + (ii // cr) * sf * sc + weights = np.zeros(numpts) + for ir in range(numpts): + i = idxmin + ir + if ii == iimin or ii == iimax: + weights[ir] = (1 - sf) * (1 - sc) + ( + (abs(cr - abs(ii - i * cr))) / (cr) + (cr / 2 - 0.5) + ) * sf * sc + else: + weights[ir] = (1 - sf) * (1 - sc) + ( + (abs(cr - abs(ii - i * cr))) / (cr) + ) * sf * sc + return [numpts, idxmin, weights] + + +## TODO Coarsening for IO: interpolation points and weights +# def coarsening_points_and_weights_for_IO( i, sf, sc, cr ): +# if ( cr==1 ): +# numpts = 1+abs(sf-sc) +# idxmin = i-sc*(1-sf) +# elif ( cr>=2 ): +# numpts = 2-sf +# idxmin = i*cr+cr//2*(1-sc)-(1-sf) +# weights = np.zeros( numpts ) +# for ir in range( numpts ): +# weights[ir] = (1/numpts)*(1-sf)*(1-sc)+(1/numpts)*sf*sc +# return [ numpts, idxmin, weights ] + +# ------------------------------------------------------------------------------- +# Main +# ------------------------------------------------------------------------------- + +# 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)) + +# Loop over possible staggering of coarse and fine grid (cell-centered or nodal) +for sc in [0, 1]: + for sf in [0, 1]: + print("\n **************************************************") + print(" * Staggering of coarse grid: sc={}".format(sc), end="") + if sc == 0: + print(" cell-centered *") + elif sc == 1: + print(" nodal *") + print(" * Staggering of fine grid: sf={}".format(sf), end="") + if sf == 0: + print(" cell-centered *") + elif sf == 1: + print(" nodal *") + print(" **************************************************") + + iimin, iimax = fine_grid_limits(sf) + imin, imax = coarse_grid_limits(sc, sf, iimin, iimax) + + 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) + ) + + # 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 + 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 + ) + ) + for ir in range(numpts): # interpolation points and weights + ii = idxmin + ir + print(" ({},{})".format(ii, weights[ir]), end="") + if not (ir == numpts - 1): + print(" ", end="") + print() + + # Coarsening for MR: check conservation properties + for ii in range(nf): # 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 + 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)) + + print("\n Refinement for MR: check interpolation points and weights") + print(" ---------------------------------------------------------") + + # 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) + 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 + print(" ({},{})".format(i, weights[ir]), end="") + if not (ir == numpts - 1): + print(" ", end="") + print() + + # Refinement for MR: check conservation properties + for i in range(nc): # 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 + ws += weights[ir] + if ws != cr: + print("\n ERROR: sum of weights ws={} should be cr".format(ws)) diff --git a/Source/Utils/check_ipaw_generate_output.sh b/Source/Utils/check_ipaw_generate_output.sh new file mode 100644 index 00000000000..2137812fd55 --- /dev/null +++ b/Source/Utils/check_ipaw_generate_output.sh @@ -0,0 +1,13 @@ +#!/usr/bin/bash +cd "/home/basso9/src/WarpX/Source/Utils/" + +CR="2" + +OUT_DEV="cipaw_dev.out" +OUT_NEW="cipaw.out" + +echo "Coarsening ratio cr=$CR" >"$OUT_DEV" +echo "Coarsening ratio cr=$CR" >"$OUT_NEW" + +/bin/python3 "check_interp_points_and_weights_dev.py" <<< "$CR" >>"$OUT_DEV" 2>&1 +/bin/python3 "check_interp_points_and_weights.py" <<< "$CR" >>"$OUT_NEW" 2>&1 diff --git a/Source/Utils/cipaw.out b/Source/Utils/cipaw.out new file mode 100644 index 00000000000..f70f8af36c8 --- /dev/null +++ b/Source/Utils/cipaw.out @@ -0,0 +1,183 @@ +Coarsening ratio cr=2 + + Select coarsening ratio (cr>=1): cr= + ************************************************** + * Staggering of coarse grid: sc=0 cell-centered * + * Staggering of fine grid: sf=0 cell-centered * + ************************************************** + + Min and max index on coarse grid: imin=0 imax=3 + Min and max index on fine grid: iimin=0 iimax=7 + + Coarsening for MR: check interpolation points and weights + --------------------------------------------------------- + + Find value at i=0 by interpolating over the following points and weights: + (0,0.5) (1,0.5) + + Find value at i=1 by interpolating over the following points and weights: + (2,0.5) (3,0.5) + + Find value at i=2 by interpolating over the following points and weights: + (4,0.5) (5,0.5) + + Find value at i=3 by interpolating over the following points and weights: + (6,0.5) (7,0.5) + + Refinement for MR: check interpolation points and weights + --------------------------------------------------------- + + Min and max index on coarse grid: imin=0 imax=3 + Min and max index on fine grid: iimin=0 iimax=7 + + Find value at ii=0 by interpolating over the following points and weights: + (0,1.0) + + Find value at ii=1 by interpolating over the following points and weights: + (0,1.0) + + Find value at ii=2 by interpolating over the following points and weights: + (1,1.0) + + Find value at ii=3 by interpolating over the following points and weights: + (1,1.0) + + Find value at ii=4 by interpolating over the following points and weights: + (2,1.0) + + Find value at ii=5 by interpolating over the following points and weights: + (2,1.0) + + Find value at ii=6 by interpolating over the following points and weights: + (3,1.0) + + Find value at ii=7 by interpolating over the following points and weights: + (3,1.0) + + ************************************************** + * Staggering of coarse grid: sc=0 cell-centered * + * Staggering of fine grid: sf=1 nodal * + ************************************************** + + Min and max index on coarse grid: imin=-1 imax=4 + Min and max index on fine grid: iimin=0 iimax=8 + + Coarsening for MR: check interpolation points and weights + --------------------------------------------------------- + + Find value at i=-1 by interpolating over the following points and weights: + (-2,0.25) (-1,0.5) (0,0.25) + + Find value at i=0 by interpolating over the following points and weights: + (0,0.25) (1,0.5) (2,0.25) + + Find value at i=1 by interpolating over the following points and weights: + (2,0.25) (3,0.5) (4,0.25) + + Find value at i=2 by interpolating over the following points and weights: + (4,0.25) (5,0.5) (6,0.25) + + Find value at i=3 by interpolating over the following points and weights: + (6,0.25) (7,0.5) (8,0.25) + + Find value at i=4 by interpolating over the following points and weights: + (8,0.25) (9,0.5) (10,0.25) + + Refinement for MR: check interpolation points and weights + --------------------------------------------------------- + + WARNING: sc=0 not equal to sf=1, not implemented for Refinement for MR, continue ... + + ************************************************** + * Staggering of coarse grid: sc=1 nodal * + * Staggering of fine grid: sf=0 cell-centered * + ************************************************** + + Min and max index on coarse grid: imin=0 imax=4 + Min and max index on fine grid: iimin=0 iimax=7 + + Coarsening for MR: check interpolation points and weights + --------------------------------------------------------- + + Find value at i=0 by interpolating over the following points and weights: + (-1,0.5) (0,0.5) + + Find value at i=1 by interpolating over the following points and weights: + (1,0.5) (2,0.5) + + Find value at i=2 by interpolating over the following points and weights: + (3,0.5) (4,0.5) + + Find value at i=3 by interpolating over the following points and weights: + (5,0.5) (6,0.5) + + Find value at i=4 by interpolating over the following points and weights: + (7,0.5) (8,0.5) + + Refinement for MR: check interpolation points and weights + --------------------------------------------------------- + + WARNING: sc=1 not equal to sf=0, not implemented for Refinement for MR, continue ... + + ************************************************** + * Staggering of coarse grid: sc=1 nodal * + * Staggering of fine grid: sf=1 nodal * + ************************************************** + + Min and max index on coarse grid: imin=0 imax=4 + Min and max index on fine grid: iimin=0 iimax=8 + + Coarsening for MR: check interpolation points and weights + --------------------------------------------------------- + + Find value at i=0 by interpolating over the following points and weights: + (-1,0.25) (0,0.5) (1,0.25) + + Find value at i=1 by interpolating over the following points and weights: + (1,0.25) (2,0.5) (3,0.25) + + Find value at i=2 by interpolating over the following points and weights: + (3,0.25) (4,0.5) (5,0.25) + + Find value at i=3 by interpolating over the following points and weights: + (5,0.25) (6,0.5) (7,0.25) + + Find value at i=4 by interpolating over the following points and weights: + (7,0.25) (8,0.5) (9,0.25) + + Refinement for MR: check interpolation points and weights + --------------------------------------------------------- + + Min and max index on coarse grid: imin=0 imax=3 + Min and max index on fine grid: iimin=-1 iimax=7 + + Find value at ii=-1 by interpolating over the following points and weights: + (-1,1.0) (0,1.0) + + Find value at ii=0 by interpolating over the following points and weights: + (0,1.0) + + Find value at ii=1 by interpolating over the following points and weights: + (0,0.5) (1,0.5) + + Find value at ii=2 by interpolating over the following points and weights: + (1,1.0) + + Find value at ii=3 by interpolating over the following points and weights: + (1,0.5) (2,0.5) + + Find value at ii=4 by interpolating over the following points and weights: + (2,1.0) + + Find value at ii=5 by interpolating over the following points and weights: + (2,0.5) (3,0.5) + + Find value at ii=6 by interpolating over the following points and weights: + (3,1.0) + + Find value at ii=7 by interpolating over the following points and weights: + (3,1.0) (4,1.0) + + ERROR: sum of weights ws=2.5 should be cr=2 for i=0 + + ERROR: sum of weights ws=2.5 should be cr=2 for i=3 diff --git a/Source/Utils/cipaw_dev.out b/Source/Utils/cipaw_dev.out new file mode 100644 index 00000000000..68e98c60880 --- /dev/null +++ b/Source/Utils/cipaw_dev.out @@ -0,0 +1,140 @@ +Coarsening ratio cr=2 + + Select coarsening ratio (cr=1,2,4): cr= + ************************************************** + * Staggering of coarse grid: sc=0 cell-centered * + * Staggering of fine grid: sf=0 cell-centered * + ************************************************** + + Min and max index on coarse grid: imin=0 imax=3 + Min and max index on fine grid: iimin=0 iimax=7 + + Number of points on coarse grid: nc=4 + Number of points on fine grid: nf=8 + + Coarsening for MR: check interpolation points and weights + --------------------------------------------------------- + + Find value at i=0 by interpolating over the following points and weights: + (0,0.5) (1,0.5) + + Find value at i=1 by interpolating over the following points and weights: + (2,0.5) (3,0.5) + + Find value at i=2 by interpolating over the following points and weights: + (4,0.5) (5,0.5) + + Find value at i=3 by interpolating over the following points and weights: + (6,0.5) (7,0.5) + + Refinement for MR: check interpolation points and weights + --------------------------------------------------------- + + Find value at ii=0 by interpolating over the following points and weights: + (0,1.0) + + Find value at ii=1 by interpolating over the following points and weights: + (0,1.0) + + Find value at ii=2 by interpolating over the following points and weights: + (1,1.0) + + Find value at ii=3 by interpolating over the following points and weights: + (1,1.0) + + Find value at ii=4 by interpolating over the following points and weights: + (2,1.0) + + Find value at ii=5 by interpolating over the following points and weights: + (2,1.0) + + Find value at ii=6 by interpolating over the following points and weights: + (3,1.0) + + Find value at ii=7 by interpolating over the following points and weights: + (3,1.0) + + ************************************************** + * Staggering of coarse grid: sc=0 cell-centered * + * Staggering of fine grid: sf=1 nodal * + ************************************************** + + Min and max index on coarse grid: imin=0 imax=3 + Min and max index on fine grid: iimin=0 iimax=8 + + Number of points on coarse grid: nc=4 + Number of points on fine grid: nf=9 + + WARNING: sc=0 not equal to sf=1, not implemented for MR, continue ... + + ************************************************** + * Staggering of coarse grid: sc=1 nodal * + * Staggering of fine grid: sf=0 cell-centered * + ************************************************** + + Min and max index on coarse grid: imin=0 imax=4 + Min and max index on fine grid: iimin=0 iimax=7 + + Number of points on coarse grid: nc=5 + Number of points on fine grid: nf=8 + + WARNING: sc=1 not equal to sf=0, not implemented for MR, continue ... + + ************************************************** + * Staggering of coarse grid: sc=1 nodal * + * Staggering of fine grid: sf=1 nodal * + ************************************************** + + Min and max index on coarse grid: imin=0 imax=4 + Min and max index on fine grid: iimin=0 iimax=8 + + Number of points on coarse grid: nc=5 + Number of points on fine grid: nf=9 + + Coarsening for MR: check interpolation points and weights + --------------------------------------------------------- + + Find value at i=0 by interpolating over the following points and weights: + (-1,0.25) (0,0.5) (1,0.25) + + Find value at i=1 by interpolating over the following points and weights: + (1,0.25) (2,0.5) (3,0.25) + + Find value at i=2 by interpolating over the following points and weights: + (3,0.25) (4,0.5) (5,0.25) + + Find value at i=3 by interpolating over the following points and weights: + (5,0.25) (6,0.5) (7,0.25) + + Find value at i=4 by interpolating over the following points and weights: + (7,0.25) (8,0.5) (9,0.25) + + Refinement for MR: check interpolation points and weights + --------------------------------------------------------- + + Find value at ii=0 by interpolating over the following points and weights: + (0,1.5) + + Find value at ii=1 by interpolating over the following points and weights: + (0,0.5) (1,0.5) + + Find value at ii=2 by interpolating over the following points and weights: + (1,1.0) + + Find value at ii=3 by interpolating over the following points and weights: + (1,0.5) (2,0.5) + + Find value at ii=4 by interpolating over the following points and weights: + (2,1.0) + + Find value at ii=5 by interpolating over the following points and weights: + (2,0.5) (3,0.5) + + Find value at ii=6 by interpolating over the following points and weights: + (3,1.0) + + Find value at ii=7 by interpolating over the following points and weights: + (3,0.5) (4,0.5) + + Find value at ii=8 by interpolating over the following points and weights: + (4,1.5) diff --git a/Source/ablastr/coarsen/average.H b/Source/ablastr/coarsen/average.H index 621b4cb54e2..4f7779470af 100644 --- a/Source/ablastr/coarsen/average.H +++ b/Source/ablastr/coarsen/average.H @@ -25,87 +25,57 @@ */ namespace ablastr::coarsen::average { - /** - * \brief Interpolates the floating point data contained in the source Array4 - * \c arr_src, extracted from a fine MultiFab, with weights defined in - * such a way that the total charge is preserved. - * - * The input (sf) and output (sc) staggering need to be the same. - * - * \param[in] arr_src floating point data to be interpolated - * \param[in] sf staggering of the source fine MultiFab - * \param[in] sc staggering of the destination coarsened MultiFab - * \param[in] cr coarsening ratio along each spatial direction - * \param[in] i index along x of the coarsened Array4 to be filled - * \param[in] j index along y of the coarsened Array4 to be filled - * \param[in] k index along z of the coarsened Array4 to be filled - * \param[in] comp index along the fourth component of the Array4 \c arr_src - * containing the data to be interpolated - * - * \return interpolated field at cell (i,j,k) of a coarsened Array4 - */ - AMREX_GPU_DEVICE AMREX_FORCE_INLINE - amrex::Real - Interp ( - amrex::Array4 const &arr_src, - amrex::GpuArray const &sf, - amrex::GpuArray const &sc, - amrex::GpuArray const &cr, - int const i, - int const j, - int const k, - int const comp - ) + void + CalculateCoarseningData ( + amrex::GpuArray &idx_min_src, + amrex::GpuArray &np_src, + amrex::GpuArray &use_half_weight, // Use int instead bool? + amrex::Real &crx_cry_crz_inv, + amrex::GpuArray const &stag_src, + amrex::GpuArray const &stag_des, + amrex::GpuArray const &crse_ratio) { using namespace amrex::literals; - AMREX_ASSERT_WITH_MESSAGE(sf[0] == sc[0], "Interp: Staggering for component 0 does not match!"); - AMREX_ASSERT_WITH_MESSAGE(sf[1] == sc[1], "Interp: Staggering for component 1 does not match!"); - AMREX_ASSERT_WITH_MESSAGE(sf[2] == sc[2], "Interp: Staggering for component 2 does not match!"); - - // Indices of destination array (coarse) - int const ic[3] = {i, j, k}; - - // Number of points and starting indices of source array (fine) - int np[3], idx_min[3]; - - // Compute number of points for (int l = 0; l < 3; ++l) { - if (cr[l] == 1) { - np[l] = 1; // no coarsening + int stag_src_minus_cr_stag_des = stag_src[l] - crse_ratio[l]*stag_des[l]; + if (stag_src_minus_cr_stag_des & 1) { + // stag_src - crse_ratio*stag_des is odd + idx_min_src[l] = (stag_src_minus_cr_stag_des - 1)/2; + np_src[l] = crse_ratio[l] + 1; + use_half_weight[l] = true; // Use half weight at end points } else { - np[l] = cr[l] * (1 - sf[l]) * (1 - sc[l]) // cell-centered - + (2 * (cr[l] - 1) + 1) * sf[l] * sc[l]; // nodal + // stag_src - crse_ratio*stag_des is even + idx_min_src[l] = stag_src_minus_cr_stag_des/2; + np_src[l] = crse_ratio[l]; + use_half_weight[l] = false; // Don't use half weight at end points } } - // Compute starting indices of source array (fine) - for (int l = 0; l < 3; ++l) { - if (cr[l] == 1) { - idx_min[l] = ic[l]; // no coarsening - } else { - idx_min[l] = ic[l] * cr[l] * (1 - sf[l]) * (1 - sc[l]) // cell-centered - + (ic[l] * cr[l] - cr[l] + 1) * sf[l] * sc[l]; // nodal - } - } + crx_cry_crz_inv = 1.0_rt/static_cast(crse_ratio[0]*crse_ratio[1]*crse_ratio[2]); + } + + AMREX_GPU_DEVICE AMREX_FORCE_INLINE + amrex::Real + InterpWithCoarseningData ( + amrex::Array4 const &arr_src, + amrex::GpuArray const &idx_min_src, + amrex::GpuArray const &np_src, + amrex::GpuArray const &use_half_weight, + amrex::Real crx_cry_crz_inv, + amrex::GpuArray const &crse_ratio, + int i, int j, int k, int comp) + { + using namespace amrex::literals; // Auxiliary integer variables - int const numx = np[0]; - int const numy = np[1]; - int const numz = np[2]; - int const imin = idx_min[0]; - int const jmin = idx_min[1]; - int const kmin = idx_min[2]; - int const sfx = sf[0]; - int const sfy = sf[1]; - int const sfz = sf[2]; - int const scx = sc[0]; - int const scy = sc[1]; - int const scz = sc[2]; - int const crx = cr[0]; - int const cry = cr[1]; - int const crz = cr[2]; + int const krefmax = np_src[0] - 1; + int const jrefmax = np_src[1] - 1; + int const irefmax = np_src[2] - 1; + int const iimin = i * crse_ratio[0] + idx_min_src[0]; + int const jjmin = j * crse_ratio[1] + idx_min_src[1]; + int const kkmin = k * crse_ratio[2] + idx_min_src[2]; // Add neutral elements (=0) beyond guard cells in source array (fine) auto const arr_src_safe = [arr_src] @@ -114,36 +84,97 @@ namespace ablastr::coarsen::average }; // Interpolate over points computed above. Weights are computed in order - // to guarantee total charge conservation for both cell-centered data - // (equal weights) and nodal data (weights depend on distance between - // points on fine and coarse grids). Terms multiplied by (1-sf)*(1-sc) - // are ON for cell-centered data and OFF for nodal data, while terms - // multiplied by sf*sc are ON for nodal data and OFF for cell-centered data. + // to guarantee total charge conservation for any staggering. // Python script Source/Utils/check_interp_points_and_weights.py can be // used to check interpolation points and weights in 1D. amrex::Real c = 0.0_rt; - for (int kref = 0; kref < numz; ++kref) { - for (int jref = 0; jref < numy; ++jref) { - for (int iref = 0; iref < numx; ++iref) { - const int ii = imin + iref; - const int jj = jmin + jref; - const int kk = kmin + kref; - const amrex::Real wx = (1.0_rt / static_cast(numx)) * (1 - sfx) * (1 - scx) // if cell-centered - + ((amrex::Math::abs(crx - amrex::Math::abs(ii - i * crx))) / - static_cast(crx * crx)) * sfx * scx; // if nodal - const amrex::Real wy = (1.0_rt / static_cast(numy)) * (1 - sfy) * (1 - scy) // if cell-centered - + ((amrex::Math::abs(cry - amrex::Math::abs(jj - j * cry))) / - static_cast(cry * cry)) * sfy * scy; // if nodal - const amrex::Real wz = (1.0_rt / static_cast(numz)) * (1 - sfz) * (1 - scz) // if cell-centered - + ((amrex::Math::abs(crz - amrex::Math::abs(kk - k * crz))) / - static_cast(crz * crz)) * sfz * scz; // if nodal - c += wx * wy * wz * arr_src_safe(ii, jj, kk, comp); + for (int kref = 0; kref <= irefmax; ++kref) { + int kk = kkmin + kref; + bool use_half_kk = use_half_weight[2] && (kref == 0 || kref == irefmax); + + for (int jref = 0; jref <= jrefmax; ++jref) { + int jj = jjmin + jref; + bool use_half_jj = use_half_weight[1] && (jref == 0 || jref == jrefmax); + + for (int iref = 0; iref <= krefmax; ++iref) { + int ii = iimin + iref; + bool use_half_ii = use_half_weight[0] && (iref == 0 || iref == krefmax); + + amrex::Real c_iijjkk = arr_src_safe(ii, jj, kk, comp); + if (use_half_kk) { c_iijjkk *= 0.5_rt; } + if (use_half_jj) { c_iijjkk *= 0.5_rt; } + if (use_half_ii) { c_iijjkk *= 0.5_rt; } + c += c_iijjkk; } } } + + c *= crx_cry_crz_inv; + return c; } + /** + * \brief Interpolates the floating point data contained in the source Array4 + * \c arr_src, extracted from a fine MultiFab, with weights defined in + * such a way that the total charge is preserved. + * + * Works for any staggering of the fine source and coarse destination. + * + * \param[in] arr_src source fine floating point data to be interpolated + * \param[in] stag_src staggering of the source fine MultiFab + * \param[in] stag_des staggering of the destination coarsened MultiFab + * \param[in] crse_ratio coarsening ratio along each spatial direction + * \param[in] i index along x of the coarsened Array4 to be filled + * \param[in] j index along y of the coarsened Array4 to be filled + * \param[in] k index along z of the coarsened Array4 to be filled + * \param[in] comp index along the fourth component of the Array4 \c arr_src + * containing the data to be interpolated + * + * \return interpolated field at cell (i,j,k) of a coarsened Array4 + */ + AMREX_GPU_DEVICE AMREX_FORCE_INLINE + amrex::Real + Interp ( + amrex::Array4 const &arr_src, + amrex::GpuArray const &stag_src, + amrex::GpuArray const &stag_des, + amrex::GpuArray const &crse_ratio, + int const i, int const j, int const k, int const comp) + { + using namespace amrex::literals; + + // Number of points and starting indices of source array (fine) + amrex::GpuArray np_src, idx_min_src; + amrex::GpuArray use_half_weight; + amrex::Real crx_cry_crz_inv; + + // The following is intended to be an exact copy of CalculateCoarseningData + // This is needed to calculate the coarsening data on the GPU device instead of the host + + for (int l = 0; l < 3; ++l) { + int stag_src_minus_cr_stag_des = stag_src[l] - crse_ratio[l]*stag_des[l]; + if (stag_src_minus_cr_stag_des & 1) { + // stag_src - crse_ratio*stag_des is odd + idx_min_src[l] = (stag_src_minus_cr_stag_des - 1)/2; + np_src[l] = crse_ratio[l] + 1; + use_half_weight[l] = true; // Use half weight at end points + } else { + // stag_src - crse_ratio*stag_des is even + idx_min_src[l] = stag_src_minus_cr_stag_des/2; + np_src[l] = crse_ratio[l]; + use_half_weight[l] = false; // Don't use half weight at end points + } + } + + crx_cry_crz_inv = 1.0_rt/static_cast(crse_ratio[0]*crse_ratio[1]*crse_ratio[2]); + + // End copy of CalculateCoarseningData + + return InterpWithCoarseningData( + arr_src,idx_min_src,np_src,use_half_weight,crx_cry_crz_inv,crse_ratio,i,j,k,comp); + } + /** * \brief Loops over the boxes of the coarsened MultiFab \c mf_dst and fills * them by interpolating the data contained in the fine MultiFab \c mf_src. diff --git a/Source/ablastr/coarsen/average.cpp b/Source/ablastr/coarsen/average.cpp index a5dde3d8aa4..f24a17ed2a1 100644 --- a/Source/ablastr/coarsen/average.cpp +++ b/Source/ablastr/coarsen/average.cpp @@ -46,6 +46,13 @@ namespace ablastr::coarsen::average cr[i] = crse_ratio[i]; } + amrex::GpuArray idx_min_src, np_src; + amrex::GpuArray use_half_weight; + amrex::Real crx_cry_crz_inv; + + CalculateCoarseningData( + idx_min_src, np_src, use_half_weight, crx_cry_crz_inv, sf, sc, cr); + #ifdef AMREX_USE_OMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif @@ -55,11 +62,11 @@ namespace ablastr::coarsen::average amrex::Box const & bx = mfi.growntilebox(ngrow); amrex::Array4 const &arr_dst = mf_dst.array(mfi); amrex::Array4 const &arr_src = mf_src.const_array(mfi); - ParallelFor(bx, ncomp, - [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) { - arr_dst(i, j, k, n) = Interp( - arr_src, sf, sc, cr, i, j, k, n); - }); + amrex::ParallelFor(bx, ncomp, + [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) { + arr_dst(i, j, k, n) = InterpWithCoarseningData( + arr_src, idx_min_src, np_src, use_half_weight, crx_cry_crz_inv, cr, i, j, k, n); + }); } } diff --git a/Source/ablastr/coarsen/sample.H b/Source/ablastr/coarsen/sample.H index 80eac14e833..8f97d0a3b70 100644 --- a/Source/ablastr/coarsen/sample.H +++ b/Source/ablastr/coarsen/sample.H @@ -84,9 +84,8 @@ namespace ablastr::coarsen::sample const int imin = idx_min[0]; const int jmin = idx_min[1]; const int kmin = idx_min[2]; - amrex::Real const wx = 1.0_rt / static_cast(numx); - amrex::Real const wy = 1.0_rt / static_cast(numy); - amrex::Real const wz = 1.0_rt / static_cast(numz); + + amrex::Real const wxwywz = 1.0_rt / static_cast(numx*numy*numz); // Interpolate over points computed above amrex::Real c = 0.0_rt; @@ -96,10 +95,11 @@ namespace ablastr::coarsen::sample const int ii = imin+iref; const int jj = jmin+jref; const int kk = kmin+kref; - c += wx*wy*wz*arr_src(ii,jj,kk,comp); + c += arr_src(ii,jj,kk,comp); } } } + c *= wxwywz; return c; }