diff --git a/src/chemistry/bmx_chem_K.H b/src/chemistry/bmx_chem_K.H index 7b0176c..148304e 100644 --- a/src/chemistry/bmx_chem_K.H +++ b/src/chemistry/bmx_chem_K.H @@ -82,8 +82,9 @@ checkSplit(Real *p_par, int *p_ipar, Real max_vol, Real max_len, AMREX_GPU_HOST_DEVICE AMREX_INLINE void checkTipFusion(Real *rij, Real *i_par, int *i_ipar, Real *j_par, - int *j_ipar, Real* fpar, int me, amrex::RandomEngine const& engine) + int *j_ipar, Real* fpar, int me, int *fusing, amrex::RandomEngine const& engine) { + *fusing = 0; if (i_ipar[intIdx::cell_type] == cellType::FUNGI && i_ipar[intIdx::position] == siteLocation::TIP && i_ipar[intIdx::fuse_flag] == 0 && @@ -245,6 +246,7 @@ checkTipFusion(Real *rij, Real *i_par, int *i_ipar, Real *j_par, i_ipar[intIdx::seg3_id2] = j_ipar[intIdx::cpu]; } i_ipar[intIdx::n_bnds] = i_bnds+1; + *fusing = 1; } } } @@ -1363,6 +1365,7 @@ void xferMeshToParticleAndUpdateChem(Real grid_vol, Real npart, length_max = chempar[9]; // check to see if particle can grow. Set growth rates to zero // if it cannot +#if 0 if (cell_type == cellType::FUNGI) { if (cell_ipar[intIdx::position] != siteLocation::TIP) { if ((cell_par[realIdx::radius] >= radius_max && @@ -1373,6 +1376,15 @@ void xferMeshToParticleAndUpdateChem(Real grid_vol, Real npart, } } } +#else + if (cell_type == cellType::FUNGI) { + if (cell_par[realIdx::radius] >= radius_max && + cell_par[realIdx::c_length] >= length_max) { + kg = 0.0; + kv = 0.0; + } + } +#endif #if 0 // Stop growth on interior segments that are maximum sized if (cell_par[realIdx::c_length] == length_max && diff --git a/src/des/bmx_pc.H b/src/des/bmx_pc.H index 8ed1065..d0b8eb7 100644 --- a/src/des/bmx_pc.H +++ b/src/des/bmx_pc.H @@ -132,7 +132,7 @@ public: std::string& knapsack_weight_type, int& nsubsteps); - void EvaluateTipFusion(const amrex::Vector cost, + bool EvaluateTipFusion(const amrex::Vector cost, std::string& knapsack_weight_type); void EvaluateInteriorFusion(const amrex::Vector cost, diff --git a/src/des/bmx_pc_fusion.cpp b/src/des/bmx_pc_fusion.cpp index 03f783f..e928123 100644 --- a/src/des/bmx_pc_fusion.cpp +++ b/src/des/bmx_pc_fusion.cpp @@ -14,9 +14,10 @@ using namespace amrex; /******************************************************************************* * Check all tips to see if they fuse with a neighbor. * ******************************************************************************/ -void BMXParticleContainer::EvaluateTipFusion (const Vector cost, +bool BMXParticleContainer::EvaluateTipFusion (const Vector cost, std::string& knapsack_weight_type) { + bool ret = false; BL_PROFILE_REGION_START("bmx_dem::EvaluateTipFusion()"); BL_PROFILE("bmx_dem::EvaluateTipFusion()"); @@ -80,6 +81,9 @@ void BMXParticleContainer::EvaluateTipFusion (const Vector cost, auto& aos = ptile.GetArrayOfStructs(); ParticleType* pstruct = aos().dataPtr(); + Gpu::DeviceScalar fused_gpu(0); + int* fused = fused_gpu.dataPtr(); + const int nrp = GetParticles(lev)[index].numRealParticles(); // Number of particles including neighbor particles @@ -98,7 +102,7 @@ void BMXParticleContainer::EvaluateTipFusion (const Vector cost, // now we loop over the neighbor list and compute the forces int me = ParallelDescriptor::MyProc(); amrex::ParallelForRNG(nrp, - [nrp,pstruct,nbor_data,fpar,xpar,me] + [nrp,pstruct,nbor_data,fpar,xpar,fused,me] AMREX_GPU_DEVICE (int i, amrex::RandomEngine const& engine) noexcept // [=] AMREX_GPU_DEVICE (int i, amrex::RandomEngine const& engine) noexcept { @@ -119,6 +123,7 @@ void BMXParticleContainer::EvaluateTipFusion (const Vector cost, { auto p2 = *mit; const int j = mit.index(); + int fusing; Real dist_x = pos1[0] - p2.pos(0); Real dist_y = pos1[1] - p2.pos(1); @@ -145,13 +150,15 @@ void BMXParticleContainer::EvaluateTipFusion (const Vector cost, // whether it is fusing to p2. If fusion occurs, mark particle // as being fused to p2. checkTipFusion(&diff[0], &particle.rdata(0), &particle.idata(0), - &p2.rdata(0), &p2.idata(0), fpar, me, engine); + &p2.rdata(0), &p2.idata(0), fpar, me, &fusing, engine); + Gpu::Atomic::Max(fused, fusing); // TODO: Do we need an OPENMP pragma here? } } // end of neighbor loop }); // end of loop over particles + if (fused_gpu.dataValue() == 1) ret = true; amrex::Gpu::Device::synchronize(); @@ -192,6 +199,8 @@ void BMXParticleContainer::EvaluateTipFusion (const Vector cost, #endif BL_PROFILE_REGION_STOP("bmx_dem::EvaluateTipFusion()"); + ParallelDescriptor::ReduceBoolOr(ret); + return ret; } /******************************************************************************* diff --git a/src/timestepping/bmx_evolve.cpp b/src/timestepping/bmx_evolve.cpp index 621395a..7125202 100644 --- a/src/timestepping/bmx_evolve.cpp +++ b/src/timestepping/bmx_evolve.cpp @@ -93,9 +93,10 @@ bmx::Evolve (int nstep, Real & dt, Real & prev_dt, Real time, Real stop_time) pc->EvolveParticles(dt, particle_cost, knapsack_weight_type, nsubsteps); pc->split_particles(time); pc->ParticleExchange(dt, particle_cost, knapsack_weight_type, nsubsteps); - pc->EvaluateTipFusion(particle_cost,knapsack_weight_type); - pc->EvaluateInteriorFusion(particle_cost,knapsack_weight_type); - pc->CleanupFusion(particle_cost,knapsack_weight_type); + if (pc->EvaluateTipFusion(particle_cost,knapsack_weight_type)) { + pc->EvaluateInteriorFusion(particle_cost,knapsack_weight_type); + pc->CleanupFusion(particle_cost,knapsack_weight_type); + } } BL_PROFILE_VAR_STOP(particlesSolve);