From ad4f13a2ed5004533620e0028d63c6ca5eeeef4f Mon Sep 17 00:00:00 2001 From: S-Explorer <1246206018@qq.com> Date: Tue, 9 Apr 2024 10:58:39 +0800 Subject: [PATCH 1/5] feat:fix bug with initialLargrangianPoint and velocityInterpolation [FUNC] --- Source/DiffusedIB.H | 6 +-- Source/DiffusedIB.cpp | 93 +++++++++++++++++++++++-------------------- 2 files changed, 53 insertions(+), 46 deletions(-) diff --git a/Source/DiffusedIB.H b/Source/DiffusedIB.H index b1359c67..07d68259 100644 --- a/Source/DiffusedIB.H +++ b/Source/DiffusedIB.H @@ -49,7 +49,7 @@ struct kernel{ RealVect varphi; Real radious; Real rho; - Real ml; + int ml; Real dv; }; @@ -118,7 +118,7 @@ public: void InitialWithLargrangianPoints(const kernel& current_kernel); - void VelocityInterpolation(const amrex::MultiFab &Euler, DELTA_FUNCTION_TYPE type); + void VelocityInterpolation(amrex::MultiFab &Euler, DELTA_FUNCTION_TYPE type); void ForceSpreading(amrex::MultiFab &Euler, DELTA_FUNCTION_TYPE type); @@ -126,7 +126,7 @@ public: void ComputeLagrangianForce(Real dt, const kernel& kernel); - int verbose = 0; + int verbose = 1; }; class Particles{ diff --git a/Source/DiffusedIB.cpp b/Source/DiffusedIB.cpp index 4eb2fdec..480c39e8 100644 --- a/Source/DiffusedIB.cpp +++ b/Source/DiffusedIB.cpp @@ -97,20 +97,18 @@ void mParticle::InteractWithEuler(MultiFab &EulerVel, MultiFab &EulerForce, int if (verbose) amrex::Print() << "mParticle::InteractWithEuler " << std::endl; for(kernel& kernel : particle_kernels){ - InitialWithLargrangianPoints(kernel); // Initialize markers for a specific particle - + Redistribute(); //UpdateParticles(Euler, kernel, dt, alpha_k); - const int EulerForceIndex = euler_force_index; //for 1 -> Ns while(loop_time > 0){ - EulerForce.setVal(0.0, EulerForceIndex, 3, EulerForce.nGrow()); //clear Euler force + EulerForce.setVal(0.0, euler_force_index, 3, EulerForce.nGrow()); //clear Euler force VelocityInterpolation(EulerVel, type); ComputeLagrangianForce(dt, kernel); ForceSpreading(EulerForce, type); MultiFab::Saxpy(EulerVel, dt, EulerForce, euler_force_index, euler_velocity_index, 3, 0); //VelocityCorrection loop_time--; - }; + } } } @@ -143,7 +141,6 @@ void mParticle::InitParticles(const Vector& x, return; } //all the particles have same radious - Real phiK = 0; Real h = m_gdb->Geom(euler_finest_level).CellSizeArray()[0]; int Ml = static_cast( Math::pi() / 3 * (12 * Math::powi<2>(radious / h))); Real dv = Math::pi() * h / 3 / Ml * (12 * radious * radious + h * h); @@ -186,9 +183,9 @@ void mParticle::InitParticles(const Vector& x, ParticleType markerP; markerP.id() = ParticleType::NextID(); markerP.cpu() = ParallelDescriptor::MyProc(); - markerP.pos(0) = 0; - markerP.pos(1) = 0; - markerP.pos(2) = 0; + markerP.pos(0) = 0.01; + markerP.pos(1) = 0.01; + markerP.pos(2) = 0.01; std::array Marker_attr; Marker_attr[U_Marker] = 0.0; @@ -202,36 +199,51 @@ void mParticle::InitParticles(const Vector& x, particleTileTmp.push_back_real(Marker_attr); } } - //Redistribute(); // No need to redistribute here! + // Redistribute(); // No need to redistribute here! if (verbose) WriteAsciiFile(amrex::Concatenate("particle", 0)); } +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void InitalLargrangianPointsLoc (amrex::ParticleReal& x, + amrex::ParticleReal& y, + amrex::ParticleReal& z, + const amrex::Long id, + const amrex::RealVect location, + const amrex::Real radious, + const int ml) noexcept{ + Real Hk = -1.0 + 2.0 * (id - 1) / ( ml - 1.0); + Real thetaK = std::acos(Hk); + Real phiK = 0; + if(id == 1 || id == ml){ + phiK = 0; + }else { + phiK = std::fmod( phiK + 3.809 / std::sqrt(ml) / std::sqrt( 1 - Math::powi<2>(Hk)) , 2 * Math::pi()); + } + // update LargrangianPoint position with particle position + x = location[0] + radious * std::sin(thetaK) * std::cos(phiK); + y = location[1] + radious * std::sin(thetaK) * std::sin(phiK); + z = location[2] + radious * std::cos(thetaK); +} + void mParticle::InitialWithLargrangianPoints(const kernel& current_kernel){ if (verbose) amrex::Print() << "mParticle::InitialWithLargrangianPoints " << std::endl; - // Update the markers' locations - if ( ParallelDescriptor::MyProc() == ParallelDescriptor::IOProcessorNumber() ) { - mParIter pti(*this, euler_finest_level); - auto *particles = pti.GetArrayOfStructs().data(); + for(mParIter pti(*this, euler_finest_level); pti.isValid(); ++pti){ - Real phiK = 0; - for(int index = 0; index < current_kernel.ml; index++){ - Real Hk = -1.0 + 2.0 * (index) / ( current_kernel.ml - 1.0); - Real thetaK = std::acos(Hk); - if(index == 0 || index == ( current_kernel.ml - 1)){ - phiK = 0; - }else { - phiK = std::fmod( phiK + 3.809 / std::sqrt(current_kernel.ml) / std::sqrt( 1 - Math::powi<2>(Hk)) , 2 * Math::pi()); + auto *particles = pti.GetArrayOfStructs().data(); + auto np = pti.numParticles(); + amrex::ParallelFor( np, [=] + AMREX_GPU_DEVICE (int i) noexcept { + InitalLargrangianPointsLoc(particles[i].pos(0),particles[i].pos(1),particles[i].pos(2), + particles[i].id(), + current_kernel.location, + current_kernel.radious, + current_kernel.ml); } - // update LargrangianPoint position with particle position - particles[index].pos(0) = current_kernel.location[0] + current_kernel.radious * std::sin(thetaK) * std::cos(phiK); - particles[index].pos(1) = current_kernel.location[1] + current_kernel.radious * std::sin(thetaK) * std::sin(phiK); - particles[index].pos(2) = current_kernel.location[2] + current_kernel.radious * std::cos(thetaK); - } + ); } // Redistribute the markers after updating their locations - Redistribute(); if (verbose) WriteAsciiFile(amrex::Concatenate("particle", 1)); } @@ -279,26 +291,28 @@ void VelocityInterpolation_cir(int p_iter, P const& p, Real& Up, Real& Vp, Real& deltaFunction( p.pos(1), yj, dx[1], tV, type); deltaFunction( p.pos(2), kz, dx[2], tW, type); const Real delta_value = tU * tV * tW; - Gpu::Atomic::AddNoRet( &Up, delta_value * E(i + ii, j + jj, k + kk, EulerVIndex ) * d ); - Gpu::Atomic::AddNoRet( &Vp, delta_value * E(i + ii, j + jj, k + kk, EulerVIndex+1) * d ); - Gpu::Atomic::AddNoRet( &Wp, delta_value * E(i + ii, j + jj, k + kk, EulerVIndex+2) * d ); + Up += delta_value * E(i + ii, j + jj, k + kk, EulerVIndex ) * d; + Vp += delta_value * E(i + ii, j + jj, k + kk, EulerVIndex + 1) * d; + Wp += delta_value * E(i + ii, j + jj, k + kk, EulerVIndex + 2) * d; } } } } -void mParticle::VelocityInterpolation(const MultiFab &EulerVel, +void mParticle::VelocityInterpolation(MultiFab &EulerVel, DELTA_FUNCTION_TYPE type)// { if (verbose) amrex::Print() << "mParticle::VelocityInterpolation " << std::endl; //amrex::Print() << "euler_finest_level " << euler_finest_level << std::endl; - const auto& gm = m_gdb->Geom(euler_finest_level); auto plo = gm.ProbLoArray(); - auto phi = gm.ProbHiArray(); auto dx = gm.CellSizeArray(); + // attention + // velocity ghost cells will be up-to-date + EulerVel.FillBoundary(euler_velocity_index, 3, gm.periodicity()); + const int EulerVelocityIndex = euler_velocity_index; // std::cout << "plo: "; @@ -309,7 +323,6 @@ void mParticle::VelocityInterpolation(const MultiFab &EulerVel, // for (const auto& val : dx) std::cout << val << " "; // std::cout << "\nEulerVelocityIndex: " << EulerVelocityIndex << std::endl; - auto ba = EulerVel.boxArray(); // amrex::Print() << "ba " << ba << std::endl; for(mParIter pti(*this, euler_finest_level); pti.isValid(); ++pti){ @@ -390,7 +403,6 @@ void mParticle::ForceSpreading(MultiFab & EulerForce, if (verbose) amrex::Print() << "mParticle::ForceSpreading " << std::endl; - int index = 0; const auto& gm = m_gdb->Geom(euler_finest_level); auto plo = gm.ProbLoArray(); auto dxi = gm.CellSizeArray(); @@ -453,9 +465,6 @@ void mParticle::UpdateParticles(const amrex::MultiFab& Euler, kernel& kernel, Re if (verbose) amrex::Print() << "mParticle::UpdateParticles " << std::endl; - const auto& gm = m_gdb->Geom(euler_finest_level); - auto plo = gm.ProbLoArray(); - auto dxi = gm.InvCellSizeArray(); //update the kernel's infomation and cal body force for(mParIter pti(*this, euler_finest_level); pti.isValid(); ++pti){ auto &particles = pti.GetArrayOfStructs(); @@ -498,7 +507,6 @@ void mParticle::UpdateParticles(const amrex::MultiFab& Euler, kernel& kernel, Re *location_ptr = *location_ptr + deltaX; *varphi_ptr = *varphi_ptr + alpha_k * dt * (*omega_ptr + oldOmega); //sum - auto Uarray = Euler[pti].array(); amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE (int i) noexcept{ //calculate the force @@ -528,17 +536,16 @@ void mParticle::ComputeLagrangianForce(Real dt, const kernel& kernel) Real Wb = kernel.velocity[2]; for(mParIter pti(*this, euler_finest_level); pti.isValid(); ++pti){ - auto& particles = pti.GetArrayOfStructs(); - auto *p_ptr = particles.data(); const Long np = pti.numParticles(); - auto& attri = pti.GetAttribs(); + auto* Up = attri[P_ATTR::U_Marker].data(); auto* Vp = attri[P_ATTR::V_Marker].data(); auto* Wp = attri[P_ATTR::W_Marker].data(); auto *FxP = attri[P_ATTR::Fx_Marker].data(); auto *FyP = attri[P_ATTR::Fy_Marker].data(); auto *FzP = attri[P_ATTR::Fz_Marker].data(); + amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE (int i) noexcept{ FxP[i] = (Ub - Up[i])/dt; // From 936dae8054a8198f86be4c7a320afc03a05b91ab Mon Sep 17 00:00:00 2001 From: S-Explorer <1246206018@qq.com> Date: Tue, 9 Apr 2024 14:14:02 +0800 Subject: [PATCH 2/5] fix multi particle calculate & deltaFunction --- Source/DiffusedIB.H | 2 +- Source/DiffusedIB.cpp | 11 ++++++----- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/Source/DiffusedIB.H b/Source/DiffusedIB.H index 07d68259..a625e59c 100644 --- a/Source/DiffusedIB.H +++ b/Source/DiffusedIB.H @@ -126,7 +126,7 @@ public: void ComputeLagrangianForce(Real dt, const kernel& kernel); - int verbose = 1; + int verbose = 0; }; class Particles{ diff --git a/Source/DiffusedIB.cpp b/Source/DiffusedIB.cpp index 480c39e8..9529a0a9 100644 --- a/Source/DiffusedIB.cpp +++ b/Source/DiffusedIB.cpp @@ -66,7 +66,7 @@ void deltaFunction(Real xf, Real xp, Real h, Real& value, DELTA_FUNCTION_TYPE ty switch (type) { case DELTA_FUNCTION_TYPE::FOUR_POINT_IB: - if(rr >= 0 && rr < 0.5 ){ + if(rr >= 0 && rr < 1 ){ value = 1.0 / 8.0 * ( 3.0 - 2.0 * rr + std::sqrt( 1.0 + 4.0 * rr - 4 * Math::powi<2>(rr))) / h; }else if (rr >= 1 && rr < 2) { value = 1.0 / 8.0 * ( 5.0 - 2.0 * rr - std::sqrt( -7.0 + 12.0 * rr - 4 * Math::powi<2>(rr))) / h; @@ -75,9 +75,9 @@ void deltaFunction(Real xf, Real xp, Real h, Real& value, DELTA_FUNCTION_TYPE ty } break; case DELTA_FUNCTION_TYPE::THREE_POINT_IB: - if(rr >= 0 && rr < 1){ + if(rr >= 0 && rr < 0.5){ value = 1.0 / 6.0 * ( 5.0 - 3.0 * rr + std::sqrt( - 3.0 * ( 1 - Math::powi<2>(rr)) + 1.0 )) / h; - }else if (rr >= 1 && rr < 2) { + }else if (rr >= 0.5 && rr < 1.5) { value = 1.0 / 3.0 * ( 1.0 + std::sqrt( 1.0 - 3 * Math::powi<2>(rr))) / h; }else { value = 0; @@ -101,13 +101,14 @@ void mParticle::InteractWithEuler(MultiFab &EulerVel, MultiFab &EulerForce, int Redistribute(); //UpdateParticles(Euler, kernel, dt, alpha_k); //for 1 -> Ns - while(loop_time > 0){ + int loop = loop_time; + while(loop > 0){ EulerForce.setVal(0.0, euler_force_index, 3, EulerForce.nGrow()); //clear Euler force VelocityInterpolation(EulerVel, type); ComputeLagrangianForce(dt, kernel); ForceSpreading(EulerForce, type); MultiFab::Saxpy(EulerVel, dt, EulerForce, euler_force_index, euler_velocity_index, 3, 0); //VelocityCorrection - loop_time--; + loop--; } } } From a4c9615f73bd098aa408e8cfcb43f1bb09e85a74 Mon Sep 17 00:00:00 2001 From: S-Explorer <1246206018@qq.com> Date: Thu, 11 Apr 2024 09:42:32 +0800 Subject: [PATCH 3/5] fix : largrangian marker's position and velocity corrected with particle velocity & deltafunction --- Source/DiffusedIB.H | 4 ++- Source/DiffusedIB.cpp | 62 ++++++++++++++++++++++++------------------- 2 files changed, 37 insertions(+), 29 deletions(-) diff --git a/Source/DiffusedIB.H b/Source/DiffusedIB.H index a625e59c..57dab832 100644 --- a/Source/DiffusedIB.H +++ b/Source/DiffusedIB.H @@ -51,6 +51,8 @@ struct kernel{ Real rho; int ml; Real dv; + std::vector thetaK; + std::vector phiK; }; class mParIter : public ParIter<0, 0, numAttri>{ @@ -120,7 +122,7 @@ public: void VelocityInterpolation(amrex::MultiFab &Euler, DELTA_FUNCTION_TYPE type); - void ForceSpreading(amrex::MultiFab &Euler, DELTA_FUNCTION_TYPE type); + void ForceSpreading(amrex::MultiFab &Euler, Real dv, DELTA_FUNCTION_TYPE type); void UpdateParticles(const amrex::MultiFab& Euler, kernel& kernel, Real dt, Real alpha_k ); diff --git a/Source/DiffusedIB.cpp b/Source/DiffusedIB.cpp index 9529a0a9..fa9a9292 100644 --- a/Source/DiffusedIB.cpp +++ b/Source/DiffusedIB.cpp @@ -75,9 +75,9 @@ void deltaFunction(Real xf, Real xp, Real h, Real& value, DELTA_FUNCTION_TYPE ty } break; case DELTA_FUNCTION_TYPE::THREE_POINT_IB: - if(rr >= 0 && rr < 0.5){ - value = 1.0 / 6.0 * ( 5.0 - 3.0 * rr + std::sqrt( - 3.0 * ( 1 - Math::powi<2>(rr)) + 1.0 )) / h; - }else if (rr >= 0.5 && rr < 1.5) { + if(rr >= 0.5 && rr < 1.5){ + value = 1.0 / 6.0 * ( 5.0 - 3.0 * rr - std::sqrt( - 3.0 * Math::powi<2>( 1 - rr) + 1.0 )) / h; + }else if (rr >= 0 && rr < 0.5) { value = 1.0 / 3.0 * ( 1.0 + std::sqrt( 1.0 - 3 * Math::powi<2>(rr))) / h; }else { value = 0; @@ -106,8 +106,8 @@ void mParticle::InteractWithEuler(MultiFab &EulerVel, MultiFab &EulerForce, int EulerForce.setVal(0.0, euler_force_index, 3, EulerForce.nGrow()); //clear Euler force VelocityInterpolation(EulerVel, type); ComputeLagrangianForce(dt, kernel); - ForceSpreading(EulerForce, type); - MultiFab::Saxpy(EulerVel, dt, EulerForce, euler_force_index, euler_velocity_index, 3, 0); //VelocityCorrection + ForceSpreading(EulerForce, kernel.dv, type); + MultiFab::Saxpy(EulerVel, dt, EulerForce, euler_force_index, euler_velocity_index, 3, EulerForce.nGrow()); //VelocityCorrection loop--; } } @@ -166,7 +166,20 @@ void mParticle::InitParticles(const Vector& x, mKernel.ml = Ml; mKernel.dv = dv; mKernel.rho = rho_s; - particle_kernels.push_back(mKernel); + + Real phiK = 0; + for(int id = 1; id <= Ml; id++){ + Real Hk = -1.0 + 2.0 * (id - 1) / ( Ml - 1.0); + Real thetaK = std::acos(Hk); + if(id == 1 || id == Ml){ + phiK = 0; + }else { + phiK = std::fmod( phiK + 3.809 / std::sqrt(Ml) / std::sqrt( 1 - Math::powi<2>(Hk)) , 2 * Math::pi()); + } + mKernel.thetaK.emplace_back(thetaK); + mKernel.phiK.emplace_back(phiK); + } + particle_kernels.emplace_back(mKernel); if (verbose) amrex::Print() << "Kernel " << index << ": Location (" << x[index] << ", " << y[index] << ", " << z[index] << "), Radius: " << radious << ", Ml: " << Ml << ", dv: " << dv << ", Rho: " << rho_s << "\n"; @@ -195,7 +208,7 @@ void mParticle::InitParticles(const Vector& x, Marker_attr[Fx_Marker] = 0.0; Marker_attr[Fy_Marker] = 0.0; Marker_attr[Fz_Marker] = 0.0; - // attr[V] = 10.0; + particleTileTmp.push_back(markerP); particleTileTmp.push_back_real(Marker_attr); } @@ -208,18 +221,11 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void InitalLargrangianPointsLoc (amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::ParticleReal& z, - const amrex::Long id, + amrex::Real thetaK, + amrex::Real phiK, const amrex::RealVect location, const amrex::Real radious, const int ml) noexcept{ - Real Hk = -1.0 + 2.0 * (id - 1) / ( ml - 1.0); - Real thetaK = std::acos(Hk); - Real phiK = 0; - if(id == 1 || id == ml){ - phiK = 0; - }else { - phiK = std::fmod( phiK + 3.809 / std::sqrt(ml) / std::sqrt( 1 - Math::powi<2>(Hk)) , 2 * Math::pi()); - } // update LargrangianPoint position with particle position x = location[0] + radious * std::sin(thetaK) * std::cos(phiK); y = location[1] + radious * std::sin(thetaK) * std::sin(phiK); @@ -237,7 +243,8 @@ void mParticle::InitialWithLargrangianPoints(const kernel& current_kernel){ amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (int i) noexcept { InitalLargrangianPointsLoc(particles[i].pos(0),particles[i].pos(1),particles[i].pos(2), - particles[i].id(), + current_kernel.thetaK.at(i), + current_kernel.phiK.at(i), current_kernel.location, current_kernel.radious, current_kernel.ml); @@ -365,6 +372,7 @@ void ForceSpreading_cic (P const& p, ParticleReal fzP, Array4 const& E, int EulerForceIndex, + Real dv, GpuArray const& plo, GpuArray const& dx, DELTA_FUNCTION_TYPE type) @@ -391,15 +399,16 @@ void ForceSpreading_cic (P const& p, deltaFunction( p.pos(1), yj, dx[1], tV, type); deltaFunction( p.pos(2), kz, dx[2], tW, type); Real delta_value = tU * tV * tW; - Gpu::Atomic::AddNoRet(&E(i + ii, j + jj, k + kk, EulerForceIndex ), delta_value * fxP * d); - Gpu::Atomic::AddNoRet(&E(i + ii, j + jj, k + kk, EulerForceIndex+1), delta_value * fyP * d); - Gpu::Atomic::AddNoRet(&E(i + ii, j + jj, k + kk, EulerForceIndex+2), delta_value * fzP * d); + Gpu::Atomic::AddNoRet(&E(i + ii, j + jj, k + kk, EulerForceIndex ), delta_value * fxP * dv); + Gpu::Atomic::AddNoRet(&E(i + ii, j + jj, k + kk, EulerForceIndex+1), delta_value * fyP * dv); + Gpu::Atomic::AddNoRet(&E(i + ii, j + jj, k + kk, EulerForceIndex+2), delta_value * fzP * dv); } } } } void mParticle::ForceSpreading(MultiFab & EulerForce, + Real dv, DELTA_FUNCTION_TYPE type){ if (verbose) amrex::Print() << "mParticle::ForceSpreading " << std::endl; @@ -410,20 +419,17 @@ void mParticle::ForceSpreading(MultiFab & EulerForce, const int EulerForceIndex = euler_force_index; for(mParIter pti(*this, euler_finest_level); pti.isValid(); ++pti){ const Long np = pti.numParticles(); - const auto& fxP = pti.GetStructOfArrays().GetRealData(P_ATTR::U_Marker);//Fx_Marker - const auto& fyP = pti.GetStructOfArrays().GetRealData(P_ATTR::V_Marker);//Fy_Marker - const auto& fzP = pti.GetStructOfArrays().GetRealData(P_ATTR::W_Marker);//Fz_Marker const auto& particles = pti.GetArrayOfStructs(); - auto Uarray = EulerForce[pti].array(); + auto& attri = pti.GetAttribs(); - const auto& fxP_ptr = fxP.data(); - const auto& fyP_ptr = fyP.data(); - const auto& fzP_ptr = fzP.data(); + const auto& fxP_ptr = attri[P_ATTR::Fx_Marker].data(); + const auto& fyP_ptr = attri[P_ATTR::Fy_Marker].data(); + const auto& fzP_ptr = attri[P_ATTR::Fz_Marker].data(); const auto& p_ptr = particles().data(); amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE (int i) noexcept{ - ForceSpreading_cic(p_ptr[i], fxP_ptr[i], fyP_ptr[i], fzP_ptr[i], Uarray, EulerForceIndex, plo, dxi, type); + ForceSpreading_cic(p_ptr[i], fxP_ptr[i], fyP_ptr[i], fzP_ptr[i], Uarray, EulerForceIndex, dv, plo, dxi, type); }); } EulerForce.SumBoundary(EulerForceIndex, 3, gm.periodicity()); From df85e486da9f6f2d85925cfa8ff0fab4ee6fb557 Mon Sep 17 00:00:00 2001 From: S-Explorer <1246206018@qq.com> Date: Fri, 12 Apr 2024 09:38:42 +0800 Subject: [PATCH 4/5] add a function for export IB force and fix bug (calc the marker's position by id) --- Source/DiffusedIB.H | 11 ++++++++--- Source/DiffusedIB.cpp | 44 +++++++++++++++++++++++++++++++++++++------ 2 files changed, 46 insertions(+), 9 deletions(-) diff --git a/Source/DiffusedIB.H b/Source/DiffusedIB.H index 57dab832..2d49a604 100644 --- a/Source/DiffusedIB.H +++ b/Source/DiffusedIB.H @@ -51,6 +51,7 @@ struct kernel{ Real rho; int ml; Real dv; + RealVect ib_forece; std::vector thetaK; std::vector phiK; }; @@ -103,11 +104,11 @@ public: int velocity_index, int finest_level); - void InteractWithEuler(MultiFab &EulerVel, MultiFab &EulerForce, int loop_time = 20, Real dt = 0.1, Real alpha_k = 0.5, DELTA_FUNCTION_TYPE type = FOUR_POINT_IB); + void InteractWithEuler(amrex::Real time, MultiFab &EulerVel, MultiFab &EulerForce, int loop_time = 20, Real dt = 0.1, Real alpha_k = 0.5, DELTA_FUNCTION_TYPE type = FOUR_POINT_IB); void WriteParticleFile(int index); -// private: +private: int euler_finest_level; int euler_force_index; @@ -122,12 +123,16 @@ public: void VelocityInterpolation(amrex::MultiFab &Euler, DELTA_FUNCTION_TYPE type); - void ForceSpreading(amrex::MultiFab &Euler, Real dv, DELTA_FUNCTION_TYPE type); + void ForceSpreading(amrex::MultiFab &Euler, RealVect& ib_forece, Real dv, DELTA_FUNCTION_TYPE type); void UpdateParticles(const amrex::MultiFab& Euler, kernel& kernel, Real dt, Real alpha_k ); void ComputeLagrangianForce(Real dt, const kernel& kernel); + void WriteIBForce(amrex::Real time); + + uint32_t ib_forece_file_index = 0; + int verbose = 0; }; diff --git a/Source/DiffusedIB.cpp b/Source/DiffusedIB.cpp index fa9a9292..cc5ef37e 100644 --- a/Source/DiffusedIB.cpp +++ b/Source/DiffusedIB.cpp @@ -92,7 +92,7 @@ void deltaFunction(Real xf, Real xp, Real h, Real& value, DELTA_FUNCTION_TYPE ty /* mParticle member function */ /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ //loop all particels -void mParticle::InteractWithEuler(MultiFab &EulerVel, MultiFab &EulerForce, int loop_time, Real dt, Real alpha_k, DELTA_FUNCTION_TYPE type){ +void mParticle::InteractWithEuler(amrex::Real time, MultiFab &EulerVel, MultiFab &EulerForce, int loop_time, Real dt, Real alpha_k, DELTA_FUNCTION_TYPE type){ if (verbose) amrex::Print() << "mParticle::InteractWithEuler " << std::endl; @@ -106,11 +106,13 @@ void mParticle::InteractWithEuler(MultiFab &EulerVel, MultiFab &EulerForce, int EulerForce.setVal(0.0, euler_force_index, 3, EulerForce.nGrow()); //clear Euler force VelocityInterpolation(EulerVel, type); ComputeLagrangianForce(dt, kernel); - ForceSpreading(EulerForce, kernel.dv, type); + ForceSpreading(EulerForce, kernel.ib_forece, kernel.dv, type); MultiFab::Saxpy(EulerVel, dt, EulerForce, euler_force_index, euler_velocity_index, 3, EulerForce.nGrow()); //VelocityCorrection loop--; } } + + WriteIBForce(time); } void mParticle::InitParticles(const Vector& x, @@ -242,9 +244,10 @@ void mParticle::InitialWithLargrangianPoints(const kernel& current_kernel){ auto np = pti.numParticles(); amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (int i) noexcept { + int id = static_cast(particles[i].id() - 1); InitalLargrangianPointsLoc(particles[i].pos(0),particles[i].pos(1),particles[i].pos(2), - current_kernel.thetaK.at(i), - current_kernel.phiK.at(i), + current_kernel.thetaK.at(id), + current_kernel.phiK.at(id), current_kernel.location, current_kernel.radious, current_kernel.ml); @@ -408,11 +411,16 @@ void ForceSpreading_cic (P const& p, } void mParticle::ForceSpreading(MultiFab & EulerForce, + RealVect& ib_forece, Real dv, DELTA_FUNCTION_TYPE type){ if (verbose) amrex::Print() << "mParticle::ForceSpreading " << std::endl; + Real fx = 0; + Real fy = 0; + Real fz = 0; + const auto& gm = m_gdb->Geom(euler_finest_level); auto plo = gm.ProbLoArray(); auto dxi = gm.CellSizeArray(); @@ -427,13 +435,18 @@ void mParticle::ForceSpreading(MultiFab & EulerForce, const auto& fyP_ptr = attri[P_ATTR::Fy_Marker].data(); const auto& fzP_ptr = attri[P_ATTR::Fz_Marker].data(); const auto& p_ptr = particles().data(); - amrex::ParallelFor(np, [=] + amrex::ParallelFor(np, [=, &fx, &fy, &fz] AMREX_GPU_DEVICE (int i) noexcept{ + fx += fxP_ptr[i] * dv; + fy += fyP_ptr[i] * dv; + fz += fzP_ptr[i] * dv; ForceSpreading_cic(p_ptr[i], fxP_ptr[i], fyP_ptr[i], fzP_ptr[i], Uarray, EulerForceIndex, dv, plo, dxi, type); }); } EulerForce.SumBoundary(EulerForceIndex, 3, gm.periodicity()); - + + ib_forece = RealVect(fx, fy, fz); + if (verbose) { // Check the Multifab // Open a file stream for output @@ -568,6 +581,25 @@ void mParticle::WriteParticleFile(int index) WriteAsciiFile(amrex::Concatenate("particle", index)); } +void mParticle::WriteIBForce(amrex::Real time) +{ + // + std::ofstream out_ib_force("IB_Force_" + std::to_string(++ib_forece_file_index) + ".csv"); + int index = 1; + out_ib_force << "# This document collects all IB force information\n" + << "# you can deal this with csv format\n" + << "# time " << std::to_string(time) << "\n\n"; + out_ib_force << "particle_num,X,Y,Z,Vx,Vy,Vz,Fx,Fy,Fz\n"; + for(auto const& kernel : particle_kernels){ + out_ib_force << "particle" << index++ << "," + << std::to_string(kernel.location[0]) << "," << std::to_string(kernel.location[1]) << "," << std::to_string(kernel.location[2]) << "," + << std::to_string(kernel.velocity[0]) << "," << std::to_string(kernel.velocity[1]) << "," << std::to_string(kernel.velocity[2]) << "," + << std::to_string(kernel.ib_forece[0]) << "," << std::to_string(kernel.ib_forece[1]) << "," << std::to_string(kernel.ib_forece[2]) + << "\n"; + } + out_ib_force.close(); +} + /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ /* Particles member function */ /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ From 76fe121ff504431066a23f5c84df0849ae512026 Mon Sep 17 00:00:00 2001 From: S-Explorer <1246206018@qq.com> Date: Fri, 12 Apr 2024 15:26:51 +0800 Subject: [PATCH 5/5] add Parallel Reduce to export IB force --- Source/DiffusedIB.H | 5 ++- Source/DiffusedIB.cpp | 88 +++++++++++++++++++++++++---------------- Source/NavierStokes.cpp | 2 +- 3 files changed, 59 insertions(+), 36 deletions(-) diff --git a/Source/DiffusedIB.H b/Source/DiffusedIB.H index 2d49a604..1018863e 100644 --- a/Source/DiffusedIB.H +++ b/Source/DiffusedIB.H @@ -43,6 +43,7 @@ enum DELTA_FUNCTION_TYPE{ }; struct kernel{ + int id; RealVect velocity; RealVect location; RealVect omega; @@ -104,7 +105,7 @@ public: int velocity_index, int finest_level); - void InteractWithEuler(amrex::Real time, MultiFab &EulerVel, MultiFab &EulerForce, int loop_time = 20, Real dt = 0.1, Real alpha_k = 0.5, DELTA_FUNCTION_TYPE type = FOUR_POINT_IB); + void InteractWithEuler(int iStep, amrex::Real time, MultiFab &EulerVel, MultiFab &EulerForce, int loop_time = 20, Real dt = 0.1, Real alpha_k = 0.5, DELTA_FUNCTION_TYPE type = FOUR_POINT_IB); void WriteParticleFile(int index); @@ -129,7 +130,7 @@ private: void ComputeLagrangianForce(Real dt, const kernel& kernel); - void WriteIBForce(amrex::Real time); + static void WriteIBForce(int step, amrex::Real time, kernel& current_kernel); uint32_t ib_forece_file_index = 0; diff --git a/Source/DiffusedIB.cpp b/Source/DiffusedIB.cpp index cc5ef37e..c4bb6e20 100644 --- a/Source/DiffusedIB.cpp +++ b/Source/DiffusedIB.cpp @@ -9,6 +9,9 @@ #include #include +#include +namespace fs = std::filesystem; + using namespace amrex; void nodal_phi_to_pvf(MultiFab& pvf, const MultiFab& phi_nodal) @@ -92,7 +95,7 @@ void deltaFunction(Real xf, Real xp, Real h, Real& value, DELTA_FUNCTION_TYPE ty /* mParticle member function */ /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ //loop all particels -void mParticle::InteractWithEuler(amrex::Real time, MultiFab &EulerVel, MultiFab &EulerForce, int loop_time, Real dt, Real alpha_k, DELTA_FUNCTION_TYPE type){ +void mParticle::InteractWithEuler(int iStep, amrex::Real time, MultiFab &EulerVel, MultiFab &EulerForce, int loop_time, Real dt, Real alpha_k, DELTA_FUNCTION_TYPE type){ if (verbose) amrex::Print() << "mParticle::InteractWithEuler " << std::endl; @@ -110,9 +113,8 @@ void mParticle::InteractWithEuler(amrex::Real time, MultiFab &EulerVel, MultiFab MultiFab::Saxpy(EulerVel, dt, EulerForce, euler_force_index, euler_velocity_index, 3, EulerForce.nGrow()); //VelocityCorrection loop--; } + WriteIBForce(iStep, time, kernel); } - - WriteIBForce(time); } void mParticle::InitParticles(const Vector& x, @@ -152,6 +154,7 @@ void mParticle::InitParticles(const Vector& x, for(int index = 0; index < x.size(); index++){ kernel mKernel; + mKernel.id = index + 1; mKernel.location[0] = x[index]; mKernel.location[1] = y[index]; mKernel.location[2] = z[index]; @@ -373,6 +376,9 @@ void ForceSpreading_cic (P const& p, ParticleReal fxP, ParticleReal fyP, ParticleReal fzP, + Real & ib_fx, + Real & ib_fy, + Real & ib_fz, Array4 const& E, int EulerForceIndex, Real dv, @@ -389,6 +395,12 @@ void ForceSpreading_cic (P const& p, int i = static_cast(Math::floor(lx)); int j = static_cast(Math::floor(ly)); int k = static_cast(Math::floor(lz)); + Real Fx = fxP * dv; + Real Fy = fyP * dv; + Real Fz = fzP * dv; + ib_fx += Fx; + ib_fy += Fy; + ib_fz += Fz; // calc_delta(i, j, k, dxi, rho); //lagrangian to Euler for(int ii = -2; ii < +3; ii++){ @@ -402,9 +414,9 @@ void ForceSpreading_cic (P const& p, deltaFunction( p.pos(1), yj, dx[1], tV, type); deltaFunction( p.pos(2), kz, dx[2], tW, type); Real delta_value = tU * tV * tW; - Gpu::Atomic::AddNoRet(&E(i + ii, j + jj, k + kk, EulerForceIndex ), delta_value * fxP * dv); - Gpu::Atomic::AddNoRet(&E(i + ii, j + jj, k + kk, EulerForceIndex+1), delta_value * fyP * dv); - Gpu::Atomic::AddNoRet(&E(i + ii, j + jj, k + kk, EulerForceIndex+2), delta_value * fzP * dv); + Gpu::Atomic::AddNoRet(&E(i + ii, j + jj, k + kk, EulerForceIndex ), delta_value * Fx); + Gpu::Atomic::AddNoRet(&E(i + ii, j + jj, k + kk, EulerForceIndex+1), delta_value * Fy); + Gpu::Atomic::AddNoRet(&E(i + ii, j + jj, k + kk, EulerForceIndex+2), delta_value * Fz); } } } @@ -416,11 +428,7 @@ void mParticle::ForceSpreading(MultiFab & EulerForce, DELTA_FUNCTION_TYPE type){ if (verbose) amrex::Print() << "mParticle::ForceSpreading " << std::endl; - - Real fx = 0; - Real fy = 0; - Real fz = 0; - + ib_forece.scale(0); const auto& gm = m_gdb->Geom(euler_finest_level); auto plo = gm.ProbLoArray(); auto dxi = gm.CellSizeArray(); @@ -435,18 +443,13 @@ void mParticle::ForceSpreading(MultiFab & EulerForce, const auto& fyP_ptr = attri[P_ATTR::Fy_Marker].data(); const auto& fzP_ptr = attri[P_ATTR::Fz_Marker].data(); const auto& p_ptr = particles().data(); - amrex::ParallelFor(np, [=, &fx, &fy, &fz] - AMREX_GPU_DEVICE (int i) noexcept{ - fx += fxP_ptr[i] * dv; - fy += fyP_ptr[i] * dv; - fz += fzP_ptr[i] * dv; - ForceSpreading_cic(p_ptr[i], fxP_ptr[i], fyP_ptr[i], fzP_ptr[i], Uarray, EulerForceIndex, dv, plo, dxi, type); + amrex::ParallelFor(np, [=, &ib_forece] + AMREX_GPU_DEVICE (int i) noexcept{ + ForceSpreading_cic(p_ptr[i], fxP_ptr[i], fyP_ptr[i], fzP_ptr[i], ib_forece[0], ib_forece[1], ib_forece[2], Uarray, EulerForceIndex, dv, plo, dxi, type); }); } EulerForce.SumBoundary(EulerForceIndex, 3, gm.periodicity()); - ib_forece = RealVect(fx, fy, fz); - if (verbose) { // Check the Multifab // Open a file stream for output @@ -581,21 +584,40 @@ void mParticle::WriteParticleFile(int index) WriteAsciiFile(amrex::Concatenate("particle", index)); } -void mParticle::WriteIBForce(amrex::Real time) +void mParticle::WriteIBForce(int step, amrex::Real time, kernel& current_kernel) { - // - std::ofstream out_ib_force("IB_Force_" + std::to_string(++ib_forece_file_index) + ".csv"); - int index = 1; - out_ib_force << "# This document collects all IB force information\n" - << "# you can deal this with csv format\n" - << "# time " << std::to_string(time) << "\n\n"; - out_ib_force << "particle_num,X,Y,Z,Vx,Vy,Vz,Fx,Fy,Fz\n"; - for(auto const& kernel : particle_kernels){ - out_ib_force << "particle" << index++ << "," - << std::to_string(kernel.location[0]) << "," << std::to_string(kernel.location[1]) << "," << std::to_string(kernel.location[2]) << "," - << std::to_string(kernel.velocity[0]) << "," << std::to_string(kernel.velocity[1]) << "," << std::to_string(kernel.velocity[2]) << "," - << std::to_string(kernel.ib_forece[0]) << "," << std::to_string(kernel.ib_forece[1]) << "," << std::to_string(kernel.ib_forece[2]) - << "\n"; + {//reduce all kernel data + amrex::ParallelAllReduce::Sum(current_kernel.velocity[0], ParallelDescriptor::Communicator()); + amrex::ParallelAllReduce::Sum(current_kernel.velocity[1], ParallelDescriptor::Communicator()); + amrex::ParallelAllReduce::Sum(current_kernel.velocity[2], ParallelDescriptor::Communicator()); + amrex::ParallelAllReduce::Sum(current_kernel.omega[0], ParallelDescriptor::Communicator()); + amrex::ParallelAllReduce::Sum(current_kernel.omega[1], ParallelDescriptor::Communicator()); + amrex::ParallelAllReduce::Sum(current_kernel.omega[2], ParallelDescriptor::Communicator()); + amrex::ParallelAllReduce::Sum(current_kernel.ib_forece[0], ParallelDescriptor::Communicator()); + amrex::ParallelAllReduce::Sum(current_kernel.ib_forece[1], ParallelDescriptor::Communicator()); + amrex::ParallelAllReduce::Sum(current_kernel.ib_forece[2], ParallelDescriptor::Communicator()); + } + if(amrex::ParallelDescriptor::MyProc() != ParallelDescriptor::IOProcessorNumber()) return; + + std::string file("IB_Force_Particle_" + std::to_string(current_kernel.id) + ".csv"); + std::ofstream out_ib_force; + + std::string head; + if(!fs::exists(file)){ + head = "iStep,time,X,Y,Z,Vx,Vy,Vz,OmegaX,OmegaY,OmegaZ,Fx,Fy,Fz\n"; + }else{ + head = ""; + } + + out_ib_force.open(file, std::ios::app); + if(!out_ib_force.is_open()){ + amrex::Print() << "[diffused IB output] write particle file error , step: " << step; + }else{ + out_ib_force << head << step << "," << time << "," + << current_kernel.location[0] << "," << current_kernel.location[1] << "," << current_kernel.location[2] << "," + << current_kernel.velocity[0] << "," << current_kernel.velocity[1] << "," << current_kernel.velocity[2] << "," + << current_kernel.omega[0] << "," << current_kernel.omega[1] << "," << current_kernel.omega[2] << "," + << current_kernel.ib_forece[0] << "," << current_kernel.ib_forece[1] << "," << current_kernel.ib_forece[2] << "\n"; } out_ib_force.close(); } diff --git a/Source/NavierStokes.cpp b/Source/NavierStokes.cpp index 79927c28..7ff1f002 100644 --- a/Source/NavierStokes.cpp +++ b/Source/NavierStokes.cpp @@ -2651,7 +2651,7 @@ NavierStokes::advance_semistaggered_fsi_diffusedib (Real time, // S_new.setVal(2.0, 1, 1, S_new.nGrow()); // v = 2 // S_new.setVal(3.0, 2, 1, S_new.nGrow()); // w = 3 MultiFab EulerForce(S_new.boxArray(), S_new.DistributionMap(), 3, S_new.nGrow()); - Particles::get_particles()->InteractWithEuler(S_new, EulerForce, 1, dt, 0.5); + Particles::get_particles()->InteractWithEuler(parent->levelSteps(0), time, S_new, EulerForce, 10, dt, 0.5); } //amrex::Abort("Stop here!"); #endif