From f8911d12c3066ade965e73f712e31e9119cfbbbf Mon Sep 17 00:00:00 2001 From: Hua Tan Date: Sat, 17 Aug 2024 02:13:01 -0700 Subject: [PATCH] light modification --- src/vof/VolumeOfFluid.cpp | 908 +++++++++++++++++++------------------- 1 file changed, 454 insertions(+), 454 deletions(-) diff --git a/src/vof/VolumeOfFluid.cpp b/src/vof/VolumeOfFluid.cpp index 94115edf..9474c006 100644 --- a/src/vof/VolumeOfFluid.cpp +++ b/src/vof/VolumeOfFluid.cpp @@ -54,7 +54,7 @@ void range_add_value (VofRange & r, Real val) /** * range_update: * @r: a #VofRange. - * + * * Updates the fields of @r. */ void range_update (VofRange & r) @@ -62,12 +62,12 @@ void range_update (VofRange & r) if (r.n > 0) { if (r.sum2 - r.sum*r.sum/(Real) r.n >= 0.) r.stddev = sqrt ((r.sum2 - r.sum*r.sum/(Real) r.n) - /(Real) r.n); + /(Real) r.n); else r.stddev = 0.; r.mean = r.sum/(Real) r.n; } - else + else r.min = r.max = r.mean = r.stddev = 0.; } @@ -89,7 +89,7 @@ static void domain_range_reduce ( VofRange & s) { double in[5]; - double out[5] = { std::numeric_limits::max(), + double out[5] = { std::numeric_limits::max(), std::numeric_limits::lowest(), 0., 0., 0. }; MPI_Op op; @@ -119,13 +119,13 @@ VolumeOfFluid::VolumeOfFluid (incflo* a_incflo) : v_incflo(a_incflo) for (int lev = 0; lev <= finest_level; ++lev){ normal.emplace_back(v_incflo->grids[lev], v_incflo->dmap[lev], AMREX_SPACEDIM, v_incflo->nghost_state(), MFInfo(), v_incflo->Factory(lev)); alpha.emplace_back(v_incflo->grids[lev], v_incflo->dmap[lev], 1, v_incflo->nghost_state(), MFInfo(), v_incflo->Factory(lev)); - tag.emplace_back(v_incflo->grids[lev], v_incflo->dmap[lev], 1, v_incflo->nghost_state(), MFInfo(), v_incflo->Factory(lev)); - Array new_height={ + tag.emplace_back(v_incflo->grids[lev], v_incflo->dmap[lev], 1, v_incflo->nghost_state(), MFInfo(), v_incflo->Factory(lev)); + Array new_height={ MultiFab(v_incflo->grids[lev], v_incflo->dmap[lev], AMREX_SPACEDIM, v_incflo->nghost_state(), MFInfo(), v_incflo->Factory(lev)), MultiFab(v_incflo->grids[lev], v_incflo->dmap[lev], AMREX_SPACEDIM, v_incflo->nghost_state(), MFInfo(), v_incflo->Factory(lev)) }; - height.emplace_back(std::move(new_height)); - kappa.emplace_back(v_incflo->grids[lev], v_incflo->dmap[lev], 1, v_incflo->nghost_state(), MFInfo(), v_incflo->Factory(lev)); + height.emplace_back(std::move(new_height)); + kappa.emplace_back(v_incflo->grids[lev], v_incflo->dmap[lev], 1, v_incflo->nghost_state(), MFInfo(), v_incflo->Factory(lev)); } } @@ -820,7 +820,7 @@ bool interface_cell (int const i, int const j, int const k, static int half_height (Array cell, Array4 const & fv, int d, - Real & H, int & n, Array range) + Real & H, int & n, Array range) { int s = 0, dim=d/2; n = 0; @@ -828,12 +828,12 @@ static int half_height (Array cell, Array4 const & fv, int while (n < HMAX && !s) { Real f = fv (cell[0],cell[1],cell[2],0); if (!CELL_IS_FULL(f)) { /* interfacial cell */ - // if (f > EPS && f < 1. - EPS) { /* interfacial cell */ + // if (f > EPS && f < 1. - EPS) { /* interfacial cell */ //hit the boundary - if (cell[dim]range[1]) - return 2; - H += f; - n++; + if (cell[dim]range[1]) + return 2; + H += f; + n++; } else /* full or empty cell */ s = (f - 0.5)>0.? 1.: -1; @@ -846,36 +846,36 @@ static int half_height (Array cell, Array4 const & fv, int static void height_propagation (Array cell, int dim, Array4 const & fv, Array4 const & hght, Array range, Real orientation) -{ +{ for (int d = 1; d >= -1; d-=2, orientation = - orientation) { - Array neighbor=cell; + Array neighbor=cell; Real H = hght(cell[0],cell[1],cell[2],dim); - neighbor[dim]+=d; + neighbor[dim]+=d; bool interface = !CELL_IS_FULL(fv(neighbor[0],neighbor[1],neighbor[2],0));//false; - while (fabs (H) < DMAX - 1.&& !interface && - neighbor[dim]>=range[0]&& neighbor[dim]<=range[1]) { + while (fabs (H) < DMAX - 1.&& !interface && + neighbor[dim]>=range[0]&& neighbor[dim]<=range[1]) { H -= orientation; hght(neighbor[0],neighbor[1],neighbor[2],dim) = H; - auto fvol = fv(neighbor[0],neighbor[1],neighbor[2],0); + auto fvol = fv(neighbor[0],neighbor[1],neighbor[2],0); interface = !CELL_IS_FULL(fvol); neighbor[dim]+=d; } } } -void calculate_height(int i, int j, int k, int dim, Array4 const & vof, +void calculate_height(int i, int j, int k, int dim, Array4 const & vof, Array4 const & hb, Array4 const & ht, Array range) { Real H = vof(i,j,k,0); - Array cell={i,j,k}; - // top part of the column + Array cell={i,j,k}; + // top part of the column int nt, st = half_height (cell, vof, 2*dim, H, nt, range); if (!st) /* still an interfacial cell */ - return; - // bottom part of the column + return; + // bottom part of the column int nb, sb = half_height (cell, vof, 2*dim + 1, H, nb, range); if (!sb) /* still an interfacial cell */ - return; + return; if (sb != 2 && st != 2) { if (st*sb > 0) /* the column does not cross the interface */ return; @@ -897,7 +897,7 @@ void calculate_height(int i, int j, int k, int dim, Array4 const & v } } -static Array4 const * boundary_hit (int i,int j,int k, int d, Array4 const & hb, +static Array4 const * boundary_hit (int i,int j,int k, int d, Array4 const & hb, Array4 const & ht) { if (hb(i,j,k,d)!= VOF_NODATA && hb(i,j,k,d)> BOUNDARY_HIT/2.) @@ -916,32 +916,32 @@ static void height_propagation_from_boundary (Array cell, int dim, int cell[dim]+=(d % 2 ? 1 : -1); Real H0=hght(cell[0],cell[1],cell[2],dim); while ( H0!=VOF_NODATA && H0 > BOUNDARY_HIT/2. && - cell[dim]>=range[0]&&cell[dim]<=range[1]) { + cell[dim]>=range[0]&&cell[dim]<=range[1]) { H += orientation; hght(cell[0],cell[1],cell[2],dim) = H; cell[dim]+=(d % 2 ? 1 : -1); - H0=hght(cell[0],cell[1],cell[2],dim); + H0=hght(cell[0],cell[1],cell[2],dim); } /* propagate to non-interfacial cells up to DMAX */ - auto fvol = fv(cell[0],cell[1],cell[2],0); + auto fvol = fv(cell[0],cell[1],cell[2],0); bool interface = !CELL_IS_FULL(fvol); - while (fabs (H) < DMAX - 1. && !interface && - cell[dim]>=range[0]&&cell[dim]<=range[1]) { + while (fabs (H) < DMAX - 1. && !interface && + cell[dim]>=range[0]&&cell[dim]<=range[1]) { H += orientation; hght(cell[0],cell[1],cell[2],dim) = H; cell[dim]+=(d % 2 ? 1 : -1); } } -Array4 const * closest_height (int i,int j,int k, int d, Array4 const & hb, +Array4 const * closest_height (int i,int j,int k, int d, Array4 const & hb, Array4 const & ht, Real * orientation) { Array4 const * hv = nullptr; Real o = 0.; if (hb(i,j,k,d)!=VOF_NODATA) { hv = &hb; o = 1.; - if (ht(i,j,k,d)!=VOF_NODATA && - fabs (ht(i,j,k,d)) < fabs (hb(i,j,k,d))) { + if (ht(i,j,k,d)!=VOF_NODATA && + fabs (ht(i,j,k,d)) < fabs (hb(i,j,k,d))) { hv = & ht; o = -1.; } } @@ -954,8 +954,8 @@ Array4 const * closest_height (int i,int j,int k, int d, Array4 const * h , int d, Real * x) +static Real neighboring_column (int i,int j,int k, int c, + Array4 const * h , int d, Real * x) { Array neighbor={i,j,k}; neighbor[d/2]+=d%2?-1:1; @@ -966,14 +966,14 @@ static Real neighboring_column (int i,int j,int k, int c, } return VOF_NODATA; } -/* +/* The function is similar to neighboring_column(). The difference is that neighboring_column_corner() returns height @h of the neighboring column in direction @(d[0], d[1]). kind of corner neighbors */ -static Real neighboring_column_corner (int i,int j,int k, int c, - Array4 const * h, int * d, Real (*x)[2]) +static Real neighboring_column_corner (int i,int j,int k, int c, + Array4 const * h, int * d, Real (*x)[2]) { Array neighbor={i,j,k}; neighbor[d[0]/2]+=d[0]%2?-1:1; @@ -981,23 +981,23 @@ static Real neighboring_column_corner (int i,int j,int k, int c, Real height=(*h)(neighbor[0],neighbor[1],neighbor[2],c); if (height!=VOF_NODATA) { (*x)[0] = d[0] % 2 ? -1. : 1.; - (*x)[1] = d[1] % 2 ? -1. : 1.; - return height; + (*x)[1] = d[1] % 2 ? -1. : 1.; + return height; } - else + else return VOF_NODATA; } -static bool height_normal (int i,int j,int k, Array4 const & hb, +static bool height_normal (int i,int j,int k, Array4 const & hb, Array4 const & ht, XDim3 & m ) { Real slope = VOF_NODATA; static int oc[3][2] = { { 1, 2 }, { 0, 2 }, { 0, 1 } }; - for (int d = 0; d < AMREX_SPACEDIM; d++){ + for (int d = 0; d < AMREX_SPACEDIM; d++){ Real orientation; - Array4 const * hv = closest_height (i,j,k,d,hb,ht,&orientation); - if (hv != nullptr && fabs ((*hv)(i,j,k,d) <= 1.)) { + Array4 const * hv = closest_height (i,j,k,d,hb,ht,&orientation); + if (hv != nullptr && fabs ((*hv)(i,j,k,d) <= 1.)) { Real H = (*hv)(i,j,k,d); Real x[2], h[2][2], hd[2]; for (int nd = 0; nd < 2; nd++) { @@ -1010,7 +1010,7 @@ static bool height_normal (int i,int j,int k, Array4 const & hb, x[1] = - x[1]; Real det = x[0]*x[1]*(x[0] - x[1]), a = x[1]*(h[nd][0] - H), b = x[0]*(h[nd][1] - H); hd[nd] = (x[0]*b - x[1]*a)/det; - } + } if (h[0][0] == VOF_NODATA || h[0][1] == VOF_NODATA || h[1][0] == VOF_NODATA || h[1][1] == VOF_NODATA) continue; @@ -1019,13 +1019,13 @@ static bool height_normal (int i,int j,int k, Array4 const & hb, (&m.x)[d] = orientation; (&m.x)[oc[d][0]] = - hd[0]; (&m.x)[oc[d][1]] = - hd[1]; - } - - } + } + + } } //Print()<<"-------slope---"< dx, - Array4 const & hb, - Array4 const & ht, - Real & kappa) +bool curvature_along_direction (int i,int j,int k, int d, GpuArray dx, + Array4 const & hb, + Array4 const & ht, + Real & kappa) { Real x[9], h[9]; - Real orientation; + Real orientation; static int oc[3][2] = { { 1, 2 }, { 0, 2 }, { 0, 1 } }; Array4 const * hv = closest_height (i,j,k,d,hb,ht,&orientation); if (!hv) { - bool loop=true; + bool loop=true; /* no data for either directions, look four neighbors to collect potential interface positions */ - for (int nd = 0; nd < 2 && loop; nd++) - for (int ndd = 0; ndd < 2; ndd++) { + for (int nd = 0; nd < 2 && loop; nd++) + for (int ndd = 0; ndd < 2; ndd++) { Array neighbor={i,j,k}; - neighbor[oc[d][nd]]+=ndd%2?-1:1; + neighbor[oc[d][nd]]+=ndd%2?-1:1; hv = closest_height (neighbor[0],neighbor[1],neighbor[2],d,hb,ht,&orientation); - if (hv){ - loop = false; - break; - } - } - if (!hv) /* give up */ + if (hv){ + loop = false; + break; + } + } + if (!hv) /* give up */ return false; } else if (fabs((*hv)(i,j,k,d))>1.) - return false; + return false; int n=0; for (int nd = 0; nd < 2; nd++) { h[n] = neighboring_column (i, j, k, d, hv, 2*oc[d][nd], &x[n]); @@ -1146,7 +1146,7 @@ bool curvature_along_direction (int i,int j,int k, int d, GpuArray dx, - Array4 const & hb, - Array4 const & ht, - Real & kappa, Vector &interface) +bool curvature_along_direction_new (int i,int j,int k, int d, GpuArray dx, + Array4 const & hb, + Array4 const & ht, + Real & kappa, Vector &interface) { Real x[9], h[9]; - Real orientation; + Real orientation; static int oc[3][2] = { { 1, 2 }, { 0, 2 }, { 0, 1 } }; Array4 const * hv = closest_height (i,j,k,d,hb,ht,&orientation); if (!hv) { - bool loop=true; + bool loop=true; /* no data for either directions, look four neighbors to collect potential interface positions */ - for (int nd = 0; nd < 2 && loop; nd++) - for (int ndd = 0; ndd < 2; ndd++) { + for (int nd = 0; nd < 2 && loop; nd++) + for (int ndd = 0; ndd < 2; ndd++) { Array neighbor={i,j,k}; - neighbor[oc[d][nd]]+=ndd%2?-1:1; + neighbor[oc[d][nd]]+=ndd%2?-1:1; hv = closest_height (neighbor[0],neighbor[1],neighbor[2],d,hb,ht,&orientation); - if (hv){ - loop = false; - break; - } - } - if (!hv) /* give up */ + if (hv){ + loop = false; + break; + } + } + if (!hv) /* give up */ return false; } else if (fabs((*hv)(i,j,k,d))>1.) - return false; + return false; int n=0; for (int nd = 0; nd < 2; nd++) { h[n] = neighboring_column (i, j, k, d, hv, 2*oc[d][nd], &x[n]); @@ -1220,15 +1220,15 @@ bool curvature_along_direction_new (int i,int j,int k, int d, GpuArray &c) @@ -1259,9 +1259,9 @@ static void orientation (VofVector m, Array &c) for (i = 0; i < AMREX_SPACEDIM - 1; i++) for (j = 0; j < AMREX_SPACEDIM - 1 - i; j++) if (fabs (m[c[j + 1]]) > fabs (m[c[j]])) { - int tmp = c[j]; - c[j] = c[j + 1]; - c[j + 1] = tmp; + int tmp = c[j]; + c[j] = c[j + 1]; + c[j + 1] = tmp; } } @@ -1277,7 +1277,7 @@ static int independent_positions (Vector &interface) for (i = 0; i < j && !depends; i++) { Real d2 = 0.; for (int c = 0; c < AMREX_SPACEDIM; c++) - d2 += (interface[i][c] - interface[j][c])*(interface[i][c] - interface[j][c]); + d2 += (interface[i][c] - interface[j][c])*(interface[i][c] - interface[j][c]); depends = d2 < 0.5*0.5; } ni += !depends; @@ -1306,43 +1306,43 @@ static int independent_positions (Vector &interface) * contained in @(i,j,k), or %VOF_NODATA if the HF method could not * compute a consistent curvature. */ -Real height_curvature_combined (int i,int j,int k, GpuArray dx, - Array4 const & hb, +Real height_curvature_combined (int i,int j,int k, GpuArray dx, + Array4 const & hb, Array4 const & ht, - Array4 const & mv, - Array4 const & alpha) + Array4 const & mv, + Array4 const & alpha) { - VofVector m; + VofVector m; Array try_dir; for (int d = 0; d < AMREX_SPACEDIM; d++) m[d] = mv(i,j,k,d); - orientation (m, try_dir); /* sort directions according to normal */ + orientation (m, try_dir); /* sort directions according to normal */ Real kappa = 0.; Vector interface; for (int d = 0; d < AMREX_SPACEDIM; d++) /* try each direction */ if (curvature_along_direction_new (i, j, k, try_dir[d], dx, hb, ht, kappa, interface)) - return kappa; + return kappa; /* Could not compute curvature from the simple algorithm along any direction: * Try parabola fitting of the collected interface positions */ if (independent_positions (interface) < 3*(AMREX_SPACEDIM - 1)) - return VOF_NODATA; + return VOF_NODATA; ParabolaFit fit; XDim3 mx={AMREX_D_DECL(m[0],m[1],m[2])},p; - + Real area=plane_area_center (mx, alpha(i,j,k,0),p); - //shift the coordinates of the center of the interfacial segment - //by using the center of the cube. plane_area_center() gives the + //shift the coordinates of the center of the interfacial segment + //by using the center of the cube. plane_area_center() gives the //coordinates of area center with the coordinate origin as (0.,0.,0.) - //After shifting, the origin becomes cell center. + //After shifting, the origin becomes cell center. for (int c = 0; c < AMREX_SPACEDIM; c++) (&p.x)[c] -= 0.5; // initialize the parameters for parabola fit parabola_fit_init (fit, p, mx); - + ////#if AMREX_SPACEDIM==2 //// parabola_fit_add (&fit, &fc.x, PARABOLA_FIT_CENTER_WEIGHT); ////#elif !PARABOLA_SIMPLER @@ -1357,11 +1357,11 @@ Real height_curvature_combined (int i,int j,int k, GpuArray dx, +Real curvature_fit (int i,int j,int k, GpuArray dx, Array4 const & vof, - Array4 const & mv, - Array4 const & alpha) + Array4 const & mv, + Array4 const & alpha) { - VofVector m; + VofVector m; for (int d = 0; d < AMREX_SPACEDIM; d++) - m[d] = mv(i,j,k,d); + m[d] = mv(i,j,k,d); ParabolaFit fit; XDim3 mx={AMREX_D_DECL(m[0],m[1],m[2])},p; - + Real area=plane_area_center (mx, alpha(i,j,k,0),p); - //shift the coordinates of the center of the interfacial segment - //by using the center of the cube. plane_area_center() gives the + //shift the coordinates of the center of the interfacial segment + //by using the center of the cube. plane_area_center() gives the //coordinates of area center with the coordinate origin as (0.,0.,0.) - //After shifting, the origin becomes cell center. + //After shifting, the origin becomes cell center. for (int c = 0; c < AMREX_SPACEDIM; c++) (&p.x)[c] -= 0.5; // initialize the parameters for parabola fit parabola_fit_init (fit, p, mx); - // add the center of the segment with the area of the segment as weight + // add the center of the segment with the area of the segment as weight parabola_fit_add (fit, {AMREX_D_DECL(p.x,p.y,p.z)}, area); int di=0,dj=0,dk=0; #if AMREX_SPACEDIM==3 @@ -1408,17 +1408,17 @@ Real curvature_fit (int i,int j,int k, GpuArray dx, #endif for (dj = -2; dj <= 2; dj++) for (di = -2; di <= 2; di++) - if (di != 0|| dj != 0|| dk != 0) { - int ni=i+di,nj=j+dj,nk=k+dk; + if (di != 0|| dj != 0|| dk != 0) { + int ni=i+di,nj=j+dj,nk=k+dk; Real fvol=vof(ni,nj,nk,0); - if (!CELL_IS_FULL(fvol)){ - mx={AMREX_D_DECL(mv(ni,nj,nk,0),mv(ni,nj,nk,1),mv(ni,nj,nk,2))}; - area=plane_area_center (mx, alpha(ni,nj,nk,0),p); - for (int c = 0; c < AMREX_SPACEDIM; c++) - (&p.x)[c] += c==0?di:c==1?dj:dk - 0.5; - parabola_fit_add (fit, {p.x,p.y,p.z}, area); - } - } + if (!CELL_IS_FULL(fvol)){ + mx={AMREX_D_DECL(mv(ni,nj,nk,0),mv(ni,nj,nk,1),mv(ni,nj,nk,2))}; + area=plane_area_center (mx, alpha(ni,nj,nk,0),p); + for (int c = 0; c < AMREX_SPACEDIM; c++) + (&p.x)[c] += c==0?di:c==1?dj:dk - 0.5; + parabola_fit_add (fit, {p.x,p.y,p.z}, area); + } + } parabola_fit_solve (fit); Real kappa = parabola_fit_curvature (fit, 2.)/dx[0]; # if PARABOLA_SIMPLER @@ -1426,11 +1426,11 @@ Real curvature_fit (int i,int j,int k, GpuArray dx, # else int nn=6; # endif - for (int c = 0; c < nn; c++) + for (int c = 0; c < nn; c++) delete[] fit.M[c]; // Delete each row delete[] fit.M; // Delete the array of pointers - return kappa; -} + return kappa; +} ////////////////////////////////////////////////////////////////////////////////////////////////// /////// /////// Update VOF properties including height values and normal direction @@ -1444,146 +1444,146 @@ VolumeOfFluid::tracer_vof_update (int lev, MultiFab & vof_mf, Array auto const& dx = geom.CellSizeArray(); auto const& problo = geom.ProbLoArray(); auto const& probhi = geom.ProbHiArray(); -/////////////////////////////////////////////////// -// update height using vof field +/////////////////////////////////////////////////// +// update height using vof field /////////////////////////////////////////////////// for (int dim = 0; dim < AMREX_SPACEDIM; dim++){ - + height[0].setVal(VOF_NODATA,dim,1); - height[1].setVal(VOF_NODATA,dim,1); - //fix me: have not thought of a way to deal with the MFIter with tiling - //an option is to use similar way as MPI's implementation. - for (MFIter mfi(vof_mf); mfi.isValid(); ++mfi) { + height[1].setVal(VOF_NODATA,dim,1); + //fix me: have not thought of a way to deal with the MFIter with tiling + //an option is to use similar way as MPI's implementation. + for (MFIter mfi(vof_mf); mfi.isValid(); ++mfi) { Box const& bx = mfi.validbox(); - Array range ={bx.smallEnd()[dim], bx.bigEnd()[dim]}; + Array range ={bx.smallEnd()[dim], bx.bigEnd()[dim]}; Array4 const& vof_arr = vof_mf.const_array(mfi); Array4 const& hb_arr = height[0].array(mfi); Array4 const& ht_arr = height[1].array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - auto fvol = vof_arr(i,j,k,0); - if (!CELL_IS_FULL(fvol)){ - calculate_height(i, j, k, dim, vof_arr, hb_arr, ht_arr, range); - }// end if + { + auto fvol = vof_arr(i,j,k,0); + if (!CELL_IS_FULL(fvol)){ + calculate_height(i, j, k, dim, vof_arr, hb_arr, ht_arr, range); + }// end if }); //end ParallelFor } //end MFIter - //fix me: temporary solution for MPI boundaries - height[0].FillBoundary(geom.periodicity()); - height[1].FillBoundary(geom.periodicity()); - + //fix me: temporary solution for MPI boundaries + height[0].FillBoundary(geom.periodicity()); + height[1].FillBoundary(geom.periodicity()); + //deal with the situation where interface goes across the MPI or periodic boundaries. if(1){ - for (MFIter mfi(vof_mf); mfi.isValid(); ++mfi) { /*fix me: no titling*/ + for (MFIter mfi(vof_mf); mfi.isValid(); ++mfi) { /*fix me: no titling*/ Box const& bx = mfi.validbox(); - Array face_min_max; + Array face_min_max; Array4 const& hb_arr = height[0].array(mfi); - Array4 const& ht_arr = height[1].array(mfi); - Array4 const& vof_arr = vof_mf.const_array(mfi); + Array4 const& ht_arr = height[1].array(mfi); + Array4 const& vof_arr = vof_mf.const_array(mfi); //search the cells on each boundary of the validbox - //we do it by creating a new indexing space (i.e., bbx) with a constant - //value for one coordinate direction. i.e., for +X face of the box, we can - // set i=imax and just vary j and k index. - Array range = {bx.smallEnd()[dim],bx.bigEnd()[dim]}; + //we do it by creating a new indexing space (i.e., bbx) with a constant + //value for one coordinate direction. i.e., for +X face of the box, we can + // set i=imax and just vary j and k index. + Array range = {bx.smallEnd()[dim],bx.bigEnd()[dim]}; auto ijk_min= bx.smallEnd(); auto ijk_max= bx.bigEnd(); //only loop through cells on two faces in the axis (defined by 'dim') - for (int nn = 0; nn < 2; nn++){ + for (int nn = 0; nn < 2; nn++){ //Note: we use the notation of Gerris for the direction of the Box (i.e.,FttDirection) -// FACE direction = 0,1,2,3,4,5 in 3D +// FACE direction = 0,1,2,3,4,5 in 3D // X+ (Right):0, X- (Left):1, Y+ (Top): 2, Y- (Bottom): 3, Z+ (Front): 4, Z- (Back):5 // direction%2=0 means the positive direction of a given axis direction (i.e.,int direction/2) -// direction%2=1 means the negative direction of a given axis direction (i.e.,int direction/2) +// direction%2=1 means the negative direction of a given axis direction (i.e.,int direction/2) // Axis direction = 0 (X-axis), 1(Y-axis), 2(Z-axis) // therefore, 'nn=0' here means the positive direction. - ijk_min[dim]= range[nn?0:1]; - ijk_max[dim]= range[nn?0:1]; + ijk_min[dim]= range[nn?0:1]; + ijk_max[dim]= range[nn?0:1]; Box bbx(ijk_min, ijk_max); -// loop through the cells on the face of the box ('bbx') - ParallelFor(bbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Array cell={i,j,k}, ghost=cell; - ghost[dim]+=nn%2?-1:1; - Array4 const * h=boundary_hit (i,j,k, dim, hb_arr,ht_arr); - /*if (i==7 && j==0 && k==5){ - AllPrint()<<"test_height_function "<<"hb "< const *hn=boundary_hit (ghost[0],ghost[1],ghost[2], dim, hb_arr,ht_arr); - if(h==hn){ - // the column crosses the interface - // propagate column height correction from one side (or PE) to the other - Real orientation = (nn%2 ? -1:1)*(h == &hb_arr ? 1 : -1); - Real h_ghost=(*h)(ghost[0],ghost[1],ghost[2],dim); +// loop through the cells on the face of the box ('bbx') + ParallelFor(bbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Array cell={i,j,k}, ghost=cell; + ghost[dim]+=nn%2?-1:1; + Array4 const * h=boundary_hit (i,j,k, dim, hb_arr,ht_arr); + /*if (i==7 && j==0 && k==5){ + AllPrint()<<"test_height_function "<<"hb "< const *hn=boundary_hit (ghost[0],ghost[1],ghost[2], dim, hb_arr,ht_arr); + if(h==hn){ + // the column crosses the interface + // propagate column height correction from one side (or PE) to the other + Real orientation = (nn%2 ? -1:1)*(h == &hb_arr ? 1 : -1); + Real h_ghost=(*h)(ghost[0],ghost[1],ghost[2],dim); Real Hn = h_ghost + 0.5 + (orientation - 1.)/2. - 2.*BOUNDARY_HIT; (*h)(i,j,k,dim) += Hn; - height_propagation_from_boundary (cell, dim, 2*dim+nn, vof_arr, *h, range, h == &hb_arr ? 1 : -1); - } - else{ - // the column does not cross the interface - Real hgh=(*h)(cell[0],cell[1],cell[2],dim); - while (!CELL_IS_BOUNDARY(cell,bx.smallEnd(),bx.bigEnd()) && - hgh!= VOF_NODATA && hgh> BOUNDARY_HIT/2.) { - (*h)(cell[0],cell[1],cell[2],dim) = VOF_NODATA; - cell[dim]+=nn%2?1:-1; - } - } - } - else{ + height_propagation_from_boundary (cell, dim, 2*dim+nn, vof_arr, *h, range, h == &hb_arr ? 1 : -1); + } + else{ + // the column does not cross the interface + Real hgh=(*h)(cell[0],cell[1],cell[2],dim); + while (!CELL_IS_BOUNDARY(cell,bx.smallEnd(),bx.bigEnd()) && + hgh!= VOF_NODATA && hgh> BOUNDARY_HIT/2.) { + (*h)(cell[0],cell[1],cell[2],dim) = VOF_NODATA; + cell[dim]+=nn%2?1:-1; + } + } + } + else{ // column did not hit a boundary, propagate height across PE boundary */ if (hb_arr(ghost[0],ghost[1],ghost[2],dim)!= VOF_NODATA) - height_propagation (ghost, dim, vof_arr, hb_arr, range, 1.); + height_propagation (ghost, dim, vof_arr, hb_arr, range, 1.); if (ht_arr(ghost[0],ghost[1],ghost[2],dim)!= VOF_NODATA) - height_propagation (ghost, dim, vof_arr, ht_arr, range, -1.); - } - //Print()<<"face_loop "<<"i "< const& hb_arr = height[0].array(mfi); - Array4 const& ht_arr = height[1].array(mfi); + Array4 const& ht_arr = height[1].array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - if (hb_arr(i,j,k,dim)!= VOF_NODATA && hb_arr(i,j,k,dim)> BOUNDARY_HIT/2) - hb_arr(i,j,k,dim)= VOF_NODATA; - if (ht_arr(i,j,k,dim)!= VOF_NODATA && ht_arr(i,j,k,dim)> BOUNDARY_HIT/2) - ht_arr(i,j,k,dim)= VOF_NODATA; - }); - } // end MFIter - //fix me: temporary solution for MPI boundaries - height[0].FillBoundary(geom.periodicity()); - height[1].FillBoundary(geom.periodicity()); + if (hb_arr(i,j,k,dim)!= VOF_NODATA && hb_arr(i,j,k,dim)> BOUNDARY_HIT/2) + hb_arr(i,j,k,dim)= VOF_NODATA; + if (ht_arr(i,j,k,dim)!= VOF_NODATA && ht_arr(i,j,k,dim)> BOUNDARY_HIT/2) + ht_arr(i,j,k,dim)= VOF_NODATA; + }); + } // end MFIter + //fix me: temporary solution for MPI boundaries + height[0].FillBoundary(geom.periodicity()); + height[1].FillBoundary(geom.periodicity()); }//end for dim -///////////////////////////////////////////////////////////////////////////////////////// +///////////////////////////////////////////////////////////////////////////////////////// //update the normal and alpha -///////////////////////////////////////////////////////////////////////////////////////// +///////////////////////////////////////////////////////////////////////////////////////// for (MFIter mfi(vof_mf,TilingIfNotGPU()); mfi.isValid(); ++mfi) { Box const& bx = mfi.tilebox(); Array4 const& vof_arr = vof_mf.const_array(mfi); Array4 const& mv = normal[lev].array(mfi); Array4 const& al = alpha[lev].array(mfi); Array4 const& hb_arr = height[0].const_array(mfi); - Array4 const& ht_arr = height[1].const_array(mfi); + Array4 const& ht_arr = height[1].const_array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { XDim3 m={0.,0.,0.}; - auto fvol = vof_arr(i,j,k,0); + auto fvol = vof_arr(i,j,k,0); THRESHOLD(fvol); if (i==5&&j==6&&k==8){ - int dddd; - Print()<<"------------"<<"\n"; - } - if (!height_normal (i,j,k, hb_arr, ht_arr, m)){ -// if(1){ + int dddd; + Print()<<"------------"<<"\n"; + } + if (!height_normal (i,j,k, hb_arr, ht_arr, m)){ +// if(1){ if (!interface_cell (i,j,k, vof_arr, fvol)) { AMREX_D_TERM(mv(i,j,k,0) = Real(0.);, mv(i,j,k,1) = Real(0.);, @@ -1594,11 +1594,11 @@ if(1){ AMREX_D_PICK( ,Real f[3][3];, Real f[3][3][3];) stencil (i,j,k, vof_arr, f); mycs (f, &m.x); - } - } + } + } Real n = 0.; for (int d = 0; d < AMREX_SPACEDIM; d++) - n += fabs ((&m.x)[d]); + n += fabs ((&m.x)[d]); if (n > 0.) for (int d = 0; d < AMREX_SPACEDIM; d++) mv(i,j,k,d)= (&m.x)[d]/n; @@ -1635,19 +1635,19 @@ if(1){ /////////////////////////////////////////////////////////////////////////////////////////////// void VolumeOfFluid::curvature_calculation (int lev, MultiFab & vof_mf, Array & height, MultiFab & kappa) -{ +{ Geometry const& geom =v_incflo->geom[lev]; auto const& dx = geom.CellSizeArray(); auto const& problo = geom.ProbLoArray(); - auto const& probhi = geom.ProbHiArray(); + auto const& probhi = geom.ProbHiArray(); MultiFab n_max(v_incflo->grids[lev], v_incflo->dmap[lev], 1, v_incflo->nghost_state(), - MFInfo(), v_incflo->Factory(lev)); - + MFInfo(), v_incflo->Factory(lev)); + kappa.setVal(VOF_NODATA); - n_max.setVal(-1.0); -// use height function method to calculate curvature - for (int dim = 0; dim < AMREX_SPACEDIM; dim++){ + n_max.setVal(-1.0); +// use height function method to calculate curvature + for (int dim = 0; dim < AMREX_SPACEDIM; dim++){ for (MFIter mfi(vof_mf,TilingIfNotGPU()); mfi.isValid(); ++mfi) { Box const& bx = mfi.tilebox(); Array4 const& vof_arr = vof_mf.const_array(mfi); @@ -1655,42 +1655,42 @@ VolumeOfFluid::curvature_calculation (int lev, MultiFab & vof_mf, Array const& hb_arr = height[0].const_array(mfi); Array4 const& ht_arr = height[1].const_array(mfi); Array4 const& kappa_arr = kappa.array(mfi); - Array4 const& nmax_arr = n_max.array(mfi); + Array4 const& nmax_arr = n_max.array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - auto fvol = vof_arr(i,j,k,0); - Real kappa0; - if (!CELL_IS_FULL(fvol)){ + { + auto fvol = vof_arr(i,j,k,0); + Real kappa0; + if (!CELL_IS_FULL(fvol)){ /* if (i==44&&j==35&&k==31){ - int dddd; - Print()<<"------------"<<"\n"; - }*/ + int dddd; + Print()<<"------------"<<"\n"; + }*/ if (curvature_along_direction(i, j, k, dim, dx, hb_arr, ht_arr, kappa0)) { if (fabs (mv(i,j,k,dim)) > nmax_arr(i,j,k,0)) { kappa_arr(i,j,k,0) = kappa0; nmax_arr(i,j,k,0) = fabs (mv(i,j,k,dim)); } - //propagate the curvature + //propagate the curvature Real orientation; - Array4 const * hv = closest_height (i,j,k,dim,hb_arr,ht_arr,&orientation); + Array4 const * hv = closest_height (i,j,k,dim,hb_arr,ht_arr,&orientation); for (int d = 0; d <= 1; d++) { - Array neighbor={i,j,k}; + Array neighbor={i,j,k}; neighbor[dim]+=d?-1:1; - int *np=&neighbor[0]; - while (!CELL_IS_BOUNDARY(neighbor,bx.smallEnd(),bx.bigEnd()) && - !CELL_IS_FULL(vof_arr(*np,*(np+1),*(np+2),0)) && - closest_height (*np,*(np+1),*(np+2),dim,hb_arr,ht_arr,&orientation) == hv) { - if (fabs (mv(*np,*(np+1),*(np+2),dim)) > nmax_arr(*np,*(np+1),*(np+2),0)) { + int *np=&neighbor[0]; + while (!CELL_IS_BOUNDARY(neighbor,bx.smallEnd(),bx.bigEnd()) && + !CELL_IS_FULL(vof_arr(*np,*(np+1),*(np+2),0)) && + closest_height (*np,*(np+1),*(np+2),dim,hb_arr,ht_arr,&orientation) == hv) { + if (fabs (mv(*np,*(np+1),*(np+2),dim)) > nmax_arr(*np,*(np+1),*(np+2),0)) { kappa_arr(*np,*(np+1),*(np+2),0) = kappa0; nmax_arr(*np,*(np+1),*(np+2),0) = fabs (mv(*np,*(np+1),*(np+2),dim)); } neighbor[dim]+=d?-1:1; } } - } - } - }); // ParallelFor - }//end MFIter + } + } + }); // ParallelFor + }//end MFIter } //remaining_curvatures @@ -1698,34 +1698,34 @@ VolumeOfFluid::curvature_calculation (int lev, MultiFab & vof_mf, Array const& vof_arr = vof_mf.const_array(mfi); Array4 const& mv = normal[lev].const_array(mfi); - Array4 const& alpha_arr = alpha[lev].const_array(mfi); + Array4 const& alpha_arr = alpha[lev].const_array(mfi); Array4 const& hb_arr = height[0].const_array(mfi); Array4 const& ht_arr = height[1].const_array(mfi); Array4 const& kappa_arr = kappa.array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - auto fvol = vof_arr(i,j,k,0); + { + auto fvol = vof_arr(i,j,k,0); /*if ((i==10&&j==7&&k==6)||(i==8&&j==5&&k==6)){ - int dddd; - Print()<<"------------"<<"\n"; - }*/ - if (!CELL_IS_FULL(fvol)){ - if (kappa_arr(i,j,k,0)==VOF_NODATA){ - // try height function and paraboloid fitting + int dddd; + Print()<<"------------"<<"\n"; + }*/ + if (!CELL_IS_FULL(fvol)){ + if (kappa_arr(i,j,k,0)==VOF_NODATA){ + // try height function and paraboloid fitting Real kappa0= height_curvature_combined (i,j,k, dx,hb_arr,ht_arr,mv,alpha_arr); - if (kappa0!=VOF_NODATA) - kappa_arr(i,j,k,0)=kappa0; - //else - // try particle method (defined in partstr.H) - //kappa_arr(i,j,k,0)= partstr_curvature (i,j,k,dx,problo,vof_arr,mv,alpha_arr); - } - } - }); // ParallelFor - }//end MFIter + if (kappa0!=VOF_NODATA) + kappa_arr(i,j,k,0)=kappa0; + //else + // try particle method (defined in partstr.H) + //kappa_arr(i,j,k,0)= partstr_curvature (i,j,k,dx,problo,vof_arr,mv,alpha_arr); + } + } + }); // ParallelFor + }//end MFIter //!!!!!!!!fix me: a temporary solution for the curvature!!!!!!!!!!! // fill value of ghost cells (BCs, MPI info.) - kappa.FillBoundary(geom.periodicity()); + kappa.FillBoundary(geom.periodicity()); // diffuse curvatures int iter = 1; MultiFab temp_K(kappa.boxArray(), kappa.DistributionMap(),1,kappa.nGrow()); @@ -1735,43 +1735,43 @@ VolumeOfFluid::curvature_calculation (int lev, MultiFab & vof_mf, Array const& temp_arr = temp_K.array(mfi); Array4 const& kappa_arr = kappa.array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - if (kappa_arr(i,j,k,0)!=VOF_NODATA) - temp_arr(i,j,k,0)=kappa_arr(i,j,k,0); - else{ - Real sa=0., s=0.; + { + if (kappa_arr(i,j,k,0)!=VOF_NODATA) + temp_arr(i,j,k,0)=kappa_arr(i,j,k,0); + else{ + Real sa=0., s=0.; /*#if AMREX_SPACEDIM==3 for (int dk = -1; dk <= 1; dk++) #endif for (int dj = -1; dj <= 1; dj++) for (int di = -1; di <= 1; di++) - if (di != 0|| dj != 0|| dk != 0) { - int ni=i+di,nj=j+dj,nk=k+dk; - if (kappa_arr(ni,nj,nk,0)!=VOF_NODATA){ - s += kappa_arr(ni,nj,nk,0); - sa += 1.; - } - }*/ + if (di != 0|| dj != 0|| dk != 0) { + int ni=i+di,nj=j+dj,nk=k+dk; + if (kappa_arr(ni,nj,nk,0)!=VOF_NODATA){ + s += kappa_arr(ni,nj,nk,0); + sa += 1.; + } + }*/ Arraynei; - for (int c = 0; c < AMREX_SPACEDIM; c++){ - nei[0]=i,nei[1]=j,nei[2]=k; - for (int di = -1; di <= 1; di+=2){ - nei[c]=(c==0?i:c==1?j:k)+di; - if (kappa_arr(nei[0],nei[1],nei[2],0)!=VOF_NODATA){ - s += kappa_arr(nei[0],nei[1],nei[2],0); - sa += 1.; - } - } - } - if (sa > 0.) + for (int c = 0; c < AMREX_SPACEDIM; c++){ + nei[0]=i,nei[1]=j,nei[2]=k; + for (int di = -1; di <= 1; di+=2){ + nei[c]=(c==0?i:c==1?j:k)+di; + if (kappa_arr(nei[0],nei[1],nei[2],0)!=VOF_NODATA){ + s += kappa_arr(nei[0],nei[1],nei[2],0); + sa += 1.; + } + } + } + if (sa > 0.) temp_arr(i,j,k,0)=s/sa; - else - temp_arr(i,j,k,0)=VOF_NODATA; - } - }); // ParallelFor - }//end MFIter + else + temp_arr(i,j,k,0)=VOF_NODATA; + } + }); // ParallelFor + }//end MFIter } - MultiFab::Copy(kappa, temp_K, 0, 0, 1, kappa.nGrow()); + MultiFab::Copy(kappa, temp_K, 0, 0, 1, kappa.nGrow()); //fit_curvatures using paraboloid fitting of the centroids of the //reconstructed interface segments @@ -1779,47 +1779,47 @@ VolumeOfFluid::curvature_calculation (int lev, MultiFab & vof_mf, Array const& vof_arr = vof_mf.const_array(mfi); Array4 const& mv = normal[lev].const_array(mfi); - Array4 const& alpha_arr = alpha[lev].const_array(mfi); + Array4 const& alpha_arr = alpha[lev].const_array(mfi); Array4 const& hb_arr = height[0].const_array(mfi); Array4 const& ht_arr = height[1].const_array(mfi); Array4 const& kappa_arr = kappa.array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { + { /* if (i==4&&j==0&&k==4){ - int dddd; - Print()<<"------------"<<"\n"; - }*/ - auto fvol = vof_arr(i,j,k,0); - if (!CELL_IS_FULL(fvol)){ - if (kappa_arr(i,j,k,0)==VOF_NODATA){ - // paraboloid fitting of the centroids of the reconstructed interface segments + int dddd; + Print()<<"------------"<<"\n"; + }*/ + auto fvol = vof_arr(i,j,k,0); + if (!CELL_IS_FULL(fvol)){ + if (kappa_arr(i,j,k,0)==VOF_NODATA){ + // paraboloid fitting of the centroids of the reconstructed interface segments kappa_arr(i,j,k,0) = curvature_fit (i,j,k,dx,vof_arr,mv,alpha_arr); - } - } - }); // ParallelFor - }//end MFIter + } + } + }); // ParallelFor + }//end MFIter //!!!!!!!!fix me: a temporary solution for the curvature!!!!!!!!!!! // fill value of ghost cells (BCs, MPI info.) - kappa.FillBoundary(geom.periodicity()); - + kappa.FillBoundary(geom.periodicity()); + ///////////////////////////////////////////////////////////////////////// /// The following is used to do statistics of calculated curvature /// for numerical tests. ///////////////////////////////////////////////////////////////////////// -if (0){ +if (0){ VofRange kappa_range; range_init(kappa_range); - + struct KappaPrint{ - Real kappa, angle; + Real kappa, angle; XDim3 center; int i,j,k; // Default constructor KappaPrint() : kappa(0), angle(0), center(), i(0), j(0), k(0) {} // Constructor to initialize the KappaPrint KappaPrint(Real ka, Real a, XDim3 o,int i, int j, int k): kappa(ka),angle(a),center(o), - i(i),j(j),k(k){} + i(i),j(j),k(k){} // Copy assignment operator KappaPrint& operator=(const KappaPrint& other) { if (this != &other) { // self-assignment check @@ -1828,18 +1828,18 @@ if (0){ i = other.i; j = other.j; k = other.k; - for (int c = 0; c < AMREX_SPACEDIM; c++) - (¢er.x)[c]=(&other.center.x)[c]; + for (int c = 0; c < AMREX_SPACEDIM; c++) + (¢er.x)[c]=(&other.center.x)[c]; } return *this; - } - }; - Vector kout,removed_elements; + } + }; + Vector kout,removed_elements; Box const& domain = geom.Domain(); - IntVect half= (domain.smallEnd()+domain.bigEnd())/2; + IntVect half= (domain.smallEnd()+domain.bigEnd())/2; Array center{AMREX_D_DECL(0.5*(problo[0]+probhi[0]), 0.5*(problo[1]+probhi[1]), - 0.5*(problo[2]+probhi[2]))}; + 0.5*(problo[2]+probhi[2]))}; #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -1848,47 +1848,47 @@ if (0){ Box const& bx = mfi.tilebox(); const auto lo = lbound(bx); const auto hi = ubound(bx); - + Array4 const& vof_arr = vof_mf.const_array(mfi); - Array4 const& mv = normal[lev].const_array(mfi); - Array4 const& alpha_arr = alpha[lev].const_array(mfi); + Array4 const& mv = normal[lev].const_array(mfi); + Array4 const& alpha_arr = alpha[lev].const_array(mfi); Array4 const& kappa_arr = kappa.const_array(mfi); for (int k = lo.z; k <= hi.z; ++k) { for (int j = lo.y; j <= hi.y; ++j) { for (int i = lo.x; i <= hi.x; ++i) { /* if(i==6&&j==4&&k==7) - Print() <<" ---- "<<"("<0? 0.:180.); - Real nnxy=(p.x>0.?1.:-1.)*sqrt(p.x*p.x+p.y*p.y); - Real angle =acos(nnxy/nn)*180./PI+(p.z>0? 0.:180.); - /* Print() <<" ---- "<<"("<0? 0.:180.); + Real nnxy=(p.x>0.?1.:-1.)*sqrt(p.x*p.x+p.y*p.y); + Real angle =acos(nnxy/nn)*180./PI+(p.z>0? 0.:180.); + /* Print() <<" ---- "<<"("< 1){ @@ -1913,27 +1913,27 @@ if (0){ // Create the MPI datatype MPI_Type_create_struct(6, lengths, disp, types, &mpi_kappa_type); MPI_Type_commit(&mpi_kappa_type); - + // Gather data from all processes Vector recvcounts(nprocs); Vector displs(nprocs, 0); MPI_Gather(&nn, 1, MPI_INT, recvcounts.data(), 1, MPI_INT, 0, MPI_COMM_WORLD); - for (int i = 1; i < nprocs; ++i) + for (int i = 1; i < nprocs; ++i) displs[i] = displs[i-1] + recvcounts[i-1]; - + Vector all_data(displs[nprocs-1] + recvcounts[nprocs-1]); - int kk=sizeof(KappaPrint),tt=sizeof(XDim3), dd=sizeof(all_data); - MPI_Gatherv(kout.data(), kout.size(), mpi_kappa_type, all_data.data(), - recvcounts.data(), displs.data(), mpi_kappa_type, 0, MPI_COMM_WORLD); - kout=all_data; - } - - if (myproc==0){ + int kk=sizeof(KappaPrint),tt=sizeof(XDim3), dd=sizeof(all_data); + MPI_Gatherv(kout.data(), kout.size(), mpi_kappa_type, all_data.data(), + recvcounts.data(), displs.data(), mpi_kappa_type, 0, MPI_COMM_WORLD); + kout=all_data; + } + + if (myproc==0){ // Sort the vector by center.x from high to low std::sort(kout.begin(), kout.end(), [](const KappaPrint& a, const KappaPrint& b) { return a.center.x > b.center.x; }); - + // Use remove_if and copy elements that match center.y<0. to removed_elements auto it = std::remove_if(kout.begin(), kout.end(), [&](const KappaPrint& kp) { if (kp.center.z < 0.) { @@ -1943,31 +1943,31 @@ if (0){ return false; }); // Erase the removed elements from the original vector - kout.erase(it, kout.end()); + kout.erase(it, kout.end()); // Append the removed elements back to the original vector kout.insert(kout.end(), removed_elements.begin(), removed_elements.end()); - - + + Print()<<"# of interfacial cells"< center1{AMREX_D_DECL((problo[0]+.45), (problo[1]+1.1), - (problo[2]+.45))}; + (problo[2]+.45))}; Array center{AMREX_D_DECL(0.5*(problo[0]+probhi[0]), 0.5*(problo[1]+probhi[1]), - 0.5*(problo[2]+probhi[2]))}; + 0.5*(problo[2]+probhi[2]))}; Real radius = .2; //5.0*dx[0]; bool fluid_is_inside = true; EB2::SphereIF my_sphere(radius, center, fluid_is_inside); - EB2::SphereIF my_sphere1(radius, center1, fluid_is_inside); + EB2::SphereIF my_sphere1(radius, center1, fluid_is_inside); // Initialise cylinder parameters int direction = 2; @@ -2266,12 +2266,12 @@ VolumeOfFluid::tracer_vof_init_fraction (int lev, MultiFab& a_tracer) { struct VOFPrint{ - Real vof; + Real vof; int i,j,k; // Default constructor VOFPrint() : vof(0), i(0), j(0), k(0) {} // Constructor to initialize the VOFPrint - VOFPrint(Real ka,int i, int j, int k): vof(ka),i(i),j(j),k(k){} + VOFPrint(Real ka,int i, int j, int k): vof(ka),i(i),j(j),k(k){} // Copy assignment operator VOFPrint& operator=(const VOFPrint& other) { if (this != &other) { // self-assignment check @@ -2281,18 +2281,18 @@ VolumeOfFluid::tracer_vof_init_fraction (int lev, MultiFab& a_tracer) k = other.k; } return *this; - } - }; - Vector vout; + } + }; + Vector vout; // Define the file name - std::string filename = "vof_value-16.dat"; + std::string filename = "vof_value-16.dat"; // Open the file std::ifstream infile(filename); - + if (!infile) { std::cerr << "Unable to open file " << filename << std::endl; - exit; - } + exit; + } // Read the file line by line std::string line; while (std::getline(infile, line)) { @@ -2302,28 +2302,28 @@ VolumeOfFluid::tracer_vof_init_fraction (int lev, MultiFab& a_tracer) if (!(iss >> i >> j >> k >> value)) { std::cerr << "Error reading line: " << line << std::endl; continue; - } + } vout.emplace_back(value,i,j,k); - - } - infile.close(); + + } + infile.close(); #ifdef AMRE_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif for (MFIter mfi(a_tracer,TilingIfNotGPU()); mfi.isValid(); ++mfi) { //Box const& vbx = mfi.validbox(); - Box const& vbx = amrex::grow(mfi.tilebox(),a_tracer.nGrow()); + Box const& vbx = amrex::grow(mfi.tilebox(),a_tracer.nGrow()); auto const& tracer = a_tracer.array(mfi); - - for(int n=0;n=vbx.smallEnd()[0]&&vout[n].i<=vbx.bigEnd()[0]&& - vout[n].j>=vbx.smallEnd()[1]&&vout[n].j<=vbx.bigEnd()[1]&& - vout[n].k>=vbx.smallEnd()[2]&&vout[n].k<=vbx.bigEnd()[2]){ - tracer(vout[n].i,vout[n].j,vout[n].k,0)=vout[n].vof; - } - - /*amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + + for(int n=0;n=vbx.smallEnd()[0]&&vout[n].i<=vbx.bigEnd()[0]&& + vout[n].j>=vbx.smallEnd()[1]&&vout[n].j<=vbx.bigEnd()[1]&& + vout[n].k>=vbx.smallEnd()[2]&&vout[n].k<=vbx.bigEnd()[2]){ + tracer(vout[n].i,vout[n].j,vout[n].k,0)=vout[n].vof; + } + + /*amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { Real x = problo[0] + Real(i+0.5)*dx[0]; Real y = problo[1] + Real(j+0.5)*dx[1]; @@ -2349,14 +2349,14 @@ VolumeOfFluid::tracer_vof_init_fraction (int lev, MultiFab& a_tracer) else tracer(i,j,k) = 0.5-rs; });*/ } - + } // Once vof tracer is initialized, we calculate the normal direction and alpha of the plane segment // intersecting each interface cell. v_incflo->p_volume_of_fluid->tracer_vof_update(lev, a_tracer, height[lev]); - v_incflo->p_volume_of_fluid->curvature_calculation(lev, a_tracer, height[lev], kappa[lev]); - + v_incflo->p_volume_of_fluid->curvature_calculation(lev, a_tracer, height[lev], kappa[lev]); + } @@ -2404,7 +2404,7 @@ VolumeOfFluid::write_tecplot_surface(Real time, int nstep) Array4 const& mv = normal[lev].const_array(mfi); Array4 const& al = alpha[lev].const_array(mfi); Array4 const& tag_arr = tag[lev].const_array(mfi); - Array4 const& kappa_arr = kappa[lev].const_array(mfi); + Array4 const& kappa_arr = kappa[lev].const_array(mfi); Vector segments; int totalnodes = 0; for (int k = lo.z; k <= hi.z; ++k) { @@ -2426,7 +2426,7 @@ VolumeOfFluid::write_tecplot_surface(Real time, int nstep) (&p.x)[dim] = problo[dim] + dx[dim]*((dim<1?i:dim<2?j:k)+(&p.x)[dim]); (¢er.x)[dim] = problo[dim] + dx[dim]*((dim<1?i:dim<2?j:k)+Real(0.5)); } - Array vars={fvol,tag_arr(i,j,k),kappa_arr(i,j,k)}; + Array vars={fvol,tag_arr(i,j,k),kappa_arr(i,j,k)}; /* Print() << " ijk index " <<"("<m_leveldata[lev]; @@ -2532,7 +2532,7 @@ void VolumeOfFluid::WriteTecPlotFile(Real time, int nstep) auto const& ijk_min= bx.smallEnd(); auto const& ijk_max= bx.bigEnd(); std::string zonetitle=("Level_"+std::to_string(lev)+"_Box_"+std::to_string(mfi.index()) - +"_Proc_"+std::to_string(myproc)+"_step_"+std::to_string(nstep)); + +"_Proc_"+std::to_string(myproc)+"_step_"+std::to_string(nstep)); TecplotFile <<(std::string("ZONE T=")+zonetitle); for (int dim = 0; dim < AMREX_SPACEDIM; ++dim) TecplotFile <<", "<<(IJK[dim]+std::string("="))<<(ijk_max[dim]-ijk_min[dim]+2); @@ -2545,9 +2545,9 @@ void VolumeOfFluid::WriteTecPlotFile(Real time, int nstep) Array4 const& mv = normal[lev].const_array(mfi); Array4 const& al = alpha[lev].const_array(mfi); Array4 const& tag_arr = tag[lev].const_array(mfi); - Array4 const& hb_arr = height[lev][0].const_array(mfi); - Array4 const& ht_arr = height[lev][1].const_array(mfi); - Array4 const& kappa_arr = kappa[lev].const_array(mfi); + Array4 const& hb_arr = height[lev][0].const_array(mfi); + Array4 const& ht_arr = height[lev][1].const_array(mfi); + Array4 const& kappa_arr = kappa[lev].const_array(mfi); int nn=0; //write coordinate variables for (int dim = 0; dim < AMREX_SPACEDIM; ++dim) { @@ -2641,7 +2641,7 @@ void VolumeOfFluid::WriteTecPlotFile(Real time, int nstep) for (int dim = 0; dim < AMREX_SPACEDIM; ++dim) { for (int k = lo.z; k <= hi.z; ++k) { for (int j = lo.y; j <= hi.y; ++j) { - for (int i = lo.x; i <= hi.x; ++i) { + for (int i = lo.x; i <= hi.x; ++i) { TecplotFile << hb_arr(i,j,k,dim)<<" "; ++nn; if (nn > 100) { @@ -2651,7 +2651,7 @@ void VolumeOfFluid::WriteTecPlotFile(Real time, int nstep) } } } - }// + }// for (int dim = 0; dim < AMREX_SPACEDIM; ++dim) { for (int k = lo.z; k <= hi.z; ++k) { @@ -2666,7 +2666,7 @@ void VolumeOfFluid::WriteTecPlotFile(Real time, int nstep) } } } - }// + }// //write curvature for (int k = lo.z; k <= hi.z; ++k) { for (int j = lo.y; j <= hi.y; ++j) { @@ -2679,8 +2679,8 @@ void VolumeOfFluid::WriteTecPlotFile(Real time, int nstep) } } } - } - + } + TecplotFile <<"\n"; } // end MFIter @@ -2974,7 +2974,7 @@ void VolumeOfFluid::output_droplet (Real time, int nstep) Real vtop[n_tag], range[AMREX_SPACEDIM][2][n_tag]; for (int n = 0; n < n_tag; n++){ ncell[n]=0; vols[n] = 0.; vels[n] = 0.; surfA[n]=0.; vtop[n] = 0.; - range_init (kappa_range[n]); + range_init (kappa_range[n]); for(int d = 0; d < AMREX_SPACEDIM; d++) { mcent[d][n]=0.; range_init (s[d][n]); @@ -2995,7 +2995,7 @@ void VolumeOfFluid::output_droplet (Real time, int nstep) Array4 const& vel_arr = ld.velocity.const_array(mfi); Array4 const& mv = normal[lev].const_array(mfi); Array4 const& al = alpha[lev].const_array(mfi); - Array4 const& ka = kappa[lev].const_array(mfi); + Array4 const& ka = kappa[lev].const_array(mfi); //fix me: not compatible with GPUs // ParallelFor(bx, [&] AMREX_GPU_DEVICE (int i, int j, int k) noexcept // { @@ -3035,15 +3035,15 @@ void VolumeOfFluid::output_droplet (Real time, int nstep) (&p.x)[d] = problo[d] + dx[d]*((d<1?i:d<2?j:k)+(&p.x)[d]); for(int d = 0; d < AMREX_SPACEDIM; d++) range_add_value (s[d][itag-1], (&p.x)[d]); - // do statistics of the curvature data - range_add_value (kappa_range[itag-1], ka(i,j,k,0)); + // do statistics of the curvature data + range_add_value (kappa_range[itag-1], ka(i,j,k,0)); } } // }); }}} //end of the ijk-loop }//end MFIter - for (int n = 0; n < n_tag; n++) + for (int n = 0; n < n_tag; n++) range_update (kappa_range[n]); // the rest of the algorithm deals with parallel BCs if (ParallelDescriptor::NProcs()> 1){ @@ -3064,11 +3064,11 @@ void VolumeOfFluid::output_droplet (Real time, int nstep) for (int n = 0; n < n_tag; n++) domain_range_reduce(s[d][n]); } - // sum curvature info. - for (int n = 0; n < n_tag; n++){ + // sum curvature info. + for (int n = 0; n < n_tag; n++){ domain_range_reduce(kappa_range[n]); - range_update (kappa_range[n]); - } + range_update (kappa_range[n]); + } } //////////////////////////////////////////////////////////////////// ////// @@ -3084,13 +3084,13 @@ if (0){ cube_min,cube_max; //theoretical centroid of regid body movement for (int d = 0; d < AMREX_SPACEDIM; d++){ - o[d] = o0[d]+1.0*time; + o[d] = o0[d]+1.0*time; cube_min[d]=o[d]-lencube*.5; cube_max[d]=o[d]+lencube*.5; int np=cube_min[d]/(probhi[d]-problo[d]); cube_min[d]-=np*(probhi[d]-problo[d]); np=cube_max[d]/(probhi[d]-problo[d]); - cube_max[d]-=np*(probhi[d]-problo[d]); + cube_max[d]-=np*(probhi[d]-problo[d]); } Print()<<"cube center"<