diff --git a/Source/Parallelization/WarpXRegrid.cpp b/Source/Parallelization/WarpXRegrid.cpp index dcea3cab58a..1bd05a2d115 100644 --- a/Source/Parallelization/WarpXRegrid.cpp +++ b/Source/Parallelization/WarpXRegrid.cpp @@ -334,6 +334,7 @@ WarpX::RemakeLevel (int lev, Real /*time*/, const BoxArray& ba, const Distributi // we can avoid redistributing these since we immediately re-build the values via BuildBufferMasks() RemakeMultiFab(current_buffer_masks[lev], dm, false ,lev); RemakeMultiFab(gather_buffer_masks[lev], dm, false ,lev); + RemakeMultiFab(interp_weight_gbuffer[lev], dm, false, lev); if (current_buffer_masks[lev] || gather_buffer_masks[lev]) { BuildBufferMasks(); diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index b77c1147a15..747d251fdcc 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -326,6 +326,28 @@ public: amrex::FArrayBox const * & ez_ptr, amrex::FArrayBox const * & bx_ptr, amrex::FArrayBox const * & by_ptr, amrex::FArrayBox const * & bz_ptr); + void InterpolateFieldsInGatherBuffer (const amrex::Box& box, + amrex::FArrayBox& bufferEx, amrex::FArrayBox& bufferEy, amrex::FArrayBox& bufferEz, + amrex::FArrayBox& bufferBx, amrex::FArrayBox& bufferBy, amrex::FArrayBox& bufferBz, + amrex::Elixir& buf_exeli, amrex::Elixir& buf_eyeli, amrex::Elixir& buf_ezeli, + amrex::Elixir& buf_bxeli, amrex::Elixir& buf_byeli, amrex::Elixir& buf_bzeli, + amrex::FArrayBox const *& ex_ptr, amrex::FArrayBox const *& ey_ptr, amrex::FArrayBox const* & ez_ptr, + amrex::FArrayBox const *& bx_ptr, amrex::FArrayBox const *& by_ptr, amrex::FArrayBox const* & bz_ptr, + amrex::FArrayBox const *& cex_ptr, amrex::FArrayBox const *& cey_ptr, amrex::FArrayBox const* & cez_ptr, + amrex::FArrayBox const *& cbx_ptr, amrex::FArrayBox const *& cby_ptr, amrex::FArrayBox const* & cbz_ptr, + const amrex::IArrayBox& gather_mask, const amrex::FArrayBox& weight_gbuffer, + const amrex::IntVect ref_ratio); + + void WeightedSumInGatherBuffer (const amrex::Box& tmp_xbox, const amrex::Box& tmp_ybox, const amrex::Box& tmp_zbox, + amrex::Array4 const & xbuf_field_arr, amrex::Array4 const& ybuf_field_arr, + amrex::Array4 const & zbuf_field_arr, + amrex::Array4 cfx_arr, amrex::Array4 cfy_arr, + amrex::Array4 cfz_arr, + amrex::Array4 fx_arr, amrex::Array4 fy_arr, + amrex::Array4 fz_arr, + amrex::Array4 gm_arr, amrex::Array4 wt_arr, + const amrex::IntVect ref_ratio); + /** * \brief This function determines if resampling should be done for the current species, and * if so, performs the resampling. diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 3d5cbe0007d..9bad2b0d558 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2016,12 +2016,11 @@ PhysicalParticleContainer::Evolve (int lev, WARPX_PROFILE_VAR_NS("PhysicalParticleContainer::Evolve::GatherAndPush", blp_fg); BL_ASSERT(OnSameGrids(lev,jx)); - amrex::LayoutData* cost = WarpX::getCosts(lev); const iMultiFab* current_masks = WarpX::CurrentBufferMasks(lev); const iMultiFab* gather_masks = WarpX::GatherBufferMasks(lev); - + const amrex::MultiFab* weight_gbuffer = WarpX::GetInstance().getInterpWeightGBuffer(lev); const bool has_buffer = cEx || cjx; if (m_do_back_transformed_particles) @@ -2050,7 +2049,8 @@ PhysicalParticleContainer::Evolve (int lev, FArrayBox filtered_Ex, filtered_Ey, filtered_Ez; FArrayBox filtered_Bx, filtered_By, filtered_Bz; - + FArrayBox bufferEx, bufferEy, bufferEz; + FArrayBox bufferBx, bufferBy, bufferBz; for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers) @@ -2079,6 +2079,7 @@ PhysicalParticleContainer::Evolve (int lev, FArrayBox const* bzfab = &Bz[pti]; Elixir exeli, eyeli, ezeli, bxeli, byeli, bzeli; + Elixir buf_exeli, buf_eyeli, buf_ezeli, buf_bxeli, buf_byeli, buf_bzeli; if (WarpX::use_fdtd_nci_corr) { @@ -2120,7 +2121,7 @@ PhysicalParticleContainer::Evolve (int lev, DepositCharge(pti, wp, ion_lev, rho, 0, 0, np_current, thread_num, lev, lev); - if (has_buffer){ + if ((np-np_current)> 0 ){ DepositCharge(pti, wp, ion_lev, crho, 0, np_current, np-np_current, thread_num, lev, lev-1); } @@ -2163,8 +2164,19 @@ PhysicalParticleContainer::Evolve (int lev, FArrayBox const* cbyfab = &(*cBy)[pti]; FArrayBox const* cbzfab = &(*cBz)[pti]; + //// buffer box (Both 2D and 3D) + InterpolateFieldsInGatherBuffer(box, bufferEx, bufferEy, bufferEz, + bufferBx, bufferBy, bufferBz, + buf_exeli, buf_eyeli, buf_ezeli, + buf_bxeli, buf_byeli, buf_bzeli, + exfab, eyfab, ezfab, bxfab, byfab, bzfab, + cexfab, ceyfab, cezfab, cbxfab, cbyfab, cbzfab, + (*gather_masks)[pti], (*weight_gbuffer)[pti], + ref_ratio); + if (WarpX::use_fdtd_nci_corr) { + // should this be bufferEFields* // Filter arrays (*cEx)[pti], store the result in // filtered_Ex and update pointer cexfab so that it // points to filtered_Ex (and do the same for all @@ -2176,15 +2188,17 @@ PhysicalParticleContainer::Evolve (int lev, (*cBx)[pti], (*cBy)[pti], (*cBz)[pti], cexfab, ceyfab, cezfab, cbxfab, cbyfab, cbzfab); } - // Field gather and push for particles in gather buffers e_is_nodal = cEx->is_nodal() and cEy->is_nodal() and cEz->is_nodal(); if (push_type == PushType::Explicit) { - PushPX(pti, cexfab, ceyfab, cezfab, - cbxfab, cbyfab, cbzfab, + //PushPX(pti, cexfab, ceyfab, cezfab, + // cbxfab, cbyfab, cbzfab, + PushPX(pti, &bufferEx, &bufferEy, &bufferEz, + &bufferBx, &bufferBy, &bufferBz, cEx->nGrowVect(), e_is_nodal, nfine_gather, np-nfine_gather, - lev, lev-1, dt, ScaleFields(false), a_dt_type); + lev, gather_lev, dt, ScaleFields(false), a_dt_type); + //lev, lev-1, dt, ScaleFields(false), a_dt_type); } else if (push_type == PushType::Implicit) { ImplicitPushXP(pti, cexfab, ceyfab, cezfab, cbxfab, cbyfab, cbzfab, @@ -2193,7 +2207,6 @@ PhysicalParticleContainer::Evolve (int lev, lev, lev-1, dt, ScaleFields(false), a_dt_type); } } - WARPX_PROFILE_VAR_STOP(blp_fg); // Current Deposition @@ -2210,7 +2223,7 @@ PhysicalParticleContainer::Evolve (int lev, 0, np_current, thread_num, lev, lev, dt, relative_time, push_type); - if (has_buffer) + if ((np-np_current)>0) { // Deposit in buffers DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, cjx, cjy, cjz, @@ -2230,7 +2243,7 @@ PhysicalParticleContainer::Evolve (int lev, DepositCharge(pti, wp, ion_lev, rho, 1, 0, np_current, thread_num, lev, lev); - if (has_buffer){ + if ((np-np_current)>0 ){ DepositCharge(pti, wp, ion_lev, crho, 1, np_current, np-np_current, thread_num, lev, lev-1); } @@ -2333,6 +2346,113 @@ PhysicalParticleContainer::applyNCIFilter ( #endif } +void +PhysicalParticleContainer::InterpolateFieldsInGatherBuffer (const Box& box, + amrex::FArrayBox& bufferEx, amrex::FArrayBox& bufferEy, amrex::FArrayBox& bufferEz, + amrex::FArrayBox& bufferBx, amrex::FArrayBox& bufferBy, amrex::FArrayBox& bufferBz, + Elixir& buf_exeli, Elixir& buf_eyeli, Elixir& buf_ezeli, + Elixir& buf_bxeli, Elixir& buf_byeli, Elixir& buf_bzeli, + FArrayBox const *& ex_ptr, FArrayBox const *& ey_ptr, FArrayBox const* & ez_ptr, + FArrayBox const *& bx_ptr, FArrayBox const *& by_ptr, FArrayBox const* & bz_ptr, + FArrayBox const *& cex_ptr, FArrayBox const *& cey_ptr, FArrayBox const* & cez_ptr, + FArrayBox const *& cbx_ptr, FArrayBox const *& cby_ptr, FArrayBox const* & cbz_ptr, + const IArrayBox& gather_mask, const FArrayBox& weight_gbuffer, + const amrex::IntVect ref_ratio) +{ + amrex::Box tmp_exbox = amrex::convert(box,ex_ptr->box().ixType()); + amrex::Box tmp_eybox = amrex::convert(box,ey_ptr->box().ixType()); + amrex::Box tmp_ezbox = amrex::convert(box,ez_ptr->box().ixType()); + + bufferEx.resize(tmp_exbox); + buf_exeli = bufferEx.elixir(); + bufferEy.resize(tmp_eybox); + buf_eyeli = bufferEy.elixir(); + bufferEz.resize(tmp_ezbox); + buf_ezeli = bufferEz.elixir(); + WeightedSumInGatherBuffer(tmp_exbox, tmp_eybox, tmp_ezbox, + bufferEx.array(), bufferEy.array(), bufferEz.array(), + cex_ptr->array(), cey_ptr->array(), cez_ptr->array(), + ex_ptr->array(), ey_ptr->array(), ez_ptr->array(), + gather_mask.array(), weight_gbuffer.array(), ref_ratio); + + amrex::Box tmp_bxbox = amrex::convert(box,bx_ptr->box().ixType()); + amrex::Box tmp_bybox = amrex::convert(box,by_ptr->box().ixType()); + amrex::Box tmp_bzbox = amrex::convert(box,bz_ptr->box().ixType()); + bufferBx.resize(tmp_bxbox); + buf_bxeli = bufferBx.elixir(); + bufferBy.resize(tmp_bybox); + buf_byeli = bufferBy.elixir(); + bufferBz.resize(tmp_bzbox); + buf_bzeli = bufferBz.elixir(); + + WeightedSumInGatherBuffer(tmp_bxbox, tmp_bybox, tmp_bzbox, + bufferBx.array(), bufferBy.array(), bufferBz.array(), + cbx_ptr->array(), cby_ptr->array(), cbz_ptr->array(), + bx_ptr->array(), by_ptr->array(), bz_ptr->array(), + gather_mask.array(), weight_gbuffer.array(), ref_ratio); +} + +void +PhysicalParticleContainer::WeightedSumInGatherBuffer ( + const amrex::Box& tmp_xbox, const amrex::Box& tmp_ybox, const amrex::Box& tmp_zbox, + amrex::Array4 const & xbuf_field_arr, amrex::Array4 const& ybuf_field_arr, + amrex::Array4 const & zbuf_field_arr, + amrex::Array4 cfx_arr, amrex::Array4 cfy_arr, + amrex::Array4 cfz_arr, + amrex::Array4 fx_arr, amrex::Array4 fy_arr, + amrex::Array4 fz_arr, + amrex::Array4 gm_arr, amrex::Array4 wt_arr, + const amrex::IntVect ref_ratio) +{ + amrex::ParallelFor( tmp_xbox, tmp_ybox, tmp_zbox, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + int ii = amrex::coarsen(i,ref_ratio[0]); + int jj = j; + int kk = k; +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + jj = amrex::coarsen(j,ref_ratio[1]); +#elif defined(WARPX_DIM_3D) + kk = amrex::coarsen(k,ref_ratio[2]); +#endif + if (gm_arr(i,j,k) == 0) { + xbuf_field_arr(i,j,k) = wt_arr(i,j,k)*fx_arr(i,j,k) + (1._rt-wt_arr(i,j,k))*cfx_arr(ii,jj,kk); + } else { + xbuf_field_arr(i,j,k) = fx_arr(i,j,k); + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + int ii = amrex::coarsen(i,ref_ratio[0]); + int jj = j; + int kk = k; +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + jj = amrex::coarsen(j,ref_ratio[1]); +#elif defined(WARPX_DIM_3D) + kk = amrex::coarsen(k,ref_ratio[2]); +#endif + if (gm_arr(i,j,k) == 0) { + ybuf_field_arr(i,j,k) = wt_arr(i,j,k)*fy_arr(i,j,k) + (1._rt-wt_arr(i,j,k))*cfy_arr(ii,jj,kk); + } else { + ybuf_field_arr(i,j,k) = fy_arr(i,j,k); + } + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + int ii = amrex::coarsen(i,ref_ratio[0]); + int jj = j; + int kk = k; +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + jj = amrex::coarsen(j,ref_ratio[1]); +#elif defined(WARPX_DIM_3D) + kk = amrex::coarsen(k,ref_ratio[2]); +#endif + if (gm_arr(i,j,k) == 0) { + zbuf_field_arr(i,j,k) = wt_arr(i,j,k)*fz_arr(i,j,k) + (1._rt-wt_arr(i,j,k))*cfz_arr(ii,jj,kk); + } else { + zbuf_field_arr(i,j,k) = fz_arr(i,j,k); + } + } + ); +} + // Loop over all particles in the particle container and // split particles tagged with p.id()=DoSplitParticleID void diff --git a/Source/WarpX.H b/Source/WarpX.H index d3900c91593..0d31e0e765a 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -355,6 +355,32 @@ public: //! #n_field_gather_buffer cells of the edge of the patch, will gather the fields //! from the lower refinement level instead of the refinement patch itself static int n_field_gather_buffer; + + /** + * If true, particles located inside refinement patch but within + * #n_field_gather_buffer cells of the patch edge will gather fields + * that are a weighted sum of the fields from lower refinement + * level and the refinement patch. By default, no interpolation will be done, + * and the fields will be gathered only from the lower refinement level. + * The option can be used to allow for a smooth transition from the fine + * and coarse solutions with a tanh profile, such that, the value + * of the gather fields is close to the solution on the refinement patch + * near the gather buffer edge, and close to the solution on the lower + * refinement level close to the edge of the patch. + * #tanh_midpoint_gather_buffer can be used to set the tanh profile + */ + static bool do_fieldinterp_gather_buffer; + + /** + * With mesh-refinement, finite number of gather buffer cells, and + * do_fieldinterp_gather_buffer, the particles on the fine patch + * gather a weighted sum of fields from fine patch and coarse patch. + * A tanh profile is used to set the smoothing kernel + * #tanh_midpoint_gather_buffer must be a fraction of the buffer width + * measured from coarsefine interface, where the tanh profile is 0.5 + * The default value is 0.5 + */ + static amrex::Real tanh_midpoint_gather_buffer; //! With mesh refinement, particles located inside a refinement patch, but within //! #n_current_deposition_buffer cells of the edge of the patch, will deposit their charge //! and current onto the lower refinement level instead of the refinement patch itself @@ -878,6 +904,10 @@ public: static std::array CellSize (int lev); static amrex::RealBox getRealBox(const amrex::Box& bx, int lev); + [[nodiscard]] const amrex::MultiFab* getInterpWeightGBuffer (int lev) const + { + return interp_weight_gbuffer[lev].get(); + } /** * \brief Return the lower corner of the box in real units. * \param bx The box @@ -1391,6 +1421,7 @@ private: return gather_buffer_masks[lev].get(); } + /** * \brief Re-orders the Fornberg coefficients so that they can be used more conveniently for * finite-order centering operations. For example, for finite-order centering of order 6, @@ -1564,6 +1595,11 @@ private: amrex::Vector, 3 > > Bfield_cax; amrex::Vector > current_buffer_masks; amrex::Vector > gather_buffer_masks; + /** Multifab that stores weights for interpolating fine and coarse solutions + * in the buffer gather region, such that, the weight for fine solution is + * 0 near the coarse-fine interface, and 1 as the buffer gather region terminates in the fine patch + */ + amrex::Vector > interp_weight_gbuffer; // If charge/current deposition buffers are used amrex::Vector, 3 > > current_buf; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index b13d7aaf63b..92ea93eeb4b 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -207,6 +207,8 @@ IntVect WarpX::filter_npass_each_dir(1); int WarpX::n_field_gather_buffer = -1; int WarpX::n_current_deposition_buffer = -1; +bool WarpX::do_fieldinterp_gather_buffer = false; +amrex::Real WarpX::tanh_midpoint_gather_buffer = 0.5; short WarpX::grid_type; amrex::IntVect m_rho_nodal_flag; @@ -399,6 +401,7 @@ WarpX::WarpX () gather_buffer_masks.resize(nlevs_max); current_buf.resize(nlevs_max); charge_buf.resize(nlevs_max); + interp_weight_gbuffer.resize(nlevs_max); pml.resize(nlevs_max); #if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD) @@ -873,6 +876,10 @@ WarpX::ReadParameters () pp_warpx, "n_field_gather_buffer", n_field_gather_buffer); utils::parser::queryWithParser( pp_warpx, "n_current_deposition_buffer", n_current_deposition_buffer); + utils::parser::queryWithParser( + pp_warpx, "do_fieldinterp_gather_buffer", do_fieldinterp_gather_buffer); + utils::parser::queryWithParser( + pp_warpx, "tanh_midpoint_gather_buffer", tanh_midpoint_gather_buffer); amrex::Real quantum_xi_tmp; const auto quantum_xi_is_specified = @@ -2052,6 +2059,7 @@ WarpX::ClearLevel (int lev) current_buffer_masks[lev].reset(); gather_buffer_masks[lev].reset(); + interp_weight_gbuffer[lev].reset(); F_fp [lev].reset(); G_fp [lev].reset(); @@ -2673,20 +2681,23 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm AllocInitMultiFab(Efield_cax[lev][1], amrex::convert(cba,Ey_nodal_flag),dm,ncomps,ngEB,lev, "Efield_cax[y]"); AllocInitMultiFab(Efield_cax[lev][2], amrex::convert(cba,Ez_nodal_flag),dm,ncomps,ngEB,lev, "Efield_cax[z]"); } - - AllocInitMultiFab(gather_buffer_masks[lev], ba, dm, ncomps, amrex::IntVect(1), lev, "gather_buffer_masks"); + amrex::Print() << " ba for allocating gahter buffer mask " << "\n"; // Gather buffer masks have 1 ghost cell, because of the fact // that particles may move by more than one cell when using subcycling. + AllocInitMultiFab(gather_buffer_masks[lev], ba, dm, ncomps, amrex::IntVect(1), lev, "gather_buffer_masks"); + AllocInitMultiFab(interp_weight_gbuffer[lev], amrex::convert(ba, IntVect::TheNodeVector()), dm, ncomps, ngEB, lev, "interp_weight_gbuffer", 0.0_rt); } if (n_current_deposition_buffer > 0) { + amrex::Print() << " current dep buffer : " << n_current_deposition_buffer << "\n"; AllocInitMultiFab(current_buf[lev][0], amrex::convert(cba,jx_nodal_flag),dm,ncomps,ngJ,lev, "current_buf[x]"); AllocInitMultiFab(current_buf[lev][1], amrex::convert(cba,jy_nodal_flag),dm,ncomps,ngJ,lev, "current_buf[y]"); AllocInitMultiFab(current_buf[lev][2], amrex::convert(cba,jz_nodal_flag),dm,ncomps,ngJ,lev, "current_buf[z]"); if (rho_cp[lev]) { AllocInitMultiFab(charge_buf[lev], amrex::convert(cba,rho_nodal_flag),dm,2*ncomps,ngRho,lev, "charge_buf"); } - AllocInitMultiFab(current_buffer_masks[lev], ba, dm, ncomps, amrex::IntVect(1), lev, "current_buffer_masks"); + amrex::Print() << " allocate current buffer mask \n"; + AllocInitMultiFab(current_buffer_masks[lev], ba, dm, ncomps, ngJ, lev, "current_buffer_masks"); // Current buffer masks have 1 ghost cell, because of the fact // that particles may move by more than one cell when using subcycling. } @@ -3048,6 +3059,8 @@ WarpX::getLoadBalanceEfficiency (const int lev) void WarpX::BuildBufferMasks () { + bool do_interpolate = WarpX::do_fieldinterp_gather_buffer; + amrex::Real tanh_midpoint = WarpX::tanh_midpoint_gather_buffer; for (int lev = 1; lev <= maxLevel(); ++lev) { for (int ipass = 0; ipass < 2; ++ipass) @@ -3073,6 +3086,75 @@ WarpX::BuildBufferMasks () const Box tbx = mfi.growntilebox(); BuildBufferMasksInBox( tbx, (*bmasks)[mfi], tmp[mfi], ngbuffer ); } + if (ipass==0) continue; + amrex::MultiFab* weight_gbuffer = interp_weight_gbuffer[lev].get(); + // Using tmp to also set weights in the interp_weight_gbuffer multifab + for (MFIter mfi(*weight_gbuffer, true); mfi.isValid(); ++mfi) + { + const Box& tbx = mfi.tilebox(IntVect::TheNodeVector(),weight_gbuffer->nGrowVect()); + auto const& gmsk = tmp[mfi].const_array(); + auto const& bmsk = (*bmasks)[mfi].array(); + auto const& wtmsk = (*weight_gbuffer)[mfi].array(); + amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) { + wtmsk(i,j,k) = 0._rt; + if (bmsk(i,j,k) == 0 && do_interpolate) { + if(gmsk(i,j,k)==0) { + wtmsk(i,j,k) = 0.; + return; + } + //for (int ii = i-1; ii >= i-ngbuffer; --ii) { + // if (gmsk(ii,j,k)==0) { + // amrex::Real arg = (static_cast(i-ii)-ngbuffer*tanh_midpoint) + // / ((1.-tanh_midpoint)*(ngbuffer/3.)); + // wtmsk(i,j,k) = std::tanh(arg)*0.5 + 0.5; + // amrex::Print() << " i edge wt is " << wtmsk(i,j,k) << "\n"; + // return; + // } + //} + //for (int ii = i+1; ii <= i+ngbuffer; ++ii) { + // if (gmsk(ii,j,k)==0) { + // amrex::Real arg = (static_cast(ii-i)-ngbuffer*tanh_midpoint) + // / ((1.-tanh_midpoint)*(ngbuffer/3.)); + // wtmsk(i,j,k) = std::tanh(arg)*0.5+0.5; + // amrex::Print() << " wt is " << wtmsk(i,j,k) << "\n"; + // return; + // } + //} + for (int jj = j-1; jj >= j-ngbuffer; --jj) { + if (gmsk(i,jj,k)==0) { + amrex::Real arg = (static_cast(j-jj)-ngbuffer*tanh_midpoint) + / ((1.-tanh_midpoint)*(ngbuffer/3.)); + wtmsk(i,j,k) = std::tanh(arg)*0.5 + 0.5; + return; + } + } + for (int jj = j+1; jj <= j+ngbuffer; ++jj) { + if (gmsk(i,jj,k)==0) { + amrex::Real arg = (static_cast(jj - j)-ngbuffer*tanh_midpoint) + / ((1.-tanh_midpoint)*(ngbuffer/3.)); + wtmsk(i,j,k) = std::tanh(arg)*0.5+0.5; + return; + } + } + //for (int kk = k-1; kk >= k-ngbuffer; --kk) { + // if (gmsk(i,j,kk)==0) { + // amrex::Real arg = (static_cast(k-kk)-ngbuffer*tanh_midpoint) + // / ((1.-tanh_midpoint)*(ngbuffer/3.)); + // wtmsk(i,j,k) = std::tanh(arg)*0.5+0.5; + // return; + // } + //} + //for (int kk = k+1; kk <= k+ngbuffer; ++kk) { + // if (gmsk(i,j,kk)==0) { + // amrex::Real arg = (static_cast(kk-k)-ngbuffer*tanh_midpoint) + // / ((1.-tanh_midpoint)*(ngbuffer/3.)); + // wtmsk(i,j,k) = std::tanh(arg)*0.5+0.5; + // return; + // } + //} + } + }); + } } } }