From a78ded2b194e84dca951f16a110a62cbea033659 Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Sun, 26 Nov 2023 19:04:07 -0800 Subject: [PATCH 01/26] Modify python file --- .../Utils/check_interp_points_and_weights.py | 367 +++++++++++++----- Source/ablastr/coarsen/average.H | 22 +- Source/ablastr/coarsen/sample.H | 3 +- 3 files changed, 278 insertions(+), 114 deletions(-) diff --git a/Source/Utils/check_interp_points_and_weights.py b/Source/Utils/check_interp_points_and_weights.py index 8bf2cf08490..e63869fb468 100644 --- a/Source/Utils/check_interp_points_and_weights.py +++ b/Source/Utils/check_interp_points_and_weights.py @@ -29,56 +29,55 @@ 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 ] -# 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 ] -# 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 ] -# 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 ] +# # 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 ): @@ -93,15 +92,139 @@ def refinement_points_and_weights( ii, sc, sf, cr ): # weights[ir] = (1/numpts)*(1-sf)*(1-sc)+(1/numpts)*sf*sc # return [ numpts, idxmin, weights ] +def GetSignString(w): + if w == 0: + return " " + if w > 0: + return "+" + if w < 0: + return "-" + +def IntToString(i): + return GetSignString(i)+str(abs(i)) + +def FloatToString(f): + if (f == 0.0): + return " none " + return GetSignString(f)+str("{:.3f}".format(abs(f))) + +def coarsening_points_and_weights( i, sc, sf, cr ): + twoImin = -cr*sc + sf - 1 + if (twoImin % 2 == 0): + idxmin = i*cr + twoImin//2 + numpts = cr+1 + weights = np.zeros( numpts ) + weights[0] = 1.0/(2*cr) + weights[numpts-1] = 1.0/(2*cr) + for ir in range( 1, numpts-1 ): + weights[ir] = 1.0/cr + return [ numpts, idxmin, weights ] + else: + idxmin = i*cr + (twoImin+1)//2 + numpts = cr + weights = np.zeros( numpts ) + for ir in range( 0, numpts ): + weights[ir] = 1.0/cr + return [ numpts, idxmin, weights ] + # 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 ] + # weights = np.zeros( iimax - iimin + 1 ) + # for ii in range( iimin, iimax + 1 ): + # if (ii < idxmin or ii > idxmin + numpts ): + # weights[ii] = 0.00 + # else: + # weights[ii] = (1/cr)*(1-sf)*(1-sc)+((abs(cr-abs(ii-i*cr)))/(cr*cr))*sf*sc + # return [ numpts, idxmin, weights ] + +# def coarseing_points(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 +# return [idxmin, numpts] + +# def coarsening_weight(i, ii, sc, sf, cr): +# idxmin, numpts = coarseing_points(i, sc, sf, cr) +# if ii < idxmin or ii >= idxmin + numpts: +# return 0.0 +# return (1/cr)*(1-sf)*(1-sc)+((abs(cr-abs(ii-i*cr)))/(cr*cr))*sf*sc + +def coarsening_fine_grid_limits( sc, sf, cr ): + if ( sf == 0 ): # cell-centered + iimin = 1 + iimax = 4*cr + elif ( sf == 1 ): # nodal + iimin = 0 + iimax = 4*cr + return [ iimin, iimax ] + +def coarsening_coarse_grid_limits( sc, sf, cr, iimin, iimax ): + # imin = int( iimin/cr ) # original + # imax = int( iimax/cr )-(1-sc)*sf+(1-sf)*sc # original + imin = 10000 + imax = -10000 + for i in range(-100,101): + numpts, idxmin, weights = coarsening_points_and_weights(i, sc, sf, cr) + if iimin <= idxmin + numpts - 1: + imin = min(imin,i) + if iimax >= idxmin: + imax = max(imax,i) + return [ imin, imax ] + +# # 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 ] + +# def refinement_fine_grid_limits( sc, sf, cr ): +# if ( sf == 0 ): # cell-centered +# iimin = 0 # original +# iimax = 7 # original +# elif ( sf == 1 ): # nodal +# iimin = 0 # original +# iimax = 8 # original +# return [ iimin, iimax ] + +# def refinement_coarse_grid_limits( sc, sf, cr, iimin, iimax ): +# imin = int( iimin/cr ) # original +# imax = int( iimax/cr )-(1-sc)*sf+(1-sf)*sc # original +# return [ imin, imax ] + #------------------------------------------------------------------------------- # 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 ) ) +# 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]: @@ -110,81 +233,129 @@ def refinement_points_and_weights( ii, sc, sf, cr ): print( '\n **************************************************' ) print( ' * Staggering of coarse grid: sc={}'.format( sc ), end='' ) if ( sc == 0 ): - print( ' cell-centered *' ) + print( ' cell-centered *' ) elif ( sc == 1 ): - print( ' nodal *' ) + print( ' nodal *' ) print( ' * Staggering of fine grid: sf={}'.format( sf ), end='' ) if ( sf == 0 ): - print( ' cell-centered *' ) + print( ' cell-centered *' ) elif ( sf == 1 ): - print( ' nodal *' ) + print( ' nodal *' ) print( ' **************************************************' ) - iimin,iimax = fine_grid_limits( sf ) - imin ,imax = coarse_grid_limits( sc, sf, iimin, iimax ) + # if ( sf!=sc ): + # print( '\n WARNING: sc={} not equal to sf={}, not implemented for MR, continue ...'.format( sc, sf ) ) + # continue - 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 ) ) + print( '\n Coarsening for MR: check interpolation points and weights' ) + print( ' ---------------------------------------------------------' ) - # Number of grid points - nc = imax-imin+1 - nf = iimax-iimin+1 + iimin,iimax = coarsening_fine_grid_limits( sc, sf, cr ) + imin, imax = coarsening_coarse_grid_limits( sc, sf, cr, iimin, iimax ) - print( '\n Number of points on coarse grid: nc={}'.format( nc ) ) - print( ' Number of points on fine grid: nf={}'.format( nf ) ) + 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 ) ) + # print( ' Min and max index on fine grid: iiminL={} iimaxL={}'.format( iiminLarge, iimaxLarge ) ) - 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( ' ---------------------------------------------------------' ) + # 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 ) ) + + # coarsening_table = range(imin,imax+1) + iiminLarge = iimin + iimaxLarge = iimax # 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 ) + # coarsening_table[i-imin] = {i,idxmin,numpts,weights} 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 - print( ' ({},{})'.format( ii, weights[ir] ), end='' ) + wtotal += weights[ir] + print( ' (ii={},w={:.3f})'.format( ii, weights[ir] ), end='' ) if not ( ir == numpts-1 ): - print( ' ', end='' ) + 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)) + + # print("iiminLarge = {}, iimaxlarge = {}".format(iiminLarge,iimaxLarge)) + + # for ii in range(iiminLarge, iimaxLarge+1): + # foundimin = False + # foundimax = False + # for i in range(iminLarge, imaxLarge+1): + # numpts,idxmin,weights = coarsening_points_and_weights( i, sc, sf, cr ) + # if (ii < idxmin): + # foundimin = True + # iminTest = min(iminTest, i) + # # print("ii={}, iminTest<={}".format(ii,i)) + # if (ii >= idxmin + numpts): + # foundimax = True + # imaxTest = max(imaxTest, i) + # # print("ii={}, imaxTest>={}".format(ii,i)) + # if (foundimin and foundimax): + # iiminTest = min(iiminTest,ii) + # iimaxTest = max(iimaxTest,ii) + + # print(" iminTest = {}, imaxTest = {}".format( iminTest, imaxTest)) + # print("iiminTest = {}, iimaxTest = {}".format(iiminTest,iimaxTest)) # 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 + for i in range(imin, imax+1): # 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 ) ) + 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( ' ---------------------------------------------------------' ) + # 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() + # iimin,iimax = refinement_fine_grid_limits( sc, sf, cr ) + # imin ,imax = refinement_coarse_grid_limits( sc, sf, cr, iimin, iimax ) + # # Number of grid points + # nc = imax-imin+1 + # nf = iimax-iimin+1 - # 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 ) ) + # 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 ) ) + # print( '\n Number of points on coarse grid: nc={}'.format( nc ) ) + # print( ' Number of points on fine grid: nf={}'.format( nf ) ) + + # # 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 + # sign = '+' + # if (i == 0): + # sign = ' ' + # if (i < 0): + # sign = '-' + # print( ' (i='+sign+'{},w={:.2f})'.format( abs(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 ( abs(ws - cr) > 1E-9 ): + # print( '\n ERROR: sum of weights ws={:.2f} should be cr={} for i={}'.format( ws, cr, i ) ) diff --git a/Source/ablastr/coarsen/average.H b/Source/ablastr/coarsen/average.H index 315938acb48..c75dcad09c0 100644 --- a/Source/ablastr/coarsen/average.H +++ b/Source/ablastr/coarsen/average.H @@ -52,10 +52,7 @@ namespace ablastr::coarsen::average amrex::GpuArray const &sf, amrex::GpuArray const &sc, amrex::GpuArray const &cr, - int const i, - int const j, - int const k, - int const comp + int const i, int const j, int const k, int const comp ) { using namespace amrex::literals; @@ -70,22 +67,16 @@ namespace ablastr::coarsen::average // Number of points and starting indices of source array (fine) int np[3], idx_min[3]; - // Compute number of points + // Compute number of points (np) and starting indices of source array (fine) for (int l = 0; l < 3; ++l) { if (cr[l] == 1) np[l] = 1; // no coarsening - else - np[l] = cr[l] * (1 - sf[l]) * (1 - sc[l]) // cell-centered - + (2 * (cr[l] - 1) + 1) * sf[l] * sc[l]; // nodal - } - - // 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 + np[l] = cr[l] * (1 - sf[l]) * (1 - sc[l]) // cell-centered + + (2 * cr[l] - 1) * sf[l] * sc[l]; // nodal + 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 } // Auxiliary integer variables @@ -139,6 +130,7 @@ namespace ablastr::coarsen::average } } } + return c; } diff --git a/Source/ablastr/coarsen/sample.H b/Source/ablastr/coarsen/sample.H index 6ef0962168a..ac0f9961198 100644 --- a/Source/ablastr/coarsen/sample.H +++ b/Source/ablastr/coarsen/sample.H @@ -96,10 +96,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*=wx*wy*wz; return c; } From a2eb42c9fa1e18345e9b00861cf7b1fb01f90bcf Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Sun, 26 Nov 2023 21:11:46 -0800 Subject: [PATCH 02/26] New methods --- .../Utils/check_interp_points_and_weights.py | 212 +++++------------- Source/ablastr/coarsen/average.H | 162 ++++++++----- Source/ablastr/coarsen/average.cpp | 50 +++++ 3 files changed, 213 insertions(+), 211 deletions(-) diff --git a/Source/Utils/check_interp_points_and_weights.py b/Source/Utils/check_interp_points_and_weights.py index e63869fb468..97c49811c96 100644 --- a/Source/Utils/check_interp_points_and_weights.py +++ b/Source/Utils/check_interp_points_and_weights.py @@ -29,36 +29,6 @@ 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 ): @@ -79,87 +49,6 @@ # 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 ] - -def GetSignString(w): - if w == 0: - return " " - if w > 0: - return "+" - if w < 0: - return "-" - -def IntToString(i): - return GetSignString(i)+str(abs(i)) - -def FloatToString(f): - if (f == 0.0): - return " none " - return GetSignString(f)+str("{:.3f}".format(abs(f))) - -def coarsening_points_and_weights( i, sc, sf, cr ): - twoImin = -cr*sc + sf - 1 - if (twoImin % 2 == 0): - idxmin = i*cr + twoImin//2 - numpts = cr+1 - weights = np.zeros( numpts ) - weights[0] = 1.0/(2*cr) - weights[numpts-1] = 1.0/(2*cr) - for ir in range( 1, numpts-1 ): - weights[ir] = 1.0/cr - return [ numpts, idxmin, weights ] - else: - idxmin = i*cr + (twoImin+1)//2 - numpts = cr - weights = np.zeros( numpts ) - for ir in range( 0, numpts ): - weights[ir] = 1.0/cr - return [ numpts, idxmin, weights ] - # 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 ] - # weights = np.zeros( iimax - iimin + 1 ) - # for ii in range( iimin, iimax + 1 ): - # if (ii < idxmin or ii > idxmin + numpts ): - # weights[ii] = 0.00 - # else: - # weights[ii] = (1/cr)*(1-sf)*(1-sc)+((abs(cr-abs(ii-i*cr)))/(cr*cr))*sf*sc - # return [ numpts, idxmin, weights ] - -# def coarseing_points(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 -# return [idxmin, numpts] - -# def coarsening_weight(i, ii, sc, sf, cr): -# idxmin, numpts = coarseing_points(i, sc, sf, cr) -# if ii < idxmin or ii >= idxmin + numpts: -# return 0.0 -# return (1/cr)*(1-sf)*(1-sc)+((abs(cr-abs(ii-i*cr)))/(cr*cr))*sf*sc - def coarsening_fine_grid_limits( sc, sf, cr ): if ( sf == 0 ): # cell-centered iimin = 1 @@ -182,6 +71,47 @@ def coarsening_coarse_grid_limits( sc, sf, cr, iimin, iimax ): imax = max(imax,i) return [ imin, imax ] +def coarsening_points_and_weights( i, sc, sf, cr ): + twoImin = -cr*sc + sf - 1 + if (twoImin % 2 == 0): + idxmin = i*cr + twoImin//2 + numpts = cr+1 + weights = np.zeros( numpts ) + weights[0] = 1.0/(2*cr) + weights[numpts-1] = 1.0/(2*cr) + for ir in range( 1, numpts-1 ): + weights[ir] = 1.0/cr + else: + idxmin = i*cr + (twoImin+1)//2 + numpts = cr + weights = np.zeros( numpts ) + for ir in range( 0, numpts ): + weights[ir] = 1.0/cr + + return [ numpts, idxmin, weights ] + +# def refinement_fine_grid_limits( sc, sf, cr ): +# if ( sf == 0 ): # cell-centered +# iimin = 0 # original +# iimax = 7 # original +# elif ( sf == 1 ): # nodal +# iimin = 0 # original +# iimax = 8 # original +# return [ iimin, iimax ] + +# def refinement_coarse_grid_limits( sc, sf, cr, iimin, iimax ): +# # imin = int( iimin/cr ) # original +# # imax = int( iimax/cr )-(1-sc)*sf+(1-sf)*sc # original +# imin = 10000 +# imax = -10000 +# for i in range(-100,101): +# numpts, idxmin, weights = refinement_points_and_weights(i, sc, sf, cr) +# if iimin <= idxmin + numpts - 1: +# imin = min(imin,i) +# if iimax >= idxmin: +# imax = max(imax,i) +# return [ imin, imax ] + # # Refinement for MR: interpolation points and weights # def refinement_points_and_weights( ii, sc, sf, cr ): # if ( cr==1 ): @@ -202,19 +132,18 @@ def coarsening_coarse_grid_limits( sc, sf, cr, iimin, iimax ): # weights[ir] = (1-sf)*(1-sc)+((abs(cr-abs(ii-i*cr)))/(cr))*sf*sc # return [ numpts, idxmin, weights ] -# def refinement_fine_grid_limits( sc, sf, cr ): -# if ( sf == 0 ): # cell-centered -# iimin = 0 # original -# iimax = 7 # original -# elif ( sf == 1 ): # nodal -# iimin = 0 # original -# iimax = 8 # original -# return [ iimin, iimax ] - -# def refinement_coarse_grid_limits( sc, sf, cr, iimin, iimax ): -# imin = int( iimin/cr ) # original -# imax = int( iimax/cr )-(1-sc)*sf+(1-sf)*sc # original -# return [ imin, imax ] +## 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 @@ -255,8 +184,6 @@ def coarsening_coarse_grid_limits( sc, sf, cr, 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 ) ) - # print( ' Min and max index on fine grid: iiminL={} iimaxL={}'.format( iiminLarge, iimaxLarge ) ) - # Number of grid points # nc = imax-imin+1 @@ -264,14 +191,9 @@ def coarsening_coarse_grid_limits( sc, sf, cr, iimin, iimax ): # print( '\n Number of points on coarse grid: nc={}'.format( nc ) ) # print( ' Number of points on fine grid: nf={}'.format( nf ) ) - # coarsening_table = range(imin,imax+1) - - iiminLarge = iimin - iimaxLarge = iimax # Coarsening for MR: interpolation points and weights - for i in range(imin, imax+1): # 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 ) - # coarsening_table[i-imin] = {i,idxmin,numpts,weights} 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 @@ -284,32 +206,10 @@ def coarsening_coarse_grid_limits( sc, sf, cr, iimin, iimax ): if (abs(wtotal - 1.0) > 1E-9): print ('\n ERROR: total weight wtotal={} should be 1 for coarse index i={}'.format(wtotal,i)) - # print("iiminLarge = {}, iimaxlarge = {}".format(iiminLarge,iimaxLarge)) - - # for ii in range(iiminLarge, iimaxLarge+1): - # foundimin = False - # foundimax = False - # for i in range(iminLarge, imaxLarge+1): - # numpts,idxmin,weights = coarsening_points_and_weights( i, sc, sf, cr ) - # if (ii < idxmin): - # foundimin = True - # iminTest = min(iminTest, i) - # # print("ii={}, iminTest<={}".format(ii,i)) - # if (ii >= idxmin + numpts): - # foundimax = True - # imaxTest = max(imaxTest, i) - # # print("ii={}, imaxTest>={}".format(ii,i)) - # if (foundimin and foundimax): - # iiminTest = min(iiminTest,ii) - # iimaxTest = max(iimaxTest,ii) - - # print(" iminTest = {}, imaxTest = {}".format( iminTest, imaxTest)) - # print("iiminTest = {}, iimaxTest = {}".format(iiminTest,iimaxTest)) - # Coarsening for MR: check conservation properties - for ii in range(iimin, iimax+1): # index on fine grid + for ii in range( iimin, iimax+1 ): # index on fine grid ws = 0.0 - for i in range(imin, imax+1): # 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 ) for ir in range( numpts ): # interpolation points and weights jj = idxmin+ir diff --git a/Source/ablastr/coarsen/average.H b/Source/ablastr/coarsen/average.H index c75dcad09c0..da8c5fb6517 100644 --- a/Source/ablastr/coarsen/average.H +++ b/Source/ablastr/coarsen/average.H @@ -44,63 +44,47 @@ namespace ablastr::coarsen::average * * \return interpolated field at cell (i,j,k) of a coarsened Array4 */ - AMREX_GPU_DEVICE - AMREX_FORCE_INLINE + 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 - ) + int const i, int const j, int const k, int const comp) { 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]; + bool useHalf[3]; - // Compute number of points (np) and starting indices of source array (fine) for (int l = 0; l < 3; ++l) { - if (cr[l] == 1) - np[l] = 1; // no coarsening - idx_min[l] = ic[l]; // no coarsening - else - np[l] = cr[l] * (1 - sf[l]) * (1 - sc[l]) // cell-centered - + (2 * cr[l] - 1) * sf[l] * sc[l]; // nodal - 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 + const int twoImin = -cr[l]*sc[l] + sf[l] - 1; + if (twoImin % 2 == 0) { + idxmin[l] = twoImin/2; + np[l] = cr[l]+1; + useHalf[l] = true; + } else { + idxmin[l] = (twoImin+1)/2; + np[l] = cr[l]; + useHalf[l] = false; + } } // 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 iimin = i * cr[0] + idx_min[0]; + int const jjmin = j * cr[1] + idx_min[1]; + int const kkmin = k * cr[2] + idx_min[2]; + // int const iimax = iimin + np[0] - 1; + // int const jjmax = jjmin + np[1] - 1; + // int const kkmax = kkmin + np[2] - 1; // Add neutral elements (=0) beyond guard cells in source array (fine) - auto const arr_src_safe = [arr_src] - AMREX_GPU_DEVICE(int const ix, int const iy, int const iz, int const n) noexcept { - return arr_src.contains(ix, iy, iz) ? arr_src(ix, iy, iz, n) : 0.0_rt; - }; + // auto const arr_src_safe = [arr_src] + // AMREX_GPU_DEVICE(int const ix, int const iy, int const iz, int const n) noexcept { + // return arr_src.contains(ix, iy, iz) ? arr_src(ix, iy, iz, n) : 0.0_rt; + // }; // Interpolate over points computed above. Weights are computed in order // to guarantee total charge conservation for both cell-centered data @@ -111,26 +95,73 @@ namespace ablastr::coarsen::average // 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 kr = 0; kr < np[2]; ++kr) { + int kk = kkmin + kr; + kk = amrex::Clamp(kk, arr_src.begin.z, arr_src.end.z - 1); + amrex::Real kfactor = (useHalf[2] && (kr == 0 || kr == np[2] - 1)) ? 0.5_rt : 1.0_rt; + + for (int jr = 0; jr < np[1]; ++jr) { + int jj = jjmin + jr; + jj = amrex::Clamp(jj, arr_src.begin.y, arr_src.end.y - 1); + amrex::Real jfactor = (useHalf[1] && (jr == 0 || jr == np[1] - 1)) ? 0.5_rt : 1.0_rt; + + for (int ir = 0; ir < np[1]; ++ir) { + int ii = iimin + ir; + ii = amrex::Clamp(ii, arr_src.begin.x, arr_src.end.x - 1); + amrex::Real ifactor = (useHalf[1] && (ir == 0 || ir == np[1] - 1)) ? 0.5_rt : 1.0_rt; + + c += ifactor*jfactor*kfactor*arr_src(ii,jj,kk,comp); + + // amrex::Real src_iijjkk = arr_src_safe(ii, jj, kk, comp); + // amrex::Real src_iijjkk = Array4ClampedEval(arr_src, ii, jj, kk, comp); + // c += ifactor*jfactor*kfactor*src_iijjkk; } } } + c *= 1.0_rt / static_cast(cr[0]*cr[1]*cr[2]); + + return c; + } + + AMREX_GPU_DEVICE AMREX_FORCE_INLINE + amrex::Real + Interp ( + amrex::Array4 const &arr_src, + amrex::GpuArray const &coarsen_ratio, + amrex::GpuArray const &ind_min, + amrex::Vector const &weights_x, + amrex::Vector const &weights_y, + amrex::Vector const &weights_z, + int const i, int const j, int const k, int const comp + ) + { + amrex::Real c = 0.0_rt; + + const int iimin = i * coarsen_ratio[0] + ind_min[0]; + const int jjmin = j * coarsen_ratio[1] + ind_min[1]; + const int kkmin = k * coarsen_ratio[2] + ind_min[2]; + + for (int kr = 0; k < weights_z.size(); ++kr) + { + int kk = kkmin + kr; + kk = amrex::Clamp(kk, arr_src.begin.z, arr_src.end.z - 1); + + for (int jr = 0; j < weights_y.size(); ++jr) + { + int jj = jjmin + jr; + jj = amrex::Clamp(jj, arr_src.begin.y, arr_src.end.y - 1); + + for (int ir = 0; i < weights_x.size(); ++ir) + { + int ii = iimin + ir; + ii = amrex::Clamp(ii, arr_src.begin.x, arr_src.end.x - 1); + + amrex::Real src_iijjkk = arr_src(ii, jj, kk, comp); + c += weights_x[ir]*weights_y[ir]*weights_z[ir]*src_iijjkk; + } + } + } return c; } @@ -173,6 +204,27 @@ namespace ablastr::coarsen::average amrex::IntVect crse_ratio ); + void + CoarseningPointsAndWeights ( + amrex::Vector &weights_x, + amrex::Vector &weights_y, + amrex::Vector &weights_z, + amrex::GpuArray const &ind_min, + amrex::GpuArray const &coarsen_ratio, + amrex::GpuArray const &stag_fine_src, + amrex::GpuArray const & stag_crse_des); + + template ::value,int>::type = 0> + [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + U& Array4ClampedEval (amrex::Array4 const& arr, + int i, int j, int k, int n) + { + int iClamped = amrex::Clamp(i, arr.begin.x, arr.end.x - 1); + int jClamped = amrex::Clamp(j, arr.begin.y, arr.end.y - 1); + int kClamped = amrex::Clamp(k, arr.begin.x, arr.end.z - 1); + return arr(iClamped,jClamped,kClamped,n); + } + } // namespace ablastr::coarsen::average #endif // ABLASTR_COARSEN_AVERAGE_H_ diff --git a/Source/ablastr/coarsen/average.cpp b/Source/ablastr/coarsen/average.cpp index a5dde3d8aa4..0fb413ecff1 100644 --- a/Source/ablastr/coarsen/average.cpp +++ b/Source/ablastr/coarsen/average.cpp @@ -82,4 +82,54 @@ namespace ablastr::coarsen::average Loop(mf_dst, mf_src, ncomp, ngrow, crse_ratio); } + void + CoarseningPointsAndWeights ( + amrex::Vector &weights_x, + amrex::Vector &weights_y, + amrex::Vector &weights_z, + amrex::GpuArray &src_index_min, + amrex::GpuArray const &coarsen_ratio, + amrex::GpuArray const &stag_fine_src, + amrex::GpuArray const &stag_crse_des + ) + { + using namespace amrex::literals; + + int num_points[3]; + bool useHalf[3]; + + for (int l = 0; l < 3; ++l) { + int const twoImin = -coarsen_ratio[l]*stag_crse_des[l] + stag_fine_src[l] - 1; + if (twoImin % 2 == 0) { + src_index_min[l] = twoImin/2; + num_points[l] = coarsen_ratio[l]+1; + useHalf[l] = true; + } else { + src_index_min[l] = (twoImin+1)/2; + num_points[l] = coarsen_ratio[l]; + useHalf[l] = false; + } + } + + // amrex::Real const coarsen_factor = 1.0_rt / static_cast(coarsen_ratio[0]*coarsen_ratio[1]*coarsen_ratio[2]); + + weights_x.resize(num_points[0], 1.0_rt); + weights_y.resize(num_points[1], 1.0_rt); + weights_z.resize(num_points[2], 1.0_rt); + + if (useHalf[0]) { + weights_x[0] *= 0.5_rt; + weights_x[num_points[0]-1] *= 0.5_rt; + } + if (useHalf[1]) { + weights_y[0] *= 0.5_rt; + weights_y[num_points[1]-1] *= 0.5_rt; + } + if (useHalf[2]) { + weights_z[0] *= 0.5_rt; + weights_z[num_points[2]-1] *= 0.5_rt; + } + + } + } // namespace ablastr::coarsen::average From 0514c26ee945864ba811ecf2f4dffaf320da7fcf Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Mon, 27 Nov 2023 11:22:14 -0800 Subject: [PATCH 03/26] Refactor with more methods --- Source/ablastr/coarsen/average.H | 274 ++++++++++++++++++----------- Source/ablastr/coarsen/average.cpp | 50 ------ Source/ablastr/coarsen/sample.H | 7 +- 3 files changed, 172 insertions(+), 159 deletions(-) diff --git a/Source/ablastr/coarsen/average.H b/Source/ablastr/coarsen/average.H index da8c5fb6517..54a37b75dd3 100644 --- a/Source/ablastr/coarsen/average.H +++ b/Source/ablastr/coarsen/average.H @@ -30,138 +30,118 @@ namespace ablastr::coarsen::average * \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. + * Works for any staggering of the fine source and coarse d * - * \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 + * \param[in] arr_src source fine floating point data to be interpolated + * \param[in] stag_src_fine staggering of the source fine MultiFab + * \param[in] stag_des_crse 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 + * \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, + amrex::GpuArray const &stag_src_fine, + amrex::GpuArray const &stag_des_crse, + 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) - int np[3], idx_min[3]; - bool useHalf[3]; - - for (int l = 0; l < 3; ++l) { - const int twoImin = -cr[l]*sc[l] + sf[l] - 1; - if (twoImin % 2 == 0) { - idxmin[l] = twoImin/2; - np[l] = cr[l]+1; - useHalf[l] = true; - } else { - idxmin[l] = (twoImin+1)/2; - np[l] = cr[l]; - useHalf[l] = false; - } - } + amrex::GpuArray np, idx_min, use_half; + amrex::Real crx_cry_crz_inv; - // Auxiliary integer variables - int const iimin = i * cr[0] + idx_min[0]; - int const jjmin = j * cr[1] + idx_min[1]; - int const kkmin = k * cr[2] + idx_min[2]; - // int const iimax = iimin + np[0] - 1; - // int const jjmax = jjmin + np[1] - 1; - // int const kkmax = kkmin + np[2] - 1; + CalculateCoarseningData(idx_min, np, use_half, crx_cry_crz_inv, stag_des_crse, stag_src_fine, crse_ratio); - // Add neutral elements (=0) beyond guard cells in source array (fine) - // auto const arr_src_safe = [arr_src] - // AMREX_GPU_DEVICE(int const ix, int const iy, int const iz, int const n) noexcept { - // return arr_src.contains(ix, iy, iz) ? arr_src(ix, iy, iz, n) : 0.0_rt; - // }; + return InterpWithCoarseningData(arr_src, idx_min, np, use_half, crx_cry_crz_inv, crse_ratio, i, j, k, comp); + } - // 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. - // 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 kr = 0; kr < np[2]; ++kr) { - int kk = kkmin + kr; - kk = amrex::Clamp(kk, arr_src.begin.z, arr_src.end.z - 1); - amrex::Real kfactor = (useHalf[2] && (kr == 0 || kr == np[2] - 1)) ? 0.5_rt : 1.0_rt; - - for (int jr = 0; jr < np[1]; ++jr) { - int jj = jjmin + jr; - jj = amrex::Clamp(jj, arr_src.begin.y, arr_src.end.y - 1); - amrex::Real jfactor = (useHalf[1] && (jr == 0 || jr == np[1] - 1)) ? 0.5_rt : 1.0_rt; - - for (int ir = 0; ir < np[1]; ++ir) { - int ii = iimin + ir; - ii = amrex::Clamp(ii, arr_src.begin.x, arr_src.end.x - 1); - amrex::Real ifactor = (useHalf[1] && (ir == 0 || ir == np[1] - 1)) ? 0.5_rt : 1.0_rt; - - c += ifactor*jfactor*kfactor*arr_src(ii,jj,kk,comp); - - // amrex::Real src_iijjkk = arr_src_safe(ii, jj, kk, comp); - // amrex::Real src_iijjkk = Array4ClampedEval(arr_src, ii, jj, kk, comp); - // c += ifactor*jfactor*kfactor*src_iijjkk; - } + AMREX_GPU_DEVICE AMREX_FORCE_INLINE + void + CalculateCoarseningData ( + amrex::GpuArray &src_fine_index_min, + amrex::GpuArray &num_src_points, + amrex::GpuArray &use_half_weight, // Use int instead bool? + amrex::Real &crx_cry_crz_inv, + amrex::GpuArray const &stag_src_fine, + amrex::GpuArray const &stag_des_crse, + amrex::GpuArray const &crse_ratio) + { + for (int l = 0; l < 3; ++l) { + int two_times_index_min = -crse_ratio[l]*stag_des_crse[l] + stag_src_fine[l] - 1; + if (two_times_index_min % 2 == 0) { + src_fine_index_min[l] = two_times_index_min/2; + num_src_points[l] = crse_ratio[l]+1; + use_half_weight[l] = 1; // True + } else { + src_fine_index_min[l] = (two_times_index_min+1)/2; + num_src_points[l] = crse_ratio[l]; + use_half_weight[l] = 0; // False } } - - c *= 1.0_rt / static_cast(cr[0]*cr[1]*cr[2]); - - return c; + 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 - Interp ( + InterpWithCoarseningData ( amrex::Array4 const &arr_src, + amrex::GpuArray const &src_fine_index_min, + amrex::GpuArray const &num_src_points, + amrex::GpuArray const &use_half, + amrex::Real crx_cry_crz_inv, amrex::GpuArray const &coarsen_ratio, - amrex::GpuArray const &ind_min, - amrex::Vector const &weights_x, - amrex::Vector const &weights_y, - amrex::Vector const &weights_z, - int const i, int const j, int const k, int const comp - ) + int i, int j, int k, int comp) { - amrex::Real c = 0.0_rt; + using namespace amrex::literals; - const int iimin = i * coarsen_ratio[0] + ind_min[0]; - const int jjmin = j * coarsen_ratio[1] + ind_min[1]; - const int kkmin = k * coarsen_ratio[2] + ind_min[2]; + // Auxiliary integer variables + int const numx = num_src_points[0]; + int const numy = num_src_points[1]; + int const numz = num_src_points[2]; + int const iimin = i * coarsen_ratio[0] + src_fine_index_min[0]; + int const jjmin = j * coarsen_ratio[1] + src_fine_index_min[1]; + int const kkmin = k * coarsen_ratio[2] + src_fine_index_min[2]; - for (int kr = 0; k < weights_z.size(); ++kr) - { - int kk = kkmin + kr; - kk = amrex::Clamp(kk, arr_src.begin.z, arr_src.end.z - 1); + // Add neutral elements (=0) beyond guard cells in source array (fine) + auto const arr_src_safe = [arr_src] + AMREX_GPU_DEVICE(int const ix, int const iy, int const iz, int const n) noexcept { + return arr_src.contains(ix, iy, iz) ? arr_src(ix, iy, iz, n) : 0.0_rt; + }; - for (int jr = 0; j < weights_y.size(); ++jr) - { - int jj = jjmin + jr; - jj = amrex::Clamp(jj, arr_src.begin.y, arr_src.end.y - 1); + // Interpolate over points computed above. Weights are computed in order + // 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) { + int kk = kkmin + kref; + // kk = amrex::Clamp(kk, arr_src.begin.z, arr_src.end.z - 1); + amrex::Real kfactor = (use_half[2] && (kref == 0 || kref == numz-1)) ? 0.5_rt : 1.0_rt; - for (int ir = 0; i < weights_x.size(); ++ir) - { - int ii = iimin + ir; - ii = amrex::Clamp(ii, arr_src.begin.x, arr_src.end.x - 1); + for (int jref = 0; jref < numy; ++jref) { + int jj = jjmin + jref; + // jj = amrex::Clamp(jj, arr_src.begin.y, arr_src.end.y - 1); + amrex::Real jfactor = (use_half[1] && (jref == 0 || jref == numy-1)) ? 0.5_rt : 1.0_rt; - amrex::Real src_iijjkk = arr_src(ii, jj, kk, comp); - c += weights_x[ir]*weights_y[ir]*weights_z[ir]*src_iijjkk; + for (int iref = 0; iref < numx; ++iref) { + int ii = iimin + iref; + // ii = amrex::Clamp(ii, arr_src.begin.x, arr_src.end.x - 1); + amrex::Real ifactor = (use_half[0] && (iref == 0 || iref == numx-1)) ? 0.5_rt : 1.0_rt; + + c += ifactor * jfactor * kfactor * arr_src_safe(ii, jj, kk, comp); } } } + + c *= crx_cry_crz_inv; + return c; } @@ -204,15 +184,99 @@ namespace ablastr::coarsen::average amrex::IntVect crse_ratio ); + AMREX_GPU_DEVICE AMREX_FORCE_INLINE + amrex::Real + InterpWithWeights ( + amrex::Array4 const &arr_src, + amrex::Vector const &weights_x, + amrex::Vector const &weights_y, + amrex::Vector const &weights_z, + amrex::GpuArray const &src_fine_ind_min, + amrex::GpuArray const &coarsen_ratio, + int i, int j, int k, int comp) + { + const int iimin = i * coarsen_ratio[0] + src_fine_ind_min[0]; + const int jjmin = j * coarsen_ratio[1] + src_fine_ind_min[1]; + const int kkmin = k * coarsen_ratio[2] + src_fine_ind_min[2]; + + auto const arr_src_safe = [arr_src] + AMREX_GPU_DEVICE(int const ix, int const iy, int const iz, int const n) noexcept { + return arr_src.contains(ix, iy, iz) ? arr_src(ix, iy, iz, n) : 0.0_rt; + }; + + amrex::Real c = 0.0_rt; + for (int kref = 0; k < weights_z.size(); ++kref) + { + int kk = kkmin + kref; + // kk = amrex::Clamp(kk, arr_src.begin.z, arr_src.end.z - 1); + amrex::Real wz = weights_z[kref]; + + for (int jref = 0; j < weights_y.size(); ++jref) + { + int jj = jjmin + jref; + // jj = amrex::Clamp(jj, arr_src.begin.y, arr_src.end.y - 1); + amrex::Real wy = weights_y[jref]; + + for (int iref = 0; i < weights_x.size(); ++iref) + { + int ii = iimin + iref; + // ii = amrex::Clamp(ii, arr_src.begin.x, arr_src.end.x - 1); + amrex::Real wx = weights_x[iref]; + + c += wx*wy*wz*arr_src_safe(ii, jj, kk, comp); + } + } + } + return c; + } + + AMREX_GPU_DEVICE AMREX_FORCE_INLINE void CoarseningPointsAndWeights ( + amrex::GpuArray &src_fine_index_min, amrex::Vector &weights_x, amrex::Vector &weights_y, amrex::Vector &weights_z, - amrex::GpuArray const &ind_min, amrex::GpuArray const &coarsen_ratio, - amrex::GpuArray const &stag_fine_src, - amrex::GpuArray const & stag_crse_des); + amrex::GpuArray const &stag_src_fine, + amrex::GpuArray const &stag_des_crse) + { + using namespace amrex::literals; + + int num_points[3]; + bool useHalf[3]; + + for (int l = 0; l < 3; ++l) { + int const twoImin = -coarsen_ratio[l]*stag_des_crse[l] + stag_src_fine[l] - 1; + if (twoImin % 2 == 0) { + src_fine_index_min[l] = twoImin/2; + num_points[l] = coarsen_ratio[l]+1; + useHalf[l] = true; + } else { + src_fine_index_min[l] = (twoImin+1)/2; + num_points[l] = coarsen_ratio[l]; + useHalf[l] = false; + } + } + + weights_x.resize(num_points[0], 1.0_rt); + weights_y.resize(num_points[1], 1.0_rt); + weights_z.resize(num_points[2], 1.0_rt); + + if (useHalf[0]) { + weights_x[0] *= 0.5_rt; + weights_x[num_points[0]-1] *= 0.5_rt; + } + if (useHalf[1]) { + weights_y[0] *= 0.5_rt; + weights_y[num_points[1]-1] *= 0.5_rt; + } + if (useHalf[2]) { + weights_z[0] *= 0.5_rt; + weights_z[num_points[2]-1] *= 0.5_rt; + } + + } template ::value,int>::type = 0> [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE diff --git a/Source/ablastr/coarsen/average.cpp b/Source/ablastr/coarsen/average.cpp index 0fb413ecff1..a5dde3d8aa4 100644 --- a/Source/ablastr/coarsen/average.cpp +++ b/Source/ablastr/coarsen/average.cpp @@ -82,54 +82,4 @@ namespace ablastr::coarsen::average Loop(mf_dst, mf_src, ncomp, ngrow, crse_ratio); } - void - CoarseningPointsAndWeights ( - amrex::Vector &weights_x, - amrex::Vector &weights_y, - amrex::Vector &weights_z, - amrex::GpuArray &src_index_min, - amrex::GpuArray const &coarsen_ratio, - amrex::GpuArray const &stag_fine_src, - amrex::GpuArray const &stag_crse_des - ) - { - using namespace amrex::literals; - - int num_points[3]; - bool useHalf[3]; - - for (int l = 0; l < 3; ++l) { - int const twoImin = -coarsen_ratio[l]*stag_crse_des[l] + stag_fine_src[l] - 1; - if (twoImin % 2 == 0) { - src_index_min[l] = twoImin/2; - num_points[l] = coarsen_ratio[l]+1; - useHalf[l] = true; - } else { - src_index_min[l] = (twoImin+1)/2; - num_points[l] = coarsen_ratio[l]; - useHalf[l] = false; - } - } - - // amrex::Real const coarsen_factor = 1.0_rt / static_cast(coarsen_ratio[0]*coarsen_ratio[1]*coarsen_ratio[2]); - - weights_x.resize(num_points[0], 1.0_rt); - weights_y.resize(num_points[1], 1.0_rt); - weights_z.resize(num_points[2], 1.0_rt); - - if (useHalf[0]) { - weights_x[0] *= 0.5_rt; - weights_x[num_points[0]-1] *= 0.5_rt; - } - if (useHalf[1]) { - weights_y[0] *= 0.5_rt; - weights_y[num_points[1]-1] *= 0.5_rt; - } - if (useHalf[2]) { - weights_z[0] *= 0.5_rt; - weights_z[num_points[2]-1] *= 0.5_rt; - } - - } - } // namespace ablastr::coarsen::average diff --git a/Source/ablastr/coarsen/sample.H b/Source/ablastr/coarsen/sample.H index ac0f9961198..1015c51ff7b 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; @@ -100,7 +99,7 @@ namespace ablastr::coarsen::sample } } } - c*=wx*wy*wz; + c *= wxwywz; return c; } From 14d5a2503d07e7f7556aa4d6c8a8db9fd8499d2d Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Mon, 27 Nov 2023 12:29:36 -0800 Subject: [PATCH 04/26] Rename some variables and rewrite Loop --- Source/ablastr/coarsen/average.H | 71 +++++++++++++++++++----------- Source/ablastr/coarsen/average.cpp | 16 ++++--- 2 files changed, 56 insertions(+), 31 deletions(-) diff --git a/Source/ablastr/coarsen/average.H b/Source/ablastr/coarsen/average.H index 54a37b75dd3..8fab3f395d3 100644 --- a/Source/ablastr/coarsen/average.H +++ b/Source/ablastr/coarsen/average.H @@ -30,17 +30,17 @@ namespace ablastr::coarsen::average * \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 d + * 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_fine staggering of the source fine MultiFab - * \param[in] stag_des_crse 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 + * \param[in] arr_src source fine 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 */ @@ -48,18 +48,19 @@ namespace ablastr::coarsen::average amrex::Real Interp ( amrex::Array4 const &arr_src, - amrex::GpuArray const &stag_src_fine, - amrex::GpuArray const &stag_des_crse, - amrex::GpuArray const &crse_ratio, + amrex::GpuArray const &sf, + amrex::GpuArray const &sc, + amrex::GpuArray const &cr, int const i, int const j, int const k, int const comp) { // Number of points and starting indices of source array (fine) - amrex::GpuArray np, idx_min, use_half; - amrex::Real crx_cry_crz_inv; + amrex::GpuArray np, idx_min; + amrex::GpuArray use_half; + amrex::Real crxyz_inv; - CalculateCoarseningData(idx_min, np, use_half, crx_cry_crz_inv, stag_des_crse, stag_src_fine, crse_ratio); + CalculateCoarseningData(idx_min, np, use_half, crxyz_inv, sf, sc, cr); - return InterpWithCoarseningData(arr_src, idx_min, np, use_half, crx_cry_crz_inv, crse_ratio, i, j, k, comp); + return InterpWithCoarseningData(arr_src, idx_min, np, use_half, crxyz_inv, cr, i, j, k, comp); } AMREX_GPU_DEVICE AMREX_FORCE_INLINE @@ -67,7 +68,7 @@ namespace ablastr::coarsen::average CalculateCoarseningData ( amrex::GpuArray &src_fine_index_min, amrex::GpuArray &num_src_points, - amrex::GpuArray &use_half_weight, // Use int instead bool? + amrex::GpuArray &use_half_weight, // Use int instead bool? amrex::Real &crx_cry_crz_inv, amrex::GpuArray const &stag_src_fine, amrex::GpuArray const &stag_des_crse, @@ -75,7 +76,7 @@ namespace ablastr::coarsen::average { for (int l = 0; l < 3; ++l) { int two_times_index_min = -crse_ratio[l]*stag_des_crse[l] + stag_src_fine[l] - 1; - if (two_times_index_min % 2 == 0) { + if ((two_times_index_min % 2) == 0) { src_fine_index_min[l] = two_times_index_min/2; num_src_points[l] = crse_ratio[l]+1; use_half_weight[l] = 1; // True @@ -94,7 +95,7 @@ namespace ablastr::coarsen::average amrex::Array4 const &arr_src, amrex::GpuArray const &src_fine_index_min, amrex::GpuArray const &num_src_points, - amrex::GpuArray const &use_half, + amrex::GpuArray const &use_half, amrex::Real crx_cry_crz_inv, amrex::GpuArray const &coarsen_ratio, int i, int j, int k, int comp) @@ -122,17 +123,14 @@ namespace ablastr::coarsen::average amrex::Real c = 0.0_rt; for (int kref = 0; kref < numz; ++kref) { int kk = kkmin + kref; - // kk = amrex::Clamp(kk, arr_src.begin.z, arr_src.end.z - 1); amrex::Real kfactor = (use_half[2] && (kref == 0 || kref == numz-1)) ? 0.5_rt : 1.0_rt; for (int jref = 0; jref < numy; ++jref) { int jj = jjmin + jref; - // jj = amrex::Clamp(jj, arr_src.begin.y, arr_src.end.y - 1); amrex::Real jfactor = (use_half[1] && (jref == 0 || jref == numy-1)) ? 0.5_rt : 1.0_rt; for (int iref = 0; iref < numx; ++iref) { int ii = iimin + iref; - // ii = amrex::Clamp(ii, arr_src.begin.x, arr_src.end.x - 1); amrex::Real ifactor = (use_half[0] && (iref == 0 || iref == numx-1)) ? 0.5_rt : 1.0_rt; c += ifactor * jfactor * kfactor * arr_src_safe(ii, jj, kk, comp); @@ -167,6 +165,30 @@ namespace ablastr::coarsen::average amrex::IntVect crse_ratio ); + void + LoopOverBox ( + amrex::Box box, int ncomp, + amrex::Array4 const &arr_dst, + amrex::Array4 const &arr_src, + amrex::GpuArray const &stag_src_fine, + amrex::GpuArray const &stag_des_crse, + amrex::GpuArray const &crse_ratio) + { + // Number of points and starting indices of source array (fine) + amrex::GpuArray idx_min, np; + amrex::GpuArray use_half; + amrex::Real crx_cry_crz_inv; + + CalculateCoarseningData(idx_min, np, use_half, crx_cry_crz_inv, + stag_des_crse, stag_src_fine, crse_ratio); + + amrex::ParallelFor(box, ncomp, + [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) { + arr_dst(i, j, k, n) = InterpWithCoarseningData( + arr_src, idx_min, np, use_half, crx_cry_crz_inv, i, j, k, n); + }); + } + /** * \brief Stores in the coarsened MultiFab \c mf_dst the values obtained by * interpolating the data contained in the fine MultiFab \c mf_src. @@ -208,19 +230,16 @@ namespace ablastr::coarsen::average for (int kref = 0; k < weights_z.size(); ++kref) { int kk = kkmin + kref; - // kk = amrex::Clamp(kk, arr_src.begin.z, arr_src.end.z - 1); amrex::Real wz = weights_z[kref]; for (int jref = 0; j < weights_y.size(); ++jref) { int jj = jjmin + jref; - // jj = amrex::Clamp(jj, arr_src.begin.y, arr_src.end.y - 1); amrex::Real wy = weights_y[jref]; for (int iref = 0; i < weights_x.size(); ++iref) { int ii = iimin + iref; - // ii = amrex::Clamp(ii, arr_src.begin.x, arr_src.end.x - 1); amrex::Real wx = weights_x[iref]; c += wx*wy*wz*arr_src_safe(ii, jj, kk, comp); diff --git a/Source/ablastr/coarsen/average.cpp b/Source/ablastr/coarsen/average.cpp index a5dde3d8aa4..2b91cf2261f 100644 --- a/Source/ablastr/coarsen/average.cpp +++ b/Source/ablastr/coarsen/average.cpp @@ -46,6 +46,12 @@ namespace ablastr::coarsen::average cr[i] = crse_ratio[i]; } + amrex::GpuArray idx_min, np; + amrex::GpuArray use_half; + amrex::Real crxyz_inv; + + CalculateCoarseningData(idx_min, np, use_half, crxyz_inv, sf, sc, cr); + #ifdef AMREX_USE_OMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif @@ -55,11 +61,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, np, use_half, crxyz_inv, cr, i, j, k, n); + }); } } From e2cee58d1b43b0bfc8b6a3981c7e9b0f772301f2 Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Mon, 27 Nov 2023 14:07:14 -0800 Subject: [PATCH 05/26] Bug fixes and move definitions to average.cpp --- Source/ablastr/coarsen/average.H | 126 +++++++---------------------- Source/ablastr/coarsen/average.cpp | 100 +++++++++++++++++++++++ 2 files changed, 130 insertions(+), 96 deletions(-) diff --git a/Source/ablastr/coarsen/average.H b/Source/ablastr/coarsen/average.H index 8fab3f395d3..e8966a73e05 100644 --- a/Source/ablastr/coarsen/average.H +++ b/Source/ablastr/coarsen/average.H @@ -51,17 +51,7 @@ namespace ablastr::coarsen::average amrex::GpuArray const &sf, amrex::GpuArray const &sc, amrex::GpuArray const &cr, - int const i, int const j, int const k, int const comp) - { - // Number of points and starting indices of source array (fine) - amrex::GpuArray np, idx_min; - amrex::GpuArray use_half; - amrex::Real crxyz_inv; - - CalculateCoarseningData(idx_min, np, use_half, crxyz_inv, sf, sc, cr); - - return InterpWithCoarseningData(arr_src, idx_min, np, use_half, crxyz_inv, cr, i, j, k, comp); - } + int const i, int const j, int const k, int const comp); AMREX_GPU_DEVICE AMREX_FORCE_INLINE void @@ -72,22 +62,7 @@ namespace ablastr::coarsen::average amrex::Real &crx_cry_crz_inv, amrex::GpuArray const &stag_src_fine, amrex::GpuArray const &stag_des_crse, - amrex::GpuArray const &crse_ratio) - { - for (int l = 0; l < 3; ++l) { - int two_times_index_min = -crse_ratio[l]*stag_des_crse[l] + stag_src_fine[l] - 1; - if ((two_times_index_min % 2) == 0) { - src_fine_index_min[l] = two_times_index_min/2; - num_src_points[l] = crse_ratio[l]+1; - use_half_weight[l] = 1; // True - } else { - src_fine_index_min[l] = (two_times_index_min+1)/2; - num_src_points[l] = crse_ratio[l]; - use_half_weight[l] = 0; // False - } - } - crx_cry_crz_inv = 1.0_rt / static_cast(crse_ratio[0]*crse_ratio[1]*crse_ratio[2]); - } + amrex::GpuArray const &crse_ratio); AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real @@ -98,50 +73,7 @@ namespace ablastr::coarsen::average amrex::GpuArray const &use_half, amrex::Real crx_cry_crz_inv, amrex::GpuArray const &coarsen_ratio, - int i, int j, int k, int comp) - { - using namespace amrex::literals; - - // Auxiliary integer variables - int const numx = num_src_points[0]; - int const numy = num_src_points[1]; - int const numz = num_src_points[2]; - int const iimin = i * coarsen_ratio[0] + src_fine_index_min[0]; - int const jjmin = j * coarsen_ratio[1] + src_fine_index_min[1]; - int const kkmin = k * coarsen_ratio[2] + src_fine_index_min[2]; - - // Add neutral elements (=0) beyond guard cells in source array (fine) - auto const arr_src_safe = [arr_src] - AMREX_GPU_DEVICE(int const ix, int const iy, int const iz, int const n) noexcept { - return arr_src.contains(ix, iy, iz) ? arr_src(ix, iy, iz, n) : 0.0_rt; - }; - - // Interpolate over points computed above. Weights are computed in order - // 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) { - int kk = kkmin + kref; - amrex::Real kfactor = (use_half[2] && (kref == 0 || kref == numz-1)) ? 0.5_rt : 1.0_rt; - - for (int jref = 0; jref < numy; ++jref) { - int jj = jjmin + jref; - amrex::Real jfactor = (use_half[1] && (jref == 0 || jref == numy-1)) ? 0.5_rt : 1.0_rt; - - for (int iref = 0; iref < numx; ++iref) { - int ii = iimin + iref; - amrex::Real ifactor = (use_half[0] && (iref == 0 || iref == numx-1)) ? 0.5_rt : 1.0_rt; - - c += ifactor * jfactor * kfactor * arr_src_safe(ii, jj, kk, comp); - } - } - } - - c *= crx_cry_crz_inv; - - return c; - } + int i, int j, int k, int comp); /** * \brief Loops over the boxes of the coarsened MultiFab \c mf_dst and fills @@ -165,30 +97,6 @@ namespace ablastr::coarsen::average amrex::IntVect crse_ratio ); - void - LoopOverBox ( - amrex::Box box, int ncomp, - amrex::Array4 const &arr_dst, - amrex::Array4 const &arr_src, - amrex::GpuArray const &stag_src_fine, - amrex::GpuArray const &stag_des_crse, - amrex::GpuArray const &crse_ratio) - { - // Number of points and starting indices of source array (fine) - amrex::GpuArray idx_min, np; - amrex::GpuArray use_half; - amrex::Real crx_cry_crz_inv; - - CalculateCoarseningData(idx_min, np, use_half, crx_cry_crz_inv, - stag_des_crse, stag_src_fine, crse_ratio); - - amrex::ParallelFor(box, ncomp, - [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) { - arr_dst(i, j, k, n) = InterpWithCoarseningData( - arr_src, idx_min, np, use_half, crx_cry_crz_inv, i, j, k, n); - }); - } - /** * \brief Stores in the coarsened MultiFab \c mf_dst the values obtained by * interpolating the data contained in the fine MultiFab \c mf_src. @@ -206,7 +114,7 @@ namespace ablastr::coarsen::average amrex::IntVect crse_ratio ); - AMREX_GPU_DEVICE AMREX_FORCE_INLINE + AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real InterpWithWeights ( amrex::Array4 const &arr_src, @@ -217,6 +125,8 @@ namespace ablastr::coarsen::average amrex::GpuArray const &coarsen_ratio, int i, int j, int k, int comp) { + using namespace amrex::literals; + const int iimin = i * coarsen_ratio[0] + src_fine_ind_min[0]; const int jjmin = j * coarsen_ratio[1] + src_fine_ind_min[1]; const int kkmin = k * coarsen_ratio[2] + src_fine_ind_min[2]; @@ -308,6 +218,30 @@ namespace ablastr::coarsen::average return arr(iClamped,jClamped,kClamped,n); } + // void + // LoopOverBox ( + // amrex::Box box, int ncomp, + // amrex::Array4 const &arr_dst, + // amrex::Array4 const &arr_src, + // amrex::GpuArray const &stag_src_fine, + // amrex::GpuArray const &stag_des_crse, + // amrex::GpuArray const &crse_ratio) + // { + // // Number of points and starting indices of source array (fine) + // amrex::GpuArray idx_min, np; + // amrex::GpuArray use_half; + // amrex::Real crx_cry_crz_inv; + + // CalculateCoarseningData(idx_min, np, use_half, crx_cry_crz_inv, + // stag_des_crse, stag_src_fine, crse_ratio); + + // amrex::ParallelFor(box, ncomp, + // [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) { + // arr_dst(i, j, k, n) = InterpWithCoarseningData( + // arr_src, idx_min, np, use_half, crx_cry_crz_inv, i, j, k, n); + // }); + // } + } // namespace ablastr::coarsen::average #endif // ABLASTR_COARSEN_AVERAGE_H_ diff --git a/Source/ablastr/coarsen/average.cpp b/Source/ablastr/coarsen/average.cpp index 2b91cf2261f..c0fee687368 100644 --- a/Source/ablastr/coarsen/average.cpp +++ b/Source/ablastr/coarsen/average.cpp @@ -22,6 +22,106 @@ namespace ablastr::coarsen::average { + 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) + { + // Number of points and starting indices of source array (fine) + amrex::GpuArray np, idx_min; + amrex::GpuArray use_half; + amrex::Real crxyz_inv; + + CalculateCoarseningData(idx_min, np, use_half, crxyz_inv, sf, sc, cr); + + return InterpWithCoarseningData(arr_src, idx_min, np, use_half, crxyz_inv, cr, i, j, k, comp); + } + + AMREX_GPU_DEVICE AMREX_FORCE_INLINE + void + CalculateCoarseningData ( + amrex::GpuArray &src_fine_index_min, + amrex::GpuArray &num_src_points, + amrex::GpuArray &use_half_weight, // Use int instead bool? + amrex::Real &crx_cry_crz_inv, + amrex::GpuArray const &stag_src_fine, + amrex::GpuArray const &stag_des_crse, + amrex::GpuArray const &crse_ratio) + { + using namespace amrex::literals; + for (int l = 0; l < 3; ++l) { + int two_times_index_min = -crse_ratio[l]*stag_des_crse[l] + stag_src_fine[l] - 1; + if ((two_times_index_min % 2) == 0) { + src_fine_index_min[l] = two_times_index_min/2; + num_src_points[l] = crse_ratio[l]+1; + use_half_weight[l] = 1; // True + } else { + src_fine_index_min[l] = (two_times_index_min+1)/2; + num_src_points[l] = crse_ratio[l]; + use_half_weight[l] = 0; // False + } + } + 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 &src_fine_index_min, + amrex::GpuArray const &num_src_points, + amrex::GpuArray const &use_half, + amrex::Real crx_cry_crz_inv, + amrex::GpuArray const &coarsen_ratio, + int i, int j, int k, int comp) + { + using namespace amrex::literals; + + // Auxiliary integer variables + int const numx = num_src_points[0]; + int const numy = num_src_points[1]; + int const numz = num_src_points[2]; + int const iimin = i * coarsen_ratio[0] + src_fine_index_min[0]; + int const jjmin = j * coarsen_ratio[1] + src_fine_index_min[1]; + int const kkmin = k * coarsen_ratio[2] + src_fine_index_min[2]; + + // Add neutral elements (=0) beyond guard cells in source array (fine) + auto const arr_src_safe = [arr_src] + AMREX_GPU_DEVICE(int const ix, int const iy, int const iz, int const n) noexcept { + return arr_src.contains(ix, iy, iz) ? arr_src(ix, iy, iz, n) : 0.0_rt; + }; + + // Interpolate over points computed above. Weights are computed in order + // 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) { + int kk = kkmin + kref; + amrex::Real kfactor = (use_half[2] && (kref == 0 || kref == numz-1)) ? 0.5_rt : 1.0_rt; + + for (int jref = 0; jref < numy; ++jref) { + int jj = jjmin + jref; + amrex::Real jfactor = (use_half[1] && (jref == 0 || jref == numy-1)) ? 0.5_rt : 1.0_rt; + + for (int iref = 0; iref < numx; ++iref) { + int ii = iimin + iref; + amrex::Real ifactor = (use_half[0] && (iref == 0 || iref == numx-1)) ? 0.5_rt : 1.0_rt; + + c += ifactor * jfactor * kfactor * arr_src_safe(ii, jj, kk, comp); + } + } + } + + c *= crx_cry_crz_inv; + + return c; + } + void Loop ( amrex::MultiFab & mf_dst, From e32f4278d63781296501fc415d2e4b93cf9a4814 Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Mon, 27 Nov 2023 14:44:41 -0800 Subject: [PATCH 06/26] Small adjustments --- Source/ablastr/coarsen/average.H | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/ablastr/coarsen/average.H b/Source/ablastr/coarsen/average.H index e8966a73e05..32c5b2c672c 100644 --- a/Source/ablastr/coarsen/average.H +++ b/Source/ablastr/coarsen/average.H @@ -42,7 +42,7 @@ namespace ablastr::coarsen::average * \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 + * \return interpolated field at cell (i,j,k) of a coarsened Array4 */ AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real @@ -51,7 +51,7 @@ namespace ablastr::coarsen::average amrex::GpuArray const &sf, amrex::GpuArray const &sc, amrex::GpuArray const &cr, - int const i, int const j, int const k, int const comp); + int i, int j, int k, int comp); AMREX_GPU_DEVICE AMREX_FORCE_INLINE void From 403665a0fdbc6178845a1374ebe18467027aa597 Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Tue, 28 Nov 2023 16:08:36 -0800 Subject: [PATCH 07/26] Restore refinement to python checker --- .../Utils/check_interp_points_and_weights.py | 246 +++++++++--------- 1 file changed, 122 insertions(+), 124 deletions(-) diff --git a/Source/Utils/check_interp_points_and_weights.py b/Source/Utils/check_interp_points_and_weights.py index 97c49811c96..b1266ff2dc9 100644 --- a/Source/Utils/check_interp_points_and_weights.py +++ b/Source/Utils/check_interp_points_and_weights.py @@ -29,25 +29,23 @@ import numpy as np -# # 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 ] +# Coarsening functions + +def coarsening_coarse_grid_limits( sc, sf, cr, ii_min, ii_max ): + i_range_start = (ii_min//cr) - 100 + i_range_end = (ii_max//cr) + 1 + 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 ] def coarsening_fine_grid_limits( sc, sf, cr ): if ( sf == 0 ): # cell-centered @@ -55,82 +53,80 @@ def coarsening_fine_grid_limits( sc, sf, cr ): iimax = 4*cr elif ( sf == 1 ): # nodal iimin = 0 - iimax = 4*cr + iimax = 4*cr+1 return [ iimin, iimax ] -def coarsening_coarse_grid_limits( sc, sf, cr, iimin, iimax ): - # imin = int( iimin/cr ) # original - # imax = int( iimax/cr )-(1-sc)*sf+(1-sf)*sc # original - imin = 10000 - imax = -10000 - for i in range(-100,101): - numpts, idxmin, weights = coarsening_points_and_weights(i, sc, sf, cr) - if iimin <= idxmin + numpts - 1: - imin = min(imin,i) - if iimax >= idxmin: - imax = max(imax,i) - return [ imin, imax ] - def coarsening_points_and_weights( i, sc, sf, cr ): - twoImin = -cr*sc + sf - 1 - if (twoImin % 2 == 0): - idxmin = i*cr + twoImin//2 - numpts = cr+1 - weights = np.zeros( numpts ) + 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[numpts-1] = 1.0/(2*cr) - for ir in range( 1, numpts-1 ): + weights[num_ii_pts-1] = 1.0/(2*cr) + for ir in range( 1, num_ii_pts-1 ): weights[ir] = 1.0/cr else: - idxmin = i*cr + (twoImin+1)//2 - numpts = cr - weights = np.zeros( numpts ) - for ir in range( 0, numpts ): + 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 [ numpts, idxmin, weights ] - -# def refinement_fine_grid_limits( sc, sf, cr ): -# if ( sf == 0 ): # cell-centered -# iimin = 0 # original -# iimax = 7 # original -# elif ( sf == 1 ): # nodal -# iimin = 0 # original -# iimax = 8 # original -# return [ iimin, iimax ] - -# def refinement_coarse_grid_limits( sc, sf, cr, iimin, iimax ): -# # imin = int( iimin/cr ) # original -# # imax = int( iimax/cr )-(1-sc)*sf+(1-sf)*sc # original -# imin = 10000 -# imax = -10000 -# for i in range(-100,101): -# numpts, idxmin, weights = refinement_points_and_weights(i, sc, sf, cr) -# if iimin <= idxmin + numpts - 1: -# imin = min(imin,i) -# if iimax >= idxmin: -# imax = max(imax,i) -# return [ imin, imax ] - -# # 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 ] + return [ num_ii_pts, ii_start, weights ] + + +# Refinement functions + +def refinement_coarse_grid_limits( sc, sf, cr ): + # imin = int( iimin/cr ) # original + # imax = int( iimax/cr )-(1-sc)*sf+(1-sf)*sc # original + 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) + 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 ): + if ( cr==1 ): + num_crse_pts = 1 + i_start = ii + elif ( cr>=2 ): + if ( ii%cr==0 ): + num_crse_pts = (1-sf)*(1-sc)+sf*sc + elif ( ii%cr!=0 ): + num_crse_pts = (1-sf)*(1-sc)+2*sf*sc + i_start = (ii//cr)*(1-sf)*(1-sc)+(ii//cr)*sf*sc + weights = np.zeros( num_crse_pts ) + for ir in range( num_crse_pts ): + i = i_start+ir + if ( ii!=iimin or ii!=iimax ): + weights[ir] = (1-sf)*(1-sc)+((abs(cr-abs(ii-i*cr)))/(cr))*sf*sc + else: + weights[ir] = (1-sf)*(1-sc)+((abs(cr-abs(ii-i*cr)))/(cr)+(cr/2-0.5))*sf*sc + return [ num_crse_pts, i_start, weights ] ## TODO Coarsening for IO: interpolation points and weights #def coarsening_points_and_weights_for_IO( i, sf, sc, cr ): @@ -210,52 +206,54 @@ def coarsening_points_and_weights( i, sc, sf, cr ): for ii in range( iimin, iimax+1 ): # index on fine grid ws = 0.0 for i in range( imin, imax+1 ): # 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 + 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 (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( ' ---------------------------------------------------------' ) + print( '\n Refinement for MR: check interpolation points and weights' ) + print( ' ---------------------------------------------------------' ) - # iimin,iimax = refinement_fine_grid_limits( sc, sf, cr ) - # imin ,imax = refinement_coarse_grid_limits( sc, sf, cr, iimin, iimax ) - # # Number of grid points + imin ,imax = refinement_coarse_grid_limits( sc, sf, cr) + iimin,iimax = refinement_fine_grid_limits( sc, sf, cr, imin, imax ) + + # Number of grid points # nc = imax-imin+1 # nf = iimax-iimin+1 - # 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 ) ) + 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 ) ) # print( '\n Number of points on coarse grid: nc={}'.format( nc ) ) # print( ' Number of points on fine grid: nf={}'.format( nf ) ) - # # 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 - # sign = '+' - # if (i == 0): - # sign = ' ' - # if (i < 0): - # sign = '-' - # print( ' (i='+sign+'{},w={:.2f})'.format( abs(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 ( abs(ws - cr) > 1E-9 ): - # print( '\n ERROR: sum of weights ws={:.2f} should be cr={} for i={}'.format( ws, cr, i ) ) + # Refinement for MR: interpolation points and weights + 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 ) + print( '\n Find value at ii={} by interpolating over the following points and weights:'.format( ii ) ) + for ir in range( num_i_pts ): # interpolation points and weights + i = i_start+ir + # sign = '+' + # if (i == 0): + # sign = ' ' + # if (i < 0): + # sign = '-' + # print( ' (i='+sign+'{},w={:.2f})'.format( abs(i), weights[ir] ), end='' ) + print( ' (i={},w={:.3f})'.format( i, weights[ir] ), end='' ) + if not ( ir == num_i_pts-1 ): + print( ' ', end='' ) + print() + + # Refinement for MR: check conservation properties + for i in range( imin, imax+1 ): # index on coarse grid + ws = 0.0 + for ii in range( iimin, iimax+1 ): # index on fine grid + num_i_pts,idxmin,weights = refinement_points_and_weights( ii, sc, sf, cr ) + 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 ( abs(ws - cr) > 1E-9 ): + print( '\n ERROR: sum of weights ws={:.3f} should be cr={} for i={}'.format( ws, cr, i ) ) From 61f2a10cd11771f47dd8e031d5b663e164931423 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 30 Nov 2023 20:15:28 +0000 Subject: [PATCH 08/26] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- Source/Utils/check_interp_points_and_weights.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Source/Utils/check_interp_points_and_weights.py b/Source/Utils/check_interp_points_and_weights.py index b1266ff2dc9..0fc8f807320 100644 --- a/Source/Utils/check_interp_points_and_weights.py +++ b/Source/Utils/check_interp_points_and_weights.py @@ -25,7 +25,6 @@ # Source/ablastr/coarsen/sample.(H/.cpp) #------------------------------------------------------------------------------- -import sys import numpy as np @@ -219,7 +218,7 @@ def refinement_points_and_weights( ii, sc, sf, cr ): imin ,imax = refinement_coarse_grid_limits( sc, sf, cr) iimin,iimax = refinement_fine_grid_limits( sc, sf, cr, imin, imax ) - + # Number of grid points # nc = imax-imin+1 # nf = iimax-iimin+1 From 3aa882f38c8ec93d21eb28b9eb4107894f3f0c94 Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Thu, 30 Nov 2023 12:28:10 -0800 Subject: [PATCH 09/26] Delete old functions --- .../Utils/check_interp_points_and_weights.py | 30 +--- Source/ablastr/coarsen/average.H | 128 ------------------ 2 files changed, 4 insertions(+), 154 deletions(-) diff --git a/Source/Utils/check_interp_points_and_weights.py b/Source/Utils/check_interp_points_and_weights.py index 0fc8f807320..adfb54c314f 100644 --- a/Source/Utils/check_interp_points_and_weights.py +++ b/Source/Utils/check_interp_points_and_weights.py @@ -78,8 +78,6 @@ def coarsening_points_and_weights( i, sc, sf, cr ): # Refinement functions def refinement_coarse_grid_limits( sc, sf, cr ): - # imin = int( iimin/cr ) # original - # imax = int( iimax/cr )-(1-sc)*sf+(1-sf)*sc # original i_min = 0 i_max = 3 return [ i_min, i_max ] @@ -146,9 +144,6 @@ 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 ) ) # Loop over possible staggering of coarse and fine grid (cell-centered or nodal) for sc in [0,1]: @@ -167,10 +162,6 @@ def refinement_points_and_weights( ii, sc, sf, cr ): print( ' nodal *' ) print( ' **************************************************' ) - # 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( ' ---------------------------------------------------------' ) @@ -180,12 +171,6 @@ def refinement_points_and_weights( ii, sc, sf, cr ): 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 ) ) - # Coarsening for MR: interpolation points and weights for i in range( imin, imax+1 ): # index on coarse grid numpts,idxmin,weights = coarsening_points_and_weights( i, sc, sf, cr ) @@ -216,17 +201,16 @@ def refinement_points_and_weights( ii, sc, sf, cr ): 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 - # nc = imax-imin+1 - # nf = iimax-iimin+1 - 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 ) ) - # print( '\n Number of points on coarse grid: nc={}'.format( nc ) ) - # print( ' Number of points on fine grid: nf={}'.format( nf ) ) # Refinement for MR: interpolation points and weights for ii in range ( iimin, iimax+1): # index on fine grid @@ -234,12 +218,6 @@ def 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( num_i_pts ): # interpolation points and weights i = i_start+ir - # sign = '+' - # if (i == 0): - # sign = ' ' - # if (i < 0): - # sign = '-' - # print( ' (i='+sign+'{},w={:.2f})'.format( abs(i), weights[ir] ), end='' ) print( ' (i={},w={:.3f})'.format( i, weights[ir] ), end='' ) if not ( ir == num_i_pts-1 ): print( ' ', end='' ) diff --git a/Source/ablastr/coarsen/average.H b/Source/ablastr/coarsen/average.H index 32c5b2c672c..1fbe7b18172 100644 --- a/Source/ablastr/coarsen/average.H +++ b/Source/ablastr/coarsen/average.H @@ -114,134 +114,6 @@ namespace ablastr::coarsen::average amrex::IntVect crse_ratio ); - AMREX_GPU_DEVICE AMREX_FORCE_INLINE - amrex::Real - InterpWithWeights ( - amrex::Array4 const &arr_src, - amrex::Vector const &weights_x, - amrex::Vector const &weights_y, - amrex::Vector const &weights_z, - amrex::GpuArray const &src_fine_ind_min, - amrex::GpuArray const &coarsen_ratio, - int i, int j, int k, int comp) - { - using namespace amrex::literals; - - const int iimin = i * coarsen_ratio[0] + src_fine_ind_min[0]; - const int jjmin = j * coarsen_ratio[1] + src_fine_ind_min[1]; - const int kkmin = k * coarsen_ratio[2] + src_fine_ind_min[2]; - - auto const arr_src_safe = [arr_src] - AMREX_GPU_DEVICE(int const ix, int const iy, int const iz, int const n) noexcept { - return arr_src.contains(ix, iy, iz) ? arr_src(ix, iy, iz, n) : 0.0_rt; - }; - - amrex::Real c = 0.0_rt; - for (int kref = 0; k < weights_z.size(); ++kref) - { - int kk = kkmin + kref; - amrex::Real wz = weights_z[kref]; - - for (int jref = 0; j < weights_y.size(); ++jref) - { - int jj = jjmin + jref; - amrex::Real wy = weights_y[jref]; - - for (int iref = 0; i < weights_x.size(); ++iref) - { - int ii = iimin + iref; - amrex::Real wx = weights_x[iref]; - - c += wx*wy*wz*arr_src_safe(ii, jj, kk, comp); - } - } - } - return c; - } - - AMREX_GPU_DEVICE AMREX_FORCE_INLINE - void - CoarseningPointsAndWeights ( - amrex::GpuArray &src_fine_index_min, - amrex::Vector &weights_x, - amrex::Vector &weights_y, - amrex::Vector &weights_z, - amrex::GpuArray const &coarsen_ratio, - amrex::GpuArray const &stag_src_fine, - amrex::GpuArray const &stag_des_crse) - { - using namespace amrex::literals; - - int num_points[3]; - bool useHalf[3]; - - for (int l = 0; l < 3; ++l) { - int const twoImin = -coarsen_ratio[l]*stag_des_crse[l] + stag_src_fine[l] - 1; - if (twoImin % 2 == 0) { - src_fine_index_min[l] = twoImin/2; - num_points[l] = coarsen_ratio[l]+1; - useHalf[l] = true; - } else { - src_fine_index_min[l] = (twoImin+1)/2; - num_points[l] = coarsen_ratio[l]; - useHalf[l] = false; - } - } - - weights_x.resize(num_points[0], 1.0_rt); - weights_y.resize(num_points[1], 1.0_rt); - weights_z.resize(num_points[2], 1.0_rt); - - if (useHalf[0]) { - weights_x[0] *= 0.5_rt; - weights_x[num_points[0]-1] *= 0.5_rt; - } - if (useHalf[1]) { - weights_y[0] *= 0.5_rt; - weights_y[num_points[1]-1] *= 0.5_rt; - } - if (useHalf[2]) { - weights_z[0] *= 0.5_rt; - weights_z[num_points[2]-1] *= 0.5_rt; - } - - } - - template ::value,int>::type = 0> - [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE - U& Array4ClampedEval (amrex::Array4 const& arr, - int i, int j, int k, int n) - { - int iClamped = amrex::Clamp(i, arr.begin.x, arr.end.x - 1); - int jClamped = amrex::Clamp(j, arr.begin.y, arr.end.y - 1); - int kClamped = amrex::Clamp(k, arr.begin.x, arr.end.z - 1); - return arr(iClamped,jClamped,kClamped,n); - } - - // void - // LoopOverBox ( - // amrex::Box box, int ncomp, - // amrex::Array4 const &arr_dst, - // amrex::Array4 const &arr_src, - // amrex::GpuArray const &stag_src_fine, - // amrex::GpuArray const &stag_des_crse, - // amrex::GpuArray const &crse_ratio) - // { - // // Number of points and starting indices of source array (fine) - // amrex::GpuArray idx_min, np; - // amrex::GpuArray use_half; - // amrex::Real crx_cry_crz_inv; - - // CalculateCoarseningData(idx_min, np, use_half, crx_cry_crz_inv, - // stag_des_crse, stag_src_fine, crse_ratio); - - // amrex::ParallelFor(box, ncomp, - // [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) { - // arr_dst(i, j, k, n) = InterpWithCoarseningData( - // arr_src, idx_min, np, use_half, crx_cry_crz_inv, i, j, k, n); - // }); - // } - } // namespace ablastr::coarsen::average #endif // ABLASTR_COARSEN_AVERAGE_H_ From b2a0590edc61a93cca5f3196ccbe5fa43706eb17 Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Thu, 30 Nov 2023 12:34:50 -0800 Subject: [PATCH 10/26] Change 1,0 to true,false --- Source/ablastr/coarsen/average.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/ablastr/coarsen/average.cpp b/Source/ablastr/coarsen/average.cpp index c0fee687368..a75b867add0 100644 --- a/Source/ablastr/coarsen/average.cpp +++ b/Source/ablastr/coarsen/average.cpp @@ -58,11 +58,11 @@ namespace ablastr::coarsen::average if ((two_times_index_min % 2) == 0) { src_fine_index_min[l] = two_times_index_min/2; num_src_points[l] = crse_ratio[l]+1; - use_half_weight[l] = 1; // True + use_half_weight[l] = true; // True } else { src_fine_index_min[l] = (two_times_index_min+1)/2; num_src_points[l] = crse_ratio[l]; - use_half_weight[l] = 0; // False + use_half_weight[l] = false; // False } } crx_cry_crz_inv = 1.0_rt / static_cast(crse_ratio[0]*crse_ratio[1]*crse_ratio[2]); From a8901dae20cbbc15b46c9dab147f23a29b0f71b8 Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Thu, 30 Nov 2023 12:38:42 -0800 Subject: [PATCH 11/26] Make CodeQL happy by initializing Python variables --- Source/Utils/check_interp_points_and_weights.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Source/Utils/check_interp_points_and_weights.py b/Source/Utils/check_interp_points_and_weights.py index adfb54c314f..72f3b02bf38 100644 --- a/Source/Utils/check_interp_points_and_weights.py +++ b/Source/Utils/check_interp_points_and_weights.py @@ -107,6 +107,8 @@ def refinement_fine_grid_limits( sc, sf, cr, i_min, i_max ): # Refinement for MR: interpolation points and weights def refinement_points_and_weights( ii, sc, sf, cr ): + i_start = 0 + num_crse_pts = 0 if ( cr==1 ): num_crse_pts = 1 i_start = ii From 83ef047183bfb3dc41a18f92d5ac1f6931d21dfa Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Thu, 30 Nov 2023 12:41:39 -0800 Subject: [PATCH 12/26] Rename num_crse_pts and fix typo in Python --- .../Utils/check_interp_points_and_weights.py | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/Source/Utils/check_interp_points_and_weights.py b/Source/Utils/check_interp_points_and_weights.py index 72f3b02bf38..2ec6d14a3c4 100644 --- a/Source/Utils/check_interp_points_and_weights.py +++ b/Source/Utils/check_interp_points_and_weights.py @@ -108,24 +108,24 @@ def refinement_fine_grid_limits( sc, sf, cr, i_min, i_max ): # Refinement for MR: interpolation points and weights def refinement_points_and_weights( ii, sc, sf, cr ): i_start = 0 - num_crse_pts = 0 + num_i_pts = 0 if ( cr==1 ): - num_crse_pts = 1 + num_i_pts = 1 i_start = ii elif ( cr>=2 ): if ( ii%cr==0 ): - num_crse_pts = (1-sf)*(1-sc)+sf*sc + num_i_pts = (1-sf)*(1-sc)+sf*sc elif ( ii%cr!=0 ): - num_crse_pts = (1-sf)*(1-sc)+2*sf*sc + num_i_pts = (1-sf)*(1-sc)+2*sf*sc i_start = (ii//cr)*(1-sf)*(1-sc)+(ii//cr)*sf*sc - weights = np.zeros( num_crse_pts ) - for ir in range( num_crse_pts ): + weights = np.zeros( num_i_pts ) + for ir in range( num_i_pts ): i = i_start+ir - if ( ii!=iimin or ii!=iimax ): - weights[ir] = (1-sf)*(1-sc)+((abs(cr-abs(ii-i*cr)))/(cr))*sf*sc - else: + 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 - return [ num_crse_pts, i_start, weights ] + else: + weights[ir] = (1-sf)*(1-sc)+((abs(cr-abs(ii-i*cr)))/(cr))*sf*sc + return [ num_i_pts, i_start, weights ] ## TODO Coarsening for IO: interpolation points and weights #def coarsening_points_and_weights_for_IO( i, sf, sc, cr ): From b669389476b9689cd4f52f37afe2b7d7af6e8deb Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Thu, 30 Nov 2023 12:44:08 -0800 Subject: [PATCH 13/26] Fix indent error --- Source/Utils/check_interp_points_and_weights.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Utils/check_interp_points_and_weights.py b/Source/Utils/check_interp_points_and_weights.py index 2ec6d14a3c4..df1a5ad5bb6 100644 --- a/Source/Utils/check_interp_points_and_weights.py +++ b/Source/Utils/check_interp_points_and_weights.py @@ -203,7 +203,7 @@ def refinement_points_and_weights( ii, sc, sf, cr ): print( '\n Refinement for MR: check interpolation points and weights' ) print( ' ---------------------------------------------------------' ) - if ( sf!=sc ): + if ( sf!=sc ): print( '\n WARNING: sc={} not equal to sf={}, not implemented for Refinement for MR, continue ...'.format( sc, sf ) ) continue From 5f8485ce75f0ccccbf8ebd4c1d48524135fd0e48 Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Thu, 30 Nov 2023 17:11:48 -0800 Subject: [PATCH 14/26] Move inline functions to header file --- Source/ablastr/coarsen/average.H | 93 ++++++++++++++++++++++++++- Source/ablastr/coarsen/average.cpp | 100 ----------------------------- 2 files changed, 92 insertions(+), 101 deletions(-) diff --git a/Source/ablastr/coarsen/average.H b/Source/ablastr/coarsen/average.H index 1fbe7b18172..91301ec3b62 100644 --- a/Source/ablastr/coarsen/average.H +++ b/Source/ablastr/coarsen/average.H @@ -51,7 +51,98 @@ namespace ablastr::coarsen::average amrex::GpuArray const &sf, amrex::GpuArray const &sc, amrex::GpuArray const &cr, - int i, int j, int k, int comp); + int const i, int const j, int const k, int const comp) + { + // Number of points and starting indices of source array (fine) + amrex::GpuArray np, idx_min; + amrex::GpuArray use_half; + amrex::Real crxyz_inv; + + CalculateCoarseningData(idx_min, np, use_half, crxyz_inv, sf, sc, cr); + + return InterpWithCoarseningData(arr_src, idx_min, np, use_half, crxyz_inv, cr, i, j, k, comp); + } + + AMREX_GPU_DEVICE AMREX_FORCE_INLINE + void + CalculateCoarseningData ( + amrex::GpuArray &src_fine_index_min, + amrex::GpuArray &num_src_points, + amrex::GpuArray &use_half_weight, // Use int instead bool? + amrex::Real &crx_cry_crz_inv, + amrex::GpuArray const &stag_src_fine, + amrex::GpuArray const &stag_des_crse, + amrex::GpuArray const &crse_ratio) + { + using namespace amrex::literals; + for (int l = 0; l < 3; ++l) { + int two_times_index_min = -crse_ratio[l]*stag_des_crse[l] + stag_src_fine[l] - 1; + if ((two_times_index_min % 2) == 0) { + src_fine_index_min[l] = two_times_index_min/2; + num_src_points[l] = crse_ratio[l]+1; + use_half_weight[l] = true; // True + } else { + src_fine_index_min[l] = (two_times_index_min+1)/2; + num_src_points[l] = crse_ratio[l]; + use_half_weight[l] = false; // False + } + } + 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 &src_fine_index_min, + amrex::GpuArray const &num_src_points, + amrex::GpuArray const &use_half, + amrex::Real crx_cry_crz_inv, + amrex::GpuArray const &coarsen_ratio, + int i, int j, int k, int comp) + { + using namespace amrex::literals; + + // Auxiliary integer variables + int const numx = num_src_points[0]; + int const numy = num_src_points[1]; + int const numz = num_src_points[2]; + int const iimin = i * coarsen_ratio[0] + src_fine_index_min[0]; + int const jjmin = j * coarsen_ratio[1] + src_fine_index_min[1]; + int const kkmin = k * coarsen_ratio[2] + src_fine_index_min[2]; + + // Add neutral elements (=0) beyond guard cells in source array (fine) + auto const arr_src_safe = [arr_src] + AMREX_GPU_DEVICE(int const ix, int const iy, int const iz, int const n) noexcept { + return arr_src.contains(ix, iy, iz) ? arr_src(ix, iy, iz, n) : 0.0_rt; + }; + + // Interpolate over points computed above. Weights are computed in order + // 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) { + int kk = kkmin + kref; + amrex::Real kfactor = (use_half[2] && (kref == 0 || kref == numz-1)) ? 0.5_rt : 1.0_rt; + + for (int jref = 0; jref < numy; ++jref) { + int jj = jjmin + jref; + amrex::Real jfactor = (use_half[1] && (jref == 0 || jref == numy-1)) ? 0.5_rt : 1.0_rt; + + for (int iref = 0; iref < numx; ++iref) { + int ii = iimin + iref; + amrex::Real ifactor = (use_half[0] && (iref == 0 || iref == numx-1)) ? 0.5_rt : 1.0_rt; + + c += ifactor * jfactor * kfactor * arr_src_safe(ii, jj, kk, comp); + } + } + } + + c *= crx_cry_crz_inv; + + return c; + } AMREX_GPU_DEVICE AMREX_FORCE_INLINE void diff --git a/Source/ablastr/coarsen/average.cpp b/Source/ablastr/coarsen/average.cpp index a75b867add0..2b91cf2261f 100644 --- a/Source/ablastr/coarsen/average.cpp +++ b/Source/ablastr/coarsen/average.cpp @@ -22,106 +22,6 @@ namespace ablastr::coarsen::average { - 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) - { - // Number of points and starting indices of source array (fine) - amrex::GpuArray np, idx_min; - amrex::GpuArray use_half; - amrex::Real crxyz_inv; - - CalculateCoarseningData(idx_min, np, use_half, crxyz_inv, sf, sc, cr); - - return InterpWithCoarseningData(arr_src, idx_min, np, use_half, crxyz_inv, cr, i, j, k, comp); - } - - AMREX_GPU_DEVICE AMREX_FORCE_INLINE - void - CalculateCoarseningData ( - amrex::GpuArray &src_fine_index_min, - amrex::GpuArray &num_src_points, - amrex::GpuArray &use_half_weight, // Use int instead bool? - amrex::Real &crx_cry_crz_inv, - amrex::GpuArray const &stag_src_fine, - amrex::GpuArray const &stag_des_crse, - amrex::GpuArray const &crse_ratio) - { - using namespace amrex::literals; - for (int l = 0; l < 3; ++l) { - int two_times_index_min = -crse_ratio[l]*stag_des_crse[l] + stag_src_fine[l] - 1; - if ((two_times_index_min % 2) == 0) { - src_fine_index_min[l] = two_times_index_min/2; - num_src_points[l] = crse_ratio[l]+1; - use_half_weight[l] = true; // True - } else { - src_fine_index_min[l] = (two_times_index_min+1)/2; - num_src_points[l] = crse_ratio[l]; - use_half_weight[l] = false; // False - } - } - 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 &src_fine_index_min, - amrex::GpuArray const &num_src_points, - amrex::GpuArray const &use_half, - amrex::Real crx_cry_crz_inv, - amrex::GpuArray const &coarsen_ratio, - int i, int j, int k, int comp) - { - using namespace amrex::literals; - - // Auxiliary integer variables - int const numx = num_src_points[0]; - int const numy = num_src_points[1]; - int const numz = num_src_points[2]; - int const iimin = i * coarsen_ratio[0] + src_fine_index_min[0]; - int const jjmin = j * coarsen_ratio[1] + src_fine_index_min[1]; - int const kkmin = k * coarsen_ratio[2] + src_fine_index_min[2]; - - // Add neutral elements (=0) beyond guard cells in source array (fine) - auto const arr_src_safe = [arr_src] - AMREX_GPU_DEVICE(int const ix, int const iy, int const iz, int const n) noexcept { - return arr_src.contains(ix, iy, iz) ? arr_src(ix, iy, iz, n) : 0.0_rt; - }; - - // Interpolate over points computed above. Weights are computed in order - // 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) { - int kk = kkmin + kref; - amrex::Real kfactor = (use_half[2] && (kref == 0 || kref == numz-1)) ? 0.5_rt : 1.0_rt; - - for (int jref = 0; jref < numy; ++jref) { - int jj = jjmin + jref; - amrex::Real jfactor = (use_half[1] && (jref == 0 || jref == numy-1)) ? 0.5_rt : 1.0_rt; - - for (int iref = 0; iref < numx; ++iref) { - int ii = iimin + iref; - amrex::Real ifactor = (use_half[0] && (iref == 0 || iref == numx-1)) ? 0.5_rt : 1.0_rt; - - c += ifactor * jfactor * kfactor * arr_src_safe(ii, jj, kk, comp); - } - } - } - - c *= crx_cry_crz_inv; - - return c; - } - void Loop ( amrex::MultiFab & mf_dst, From f00e8993fcd22bfb39e1678cd80220575f62672d Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Thu, 30 Nov 2023 17:31:38 -0800 Subject: [PATCH 15/26] Python file works for any coarsening ratio, not just 1,2, or 4 --- Source/Utils/check_interp_points_and_weights.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Utils/check_interp_points_and_weights.py b/Source/Utils/check_interp_points_and_weights.py index df1a5ad5bb6..85394856478 100644 --- a/Source/Utils/check_interp_points_and_weights.py +++ b/Source/Utils/check_interp_points_and_weights.py @@ -145,7 +145,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=" ) ) +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]: From 7b15051e85bcd73d0532ddc41eb5b238bfe21c8f Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Thu, 30 Nov 2023 17:36:57 -0800 Subject: [PATCH 16/26] Reorder functions in header file --- Source/ablastr/coarsen/average.H | 90 ++++++++++++-------------------- 1 file changed, 34 insertions(+), 56 deletions(-) diff --git a/Source/ablastr/coarsen/average.H b/Source/ablastr/coarsen/average.H index 91301ec3b62..c4454a0bf2f 100644 --- a/Source/ablastr/coarsen/average.H +++ b/Source/ablastr/coarsen/average.H @@ -25,44 +25,6 @@ */ 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. - * - * 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] 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) - { - // Number of points and starting indices of source array (fine) - amrex::GpuArray np, idx_min; - amrex::GpuArray use_half; - amrex::Real crxyz_inv; - - CalculateCoarseningData(idx_min, np, use_half, crxyz_inv, sf, sc, cr); - - return InterpWithCoarseningData(arr_src, idx_min, np, use_half, crxyz_inv, cr, i, j, k, comp); - } - AMREX_GPU_DEVICE AMREX_FORCE_INLINE void CalculateCoarseningData ( @@ -144,27 +106,43 @@ namespace ablastr::coarsen::average return c; } - AMREX_GPU_DEVICE AMREX_FORCE_INLINE - void - CalculateCoarseningData ( - amrex::GpuArray &src_fine_index_min, - amrex::GpuArray &num_src_points, - amrex::GpuArray &use_half_weight, // Use int instead bool? - amrex::Real &crx_cry_crz_inv, - amrex::GpuArray const &stag_src_fine, - amrex::GpuArray const &stag_des_crse, - amrex::GpuArray const &crse_ratio); - + /** + * \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] 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 - InterpWithCoarseningData ( + Interp ( amrex::Array4 const &arr_src, - amrex::GpuArray const &src_fine_index_min, - amrex::GpuArray const &num_src_points, - amrex::GpuArray const &use_half, - amrex::Real crx_cry_crz_inv, - amrex::GpuArray const &coarsen_ratio, - int i, int j, int k, int comp); + amrex::GpuArray const &sf, + amrex::GpuArray const &sc, + amrex::GpuArray const &cr, + int const i, int const j, int const k, int const comp) + { + // Number of points and starting indices of source array (fine) + amrex::GpuArray np, idx_min; + amrex::GpuArray use_half; + amrex::Real crxyz_inv; + + CalculateCoarseningData(idx_min, np, use_half, crxyz_inv, sf, sc, cr); + + return InterpWithCoarseningData(arr_src, idx_min, np, use_half, crxyz_inv, cr, i, j, k, comp); + } /** * \brief Loops over the boxes of the coarsened MultiFab \c mf_dst and fills From 2e443336c4a05c87b6c29100b069bc4cc6867511 Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Mon, 11 Dec 2023 09:10:12 -0800 Subject: [PATCH 17/26] 2023 Dec 11 update --- .../Utils/check_interp_points_and_weights.py | 18 +++++++-------- Source/ablastr/coarsen/average.H | 22 +++++++++++-------- 2 files changed, 22 insertions(+), 18 deletions(-) diff --git a/Source/Utils/check_interp_points_and_weights.py b/Source/Utils/check_interp_points_and_weights.py index 85394856478..a61f5b625d9 100644 --- a/Source/Utils/check_interp_points_and_weights.py +++ b/Source/Utils/check_interp_points_and_weights.py @@ -32,7 +32,7 @@ def coarsening_coarse_grid_limits( sc, sf, cr, ii_min, ii_max ): i_range_start = (ii_min//cr) - 100 - i_range_end = (ii_max//cr) + 1 + 100 + i_range_end = (ii_max//cr) + 100 i_min = i_range_end i_max = i_range_start @@ -41,18 +41,18 @@ def coarsening_coarse_grid_limits( sc, sf, cr, ii_min, ii_max ): 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) + i_min = min(i_min, i) if (ii_max >= ii_start): - i_max = max(i_max,i) + i_max = max(i_max, i) return [ i_min, i_max ] def coarsening_fine_grid_limits( sc, sf, cr ): if ( sf == 0 ): # cell-centered - iimin = 1 - iimax = 4*cr + iimin = 0 + iimax = 4*cr - 1 elif ( sf == 1 ): # nodal iimin = 0 - iimax = 4*cr+1 + iimax = 4*cr return [ iimin, iimax ] def coarsening_points_and_weights( i, sc, sf, cr ): @@ -66,7 +66,7 @@ def coarsening_points_and_weights( i, sc, sf, 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 + 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 ): @@ -97,9 +97,9 @@ def refinement_fine_grid_limits( sc, sf, cr, i_min, i_max ): num_i_pts, i_start, weights = refinement_points_and_weights(ii, sc, sf, cr) i_end = i_start + num_i_pts - 1 if i_min <= i_end: - ii_min = min(ii_min,ii) + ii_min = min(ii_min, ii) if i_max >= i_start: - ii_max = max(ii_max,ii) + ii_max = max(ii_max, ii) print("After ii_min={} and ii_max={}".format(ii_min,ii_max)) diff --git a/Source/ablastr/coarsen/average.H b/Source/ablastr/coarsen/average.H index c4454a0bf2f..5e84bfe2fb1 100644 --- a/Source/ablastr/coarsen/average.H +++ b/Source/ablastr/coarsen/average.H @@ -25,7 +25,7 @@ */ namespace ablastr::coarsen::average { - AMREX_GPU_DEVICE AMREX_FORCE_INLINE + AMREX_FORCE_INLINE void CalculateCoarseningData ( amrex::GpuArray &src_fine_index_min, @@ -37,19 +37,23 @@ namespace ablastr::coarsen::average amrex::GpuArray const &crse_ratio) { using namespace amrex::literals; + for (int l = 0; l < 3; ++l) { - int two_times_index_min = -crse_ratio[l]*stag_des_crse[l] + stag_src_fine[l] - 1; - if ((two_times_index_min % 2) == 0) { - src_fine_index_min[l] = two_times_index_min/2; - num_src_points[l] = crse_ratio[l]+1; - use_half_weight[l] = true; // True + int sf_minus_cr_sc = stag_src_fine[l] - crse_ratio[l]*stag_des_crse[l]; + if (sf_minus_cr_sc & 1) { + // sf - cr*sc is odd + src_fine_index_min[l] = (sf_minus_cr_sc - 1)/2; + num_src_points[l] = crse_ratio[l] + 1; + use_half_weight[l] = true; // Use half weight at end points } else { - src_fine_index_min[l] = (two_times_index_min+1)/2; + // sf - cr*sc is even + src_fine_index_min[l] = sf_minus_cr_sc/2; num_src_points[l] = crse_ratio[l]; - use_half_weight[l] = false; // False + 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]); + + crx_cry_crz_inv = 1.0_rt/static_cast(crse_ratio[0]*crse_ratio[1]*crse_ratio[2]); } AMREX_GPU_DEVICE AMREX_FORCE_INLINE From 48ef97e88f2e0cea3a62b32da04e866ee3be1222 Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Mon, 11 Dec 2023 10:59:52 -0800 Subject: [PATCH 18/26] Rename variables and paste contents of CalculateCoarseningData into Interp --- Source/ablastr/coarsen/average.H | 118 +++++++++++++++++------------ Source/ablastr/coarsen/average.cpp | 27 +++---- 2 files changed, 85 insertions(+), 60 deletions(-) diff --git a/Source/ablastr/coarsen/average.H b/Source/ablastr/coarsen/average.H index 5e84bfe2fb1..746fa047e6e 100644 --- a/Source/ablastr/coarsen/average.H +++ b/Source/ablastr/coarsen/average.H @@ -25,30 +25,29 @@ */ namespace ablastr::coarsen::average { - AMREX_FORCE_INLINE void CalculateCoarseningData ( - amrex::GpuArray &src_fine_index_min, - amrex::GpuArray &num_src_points, + 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_fine, - amrex::GpuArray const &stag_des_crse, + amrex::GpuArray const &stag_src, + amrex::GpuArray const &stag_des, amrex::GpuArray const &crse_ratio) { using namespace amrex::literals; for (int l = 0; l < 3; ++l) { - int sf_minus_cr_sc = stag_src_fine[l] - crse_ratio[l]*stag_des_crse[l]; - if (sf_minus_cr_sc & 1) { - // sf - cr*sc is odd - src_fine_index_min[l] = (sf_minus_cr_sc - 1)/2; - num_src_points[l] = crse_ratio[l] + 1; + 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 { - // sf - cr*sc is even - src_fine_index_min[l] = sf_minus_cr_sc/2; - num_src_points[l] = crse_ratio[l]; + // 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 } } @@ -60,22 +59,22 @@ namespace ablastr::coarsen::average amrex::Real InterpWithCoarseningData ( amrex::Array4 const &arr_src, - amrex::GpuArray const &src_fine_index_min, - amrex::GpuArray const &num_src_points, - amrex::GpuArray const &use_half, + 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 &coarsen_ratio, + amrex::GpuArray const &crse_ratio, int i, int j, int k, int comp) { using namespace amrex::literals; // Auxiliary integer variables - int const numx = num_src_points[0]; - int const numy = num_src_points[1]; - int const numz = num_src_points[2]; - int const iimin = i * coarsen_ratio[0] + src_fine_index_min[0]; - int const jjmin = j * coarsen_ratio[1] + src_fine_index_min[1]; - int const kkmin = k * coarsen_ratio[2] + src_fine_index_min[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] @@ -88,19 +87,23 @@ namespace ablastr::coarsen::average // 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 kref = 0; kref <= irefmax; ++kref) { int kk = kkmin + kref; - amrex::Real kfactor = (use_half[2] && (kref == 0 || kref == numz-1)) ? 0.5_rt : 1.0_rt; + bool use_half_kk = use_half_weight[2] && (kref == 0 || kref == irefmax); - for (int jref = 0; jref < numy; ++jref) { + for (int jref = 0; jref <= jrefmax; ++jref) { int jj = jjmin + jref; - amrex::Real jfactor = (use_half[1] && (jref == 0 || jref == numy-1)) ? 0.5_rt : 1.0_rt; + bool use_half_jj = use_half_weight[1] && (jref == 0 || jref == jrefmax); - for (int iref = 0; iref < numx; ++iref) { + for (int iref = 0; iref <= krefmax; ++iref) { int ii = iimin + iref; - amrex::Real ifactor = (use_half[0] && (iref == 0 || iref == numx-1)) ? 0.5_rt : 1.0_rt; + bool use_half_ii = use_half_weight[0] && (iref == 0 || iref == krefmax); - c += ifactor * jfactor * kfactor * arr_src_safe(ii, jj, kk, comp); + 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; } } } @@ -117,15 +120,15 @@ namespace ablastr::coarsen::average * * 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] 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 + * \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 */ @@ -133,19 +136,40 @@ namespace ablastr::coarsen::average amrex::Real Interp ( amrex::Array4 const &arr_src, - amrex::GpuArray const &sf, - amrex::GpuArray const &sc, - amrex::GpuArray const &cr, + 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, idx_min; - amrex::GpuArray use_half; - amrex::Real crxyz_inv; + 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 - CalculateCoarseningData(idx_min, np, use_half, crxyz_inv, sf, sc, cr); + 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]); - return InterpWithCoarseningData(arr_src, idx_min, np, use_half, crxyz_inv, cr, i, j, k, comp); + return InterpWithCoarseningData( + arr_src,idx_min_src,np_src,use_half_weight,crx_cry_crz_inv,crse_ratio,i,j,k,comp); } /** diff --git a/Source/ablastr/coarsen/average.cpp b/Source/ablastr/coarsen/average.cpp index 2b91cf2261f..7ff3f233a42 100644 --- a/Source/ablastr/coarsen/average.cpp +++ b/Source/ablastr/coarsen/average.cpp @@ -32,25 +32,26 @@ namespace ablastr::coarsen::average ) { // Staggering of source fine MultiFab and destination coarse MultiFab - amrex::IntVect const stag_src = mf_src.boxArray().ixType().toIntVect(); - amrex::IntVect const stag_dst = mf_dst.boxArray().ixType().toIntVect(); + amrex::IntVect const ixtype_src = mf_src.boxArray().ixType().toIntVect(); + amrex::IntVect const ixtype_dst = mf_dst.boxArray().ixType().toIntVect(); // Auxiliary integer arrays (always 3D) - auto sf = amrex::GpuArray{0,0,0}; // staggering of source fine MultiFab - auto sc = amrex::GpuArray{0,0,0}; // staggering of destination coarse MultiFab - auto cr = amrex::GpuArray{1,1,1}; // coarsening ratio + auto stag_src = amrex::GpuArray{0,0,0}; // staggering of source fine MultiFab + auto stag_des = amrex::GpuArray{0,0,0}; // staggering of destination coarse MultiFab + auto crse_ratio = amrex::GpuArray{1,1,1}; // coarsening ratio for (int i=0; i idx_min, np; - amrex::GpuArray use_half; - amrex::Real crxyz_inv; + amrex::GpuArray idx_min_src, np_src; + amrex::GpuArray use_half_weight; + amrex::Real crx_cry_crz_inv; - CalculateCoarseningData(idx_min, np, use_half, crxyz_inv, sf, sc, cr); + CalculateCoarseningData( + idx_min_src, np_src, use_half_weight, crx_cry_crz_inv, stag_src, stag_des, crse_ratio); #ifdef AMREX_USE_OMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) @@ -64,7 +65,7 @@ namespace ablastr::coarsen::average 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, np, use_half, crxyz_inv, cr, i, j, k, n); + arr_src, idx_min_src, np_src, use_half_weight, crx_cry_crz_inv, crse_ratio, i, j, k, n); }); } } From ba0877df4bb3211b34b1f02153348d6ede7d730c Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Mon, 11 Dec 2023 11:00:47 -0800 Subject: [PATCH 19/26] Small comment --- Source/ablastr/coarsen/average.H | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Source/ablastr/coarsen/average.H b/Source/ablastr/coarsen/average.H index 746fa047e6e..1efa59aae55 100644 --- a/Source/ablastr/coarsen/average.H +++ b/Source/ablastr/coarsen/average.H @@ -168,6 +168,8 @@ namespace ablastr::coarsen::average 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); } From e481179e56e780184e0cacf6a1e4838b10292ab5 Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Mon, 11 Dec 2023 17:54:59 -0800 Subject: [PATCH 20/26] Rename variables and inline CalculateCoarseningData --- Source/ablastr/coarsen/average.H | 1 + Source/ablastr/coarsen/average.cpp | 20 ++++++++++---------- 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/Source/ablastr/coarsen/average.H b/Source/ablastr/coarsen/average.H index 1efa59aae55..4f7779470af 100644 --- a/Source/ablastr/coarsen/average.H +++ b/Source/ablastr/coarsen/average.H @@ -25,6 +25,7 @@ */ namespace ablastr::coarsen::average { + AMREX_FORCE_INLINE void CalculateCoarseningData ( amrex::GpuArray &idx_min_src, diff --git a/Source/ablastr/coarsen/average.cpp b/Source/ablastr/coarsen/average.cpp index 7ff3f233a42..f24a17ed2a1 100644 --- a/Source/ablastr/coarsen/average.cpp +++ b/Source/ablastr/coarsen/average.cpp @@ -32,18 +32,18 @@ namespace ablastr::coarsen::average ) { // Staggering of source fine MultiFab and destination coarse MultiFab - amrex::IntVect const ixtype_src = mf_src.boxArray().ixType().toIntVect(); - amrex::IntVect const ixtype_dst = mf_dst.boxArray().ixType().toIntVect(); + amrex::IntVect const stag_src = mf_src.boxArray().ixType().toIntVect(); + amrex::IntVect const stag_dst = mf_dst.boxArray().ixType().toIntVect(); // Auxiliary integer arrays (always 3D) - auto stag_src = amrex::GpuArray{0,0,0}; // staggering of source fine MultiFab - auto stag_des = amrex::GpuArray{0,0,0}; // staggering of destination coarse MultiFab - auto crse_ratio = amrex::GpuArray{1,1,1}; // coarsening ratio + auto sf = amrex::GpuArray{0,0,0}; // staggering of source fine MultiFab + auto sc = amrex::GpuArray{0,0,0}; // staggering of destination coarse MultiFab + auto cr = amrex::GpuArray{1,1,1}; // coarsening ratio for (int i=0; i idx_min_src, np_src; @@ -51,7 +51,7 @@ namespace ablastr::coarsen::average amrex::Real crx_cry_crz_inv; CalculateCoarseningData( - idx_min_src, np_src, use_half_weight, crx_cry_crz_inv, stag_src, stag_des, crse_ratio); + 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()) @@ -65,7 +65,7 @@ namespace ablastr::coarsen::average 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, crse_ratio, i, j, k, n); + arr_src, idx_min_src, np_src, use_half_weight, crx_cry_crz_inv, cr, i, j, k, n); }); } } From 8caede943675ac87d083523ba8be167094ab3873 Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Mon, 11 Dec 2023 18:04:13 -0800 Subject: [PATCH 21/26] Tab indent --- Source/Utils/check_interp_points_and_weights.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Utils/check_interp_points_and_weights.py b/Source/Utils/check_interp_points_and_weights.py index a61f5b625d9..c21f8f111a7 100644 --- a/Source/Utils/check_interp_points_and_weights.py +++ b/Source/Utils/check_interp_points_and_weights.py @@ -196,7 +196,7 @@ def refinement_points_and_weights( ii, 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] + ws += weights[ir] 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 ) ) From 6244867f1bbdc7f4b447029f795cd071b0d5a023 Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Wed, 20 Nov 2024 13:31:49 -0800 Subject: [PATCH 22/26] Script to compare versions of check_interp_points_and_weights.py --- .../check_interp_points_and_weights_dev.py | 215 ++++++++++++++++++ Source/Utils/check_ipaw_generate_output.sh | 13 ++ Source/Utils/cipaw.out | 187 +++++++++++++++ Source/Utils/cipaw_dev.out | 140 ++++++++++++ 4 files changed, 555 insertions(+) create mode 100644 Source/Utils/check_interp_points_and_weights_dev.py create mode 100644 Source/Utils/check_ipaw_generate_output.sh create mode 100644 Source/Utils/cipaw.out create mode 100644 Source/Utils/cipaw_dev.out 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..85057b431d7 --- /dev/null +++ b/Source/Utils/cipaw.out @@ -0,0 +1,187 @@ +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 * + ************************************************** + + Coarsening 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 i=0 by interpolating over the following points and weights: + (ii=0,w=0.500) (ii=1,w=0.500) + + Find value at i=1 by interpolating over the following points and weights: + (ii=2,w=0.500) (ii=3,w=0.500) + + Find value at i=2 by interpolating over the following points and weights: + (ii=4,w=0.500) (ii=5,w=0.500) + + Find value at i=3 by interpolating over the following points and weights: + (ii=6,w=0.500) (ii=7,w=0.500) + + Refinement for MR: check interpolation points and weights + --------------------------------------------------------- +Before ii_min=206 and ii_max=-200 +After ii_min=0 and ii_max=7 + + 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: + (i=0,w=1.000) + + Find value at ii=1 by interpolating over the following points and weights: + (i=0,w=1.000) + + Find value at ii=2 by interpolating over the following points and weights: + (i=1,w=1.000) + + Find value at ii=3 by interpolating over the following points and weights: + (i=1,w=1.000) + + Find value at ii=4 by interpolating over the following points and weights: + (i=2,w=1.000) + + Find value at ii=5 by interpolating over the following points and weights: + (i=2,w=1.000) + + Find value at ii=6 by interpolating over the following points and weights: + (i=3,w=1.000) + + Find value at ii=7 by interpolating over the following points and weights: + (i=3,w=1.000) + + ************************************************** + * Staggering of coarse grid: sc=0 cell-centered * + * Staggering of fine grid: sf=1 nodal * + ************************************************** + + Coarsening for MR: check interpolation points and weights + --------------------------------------------------------- + + Min and max index on coarse grid: imin=-1 imax=4 + Min and max index on fine grid: iimin=0 iimax=8 + + Find value at i=-1 by interpolating over the following points and weights: + (ii=-2,w=0.250) (ii=-1,w=0.500) (ii=0,w=0.250) + + Find value at i=0 by interpolating over the following points and weights: + (ii=0,w=0.250) (ii=1,w=0.500) (ii=2,w=0.250) + + Find value at i=1 by interpolating over the following points and weights: + (ii=2,w=0.250) (ii=3,w=0.500) (ii=4,w=0.250) + + Find value at i=2 by interpolating over the following points and weights: + (ii=4,w=0.250) (ii=5,w=0.500) (ii=6,w=0.250) + + Find value at i=3 by interpolating over the following points and weights: + (ii=6,w=0.250) (ii=7,w=0.500) (ii=8,w=0.250) + + Find value at i=4 by interpolating over the following points and weights: + (ii=8,w=0.250) (ii=9,w=0.500) (ii=10,w=0.250) + + 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 * + ************************************************** + + Coarsening for MR: check interpolation points and weights + --------------------------------------------------------- + + Min and max index on coarse grid: imin=0 imax=4 + Min and max index on fine grid: iimin=0 iimax=7 + + Find value at i=0 by interpolating over the following points and weights: + (ii=-1,w=0.500) (ii=0,w=0.500) + + Find value at i=1 by interpolating over the following points and weights: + (ii=1,w=0.500) (ii=2,w=0.500) + + Find value at i=2 by interpolating over the following points and weights: + (ii=3,w=0.500) (ii=4,w=0.500) + + Find value at i=3 by interpolating over the following points and weights: + (ii=5,w=0.500) (ii=6,w=0.500) + + Find value at i=4 by interpolating over the following points and weights: + (ii=7,w=0.500) (ii=8,w=0.500) + + 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 * + ************************************************** + + Coarsening for MR: check interpolation points and weights + --------------------------------------------------------- + + Min and max index on coarse grid: imin=0 imax=4 + Min and max index on fine grid: iimin=0 iimax=8 + + Find value at i=0 by interpolating over the following points and weights: + (ii=-1,w=0.250) (ii=0,w=0.500) (ii=1,w=0.250) + + Find value at i=1 by interpolating over the following points and weights: + (ii=1,w=0.250) (ii=2,w=0.500) (ii=3,w=0.250) + + Find value at i=2 by interpolating over the following points and weights: + (ii=3,w=0.250) (ii=4,w=0.500) (ii=5,w=0.250) + + Find value at i=3 by interpolating over the following points and weights: + (ii=5,w=0.250) (ii=6,w=0.500) (ii=7,w=0.250) + + Find value at i=4 by interpolating over the following points and weights: + (ii=7,w=0.250) (ii=8,w=0.500) (ii=9,w=0.250) + + Refinement for MR: check interpolation points and weights + --------------------------------------------------------- +Before ii_min=206 and ii_max=-200 +After ii_min=-1 and ii_max=7 + + 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: + (i=-1,w=1.000) (i=0,w=1.000) + + Find value at ii=0 by interpolating over the following points and weights: + (i=0,w=1.000) + + Find value at ii=1 by interpolating over the following points and weights: + (i=0,w=0.500) (i=1,w=0.500) + + Find value at ii=2 by interpolating over the following points and weights: + (i=1,w=1.000) + + Find value at ii=3 by interpolating over the following points and weights: + (i=1,w=0.500) (i=2,w=0.500) + + Find value at ii=4 by interpolating over the following points and weights: + (i=2,w=1.000) + + Find value at ii=5 by interpolating over the following points and weights: + (i=2,w=0.500) (i=3,w=0.500) + + Find value at ii=6 by interpolating over the following points and weights: + (i=3,w=1.000) + + Find value at ii=7 by interpolating over the following points and weights: + (i=3,w=1.000) (i=4,w=1.000) + + ERROR: sum of weights ws=2.500 should be cr=2 for i=0 + + ERROR: sum of weights ws=2.500 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) From 857b65a4d6a46ff88d6cdc201d71d4c8c380f170 Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Wed, 20 Nov 2024 13:57:13 -0800 Subject: [PATCH 23/26] Move some blocks around to better match current development branch --- .../Utils/check_interp_points_and_weights.py | 29 +++++++++++-------- Source/Utils/cipaw.out | 24 +++++++-------- 2 files changed, 29 insertions(+), 24 deletions(-) diff --git a/Source/Utils/check_interp_points_and_weights.py b/Source/Utils/check_interp_points_and_weights.py index 4f85e005208..28dac1c16ee 100644 --- a/Source/Utils/check_interp_points_and_weights.py +++ b/Source/Utils/check_interp_points_and_weights.py @@ -25,11 +25,23 @@ # Source/ablastr/coarsen/sample.(H/.cpp) # ------------------------------------------------------------------------------- +import sys import numpy as np -# Coarsening functions +# Fine grid limits (without ghost cells) +def coarsening_fine_grid_limits(sc, sf, cr): + if sf == 0: # cell-centered + iimin = 0 + iimax = 4 * cr - 1 + elif sf == 1: # nodal + iimin = 0 + iimax = 4 * cr + return [iimin, iimax] + + +# Coarse grid limits (without ghost cells) def coarsening_coarse_grid_limits(sc, sf, cr, ii_min, ii_max): i_range_start = (ii_min // cr) - 100 i_range_end = (ii_max // cr) + 100 @@ -46,15 +58,8 @@ def coarsening_coarse_grid_limits(sc, sf, cr, ii_min, ii_max): i_max = max(i_max, i) return [i_min, i_max] -def coarsening_fine_grid_limits(sc, sf, cr): - if sf == 0: # cell-centered - iimin = 0 - iimax = 4 * cr - 1 - elif sf == 1: # nodal - iimin = 0 - iimax = 4 * cr - return [iimin, iimax] +# Coarsening for MR: interpolation points and weights def coarsening_points_and_weights(i, sc, sf, cr): two_ii_start = -cr * sc + sf - 1 if two_ii_start % 2 == 0: @@ -167,9 +172,6 @@ def refinement_points_and_weights(ii, sc, sf, cr): print(" nodal *") print(" **************************************************") - print("\n Coarsening for MR: check interpolation points and weights") - print(" ---------------------------------------------------------") - iimin, iimax = coarsening_fine_grid_limits(sc, sf, cr) imin, imax = coarsening_coarse_grid_limits(sc, sf, cr, iimin, iimax) @@ -180,6 +182,9 @@ def refinement_points_and_weights(ii, sc, sf, cr): " Min and max index on fine grid: iimin={} iimax={}".format(iimin, iimax) ) + print("\n Coarsening for MR: check interpolation points and weights") + print(" ---------------------------------------------------------") + # Coarsening for MR: interpolation points and weights for i in range(imin, imax+1): # index on coarse grid numpts, idxmin, weights = coarsening_points_and_weights(i, sc, sf, cr) diff --git a/Source/Utils/cipaw.out b/Source/Utils/cipaw.out index 85057b431d7..a49c0038711 100644 --- a/Source/Utils/cipaw.out +++ b/Source/Utils/cipaw.out @@ -6,12 +6,12 @@ Coarsening ratio cr=2 * Staggering of fine grid: sf=0 cell-centered * ************************************************** - Coarsening 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 + Coarsening for MR: check interpolation points and weights + --------------------------------------------------------- + Find value at i=0 by interpolating over the following points and weights: (ii=0,w=0.500) (ii=1,w=0.500) @@ -61,12 +61,12 @@ After ii_min=0 and ii_max=7 * Staggering of fine grid: sf=1 nodal * ************************************************** - Coarsening for MR: check interpolation points and weights - --------------------------------------------------------- - 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: (ii=-2,w=0.250) (ii=-1,w=0.500) (ii=0,w=0.250) @@ -95,12 +95,12 @@ After ii_min=0 and ii_max=7 * Staggering of fine grid: sf=0 cell-centered * ************************************************** - Coarsening for MR: check interpolation points and weights - --------------------------------------------------------- - 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: (ii=-1,w=0.500) (ii=0,w=0.500) @@ -126,12 +126,12 @@ After ii_min=0 and ii_max=7 * Staggering of fine grid: sf=1 nodal * ************************************************** - Coarsening for MR: check interpolation points and weights - --------------------------------------------------------- - 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: (ii=-1,w=0.250) (ii=0,w=0.500) (ii=1,w=0.250) From d0ada5f5624033ca3be01b2cbe2608ab474ca65b Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Wed, 20 Nov 2024 15:31:03 -0800 Subject: [PATCH 24/26] Added comments, modified print statements, add brute force method --- .../Utils/check_interp_points_and_weights.py | 29 ++++--- Source/Utils/cipaw.out | 78 +++++++++---------- 2 files changed, 55 insertions(+), 52 deletions(-) diff --git a/Source/Utils/check_interp_points_and_weights.py b/Source/Utils/check_interp_points_and_weights.py index 28dac1c16ee..230b69bc204 100644 --- a/Source/Utils/check_interp_points_and_weights.py +++ b/Source/Utils/check_interp_points_and_weights.py @@ -29,20 +29,27 @@ import numpy as np +# Coarsening: map fine grid to coarse grid # Fine grid limits (without ghost cells) +# 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 = 4 * cr - 1 + ii_min = 0 + ii_max = 4 * cr - 1 elif sf == 1: # nodal - iimin = 0 - iimax = 4 * cr - return [iimin, iimax] + ii_min = 0 + ii_max = 4 * cr + return [ii_min, ii_max] # Coarse grid limits (without ghost cells) 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 @@ -80,7 +87,7 @@ def coarsening_points_and_weights(i, sc, sf, cr): return [num_ii_pts, ii_start, weights] -# Refinement functions +# Refinement: map coarse grid to fine grid def refinement_coarse_grid_limits(sc, sf, cr): i_min = 0 @@ -96,7 +103,7 @@ def refinement_fine_grid_limits(sc, sf, cr, i_min, i_max): ii_min = ii_range_end ii_max = ii_range_start - print("Before ii_min={} and ii_max={}".format(ii_min,ii_max)) + # 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) @@ -106,7 +113,7 @@ def refinement_fine_grid_limits(sc, sf, cr, i_min, i_max): if i_max >= i_start: ii_max = max(ii_max, ii) - print("After ii_min={} and ii_max={}".format(ii_min,ii_max)) + # print(" After ii_min={} and ii_max={}".format(ii_min,ii_max)) return [ii_min, ii_max] @@ -197,7 +204,7 @@ def refinement_points_and_weights(ii, sc, sf, cr): for ir in range(numpts): # interpolation points and weights ii = idxmin + ir wtotal += weights[ir] - print(" (ii={},w={:.3f})".format(ii, weights[ir]), end="") + print(" ({},{})".format(ii, weights[ir]), end="") if not (ir == numpts - 1): print(" ", end="") print() @@ -240,8 +247,8 @@ def refinement_points_and_weights(ii, sc, sf, cr): ) for ir in range(num_i_pts): # interpolation points and weights i = i_start + ir - print( ' (i={},w={:.3f})'.format( i, weights[ir] ), end="") - if not ( ir == num_i_pts-1 ): + print(" ({},{})".format(i, weights[ir]), end="") + if not (ir == num_i_pts - 1): print(" ", end="") print() diff --git a/Source/Utils/cipaw.out b/Source/Utils/cipaw.out index a49c0038711..2b92e8f6650 100644 --- a/Source/Utils/cipaw.out +++ b/Source/Utils/cipaw.out @@ -13,48 +13,46 @@ Coarsening ratio cr=2 --------------------------------------------------------- Find value at i=0 by interpolating over the following points and weights: - (ii=0,w=0.500) (ii=1,w=0.500) + (0,0.5) (1,0.5) Find value at i=1 by interpolating over the following points and weights: - (ii=2,w=0.500) (ii=3,w=0.500) + (2,0.5) (3,0.5) Find value at i=2 by interpolating over the following points and weights: - (ii=4,w=0.500) (ii=5,w=0.500) + (4,0.5) (5,0.5) Find value at i=3 by interpolating over the following points and weights: - (ii=6,w=0.500) (ii=7,w=0.500) + (6,0.5) (7,0.5) Refinement for MR: check interpolation points and weights --------------------------------------------------------- -Before ii_min=206 and ii_max=-200 -After ii_min=0 and ii_max=7 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: - (i=0,w=1.000) + (0,1.0) Find value at ii=1 by interpolating over the following points and weights: - (i=0,w=1.000) + (0,1.0) Find value at ii=2 by interpolating over the following points and weights: - (i=1,w=1.000) + (1,1.0) Find value at ii=3 by interpolating over the following points and weights: - (i=1,w=1.000) + (1,1.0) Find value at ii=4 by interpolating over the following points and weights: - (i=2,w=1.000) + (2,1.0) Find value at ii=5 by interpolating over the following points and weights: - (i=2,w=1.000) + (2,1.0) Find value at ii=6 by interpolating over the following points and weights: - (i=3,w=1.000) + (3,1.0) Find value at ii=7 by interpolating over the following points and weights: - (i=3,w=1.000) + (3,1.0) ************************************************** * Staggering of coarse grid: sc=0 cell-centered * @@ -68,22 +66,22 @@ After ii_min=0 and ii_max=7 --------------------------------------------------------- Find value at i=-1 by interpolating over the following points and weights: - (ii=-2,w=0.250) (ii=-1,w=0.500) (ii=0,w=0.250) + (-2,0.25) (-1,0.5) (0,0.25) Find value at i=0 by interpolating over the following points and weights: - (ii=0,w=0.250) (ii=1,w=0.500) (ii=2,w=0.250) + (0,0.25) (1,0.5) (2,0.25) Find value at i=1 by interpolating over the following points and weights: - (ii=2,w=0.250) (ii=3,w=0.500) (ii=4,w=0.250) + (2,0.25) (3,0.5) (4,0.25) Find value at i=2 by interpolating over the following points and weights: - (ii=4,w=0.250) (ii=5,w=0.500) (ii=6,w=0.250) + (4,0.25) (5,0.5) (6,0.25) Find value at i=3 by interpolating over the following points and weights: - (ii=6,w=0.250) (ii=7,w=0.500) (ii=8,w=0.250) + (6,0.25) (7,0.5) (8,0.25) Find value at i=4 by interpolating over the following points and weights: - (ii=8,w=0.250) (ii=9,w=0.500) (ii=10,w=0.250) + (8,0.25) (9,0.5) (10,0.25) Refinement for MR: check interpolation points and weights --------------------------------------------------------- @@ -102,19 +100,19 @@ After ii_min=0 and ii_max=7 --------------------------------------------------------- Find value at i=0 by interpolating over the following points and weights: - (ii=-1,w=0.500) (ii=0,w=0.500) + (-1,0.5) (0,0.5) Find value at i=1 by interpolating over the following points and weights: - (ii=1,w=0.500) (ii=2,w=0.500) + (1,0.5) (2,0.5) Find value at i=2 by interpolating over the following points and weights: - (ii=3,w=0.500) (ii=4,w=0.500) + (3,0.5) (4,0.5) Find value at i=3 by interpolating over the following points and weights: - (ii=5,w=0.500) (ii=6,w=0.500) + (5,0.5) (6,0.5) Find value at i=4 by interpolating over the following points and weights: - (ii=7,w=0.500) (ii=8,w=0.500) + (7,0.5) (8,0.5) Refinement for MR: check interpolation points and weights --------------------------------------------------------- @@ -133,54 +131,52 @@ After ii_min=0 and ii_max=7 --------------------------------------------------------- Find value at i=0 by interpolating over the following points and weights: - (ii=-1,w=0.250) (ii=0,w=0.500) (ii=1,w=0.250) + (-1,0.25) (0,0.5) (1,0.25) Find value at i=1 by interpolating over the following points and weights: - (ii=1,w=0.250) (ii=2,w=0.500) (ii=3,w=0.250) + (1,0.25) (2,0.5) (3,0.25) Find value at i=2 by interpolating over the following points and weights: - (ii=3,w=0.250) (ii=4,w=0.500) (ii=5,w=0.250) + (3,0.25) (4,0.5) (5,0.25) Find value at i=3 by interpolating over the following points and weights: - (ii=5,w=0.250) (ii=6,w=0.500) (ii=7,w=0.250) + (5,0.25) (6,0.5) (7,0.25) Find value at i=4 by interpolating over the following points and weights: - (ii=7,w=0.250) (ii=8,w=0.500) (ii=9,w=0.250) + (7,0.25) (8,0.5) (9,0.25) Refinement for MR: check interpolation points and weights --------------------------------------------------------- -Before ii_min=206 and ii_max=-200 -After ii_min=-1 and ii_max=7 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: - (i=-1,w=1.000) (i=0,w=1.000) + (-1,1.0) (0,1.0) Find value at ii=0 by interpolating over the following points and weights: - (i=0,w=1.000) + (0,1.0) Find value at ii=1 by interpolating over the following points and weights: - (i=0,w=0.500) (i=1,w=0.500) + (0,0.5) (1,0.5) Find value at ii=2 by interpolating over the following points and weights: - (i=1,w=1.000) + (1,1.0) Find value at ii=3 by interpolating over the following points and weights: - (i=1,w=0.500) (i=2,w=0.500) + (1,0.5) (2,0.5) Find value at ii=4 by interpolating over the following points and weights: - (i=2,w=1.000) + (2,1.0) Find value at ii=5 by interpolating over the following points and weights: - (i=2,w=0.500) (i=3,w=0.500) + (2,0.5) (3,0.5) Find value at ii=6 by interpolating over the following points and weights: - (i=3,w=1.000) + (3,1.0) Find value at ii=7 by interpolating over the following points and weights: - (i=3,w=1.000) (i=4,w=1.000) + (3,1.0) (4,1.0) ERROR: sum of weights ws=2.500 should be cr=2 for i=0 From 40311d42bcc06f91bea5333575d65b5868ce6ad1 Mon Sep 17 00:00:00 2001 From: Edward Basso Date: Wed, 20 Nov 2024 15:47:30 -0800 Subject: [PATCH 25/26] Add explicit iimin and iimax arguments to refinement_points_and_weights --- .../Utils/check_interp_points_and_weights.py | 43 ++++++++++--------- Source/Utils/cipaw.out | 4 +- 2 files changed, 24 insertions(+), 23 deletions(-) diff --git a/Source/Utils/check_interp_points_and_weights.py b/Source/Utils/check_interp_points_and_weights.py index 230b69bc204..769cd03414d 100644 --- a/Source/Utils/check_interp_points_and_weights.py +++ b/Source/Utils/check_interp_points_and_weights.py @@ -106,7 +106,7 @@ def refinement_fine_grid_limits(sc, sf, cr, i_min, i_max): # 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) + 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) @@ -118,21 +118,21 @@ def refinement_fine_grid_limits(sc, sf, cr, i_min, i_max): return [ii_min, ii_max] # Refinement for MR: interpolation points and weights -def refinement_points_and_weights(ii, sc, sf, cr): - i_start = 0 - num_i_pts = 0 +def refinement_points_and_weights(ii, sc, sf, cr, iimin, iimax): + idxmin = 0 + numpts = 0 if cr == 1: - num_i_pts = 1 - i_start = ii + numpts = 1 + idxmin = ii elif cr >= 2: if ii % cr == 0: - num_i_pts = (1 - sf) * (1 - sc) + sf * sc + numpts = (1 - sf) * (1 - sc) + sf * sc elif ii % cr != 0: - num_i_pts = (1 - sf) * (1 - sc) + 2 * sf * sc - i_start = (ii // cr) * (1 - sf) * (1 - sc) + (ii // cr) * sf * sc - weights = np.zeros(num_i_pts) - for ir in range(num_i_pts): - i = i_start + ir + 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) @@ -141,7 +141,8 @@ def refinement_points_and_weights(ii, sc, sf, cr): weights[ir] = (1 - sf) * (1 - sc) + ( (abs(cr - abs(ii - i * cr))) / (cr) ) * sf * sc - return [num_i_pts, i_start, weights] + return [numpts, idxmin, weights] + ## TODO Coarsening for IO: interpolation points and weights # def coarsening_points_and_weights_for_IO( i, sf, sc, cr ): @@ -230,7 +231,7 @@ def refinement_points_and_weights(ii, sc, sf, cr): 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) + imin, imax = refinement_coarse_grid_limits(sc, sf, cr) iimin, iimax = refinement_fine_grid_limits(sc, sf, cr, imin, imax) # Number of grid points @@ -238,8 +239,8 @@ def refinement_points_and_weights(ii, sc, sf, cr): print( ' Min and max index on fine grid: iimin={} iimax={}'.format( iimin, iimax ) ) # Refinement for MR: interpolation points and weights - 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 ) + 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 @@ -253,13 +254,13 @@ def refinement_points_and_weights(ii, sc, sf, cr): print() # Refinement for MR: check conservation properties - for i in range(imin, imax + 1): # index on coarse grid + for i in range(imin, imax + 1): # index on coarse grid ws = 0.0 - for ii in range(iimin, iimax + 1): # index on fine grid - num_i_pts, idxmin, weights = refinement_points_and_weights(ii, sc, sf, cr) + 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 + if j == i: # interpolation point matches point on coarse grid ws += weights[ir] if abs(ws - cr) > 1E-9: - print("\n ERROR: sum of weights ws={:.3f} should be cr={} for i={}".format(ws, cr, i)) + print("\n ERROR: sum of weights ws={} should be cr={} for i={}".format(ws, cr, i)) diff --git a/Source/Utils/cipaw.out b/Source/Utils/cipaw.out index 2b92e8f6650..f70f8af36c8 100644 --- a/Source/Utils/cipaw.out +++ b/Source/Utils/cipaw.out @@ -178,6 +178,6 @@ Coarsening ratio cr=2 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.500 should be cr=2 for i=0 + ERROR: sum of weights ws=2.5 should be cr=2 for i=0 - ERROR: sum of weights ws=2.500 should be cr=2 for i=3 + ERROR: sum of weights ws=2.5 should be cr=2 for i=3 From 0d4026f005685505b07a3820a24e220433b64c64 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 24 Feb 2025 18:47:16 +0000 Subject: [PATCH 26/26] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../Utils/check_interp_points_and_weights.py | 84 +++++++++++++------ 1 file changed, 58 insertions(+), 26 deletions(-) diff --git a/Source/Utils/check_interp_points_and_weights.py b/Source/Utils/check_interp_points_and_weights.py index 769cd03414d..1e6a62e911e 100644 --- a/Source/Utils/check_interp_points_and_weights.py +++ b/Source/Utils/check_interp_points_and_weights.py @@ -25,12 +25,12 @@ # 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) # The find grid limits are arbitrarily chosen here to act as an example def coarsening_fine_grid_limits(sc, sf, cr): @@ -47,6 +47,7 @@ def coarsening_fine_grid_limits(sc, sf, cr): 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 @@ -59,9 +60,9 @@ def coarsening_coarse_grid_limits_brute_force(sc, sf, cr, ii_min, ii_max): 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): + if ii_min <= ii_end: i_min = min(i_min, i) - if (ii_max >= ii_start): + if ii_max >= ii_start: i_max = max(i_max, i) return [i_min, i_max] @@ -75,7 +76,7 @@ def coarsening_points_and_weights(i, sc, sf, cr): 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): + for ir in range(1, num_ii_pts - 1): weights[ir] = 1.0 / cr else: ii_start = i * cr + (two_ii_start + 1) // 2 @@ -89,14 +90,16 @@ def coarsening_points_and_weights(i, sc, sf, cr): # 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 + 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)) @@ -105,8 +108,10 @@ def refinement_fine_grid_limits(sc, sf, cr, i_min, i_max): # 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) + 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) @@ -117,6 +122,7 @@ def refinement_fine_grid_limits(sc, sf, cr, i_min, i_max): return [ii_min, ii_max] + # Refinement for MR: interpolation points and weights def refinement_points_and_weights(ii, sc, sf, cr, iimin, iimax): idxmin = 0 @@ -181,7 +187,7 @@ def refinement_points_and_weights(ii, sc, sf, cr, iimin, iimax): print(" **************************************************") iimin, iimax = coarsening_fine_grid_limits(sc, sf, cr) - imin, imax = coarsening_coarse_grid_limits(sc, sf, cr, iimin, iimax) + 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) @@ -194,7 +200,7 @@ def refinement_points_and_weights(ii, sc, sf, cr, iimin, iimax): print(" ---------------------------------------------------------") # Coarsening for MR: interpolation points and weights - for i in range(imin, imax+1): # 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( @@ -209,38 +215,58 @@ def refinement_points_and_weights(ii, sc, sf, cr, iimin, iimax): 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)) + 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(iimin, iimax+1): # index on fine grid + for ii in range(iimin, iimax + 1): # index on fine grid ws = 0.0 - 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 + 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 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)) + 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)) + 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 ) ) + 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 (iimin, iimax+1): # index on fine grid - num_i_pts, i_start, weights = refinement_points_and_weights(ii, sc, sf, cr, iimin, iimax) + 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 @@ -257,10 +283,16 @@ def refinement_points_and_weights(ii, sc, sf, cr, iimin, iimax): for i in range(imin, imax + 1): # index on coarse grid ws = 0.0 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) + 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 abs(ws - cr) > 1E-9: - print("\n ERROR: sum of weights ws={} should be cr={} for i={}".format(ws, cr, i)) + if abs(ws - cr) > 1e-9: + print( + "\n ERROR: sum of weights ws={} should be cr={} for i={}".format( + ws, cr, i + ) + )