diff --git a/src/convection/incflo_compute_MAC_projected_velocities.cpp b/src/convection/incflo_compute_MAC_projected_velocities.cpp index 0d5b6961..0ba9aae1 100644 --- a/src/convection/incflo_compute_MAC_projected_velocities.cpp +++ b/src/convection/incflo_compute_MAC_projected_velocities.cpp @@ -180,49 +180,49 @@ incflo::compute_MAC_projected_velocities ( //add surface tension if(m_vof_advect_tracer && m_use_cc_proj) get_volume_of_fluid ()->velocity_face_source(lev,0.5*l_dt, AMREX_D_DECL(*u_mac[lev], *v_mac[lev], *w_mac[lev]), - AMREX_D_DECL(nullptr, nullptr, nullptr)); + AMREX_D_DECL(nullptr, nullptr, nullptr)); if(0){ //The following is only used for testing the pure advection of VOF algorithm -//Average the cell-centered velocity to face center as MAC velocity +//Average the cell-centered velocity to face center as MAC velocity #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif for (MFIter mfi(*vel[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi) - { + { // Note nodaltilebox will not include the nodal index beyond boundaries between neighboring - // titles. Therefore,if we want to use face values (i.e., face_val) immediately below (commented - // out), we must create index space for the face-centered values of the tiled region - // (i.e., surroundingNodes()). - Box const& bx = mfi.tilebox(); - AMREX_D_TERM(Box const& xbx = mfi.nodaltilebox(0);, + // titles. Therefore,if we want to use face values (i.e., face_val) immediately below (commented + // out), we must create index space for the face-centered values of the tiled region + // (i.e., surroundingNodes()). + Box const& bx = mfi.tilebox(); + AMREX_D_TERM(Box const& xbx = mfi.nodaltilebox(0);, Box const& ybx = mfi.nodaltilebox(1);, - Box const& zbx = mfi.nodaltilebox(2);); - Array4 const& velocity = vel[lev]->const_array(mfi); - AMREX_D_TERM(Array4 const& xfv = u_mac[lev]->array(mfi);, - Array4 const& yfv = v_mac[lev]->array(mfi);, - Array4 const& zfv = w_mac[lev]->array(mfi);); - AMREX_D_TERM( + Box const& zbx = mfi.nodaltilebox(2);); + Array4 const& velocity = vel[lev]->const_array(mfi); + AMREX_D_TERM(Array4 const& xfv = u_mac[lev]->array(mfi);, + Array4 const& yfv = v_mac[lev]->array(mfi);, + Array4 const& zfv = w_mac[lev]->array(mfi);); + AMREX_D_TERM( ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - xfv(i,j,k) = .5*(velocity(i,j,k,0)+velocity(i-1,j,k,0)); - });, - ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + xfv(i,j,k) = .5*(velocity(i,j,k,0)+velocity(i-1,j,k,0)); + });, + ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - yfv(i,j,k) = .5*(velocity(i,j,k,1)+velocity(i,j-1,k,1)); - });, + yfv(i,j,k) = .5*(velocity(i,j,k,1)+velocity(i,j-1,k,1)); + });, - ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - zfv(i,j,k) = .5*(velocity(i,j,k,2)+velocity(i,j,k-1,2)); + zfv(i,j,k) = .5*(velocity(i,j,k,2)+velocity(i,j,k-1,2)); }); - ) // end AMREX_D_TERM - - } - return; + ) // end AMREX_D_TERM + + } + return; }//test - } + } Vector > mac_vec(finest_level+1); for (int lev=0; lev <= finest_level; ++lev) diff --git a/src/incflo_compute_forces.cpp b/src/incflo_compute_forces.cpp index 0aab7941..ec13e68a 100644 --- a/src/incflo_compute_forces.cpp +++ b/src/incflo_compute_forces.cpp @@ -41,7 +41,7 @@ void incflo::compute_vel_forces (Vector const& vel_forces, Vector const& tracer_old, Vector const& tracer_new, bool include_pressure_gradient, - bool include_SF) + bool include_SF) { for (int lev = 0; lev <= finest_level; ++lev) compute_vel_forces_on_level (lev, *vel_forces[lev], *velocity[lev], *density[lev], @@ -55,7 +55,7 @@ void incflo::compute_vel_forces_on_level (int lev, const MultiFab& tracer_old, const MultiFab& tracer_new, bool include_pressure_gradient, - bool include_SF) + bool include_SF) { GpuArray l_gravity{m_gravity[0],m_gravity[1],m_gravity[2]}; GpuArray l_gp0{m_gp0[0], m_gp0[1], m_gp0[2]}; @@ -150,262 +150,262 @@ void incflo::compute_vel_forces_on_level (int lev, }); } } -/////////////////////////////////////////////////////////////////////////// - // add surface tension - // surface tension = sigma*kappa*grad(VOF)/rho - // sigma: surface tension coefficient - // kappa: curvature of the interface - // grad(VOF): gradient of VOF field variable - // rho: density - - //fixme: we just consider the surface tension for first tracer +/////////////////////////////////////////////////////////////////////////// + // add surface tension + // surface tension = sigma*kappa*grad(VOF)/rho + // sigma: surface tension coefficient + // kappa: curvature of the interface + // grad(VOF): gradient of VOF field variable + // rho: density + + //fixme: we just consider the surface tension for first tracer if (m_vof_advect_tracer && m_sigma[0]!=0.&&!m_use_cc_proj&&include_SF){ - - VolumeOfFluid* vof_p = get_volume_of_fluid (); - - //choice 1: The original cell-centered kappa and rho are averaged to face center. Grad(VOF) and - // surface tension (SF) are calculated at face center. Then the face-centered SF is finally averaged to cell center. - //choice 2: Similar to choice 1, SF is estimated at the face center and then averaged to the cell nodes. - // the node-centered SF is finally averaged to cell center. - //Choice 3: The original cell-centered rho and VOF are averaged to nodes. The node-centered rho is averaged + + VolumeOfFluid* vof_p = get_volume_of_fluid (); + + //choice 1: The original cell-centered kappa and rho are averaged to face center. Grad(VOF) and + // surface tension (SF) are calculated at face center. Then the face-centered SF is finally averaged to cell center. + //choice 2: Similar to choice 1, SF is estimated at the face center and then averaged to the cell nodes. + // the node-centered SF is finally averaged to cell center. + //Choice 3: The original cell-centered rho and VOF are averaged to nodes. The node-centered rho is averaged // to cell center. grad(VOF) is estimated at the center of the cell edge and then averaged to the cell center. // finally, SF is calculated using original cell-centered kappa, averaged rho, and averaged grad(VOF). - //Choice 4: SF is calculated using original cell-centered kappa, rho and center-difference for grad(VOF). - - int choice = 3; - - - const auto& ba = density.boxArray(); + //Choice 4: SF is calculated using original cell-centered kappa, rho and center-difference for grad(VOF). + + int choice = 3; + + + const auto& ba = density.boxArray(); const auto& dm = density.DistributionMap(); const auto& fact = density.Factory(); Array face_val{AMREX_D_DECL( - MultiFab(amrex::convert(ba,IntVect::TheDimensionVector(0)), + MultiFab(amrex::convert(ba,IntVect::TheDimensionVector(0)), dm, 1, 0, MFInfo(), fact), MultiFab(amrex::convert(ba,IntVect::TheDimensionVector(1)), dm, 1, 0, MFInfo(), fact), MultiFab(amrex::convert(ba,IntVect::TheDimensionVector(2)), dm, 1, 0, MFInfo(), fact))}; - // store the nodal values (the last component stores the node-centered VOF) - MultiFab node_val(amrex::convert(ba,IntVect::TheNodeVector()),dm, AMREX_SPACEDIM+1, 0 , MFInfo(), fact); + // store the nodal values (the last component stores the node-centered VOF) + MultiFab node_val(amrex::convert(ba,IntVect::TheNodeVector()),dm, AMREX_SPACEDIM+1, 0 , MFInfo(), fact); if (choice==1) { - // kappa, rho, and grad(VOF) are first averaged to the center of cell faces. - // The cell-centered force is then obtained by averaging the face-centered value. - average_cellcenter_to_face(GetArrOfPtrs(face_val), density, Geom(lev)); + // kappa, rho, and grad(VOF) are first averaged to the center of cell faces. + // The cell-centered force is then obtained by averaging the face-centered value. + average_cellcenter_to_face(GetArrOfPtrs(face_val), density, Geom(lev)); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { face_val[idim].invert(m_sigma[0], 0); } - + #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif for (MFIter mfi(vel_forces,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { + { // Note nodaltilebox will not include the nodal index beyond boundaries between neighboring - // titles. Therefore,if we want to use face values (i.e., face_val) immediately below (commented - // out), we must create index space for the face-centered values of the tiled region - // (i.e., surroundingNodes()). - Box const& bx = mfi.tilebox(); - AMREX_D_TERM( - Box const& xbx = surroundingNodes(bx,0);, - Box const& ybx = surroundingNodes(bx,1);, - Box const& zbx = surroundingNodes(bx,2);); - Array4 const& tra = tracer_new.const_array(mfi); + // titles. Therefore,if we want to use face values (i.e., face_val) immediately below (commented + // out), we must create index space for the face-centered values of the tiled region + // (i.e., surroundingNodes()). + Box const& bx = mfi.tilebox(); + AMREX_D_TERM( + Box const& xbx = surroundingNodes(bx,0);, + Box const& ybx = surroundingNodes(bx,1);, + Box const& zbx = surroundingNodes(bx,2);); + Array4 const& tra = tracer_new.const_array(mfi); Array4 const& kap = vof_p->kappa[lev].const_array(mfi); - AMREX_D_TERM(Array4 const& xfv = face_val[0].array(mfi);, - Array4 const& yfv = face_val[1].array(mfi);, - Array4 const& zfv = face_val[2].array(mfi);); - AMREX_D_TERM( + AMREX_D_TERM(Array4 const& xfv = face_val[0].array(mfi);, + Array4 const& yfv = face_val[1].array(mfi);, + Array4 const& zfv = face_val[2].array(mfi);); + AMREX_D_TERM( ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { Real kaf; - if(kap(i,j,k,0)!=VOF_NODATA && kap(i-1,j,k,0)!=VOF_NODATA) - kaf=Real(0.5)*(kap(i,j,k,0)+kap(i-1,j,k,0)); - else if (kap(i,j,k,0)!=VOF_NODATA) - kaf=kap(i,j,k); - else if (kap(i-1,j,k,0)!=VOF_NODATA) + if(kap(i,j,k,0)!=VOF_NODATA && kap(i-1,j,k,0)!=VOF_NODATA) + kaf=Real(0.5)*(kap(i,j,k,0)+kap(i-1,j,k,0)); + else if (kap(i,j,k,0)!=VOF_NODATA) + kaf=kap(i,j,k); + else if (kap(i-1,j,k,0)!=VOF_NODATA) kaf=kap(i-1,j,k); - else + else kaf=0.; - xfv(i,j,k) *= kaf*(tra(i,j,k,0)-tra(i-1,j,k,0))/dx[0]; - });, - ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + xfv(i,j,k) *= kaf*(tra(i,j,k,0)-tra(i-1,j,k,0))/dx[0]; + });, + ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { Real kaf; - if(kap(i,j,k,0)!=VOF_NODATA && kap(i,j-1,k,0)!=VOF_NODATA) - kaf=Real(0.5)*(kap(i,j,k,0)+kap(i,j-1,k,0)); - else if (kap(i,j,k,0)!=VOF_NODATA) - kaf=kap(i,j,k); - else if (kap(i,j-1,k,0)!=VOF_NODATA) + if(kap(i,j,k,0)!=VOF_NODATA && kap(i,j-1,k,0)!=VOF_NODATA) + kaf=Real(0.5)*(kap(i,j,k,0)+kap(i,j-1,k,0)); + else if (kap(i,j,k,0)!=VOF_NODATA) + kaf=kap(i,j,k); + else if (kap(i,j-1,k,0)!=VOF_NODATA) kaf=kap(i,j-1,k); - else + else kaf=0.; - yfv(i,j,k) *= kaf*(tra(i,j,k,0)-tra(i,j-1,k,0))/dx[1]; - });, + yfv(i,j,k) *= kaf*(tra(i,j,k,0)-tra(i,j-1,k,0))/dx[1]; + });, - ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { Real kaf; - if(kap(i,j,k,0)!=VOF_NODATA && kap(i,j,k-1,0)!=VOF_NODATA) - kaf=Real(0.5)*(kap(i,j,k,0)+kap(i,j,k-1,0)); - else if (kap(i,j,k,0)!=VOF_NODATA) - kaf=kap(i,j,k); - else if (kap(i,j,k-1,0)!=VOF_NODATA) + if(kap(i,j,k,0)!=VOF_NODATA && kap(i,j,k-1,0)!=VOF_NODATA) + kaf=Real(0.5)*(kap(i,j,k,0)+kap(i,j,k-1,0)); + else if (kap(i,j,k,0)!=VOF_NODATA) + kaf=kap(i,j,k); + else if (kap(i,j,k-1,0)!=VOF_NODATA) kaf=kap(i,j,k-1); - else + else kaf=0.; - zfv(i,j,k) *= kaf*(tra(i,j,k,0)-tra(i,j,k-1,0))/dx[2]; - + zfv(i,j,k) *= kaf*(tra(i,j,k,0)-tra(i,j,k-1,0))/dx[2]; + }); - ) // end AMREX_D_TERM - - // immediately calculate the cell-centered surface tension force using the face-centered value - // If we uncomment the following, we must use surroundingNodes() to define the index space (i.e., - // xbx, ybx, zbx) for cell faces in the beginning. - Array4 const& forarr = vof_p->force[lev].array(mfi); - Array4 const& vel_f = vel_forces.array(mfi); + ) // end AMREX_D_TERM + + // immediately calculate the cell-centered surface tension force using the face-centered value + // If we uncomment the following, we must use surroundingNodes() to define the index space (i.e., + // xbx, ybx, zbx) for cell faces in the beginning. + Array4 const& forarr = vof_p->force[lev].array(mfi); + Array4 const& vel_f = vel_forces.array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - - AMREX_D_TERM(vel_f(i,j,k,0) -= Real(0.5)*(xfv(i,j,k,0)+xfv(i+1,j,k,0));, - vel_f(i,j,k,1) -= Real(0.5)*(yfv(i,j,k,0)+yfv(i,j+1,k,0));, - vel_f(i,j,k,2) -= Real(0.5)*(zfv(i,j,k,0)+zfv(i,j,k+1,0));); - AMREX_D_TERM(forarr(i,j,k,0) = -Real(0.5)*(xfv(i,j,k,0)+xfv(i+1,j,k,0));, - forarr(i,j,k,1) = -Real(0.5)*(yfv(i,j,k,0)+yfv(i,j+1,k,0));, - forarr(i,j,k,2) = -Real(0.5)*(zfv(i,j,k,0)+zfv(i,j,k+1,0));); - }); - } - } + { + + AMREX_D_TERM(vel_f(i,j,k,0) -= Real(0.5)*(xfv(i,j,k,0)+xfv(i+1,j,k,0));, + vel_f(i,j,k,1) -= Real(0.5)*(yfv(i,j,k,0)+yfv(i,j+1,k,0));, + vel_f(i,j,k,2) -= Real(0.5)*(zfv(i,j,k,0)+zfv(i,j,k+1,0));); + AMREX_D_TERM(forarr(i,j,k,0) = -Real(0.5)*(xfv(i,j,k,0)+xfv(i+1,j,k,0));, + forarr(i,j,k,1) = -Real(0.5)*(yfv(i,j,k,0)+yfv(i,j+1,k,0));, + forarr(i,j,k,2) = -Real(0.5)*(zfv(i,j,k,0)+zfv(i,j,k+1,0));); + }); + } + } else if (choice ==2){ // kappa, rho, and grad(VOF) are first averaged to the center of cell faces. // The face-centered values are then averaged to the cell node. - // Finally, the nodal values are averaged to the cell center. - average_cellcenter_to_face(GetArrOfPtrs(face_val), density, Geom(lev)); + // Finally, the nodal values are averaged to the cell center. + average_cellcenter_to_face(GetArrOfPtrs(face_val), density, Geom(lev)); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { face_val[idim].invert(m_sigma[0], 0); } - + #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif for (MFIter mfi(vel_forces,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { + { // Note nodaltilebox will not include the nodal index beyond boundaries between neighboring - // titles. Therefore,if we want to use face values (i.e., face_val) immediately below (commented - // out), we must create index space for the face-centered values of the tiled region - // (i.e., surroundingNodes()).Since we currently calculate all face values in all boxes and then - // convert them to node-centered value, it is better that we use (nodaltilebox()) to avoid the - // repeated calculation of face-centered values at cell faces which are shared by two tiles. - AMREX_D_TERM(Box const& xbx = mfi.nodaltilebox(0);, + // titles. Therefore,if we want to use face values (i.e., face_val) immediately below (commented + // out), we must create index space for the face-centered values of the tiled region + // (i.e., surroundingNodes()).Since we currently calculate all face values in all boxes and then + // convert them to node-centered value, it is better that we use (nodaltilebox()) to avoid the + // repeated calculation of face-centered values at cell faces which are shared by two tiles. + AMREX_D_TERM(Box const& xbx = mfi.nodaltilebox(0);, Box const& ybx = mfi.nodaltilebox(1);, - Box const& zbx = mfi.nodaltilebox(2);); - Array4 const& tra = tracer_new.const_array(mfi); + Box const& zbx = mfi.nodaltilebox(2);); + Array4 const& tra = tracer_new.const_array(mfi); Array4 const& kap = vof_p->kappa[lev].const_array(mfi); - AMREX_D_TERM( Array4 const& xfv = face_val[0].array(mfi);, - Array4 const& yfv = face_val[1].array(mfi);, - Array4 const& zfv = face_val[2].array(mfi);); - - AMREX_D_TERM( + AMREX_D_TERM( Array4 const& xfv = face_val[0].array(mfi);, + Array4 const& yfv = face_val[1].array(mfi);, + Array4 const& zfv = face_val[2].array(mfi);); + + AMREX_D_TERM( ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { Real kaf; - if(kap(i,j,k,0)!=VOF_NODATA && kap(i-1,j,k,0)!=VOF_NODATA) - kaf=Real(0.5)*(kap(i,j,k,0)+kap(i-1,j,k,0)); - else if (kap(i,j,k,0)!=VOF_NODATA) - kaf=kap(i,j,k); - else if (kap(i-1,j,k,0)!=VOF_NODATA) + if(kap(i,j,k,0)!=VOF_NODATA && kap(i-1,j,k,0)!=VOF_NODATA) + kaf=Real(0.5)*(kap(i,j,k,0)+kap(i-1,j,k,0)); + else if (kap(i,j,k,0)!=VOF_NODATA) + kaf=kap(i,j,k); + else if (kap(i-1,j,k,0)!=VOF_NODATA) kaf=kap(i-1,j,k); - else + else kaf=0.; - xfv(i,j,k) *= kaf*(tra(i,j,k,0)-tra(i-1,j,k,0))/dx[0]; - });, - ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + xfv(i,j,k) *= kaf*(tra(i,j,k,0)-tra(i-1,j,k,0))/dx[0]; + });, + ParallelFor(ybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { Real kaf; - if(kap(i,j,k,0)!=VOF_NODATA && kap(i,j-1,k,0)!=VOF_NODATA) - kaf=Real(0.5)*(kap(i,j,k,0)+kap(i,j-1,k,0)); - else if (kap(i,j,k,0)!=VOF_NODATA) - kaf=kap(i,j,k); - else if (kap(i,j-1,k,0)!=VOF_NODATA) + if(kap(i,j,k,0)!=VOF_NODATA && kap(i,j-1,k,0)!=VOF_NODATA) + kaf=Real(0.5)*(kap(i,j,k,0)+kap(i,j-1,k,0)); + else if (kap(i,j,k,0)!=VOF_NODATA) + kaf=kap(i,j,k); + else if (kap(i,j-1,k,0)!=VOF_NODATA) kaf=kap(i,j-1,k); - else + else kaf=0.; - yfv(i,j,k) *= kaf*(tra(i,j,k,0)-tra(i,j-1,k,0))/dx[1]; - });, - ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + yfv(i,j,k) *= kaf*(tra(i,j,k,0)-tra(i,j-1,k,0))/dx[1]; + });, + ParallelFor(zbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { Real kaf; - if(kap(i,j,k,0)!=VOF_NODATA && kap(i,j,k-1,0)!=VOF_NODATA) - kaf=Real(0.5)*(kap(i,j,k,0)+kap(i,j,k-1,0)); - else if (kap(i,j,k,0)!=VOF_NODATA) - kaf=kap(i,j,k); - else if (kap(i,j,k-1,0)!=VOF_NODATA) + if(kap(i,j,k,0)!=VOF_NODATA && kap(i,j,k-1,0)!=VOF_NODATA) + kaf=Real(0.5)*(kap(i,j,k,0)+kap(i,j,k-1,0)); + else if (kap(i,j,k,0)!=VOF_NODATA) + kaf=kap(i,j,k); + else if (kap(i,j,k-1,0)!=VOF_NODATA) kaf=kap(i,j,k-1); - else + else kaf=0.; - zfv(i,j,k) *= kaf*(tra(i,j,k,0)-tra(i,j,k-1,0))/dx[2]; - + zfv(i,j,k) *= kaf*(tra(i,j,k,0)-tra(i,j,k-1,0))/dx[2]; + }); - ) // end AMREX_D_TERM - - } + ) // end AMREX_D_TERM + + } static int oct[3][2] = { { 1, 2 }, { 0, 2 }, { 0, 1 } }; #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif for (MFIter mfi(vel_forces,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box const& bx = mfi.tilebox(); - //We will immediately use the nodal values so we must use surroundingNodes()instead of nodaltilebox() - //See the previous comments for calculating face-centered value. - Box const& nbx = surroundingNodes(bx); - Box const& vbx = mfi.validbox(); - AMREX_D_TERM( - Box const& xvbx = surroundingNodes(vbx,0);, - Box const& yvbx = surroundingNodes(vbx,1);, - Box const& zvbx = surroundingNodes(vbx,2);); - AMREX_D_TERM( - Array4 const& xfv = face_val[0].array(mfi);, - Array4 const& yfv = face_val[1].array(mfi);, - Array4 const& zfv = face_val[2].array(mfi);); - Array4 const& nv = node_val.array(mfi); + { + Box const& bx = mfi.tilebox(); + //We will immediately use the nodal values so we must use surroundingNodes()instead of nodaltilebox() + //See the previous comments for calculating face-centered value. + Box const& nbx = surroundingNodes(bx); + Box const& vbx = mfi.validbox(); + AMREX_D_TERM( + Box const& xvbx = surroundingNodes(vbx,0);, + Box const& yvbx = surroundingNodes(vbx,1);, + Box const& zvbx = surroundingNodes(vbx,2);); + AMREX_D_TERM( + Array4 const& xfv = face_val[0].array(mfi);, + Array4 const& yfv = face_val[1].array(mfi);, + Array4 const& zfv = face_val[2].array(mfi);); + Array4 const& nv = node_val.array(mfi); ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - for (int dim = 0; dim < AMREX_SPACEDIM; ++dim){ - nv(i,j,k,dim)=0.; - int nt=0; - for (int detj = 0; detj < 2; ++detj) + { + for (int dim = 0; dim < AMREX_SPACEDIM; ++dim){ + nv(i,j,k,dim)=0.; + int nt=0; + for (int detj = 0; detj < 2; ++detj) for (int deti = 0; deti < 2; ++deti){ - Array in{i, j, k}; - in[oct[dim][0]]-=deti,in[oct[dim][1]]-=detj; - if (dim==0&& xvbx.contains(in[0],in[1],in[2])){ - nv(i,j,k,dim)+= xfv(in[0],in[1],in[2]); - nt++; - } - else if (dim==1&& yvbx.contains(in[0],in[1],in[2])){ - nv(i,j,k,dim)+= yfv(in[0],in[1],in[2]); - nt++; - } -#if AMREX_SPACEDIM == 3 - else if (dim==2&& zvbx.contains(in[0],in[1],in[2])){ - nv(i,j,k,dim)+= zfv(in[0],in[1],in[2]); - nt++; - } -#endif - } - if(nt>0) nv(i,j,k,dim)/= Real(nt); - } - }); - Array4 const& forarr = vof_p->force[lev].array(mfi); - Array4 const& vel_f = vel_forces.array(mfi); + Array in{i, j, k}; + in[oct[dim][0]]-=deti,in[oct[dim][1]]-=detj; + if (dim==0&& xvbx.contains(in[0],in[1],in[2])){ + nv(i,j,k,dim)+= xfv(in[0],in[1],in[2]); + nt++; + } + else if (dim==1&& yvbx.contains(in[0],in[1],in[2])){ + nv(i,j,k,dim)+= yfv(in[0],in[1],in[2]); + nt++; + } +#if AMREX_SPACEDIM == 3 + else if (dim==2&& zvbx.contains(in[0],in[1],in[2])){ + nv(i,j,k,dim)+= zfv(in[0],in[1],in[2]); + nt++; + } +#endif + } + if(nt>0) nv(i,j,k,dim)/= Real(nt); + } + }); + Array4 const& forarr = vof_p->force[lev].array(mfi); + Array4 const& vel_f = vel_forces.array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - for (int dim = 0; dim < AMREX_SPACEDIM; ++dim){ + { + for (int dim = 0; dim < AMREX_SPACEDIM; ++dim){ #if AMREX_SPACEDIM==2 /* 2D */ vel_f(i,j,k,dim) -= Real(0.25)*(nv(i,j ,k ,dim) + nv(i+1,j ,k ,dim) + nv(i,j+1,k ,dim) + nv(i+1,j+1,k ,dim)); forarr(i,j,k,dim) =-Real(0.25)*(nv(i,j ,k ,dim) + nv(i+1,j ,k ,dim) + nv(i,j+1,k ,dim) + nv(i+1,j+1,k ,dim)); -#else /* 3D */ +#else /* 3D */ vel_f(i,j,k,dim) -= Real(0.125)*(nv(i,j ,k ,dim) + nv(i+1,j ,k ,dim) + nv(i,j+1,k ,dim) + nv(i+1,j+1,k ,dim) + nv(i,j ,k+1,dim) + nv(i+1,j ,k+1,dim) @@ -414,124 +414,124 @@ void incflo::compute_vel_forces_on_level (int lev, + nv(i,j+1,k ,dim) + nv(i+1,j+1,k ,dim) + nv(i,j ,k+1,dim) + nv(i+1,j ,k+1,dim) + nv(i,j+1,k+1,dim) + nv(i+1,j+1,k+1,dim)); -#endif - } - }); - } - } - else if (choice ==3){ - MultiFab center_val(ba,dm,2,0,MFInfo(), fact); +#endif + } + }); + } + } + else if (choice ==3){ + MultiFab center_val(ba,dm,2,0,MFInfo(), fact); static int oct[3][2] = { { 1, 2 }, { 0, 2 }, { 0, 1 } }; #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif for (MFIter mfi(vel_forces,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box const& nbx = surroundingNodes(mfi.tilebox()); - Box const& vbx = mfi.validbox(); - Array4 const& nv = node_val.array(mfi); - Array4 const& tra = tracer_new.const_array(mfi); - Array4 const& kappa = vof_p->kappa[lev].const_array(mfi); - Array4 const& rho = density.const_array(mfi); + { + Box const& nbx = surroundingNodes(mfi.tilebox()); + Box const& vbx = mfi.validbox(); + Array4 const& nv = node_val.array(mfi); + Array4 const& tra = tracer_new.const_array(mfi); + Array4 const& kappa = vof_p->kappa[lev].const_array(mfi); + Array4 const& rho = density.const_array(mfi); ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - /* if(i==8&&j==8&&k==8){ - Print()<<"zbx "<<"low "< -2; --detk) -#endif - for (int detj = 0; detj > -2; --detj) + { + /* if(i==8&&j==8&&k==8){ + Print()<<"zbx "<<"low "< -2; --detk) +#endif + for (int detj = 0; detj > -2; --detj) for (int deti = 0; deti > -2; --deti) { - Array in{i+deti,j+detj,k+detk}; - //if (vbx.contains(in[0],in[1],in[2])){ - // averaging VOF to nodes - nv(i,j,k,AMREX_SPACEDIM)+= tra(in[0],in[1],in[2]); - nt++; - //averaging kappa to nodes - if(kappa(in[0],in[1],in[2],0)!=VOF_NODATA){ - nv(i,j,k,0)+= kappa(in[0],in[1],in[2],0); - nkap++; - } - //averaging density to nodes - nv(i,j,k,1)+= rho(in[0],in[1],in[2],0); - nrho++; - //} - } - if(nt>0) nv(i,j,k,AMREX_SPACEDIM)/= Real(nt); - if(nkap>0) - nv(i,j,k,0)/= Real(nkap); - else - nv(i,j,k,0)=0.; - nv(i,j,k,1)/= Real(nrho); - }); - } - average_node_to_cellcenter(center_val, 0, node_val, 0, 2); - //MultiFab::Copy(density, center_val , 1, 0, 1, density.nGrow()); + Array in{i+deti,j+detj,k+detk}; + //if (vbx.contains(in[0],in[1],in[2])){ + // averaging VOF to nodes + nv(i,j,k,AMREX_SPACEDIM)+= tra(in[0],in[1],in[2]); + nt++; + //averaging kappa to nodes + if(kappa(in[0],in[1],in[2],0)!=VOF_NODATA){ + nv(i,j,k,0)+= kappa(in[0],in[1],in[2],0); + nkap++; + } + //averaging density to nodes + nv(i,j,k,1)+= rho(in[0],in[1],in[2],0); + nrho++; + //} + } + if(nt>0) nv(i,j,k,AMREX_SPACEDIM)/= Real(nt); + if(nkap>0) + nv(i,j,k,0)/= Real(nkap); + else + nv(i,j,k,0)=0.; + nv(i,j,k,1)/= Real(nrho); + }); + } + average_node_to_cellcenter(center_val, 0, node_val, 0, 2); + //MultiFab::Copy(density, center_val , 1, 0, 1, density.nGrow()); #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif for (MFIter mfi(vel_forces,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box const& bx = mfi.tilebox(); - //We will immediately use the nodal values so we must use surroundingNodes()instead of nodaltilebox() - //See the previous comments for calculating face-centered value. - Box const& nbx = surroundingNodes(bx); - Box const& vbx = mfi.validbox(); - Array4 const& nv = node_val.array(mfi); - Array4 const& tra = tracer_new.const_array(mfi); - Array4 const& kappa = vof_p->kappa[lev].const_array(mfi); - Array4 const& rho = density.const_array(mfi); - Array4 const& forarr = vof_p->force[lev].array(mfi); - Array4 const& vel_f = vel_forces.array(mfi); - Array4 const& center = center_val.array(mfi); + { + Box const& bx = mfi.tilebox(); + //We will immediately use the nodal values so we must use surroundingNodes()instead of nodaltilebox() + //See the previous comments for calculating face-centered value. + Box const& nbx = surroundingNodes(bx); + Box const& vbx = mfi.validbox(); + Array4 const& nv = node_val.array(mfi); + Array4 const& tra = tracer_new.const_array(mfi); + Array4 const& kappa = vof_p->kappa[lev].const_array(mfi); + Array4 const& rho = density.const_array(mfi); + Array4 const& forarr = vof_p->force[lev].array(mfi); + Array4 const& vel_f = vel_forces.array(mfi); + Array4 const& center = center_val.array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Array gradVof{0.}; - for (int dim = 0; dim < AMREX_SPACEDIM; ++dim){ - int ng=0; -#if AMREX_SPACEDIM==3 - for (int detj = 0; detj < 2; ++detj) -#endif + { + Array gradVof{0.}; + for (int dim = 0; dim < AMREX_SPACEDIM; ++dim){ + int ng=0; +#if AMREX_SPACEDIM==3 + for (int detj = 0; detj < 2; ++detj) +#endif for (int deti = 0; deti < 2; ++deti){ - Array in0{i, j, k},in1{i, j, k}; - in0[oct[dim][0]]+=deti,in1[oct[dim][0]]+=deti; -#if AMREX_SPACEDIM==3 + Array in0{i, j, k},in1{i, j, k}; + in0[oct[dim][0]]+=deti,in1[oct[dim][0]]+=deti; +#if AMREX_SPACEDIM==3 in0[oct[dim][1]]+=detj,in1[oct[dim][1]]+=detj; -#endif - in1[dim] +=1; - gradVof[dim] +=(nv(in1[0],in1[1],in1[2],AMREX_SPACEDIM)- - nv(in0[0],in0[1],in0[2],AMREX_SPACEDIM))/dx[dim]; - ng++; - } - gradVof[dim]/=Real(ng); - /*Array in0{i, j, k},in1{i, j, k}; - in1[dim] +=1,in0[dim] -=1;; - gradVof[dim] = Real(0.5)*(tra(in1[0],in1[1],in1[2],0)-tra(in0[0],in0[1],in0[2],0))/dx[dim];*/ - } - for (int dim = 0; dim < AMREX_SPACEDIM; ++dim){ +#endif + in1[dim] +=1; + gradVof[dim] +=(nv(in1[0],in1[1],in1[2],AMREX_SPACEDIM)- + nv(in0[0],in0[1],in0[2],AMREX_SPACEDIM))/dx[dim]; + ng++; + } + gradVof[dim]/=Real(ng); + /*Array in0{i, j, k},in1{i, j, k}; + in1[dim] +=1,in0[dim] -=1;; + gradVof[dim] = Real(0.5)*(tra(in1[0],in1[1],in1[2],0)-tra(in0[0],in0[1],in0[2],0))/dx[dim];*/ + } + for (int dim = 0; dim < AMREX_SPACEDIM; ++dim){ if(kappa(i,j,k,0)!=VOF_NODATA){ - vel_f(i,j,k,dim) -= m_sigma[0]*kappa(i,j,k,0)/center(i,j,k,1)*gradVof[dim]; - forarr(i,j,k,dim) =-m_sigma[0]*kappa(i,j,k,0)/center(i,j,k,1)*gradVof[dim]; - } + vel_f(i,j,k,dim) -= m_sigma[0]*kappa(i,j,k,0)/center(i,j,k,1)*gradVof[dim]; + forarr(i,j,k,dim) =-m_sigma[0]*kappa(i,j,k,0)/center(i,j,k,1)*gradVof[dim]; + } else { for (int dim = 0; dim < AMREX_SPACEDIM; ++dim) - forarr(i,j,k,dim) = 0.; - } - } - }); - - - } - } + forarr(i,j,k,dim) = 0.; + } + } + }); + + + } + } else if (choice ==4) { - //cell-centered surface tension force + //cell-centered surface tension force #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif @@ -540,36 +540,36 @@ static int oct[3][2] = { { 1, 2 }, { 0, 2 }, { 0, 1 } }; Box const& bx = mfi.tilebox(); Array4 const& vel_f = vel_forces.array(mfi); Array4 const& rho = density.const_array(mfi); - Array4 const& tra = tracer_new.const_array(mfi); - Array4 const& kappa = vof_p->kappa[lev].const_array(mfi); - Array4 const& forarr = vof_p->force[lev].array(mfi); + Array4 const& tra = tracer_new.const_array(mfi); + Array4 const& kappa = vof_p->kappa[lev].const_array(mfi); + Array4 const& forarr = vof_p->force[lev].array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - if(kappa(i,j,k,0)!=VOF_NODATA){ + { + if(kappa(i,j,k,0)!=VOF_NODATA){ Real sig_kappa = m_sigma[0]*kappa(i,j,k,0)/rho(i,j,k); - //note: the minus sign is used because of the way curvature is calculated - AMREX_D_TERM( - vel_f(i,j,k,0) -= Real(0.5)*(tra(i+1,j,k,0)-tra(i-1,j,k,0))/dx[0]*sig_kappa;, - vel_f(i,j,k,1) -= Real(0.5)*(tra(i,j+1,k,0)-tra(i,j-1,k,0))/dx[1]*sig_kappa;, - vel_f(i,j,k,2) -= Real(0.5)*(tra(i,j,k+1,0)-tra(i,j,k-1,0))/dx[2]*sig_kappa;); - - AMREX_D_TERM( - forarr(i,j,k,0) = -Real(0.5)*(tra(i+1,j,k,0)-tra(i-1,j,k,0))/dx[0]*sig_kappa;, - forarr(i,j,k,1) = -Real(0.5)*(tra(i,j+1,k,0)-tra(i,j-1,k,0))/dx[1]*sig_kappa;, - forarr(i,j,k,2) = -Real(0.5)*(tra(i,j,k+1,0)-tra(i,j,k-1,0))/dx[2]*sig_kappa;); - - } + //note: the minus sign is used because of the way curvature is calculated + AMREX_D_TERM( + vel_f(i,j,k,0) -= Real(0.5)*(tra(i+1,j,k,0)-tra(i-1,j,k,0))/dx[0]*sig_kappa;, + vel_f(i,j,k,1) -= Real(0.5)*(tra(i,j+1,k,0)-tra(i,j-1,k,0))/dx[1]*sig_kappa;, + vel_f(i,j,k,2) -= Real(0.5)*(tra(i,j,k+1,0)-tra(i,j,k-1,0))/dx[2]*sig_kappa;); + + AMREX_D_TERM( + forarr(i,j,k,0) = -Real(0.5)*(tra(i+1,j,k,0)-tra(i-1,j,k,0))/dx[0]*sig_kappa;, + forarr(i,j,k,1) = -Real(0.5)*(tra(i,j+1,k,0)-tra(i,j-1,k,0))/dx[1]*sig_kappa;, + forarr(i,j,k,2) = -Real(0.5)*(tra(i,j,k+1,0)-tra(i,j,k-1,0))/dx[2]*sig_kappa;); + + } else { for (int dim = 0; dim < AMREX_SPACEDIM; ++dim) - forarr(i,j,k,dim) = 0.; - } - }); - } - - } - - - - - }// end if (m_vof_advect_tracer) + forarr(i,j,k,dim) = 0.; + } + }); + } + + } + + + + + }// end if (m_vof_advect_tracer) } diff --git a/src/incflo_update_density.cpp b/src/incflo_update_density.cpp index c3680226..7a16830a 100644 --- a/src/incflo_update_density.cpp +++ b/src/incflo_update_density.cpp @@ -67,9 +67,9 @@ void incflo::update_density (StepType step_type) } else { for (int lev = 0; lev <= finest_level; lev++) { - if (m_vof_advect_tracer){ - MultiFab::Copy(m_leveldata[lev]->density, m_leveldata[lev]->density_o, 0, 0, 1, m_leveldata[lev]->density_o.nGrow()); - } + if (m_vof_advect_tracer){ + MultiFab::Copy(m_leveldata[lev]->density, m_leveldata[lev]->density_o, 0, 0, 1, m_leveldata[lev]->density_o.nGrow()); + } MultiFab::Copy(m_leveldata[lev]->density_nph, m_leveldata[lev]->density_o, 0, 0, 1, ng); } } diff --git a/src/prob/prob_init_fluid.cpp b/src/prob/prob_init_fluid.cpp index d095aebb..09ed07bf 100644 --- a/src/prob/prob_init_fluid.cpp +++ b/src/prob/prob_init_fluid.cpp @@ -207,15 +207,15 @@ void incflo::prob_init_fluid (int lev) if (1109 == m_probtype) { get_volume_of_fluid ()->tracer_vof_init_fraction(lev, ld.tracer); MultiFab::Copy(ld.tracer_o, ld.tracer, 0, 0, 1, ld.tracer.nGrow()); - ld.tracer_o.FillBoundary(geom[lev].periodicity()); - if (m_vof_advect_tracer){ + ld.tracer_o.FillBoundary(geom[lev].periodicity()); + if (m_vof_advect_tracer){ update_vof_density (lev, get_density_new(),get_tracer_new()); MultiFab::Copy(ld.density_o, ld.density, 0, 0, 1, ld.density.nGrow()); fillpatch_density(lev, m_t_new[lev], ld.density_o, 3); - MultiFab::Copy(ld.density_nph, ld.density, 0, 0, 1, ld.density.nGrow()); - fillpatch_density(lev, m_t_new[lev], ld.density_nph, 3); - } - + MultiFab::Copy(ld.density_nph, ld.density, 0, 0, 1, ld.density.nGrow()); + fillpatch_density(lev, m_t_new[lev], ld.density_nph, 3); + } + } } diff --git a/src/vof/VolumeOfFluid.cpp b/src/vof/VolumeOfFluid.cpp index 66d2ac7c..bfd79ab7 100644 --- a/src/vof/VolumeOfFluid.cpp +++ b/src/vof/VolumeOfFluid.cpp @@ -58,7 +58,7 @@ void range_add_value (VofRange & r, Real val) /** * range_update: * @r: a #VofRange. - * + * * Updates the fields of @r. */ void range_update (VofRange & r) @@ -66,12 +66,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.; } @@ -93,7 +93,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; @@ -123,17 +123,17 @@ 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)); - force.emplace_back(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)); + force.emplace_back(v_incflo->grids[lev], v_incflo->dmap[lev], AMREX_SPACEDIM, v_incflo->nghost_state(), MFInfo(), v_incflo->Factory(lev)); //fixme - force[lev].setVal(0,0,AMREX_SPACEDIM,v_incflo->nghost_state()); - } + force[lev].setVal(0,0,AMREX_SPACEDIM,v_incflo->nghost_state()); + } } static XDim3 edge[12][2] = { @@ -290,12 +290,12 @@ struct Segment{ //vof, tag Array vars; // Constructor to initialize the Segment - Segment(int n, + Segment(int n, #if AMREX_SPACEDIM==2 /* 2D */ - Array const& nodes, + Array const& nodes, #else /* 3D */ Array const& nodes, -#endif +#endif XDim3 m, Real a, Array v, int ns=0) : nnodes(n), mv (m), alpha(a), vars(v) { for (int i = 0; i < n; ++i) @@ -310,12 +310,12 @@ static void add_segment (XDim3 const & center, GpuArray co /*Print() <<" add_segment "<< *o<<" "<<" vector "<<*m <<"vof"<<" "< nodecutface; + Array nodecutface; Real x, y, h=dx[0]; - int n=0, nnodecutface; + int n=0, nnodecutface; if (fabs (m.y) > EPS) { y = (alpha - m.x)/m.y; if (y >= 0. && y <= 1.) { @@ -339,18 +339,18 @@ static void add_segment (XDim3 const & center, GpuArray co if (x >= 0. && x <= 1.) { nodecutface[n].x = center.x + h*(x - 0.5); nodecutface[n].y = center.y - h/2.; nodecutface[n++].z = 0.; } - } + } nnodecutface = n; if (n > 2) { /*check if there are duplicated points*/ int i,j; bool ok[n]; for (i=0; i co if (!ok[i]){ if (i!=n-1) for (j=i+1; j nodecutface; @@ -461,8 +461,8 @@ static void add_segment (XDim3 const & center, GpuArray co nt += nnodecutface; segments.emplace_back(nnodecutface, nodecutface, m, alpha, vars, 3); } /* cut face must be divided into 2 quadrilateral/triangular faces */ - -#endif + +#endif } /** @@ -505,7 +505,7 @@ Real line_area (Array &m, Real alpha) a = alpha1 - n.x; if (a > 0.) v -= a*a; - + a = alpha1 - n.y; if (a > 0.) v -= a*a; @@ -522,17 +522,17 @@ Real line_area (Array &m, Real alpha) * @c: a volume fraction. * * Returns: the value @alpha such that the area of a square cell - * lying under the line defined by @m.@x = @alpha is equal to @c. + * lying under the line defined by @m.@x = @alpha is equal to @c. */ Real line_alpha (XDim3 & m, Real c) { Real alpha, m1, m2, v1; - + m1 = fabs (m.x); m2 = fabs (m.y); if (m1 > m2) { v1 = m1; m1 = m2; m2 = v1; } - + v1 = m1/2.; if (c <= v1/m2) alpha = sqrt (2.*c*m1*m2); @@ -610,7 +610,7 @@ void line_center (XDim3 const & m, Real alpha, Real a, XDim3 & p) p.y -= b*b*(alpha + 2.*n.y); p.x -= b*b*b; } - + p.x /= 6.*n.x*n.x*n.y*a; p.y /= 6.*n.x*n.y*n.y*a; @@ -1060,7 +1060,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; @@ -1068,12 +1068,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; @@ -1086,36 +1086,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; @@ -1137,7 +1137,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.) @@ -1156,32 +1156,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.; } } @@ -1194,8 +1194,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; @@ -1206,14 +1206,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; @@ -1221,23 +1221,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]={0.}, hd[2]={0.}; for (int nd = 0; nd < (AMREX_SPACEDIM==3?2:1); nd++) { @@ -1250,7 +1250,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; @@ -1258,16 +1258,16 @@ static bool height_normal (int i,int j,int k, Array4 const & hb, slope = hd[0]*hd[0] + hd[1]*hd[1]; (&m.x)[d] = orientation; (&m.x)[oc[d][0]] = - hd[0]; -#if AMREX_SPACEDIM==3 +#if AMREX_SPACEDIM==3 (&m.x)[oc[d][1]] = - hd[1]; -#endif - } - - } +#endif + } + + } } //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[AMREX_SPACEDIM==3?9:3], h[AMREX_SPACEDIM==3?9:3]; - 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 < (AMREX_SPACEDIM==3?2:1) && loop; nd++) - for (int ndd = 0; ndd < 2; ndd++) { + for (int nd = 0; nd < (AMREX_SPACEDIM==3?2:1) && 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 < (AMREX_SPACEDIM==3?2:1); nd++) { h[n] = neighboring_column (i, j, k, d, hv, 2*oc[d][nd], &x[n]); @@ -1369,13 +1369,13 @@ 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[AMREX_SPACEDIM==3?9:3], h[AMREX_SPACEDIM==3?9:3]; - 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 < (AMREX_SPACEDIM==3?2:1) && loop; nd++) - for (int ndd = 0; ndd < 2; ndd++) { + for (int nd = 0; nd < (AMREX_SPACEDIM==3?2:1) && 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 < (AMREX_SPACEDIM==3?2:1); nd++) { h[n] = neighboring_column (i, j, k, d, hv, 2*oc[d][nd], &x[n]); @@ -1460,25 +1460,25 @@ bool curvature_along_direction_new (int i,int j,int k, int d, GpuArray &c) @@ -1536,9 +1536,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; } } @@ -1554,7 +1554,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; @@ -1583,43 +1583,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 @@ -1634,11 +1634,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 @@ -1685,17 +1685,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 || AMREX_SPACEDIM==2 @@ -1703,11 +1703,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 @@ -1721,146 +1721,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,v_incflo->nghost_state()); - height[1].setVal(VOF_NODATA,dim,1,v_incflo->nghost_state()); - //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,v_incflo->nghost_state()); + //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.);, @@ -1871,11 +1871,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; @@ -1912,19 +1912,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)); - //fixme: need to change for BCs + MFInfo(), v_incflo->Factory(lev)); + //fixme: need to change for BCs kappa.setVal(VOF_NODATA,0,1,v_incflo->nghost_state()); - 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); @@ -1932,42 +1932,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==4||i==11)&&j==9&&k==0){ - 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 @@ -1975,84 +1975,84 @@ 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==4||i==11)&&j==11&&k==0){ - 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 //!!!!!!!!fixme: 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 = 0; if (iter >0){ MultiFab temp_K(kappa.boxArray(), kappa.DistributionMap(),1,kappa.nGrow()); - //fixme: need to change for BCs - temp_K.setVal(VOF_NODATA,0,1,kappa.nGrow()); + //fixme: need to change for BCs + temp_K.setVal(VOF_NODATA,0,1,kappa.nGrow()); while (iter--){ for (MFIter mfi(vof_mf,TilingIfNotGPU()); mfi.isValid(); ++mfi) { Box const& bx = mfi.tilebox(); Array4 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 @@ -2060,47 +2060,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||i==11)&&j==11&&k==0){ - int dddd; - Print()<<"------------"< 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 @@ -2129,47 +2129,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){ @@ -2194,27 +2194,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.) { @@ -2224,31 +2224,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"< const& tracer, m_total_flux[lev].setVal(0.0); vof_total_flux[lev].setVal(0.0); MultiFab const * U_MF = dir < 1? u_mac[lev]: -#if AMREX_SPACEDIM == 3 +#if AMREX_SPACEDIM == 3 dir >= 2? w_mac[lev]: -#endif - v_mac[lev]; +#endif + v_mac[lev]; #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) @@ -2464,22 +2464,22 @@ VolumeOfFluid::tracer_vof_advection(Vector const& tracer, // add surface tension force (F) to the MAC velocity at the center of cell faces. // F = dt*sigma*kappa*grad(VOF)/rho // Umac <- Uma-F, note minus sign before F because of the way which curvature being calculated. -// kappa and rho need to be estimated at the cell face center. The simple average of the cell-centered +// kappa and rho need to be estimated at the cell face center. The simple average of the cell-centered // values of two neighboring cells dilimited by the face is used to calculate the face-centered value. -// grad(VOF) is also estimated at the face center using the center-difference method for two cells +// grad(VOF) is also estimated at the face center using the center-difference method for two cells // i.e., in x-dir, grad(VOF)= (VOFcell[1]-VOFcell[0])/dx[0]. // u_mac/v_mac/w_mac stores the face-centered velocity (MAC). // // gu_mac/gv_mac/gw_mac stores the face-centered value of F/dt (i.e., sigma*kappa*grad(VOF)/rho) // gu_mac/gv_mac/gw_mac will be averaged to the cell center when correcting the cell-centered velocity -// after final cell-centered projection. -void +// after final cell-centered projection. +void VolumeOfFluid:: velocity_face_source (int lev, Real dt, AMREX_D_DECL(MultiFab& u_mac, MultiFab& v_mac, MultiFab& w_mac), AMREX_D_DECL(MultiFab* gu_mac, MultiFab* gv_mac, MultiFab* gw_mac)) { auto& ld = *v_incflo->m_leveldata[lev]; Geometry const& geom = v_incflo->Geom(lev); - auto const& dx = geom.CellSizeArray(); + auto const& dx = geom.CellSizeArray(); Real sigma = v_incflo->m_sigma[0]; #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) @@ -2489,71 +2489,71 @@ VolumeOfFluid:: velocity_face_source (int lev, Real dt, AMREX_D_DECL(MultiFab& u Box const& bx = mfi.tilebox(); Box const& xbx = mfi.nodaltilebox(0); Box const& ybx = mfi.nodaltilebox(1); - Box const& zbx = mfi.nodaltilebox(2); + Box const& zbx = mfi.nodaltilebox(2); Array4 const& rho = ld.density.const_array(mfi); - Array4 const& tra = ld.tracer.const_array(mfi); + Array4 const& tra = ld.tracer.const_array(mfi); Array4 const& kap = kappa[lev].const_array(mfi); - AMREX_D_TERM(Array4 const& umac = u_mac.array(mfi);, - Array4 const& vmac = v_mac.array(mfi);, - Array4 const& wmac = w_mac.array(mfi);); + AMREX_D_TERM(Array4 const& umac = u_mac.array(mfi);, + Array4 const& vmac = v_mac.array(mfi);, + Array4 const& wmac = w_mac.array(mfi);); AMREX_D_TERM(Array4 gumac;, Array4 gvmac;,Array4 gwmac;); AMREX_D_TERM(if(gu_mac) gumac = gu_mac->array(mfi);, - if(gv_mac) gvmac = gv_mac->array(mfi);, - if(gw_mac) gwmac = gw_mac->array(mfi);); + if(gv_mac) gvmac = gv_mac->array(mfi);, + if(gw_mac) gwmac = gw_mac->array(mfi);); ParallelFor(xbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { Real kaf; - if(kap(i,j,k,0)!=VOF_NODATA && kap(i-1,j,k,0)!=VOF_NODATA) - kaf=Real(0.5)*(kap(i,j,k,0)+kap(i-1,j,k,0)); - else if (kap(i,j,k,0)!=VOF_NODATA) - kaf=kap(i,j,k); - else if (kap(i-1,j,k,0)!=VOF_NODATA) + if(kap(i,j,k,0)!=VOF_NODATA && kap(i-1,j,k,0)!=VOF_NODATA) + kaf=Real(0.5)*(kap(i,j,k,0)+kap(i-1,j,k,0)); + else if (kap(i,j,k,0)!=VOF_NODATA) + kaf=kap(i,j,k); + else if (kap(i-1,j,k,0)!=VOF_NODATA) kaf=kap(i-1,j,k); - else + else kaf=0.; - // density estimated at the face center, time increment is half of the timestep, i.e., dt/2 - //Print()< center1{AMREX_D_DECL((problo[0]+.5), (problo[1]+.75), - (problo[2]+.35))}; + (problo[2]+.35))}; 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); - - Array radii{AMREX_D_DECL(.2, .3, .25)}; - EB2::EllipsoidIF my_ellipsoid(radii, center, fluid_is_inside); + EB2::SphereIF my_sphere1(radius, center1, fluid_is_inside); + + Array radii{AMREX_D_DECL(.2, .3, .25)}; + EB2::EllipsoidIF my_ellipsoid(radii, center, fluid_is_inside); // Initialise cylinder parameters int direction = 2; @@ -2619,7 +2619,7 @@ VolumeOfFluid::tracer_vof_init_fraction (int lev, MultiFab& a_tracer) (0.5*(problo[2]+probhi[2])-.2))}; Array high{AMREX_D_DECL((0.5*(problo[0]+probhi[0])+.2), (0.5*(problo[1]+probhi[1])+.2), - (0.5*(problo[2]+probhi[2])+.2))}; + (0.5*(problo[2]+probhi[2])+.2))}; auto my_box= EB2::BoxIF( low, high, fluid_is_inside); //auto my_box= EB2::rotate(EB2::BoxIF( low, high, fluid_is_inside), .3, 1); auto my_box1= EB2::rotate(my_box, .3, 0); @@ -2630,7 +2630,7 @@ VolumeOfFluid::tracer_vof_init_fraction (int lev, MultiFab& a_tracer) // Generate GeometryShop //auto gshop = EB2::makeShop(two); //auto gshop = EB2::makeShop(my_box); - auto gshop = EB2::makeShop(my_box); + auto gshop = EB2::makeShop(my_box); //auto gshop = EB2::makeShop(my_cyl); int max_level = v_incflo->maxLevel(); EB2::Build(gshop, v_incflo->Geom(max_level), max_level, max_level); @@ -2644,18 +2644,18 @@ VolumeOfFluid::tracer_vof_init_fraction (int lev, MultiFab& a_tracer) if (lev == v_incflo->finestLevel()) { EB2::IndexSpace::pop(); } - } - else + } + else #endif { 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 @@ -2665,18 +2665,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-32.dat"; + std::string filename = "vof_value-32.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)) { @@ -2686,37 +2686,37 @@ 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] + + 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] #if AMREX_SPACEDIM==2 &&vout[n].k==7) -#else - &&vout[n].k>=vbx.smallEnd()[2]&&vout[n].k<=vbx.bigEnd()[2]) +#else + &&vout[n].k>=vbx.smallEnd()[2]&&vout[n].k<=vbx.bigEnd()[2]) #endif - { - //if(vbx.contains(vout[n].i,vout[n].j,vout[n].k)){ -#if AMREX_SPACEDIM==2 + { + //if(vbx.contains(vout[n].i,vout[n].j,vout[n].k)){ +#if AMREX_SPACEDIM==2 tracer(vout[n].i,vout[n].j,0,0)=vout[n].vof; -#else - tracer(vout[n].i,vout[n].j,vout[n].k,0)=vout[n].vof; -#endif - } - - /*amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept +#else + tracer(vout[n].i,vout[n].j,vout[n].k,0)=vout[n].vof; +#endif + } + + /*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]; @@ -2742,15 +2742,15 @@ 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. tracer_vof_update(lev, a_tracer, height[lev]); - curvature_calculation(lev, a_tracer, height[lev], kappa[lev]); - int n_tag=domain_tag_droplets (finest_level, v_incflo->grids,v_incflo->geom, - v_incflo->get_tracer_new (),GetVecOfPtrs(tag)); + curvature_calculation(lev, a_tracer, height[lev], kappa[lev]); + int n_tag=domain_tag_droplets (finest_level, v_incflo->grids,v_incflo->geom, + v_incflo->get_tracer_new (),GetVecOfPtrs(tag)); } @@ -2798,11 +2798,11 @@ 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, k=0; -#if AMREX_SPACEDIM==3 - for (k = lo.z; k <= hi.z; ++k) +#if AMREX_SPACEDIM==3 + for (k = lo.z; k <= hi.z; ++k) #endif for (int j = lo.y; j <= hi.y; ++j) { for (int i = lo.x; i <= hi.x; ++i) { @@ -2822,7 +2822,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 " <<"("< 0) { // std::ofstream TecplotFile; @@ -2840,30 +2840,30 @@ VolumeOfFluid::write_tecplot_surface(Real time, int nstep) TecplotFile << (AMREX_SPACEDIM== 2 ? "VARIABLES = \"X\", \"Y\"":"VARIABLES = \"X\", \"Y\", \"Z\""); //output variables TecplotFile <<", \"F\""<<", \"m_x\""<<", \"m_y\"" -#if AMREX_SPACEDIM==3 - <<", \"m_z\"" +#if AMREX_SPACEDIM==3 + <<", \"m_z\"" #endif - <<", \"alpha\""<<", \"tag\""<<", \"kappa\""<<"\n"; + <<", \"alpha\""<<", \"tag\""<<", \"kappa\""<<"\n"; std::string zonetitle=("Level_"+std::to_string(lev)+ "_Box_" +std::to_string(mfi.index())+ "_Proc_"+std::to_string(myproc)+"_step_"+std::to_string(nstep)); TecplotFile <<(std::string("ZONE T=")+zonetitle); - TecplotFile <<", DATAPACKING=POINT"<<", NODES="<m_use_cc_proj; + bool m_use_cc_proj=v_incflo->m_use_cc_proj; auto& ld = *v_incflo->m_leveldata[lev]; Geometry const& geom = v_incflo->Geom(lev); Box const& domain = geom.Domain(); @@ -2954,35 +2954,35 @@ 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); TecplotFile <<", DATAPACKING=BLOCK"<<", VARLOCATION=(["<<(m_use_cc_proj?AMREX_SPACEDIM+1:AMREX_SPACEDIM+2)<<"-" - <<(AMREX_SPACEDIM==3?24:19)<<"]=CELLCENTERED)" + <<(AMREX_SPACEDIM==3?24:19)<<"]=CELLCENTERED)" <<", SOLUTIONTIME="< const& pa_nd = ld.p_nd.const_array(mfi); - Array4 const& pa_cc = ld.p_cc.const_array(mfi); - Array4 const& pa_mac = ld.mac_phi.const_array(mfi); + Array4 const& pa_cc = ld.p_cc.const_array(mfi); + Array4 const& pa_mac = ld.mac_phi.const_array(mfi); Array4 const& tracer = ld.tracer.const_array(mfi); Array4 const& vel = ld.velocity.const_array(mfi); 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& density_arr = ld.density.const_array(mfi); - Array4 const& force_arr = force[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& density_arr = ld.density.const_array(mfi); + Array4 const& force_arr = force[lev].const_array(mfi); int nn=0, k=0; //write coordinate variables for (int dim = 0; dim < AMREX_SPACEDIM; ++dim) { -#if AMREX_SPACEDIM==3 +#if AMREX_SPACEDIM==3 for (k = lo.z; k <= hi.z +1; ++k) -#endif +#endif for (int j = lo.y; j <= hi.y +1; ++j) { for (int i = lo.x; i <= hi.x +1; ++i) { TecplotFile << (problo[dim]+dx[dim]*(dim<1?i:dim<2?j:k))<<" "; @@ -2993,13 +2993,13 @@ void VolumeOfFluid::WriteTecPlotFile(Real time, int nstep) } } } - + }// //write presure - int nt=m_use_cc_proj?0:1; -#if AMREX_SPACEDIM==3 + int nt=m_use_cc_proj?0:1; +#if AMREX_SPACEDIM==3 for (k = lo.z; k <= hi.z+nt; ++k) -#endif +#endif for (int j = lo.y; j <= hi.y+nt; ++j) { for (int i = lo.x; i <= hi.x+nt; ++i) { TecplotFile << (m_use_cc_proj?pa_cc(i,j,k):pa_nd(i,j,k))<<" "; @@ -3010,11 +3010,11 @@ void VolumeOfFluid::WriteTecPlotFile(Real time, int nstep) } } } - - //write VOF -#if AMREX_SPACEDIM==3 - for (k = lo.z; k <= hi.z; ++k) -#endif + + //write VOF +#if AMREX_SPACEDIM==3 + for (k = lo.z; k <= hi.z; ++k) +#endif for (int j = lo.y; j <= hi.y; ++j) { for (int i = lo.x; i <= hi.x; ++i) { TecplotFile << tracer(i,j,k,0)<<" "; @@ -3025,12 +3025,12 @@ void VolumeOfFluid::WriteTecPlotFile(Real time, int nstep) } } } - + //write velocity for (int dim = 0; dim < AMREX_SPACEDIM; ++dim) { -#if AMREX_SPACEDIM==3 - for (k = lo.z; k <= hi.z; ++k) -#endif +#if AMREX_SPACEDIM==3 + for (k = lo.z; k <= hi.z; ++k) +#endif for (int j = lo.y; j <= hi.y; ++j) { for (int i = lo.x; i <= hi.x; ++i) { TecplotFile << vel(i,j,k,dim)<<" "; @@ -3041,14 +3041,14 @@ void VolumeOfFluid::WriteTecPlotFile(Real time, int nstep) } } } - + }// //write variables of the normal direction of the interface for (int dim = 0; dim < AMREX_SPACEDIM; ++dim) { #if AMREX_SPACEDIM==3 - for (k = lo.z; k <= hi.z; ++k) -#endif + for (k = lo.z; k <= hi.z; ++k) +#endif for (int j = lo.y; j <= hi.y; ++j) { for (int i = lo.x; i <= hi.x; ++i) { TecplotFile << mv(i,j,k,dim)<<" "; @@ -3059,13 +3059,13 @@ void VolumeOfFluid::WriteTecPlotFile(Real time, int nstep) } } } - + }// //write alpha of the interface #if AMREX_SPACEDIM==3 - for (k = lo.z; k <= hi.z; ++k) -#endif + for (k = lo.z; k <= hi.z; ++k) +#endif for (int j = lo.y; j <= hi.y; ++j) { for (int i = lo.x; i <= hi.x; ++i) { TecplotFile << al(i,j,k)<<" "; @@ -3076,12 +3076,12 @@ void VolumeOfFluid::WriteTecPlotFile(Real time, int nstep) } } } - + //write id of the droplets or bubbles -#if AMREX_SPACEDIM==3 - for (k = lo.z; k <= hi.z; ++k) -#endif +#if AMREX_SPACEDIM==3 + for (k = lo.z; k <= hi.z; ++k) +#endif for (int j = lo.y; j <= hi.y; ++j) { for (int i = lo.x; i <= hi.x; ++i) { TecplotFile << tag_arr(i,j,k)<<" "; @@ -3092,14 +3092,14 @@ void VolumeOfFluid::WriteTecPlotFile(Real time, int nstep) } } } - + //write height function values for (int dim = 0; dim < AMREX_SPACEDIM; ++dim) { -#if AMREX_SPACEDIM==3 - for (k = lo.z; k <= hi.z; ++k) +#if AMREX_SPACEDIM==3 + for (k = lo.z; k <= hi.z; ++k) #endif 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) { @@ -3108,12 +3108,12 @@ void VolumeOfFluid::WriteTecPlotFile(Real time, int nstep) } } } - }// + }// for (int dim = 0; dim < AMREX_SPACEDIM; ++dim) { -#if AMREX_SPACEDIM==3 +#if AMREX_SPACEDIM==3 for (k = lo.z; k <= hi.z; ++k) -#endif +#endif for (int j = lo.y; j <= hi.y; ++j) { for (int i = lo.x; i <= hi.x; ++i) { TecplotFile << ht_arr(i,j,k,dim)<<" "; @@ -3124,11 +3124,11 @@ void VolumeOfFluid::WriteTecPlotFile(Real time, int nstep) } } } - }// + }// //write curvature -#if AMREX_SPACEDIM==3 +#if AMREX_SPACEDIM==3 for (k = lo.z; k <= hi.z; ++k) -#endif +#endif for (int j = lo.y; j <= hi.y; ++j) { for (int i = lo.x; i <= hi.x; ++i) { TecplotFile << kappa_arr(i,j,k)<<" "; @@ -3139,11 +3139,11 @@ void VolumeOfFluid::WriteTecPlotFile(Real time, int nstep) } } } - + //write density -#if AMREX_SPACEDIM==3 +#if AMREX_SPACEDIM==3 for (k = lo.z; k <= hi.z; ++k) -#endif +#endif for (int j = lo.y; j <= hi.y; ++j) { for (int i = lo.x; i <= hi.x; ++i) { TecplotFile << density_arr(i,j,k)<<" "; @@ -3154,13 +3154,13 @@ void VolumeOfFluid::WriteTecPlotFile(Real time, int nstep) } } } - + //write force vector for (int dim = 0; dim < AMREX_SPACEDIM; ++dim) { -#if AMREX_SPACEDIM==3 +#if AMREX_SPACEDIM==3 for (k = lo.z; k <= hi.z; ++k) -#endif +#endif for (int j = lo.y; j <= hi.y; ++j) { for (int i = lo.x; i <= hi.x; ++i) { TecplotFile << force_arr(i,j,k,dim)<<" "; @@ -3171,9 +3171,9 @@ void VolumeOfFluid::WriteTecPlotFile(Real time, int nstep) } } } - + }// - + TecplotFile <<"\n"; } // end MFIter @@ -3348,7 +3348,7 @@ int VolumeOfFluid::domain_tag_droplets (int finest_level, Vector cons ,ort2=ORTHOGONAL_COMPONENT(ort1);//2nd transverse direction for (int n=0;n<2;n++){ int k0=dim_limit[n][d],gd=k0+(n==0?-1:1); - for(int i0=ijk_min[ort1];i0<=ijk_max[ort1];i0++){ + for(int i0=ijk_min[ort1];i0<=ijk_max[ort1];i0++){ #if AMREX_SPACEDIM==2 /*2D*/ Real tag_cell=(d==0?tag_arr(k0,i0,0):tag_arr(i0,k0,0)); if(tag_cell > 0){ @@ -3356,7 +3356,7 @@ int VolumeOfFluid::domain_tag_droplets (int finest_level, Vector cons if(tag_gcell > 0) touching_regions (tag_cell, tag_gcell, touch); } -#else /*3D */ +#else /*3D */ for(int j0=ijk_min[ort2];j0<=ijk_max[ort2];j0++){ Real tag_cell=(d==0?tag_arr(k0,i0,j0): d==1?tag_arr(j0,k0,i0): @@ -3369,8 +3369,8 @@ int VolumeOfFluid::domain_tag_droplets (int finest_level, Vector cons touching_regions (tag_cell, tag_gcell, touch); } }// end for-loop for searching cells in the boundaries. -#endif - } +#endif + } }// end for-loop for low and high boundary }// end for-loop for AMREX_SPACEDIM } @@ -3475,7 +3475,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]); @@ -3496,7 +3496,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 // { @@ -3536,15 +3536,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){ @@ -3565,11 +3565,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]); + } } //////////////////////////////////////////////////////////////////// ////// @@ -3585,13 +3585,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"<