From 879ac02f87940197f9dfe7948532f759a0fc9582 Mon Sep 17 00:00:00 2001 From: S-Explorer <1246206018@qq.com> Date: Mon, 15 Apr 2024 11:18:55 +0800 Subject: [PATCH 1/9] Add parameter reading functions and organize the code structure. Add configuration information for particles to the case setup file --- Source/DiffusedIB.H | 75 +-- Source/DiffusedIB.cpp | 437 +++++++++++------- Source/NavierStokes.cpp | 6 +- Source/NavierStokesBase.H | 1 - Source/NavierStokesBase.cpp | 13 +- .../FlowPastSphere/inputs.3d.flow_past_sphere | 50 ++ 6 files changed, 382 insertions(+), 200 deletions(-) diff --git a/Source/DiffusedIB.H b/Source/DiffusedIB.H index 1018863e..7f88fe2f 100644 --- a/Source/DiffusedIB.H +++ b/Source/DiffusedIB.H @@ -44,17 +44,26 @@ enum DELTA_FUNCTION_TYPE{ struct kernel{ int id; - RealVect velocity; - RealVect location; - RealVect omega; - RealVect varphi; - Real radious; - Real rho; - int ml; - Real dv; - RealVect ib_forece; + RealVect velocity{0.0,0.0,0.0}; + RealVect location{0.0,0.0,0.0}; + RealVect omega{0.0,0.0,0.0}; + RealVect varphi{0.0,0.0,0.0}; + Real radius{0.0}; + Real rho{0.0}; + int ml{0}; + Real dv{0.0}; + RealVect ib_forece{0.0,0.0,0.0}; + std::vector thetaK; std::vector phiK; + + RealVect sum_u{0.0,0.0,0.0}; + RealVect sum_t{0.0,0.0,0.0}; + + RealVect Fcp{0.0,0.0,0.0}; + RealVect Tcp{0.0,0.0,0.0}; + + int TLX{0},TLY{0},TLZ{0},RLX{0},RLY{0},RLZ{0}; }; class mParIter : public ParIter<0, 0, numAttri>{ @@ -98,14 +107,26 @@ public: void InitParticles(const Vector& x, const Vector& y, const Vector& z, - Real rho_s, - Real radious, - Real rho_f, + const Vector& rho_s, + const Vector& Vx, + const Vector& Vy, + const Vector& Vz, + const Vector& TLX, + const Vector& TLY, + const Vector& TLZ, + const Vector& RLX, + const Vector& RLY, + const Vector& RLZ, + Real radius, + Real rho_f, + Real gravity, int force_index, int velocity_index, - int finest_level); + int finest_level, + int verbose = 0, + int loop_time = 10); - 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 InteractWithEuler(int iStep, amrex::Real time, MultiFab &EulerVel, MultiFab &EulerForce, Real dt = 0.1, DELTA_FUNCTION_TYPE type = FOUR_POINT_IB); void WriteParticleFile(int index); @@ -118,22 +139,30 @@ private: Real euler_fluid_rho; + int loop_time; + Vector particle_kernels; void InitialWithLargrangianPoints(const kernel& current_kernel); + void UpdateMarkers(kernel& current_kernel, Real dt); + void VelocityInterpolation(amrex::MultiFab &Euler, DELTA_FUNCTION_TYPE type); + void ComputeLagrangianForce(Real dt, const kernel& kernel); + 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 VelocityCorrection(amrex::MultiFab &Euler, amrex::MultiFab &EulerForce, Real dt); - void ComputeLagrangianForce(Real dt, const kernel& kernel); + void UpdateParticles(const MultiFab& Euler_old, const MultiFab& Euler, const MultiFab& pvf, kernel& kernel, Real dt); static void WriteIBForce(int step, amrex::Real time, kernel& current_kernel); uint32_t ib_forece_file_index = 0; + RealVect m_gravity{0.0,0.0,0.0}; + int verbose = 0; }; @@ -144,27 +173,19 @@ public: const BoxArray & ba); static mParticle* get_particles(); - static void init_particle(); - + static void init_particle( int level, Real gravity); + static void Initialize(); static void define_para(const Vector& x, const Vector& y, const Vector& z, Real rho_s, - Real radious, + Real radius, Real rho_f, int force_index, int velocity_index, int finest_level); private: - inline static Vector _x{}, _y{}, _z{}; - inline static int euler_finest_level = 0; - inline static int euler_velocity_index = 0; - inline static int euler_force_index = 0; - inline static Real euler_fluid_rho = 0.0; - inline static Real euler_solid_rho = 0.0; - inline static Real particle_radious = 0.0; - inline static mParticle* particle = nullptr; }; diff --git a/Source/DiffusedIB.cpp b/Source/DiffusedIB.cpp index c4bb6e20..ed889c89 100644 --- a/Source/DiffusedIB.cpp +++ b/Source/DiffusedIB.cpp @@ -14,6 +14,29 @@ namespace fs = std::filesystem; using namespace amrex; +/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ +/* global variable */ +/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ + +namespace ParticleProperties{ + Vector _x{}, _y{}, _z{}, _rho{}; + Vector _Vx{}, _Vy{}, _Vz{}; + Real _radius; + Vector _TLX{},_TLY{},_TLZ{},_RLX{},_RLY{},_RLZ{}; + int euler_finest_level{0}; + int euler_velocity_index{0}; + int euler_force_index{0}; + Real euler_fluid_rho{0.0}; + int verbose{0}; + int loop_time{0}; + + GpuArray plo{0.0,0.0,0.0}, phi{0.0,0.0,0.0}, dx{0.0, 0.0, 0.0}; +} + +/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ +/* other function */ +/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ + void nodal_phi_to_pvf(MultiFab& pvf, const MultiFab& phi_nodal) { @@ -52,14 +75,64 @@ void nodal_phi_to_pvf(MultiFab& pvf, const MultiFab& phi_nodal) } -/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ -/* other function */ -/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ +void calculate_phi_nodal(MultiFab& phi_nodal, kernel& current_kernel) +{ + Real Xp2 = Math::powi<2>(current_kernel.location[0]); + Real Yp2 = Math::powi<2>(current_kernel.location[1]); + Real Zp2 = Math::powi<2>(current_kernel.location[2]); + Real Rp2 = Math::powi<2>(current_kernel.radius); + + for (MFIter mfi(phi_nodal,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + auto const& pnfab = phi_nodal.array(mfi); + amrex::ParallelFor(bx, [=, &pnfab] + AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + Real Xn = i * ParticleProperties::dx[0] + 0.5 * ParticleProperties::dx[0] + ParticleProperties::plo[0]; + Real Yn = j * ParticleProperties::dx[1] + 0.5 * ParticleProperties::dx[1] + ParticleProperties::plo[1]; + Real Zn = k * ParticleProperties::dx[2] + 0.5 * ParticleProperties::dx[2] + ParticleProperties::plo[2]; + + pnfab(i,j,k) = std::sqrt((Math::powi<2>(Xn) - Xp2)/Rp2 + + (Math::powi<2>(Yn) - Yp2)/Rp2 + + (Math::powi<2>(Zn) - Zp2)/Rp2 + ) - 1.0; + } + ); + } +} + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +void CalculateSum_cir (RealVect& sum, + const MultiFab& E_old, + const MultiFab& E, + const MultiFab& pvf, + int EulerVelIndex) +{ + const Real d = Math::powi<3>(ParticleProperties::dx[0]); + + for (MFIter mfi(E,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + auto const& uarray_old = E_old.array(mfi); + auto const& uarray = E.array(mfi); + auto const& pvffab = pvf.array(mfi); + amrex::ParallelFor(bx, [&, d] + AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + Gpu::Atomic::AddNoRet(&sum[0], (uarray(i, j, k, EulerVelIndex ) - uarray_old(i, j, k, EulerVelIndex )) * d * pvffab(i, j, k)); + Gpu::Atomic::AddNoRet(&sum[1], (uarray(i, j, k, EulerVelIndex + 1) - uarray_old(i, j, k, EulerVelIndex + 1)) * d * pvffab(i, j, k)); + Gpu::Atomic::AddNoRet(&sum[2], (uarray(i, j, k, EulerVelIndex + 2) - uarray_old(i, j, k, EulerVelIndex + 2)) * d * pvffab(i, j, k)); + } + ); + } + +} [[nodiscard]] AMREX_FORCE_INLINE -Real cal_momentum(Real rho, Real radious) +Real cal_momentum(Real rho, Real radius) { - return 8.0 * Math::pi() * rho * Math::powi<5>(radious) / 15.0; + return 8.0 * Math::pi() * rho * Math::powi<5>(radius) / 15.0; } AMREX_FORCE_INLINE @@ -95,22 +168,28 @@ void deltaFunction(Real xf, Real xp, Real h, Real& value, DELTA_FUNCTION_TYPE ty /* mParticle member function */ /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ //loop all particels -void mParticle::InteractWithEuler(int iStep, 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, + Real dt, + DELTA_FUNCTION_TYPE type){ - if (verbose) amrex::Print() << "mParticle::InteractWithEuler " << std::endl; + if (verbose) amrex::Print() << "[Particle] mParticle::InteractWithEuler " << std::endl; for(kernel& kernel : particle_kernels){ InitialWithLargrangianPoints(kernel); // Initialize markers for a specific particle Redistribute(); - //UpdateParticles(Euler, kernel, dt, alpha_k); + UpdateMarkers(kernel, dt); //for 1 -> Ns int loop = loop_time; while(loop > 0){ + if(verbose) amrex::Print() << "[Particle] Ns loop index : " << loop + 1 << "\n"; EulerForce.setVal(0.0, euler_force_index, 3, EulerForce.nGrow()); //clear Euler force VelocityInterpolation(EulerVel, type); ComputeLagrangianForce(dt, kernel); ForceSpreading(EulerForce, kernel.ib_forece, kernel.dv, type); - MultiFab::Saxpy(EulerVel, dt, EulerForce, euler_force_index, euler_velocity_index, 3, EulerForce.nGrow()); //VelocityCorrection + VelocityCorrection(EulerVel, EulerForce, dt); loop--; } WriteIBForce(iStep, time, kernel); @@ -120,19 +199,34 @@ void mParticle::InteractWithEuler(int iStep, amrex::Real time, MultiFab &EulerVe void mParticle::InitParticles(const Vector& x, const Vector& y, const Vector& z, - Real rho_s, - Real radious, - Real rho_f, + const Vector& rho_s, + const Vector& Vx, + const Vector& Vy, + const Vector& Vz, + const Vector& TLXt, + const Vector& TLYt, + const Vector& TLZt, + const Vector& RLXt, + const Vector& RLYt, + const Vector& RLZt, + Real radius, + Real rho_f, + Real gravity, int force_index, int velocity_index, - int finest_level){ - - if (verbose) amrex::Print() << "mParticle::InitParticles " << std::endl; + int finest_level, + int _verbose, + int _loop_time) +{ + verbose = _verbose; + loop_time = _loop_time; + if (verbose) amrex::Print() << "[Particle] mParticle::InitParticles " << std::endl; euler_finest_level = finest_level; euler_force_index = force_index; euler_fluid_rho = rho_f; euler_velocity_index = velocity_index; + m_gravity[2] = gravity; // Assuming the variables are defined similarly // amrex::Print() << "euler_finest_level: " << euler_finest_level << "\n" @@ -145,10 +239,10 @@ void mParticle::InitParticles(const Vector& x, Print() << "particle's position container are all different size"; return; } - //all the particles have same radious + //all the particles have same radius 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); + int Ml = static_cast( Math::pi() / 3 * (12 * Math::powi<2>(radius / h))); + Real dv = Math::pi() * h / 3 / Ml * (12 * radius * radius + h * h); if (verbose) amrex::Print() << "h: " << h << ", Ml: " << Ml << ", dv: " << dv << "\n"; @@ -158,19 +252,19 @@ void mParticle::InitParticles(const Vector& x, mKernel.location[0] = x[index]; mKernel.location[1] = y[index]; mKernel.location[2] = z[index]; - mKernel.velocity[0] = 0.0; - mKernel.velocity[1] = 0.0; - mKernel.velocity[2] = 0.0; - mKernel.omega[0] = 0.0; - mKernel.omega[1] = 0.0; - mKernel.omega[2] = 0.0; - mKernel.varphi[0] = 0.0; - mKernel.varphi[1] = 0.0; - mKernel.varphi[2] = 0.0; - mKernel.radious = radious; + mKernel.velocity[0] = Vx[index]; + mKernel.velocity[1] = Vy[index]; + mKernel.velocity[2] = Vz[index]; + mKernel.TLX = TLXt[index]; + mKernel.TLY = TLYt[index]; + mKernel.TLZ = TLZt[index]; + mKernel.RLX = RLXt[index]; + mKernel.RLY = RLYt[index]; + mKernel.RLZ = RLZt[index]; + mKernel.rho = rho_s[index]; + mKernel.radius = radius; mKernel.ml = Ml; mKernel.dv = dv; - mKernel.rho = rho_s; Real phiK = 0; for(int id = 1; id <= Ml; id++){ @@ -187,9 +281,9 @@ void mParticle::InitParticles(const Vector& x, 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"; + << "), Velocity : (" << mKernel.velocity[0] << ", " << mKernel.velocity[1] << ", "<< mKernel.velocity[2] + << "), Radius: " << radius << ", Ml: " << Ml << ", dv: " << dv << ", Rho: " << mKernel.rho << "\n"; } - //get particle tile std::pair key{0,0}; auto& particleTileTmp = GetParticles(0)[key]; @@ -224,22 +318,23 @@ void mParticle::InitParticles(const Vector& x, AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void InitalLargrangianPointsLoc (amrex::ParticleReal& x, - amrex::ParticleReal& y, - amrex::ParticleReal& z, - amrex::Real thetaK, - amrex::Real phiK, - const amrex::RealVect location, - const amrex::Real radious, - const int ml) noexcept{ + amrex::ParticleReal& y, + amrex::ParticleReal& z, + amrex::Real thetaK, + amrex::Real phiK, + const amrex::RealVect location, + const amrex::Real radius, + const int ml) noexcept +{ // 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); + x = location[0] + radius * std::sin(thetaK) * std::cos(phiK); + y = location[1] + radius * std::sin(thetaK) * std::sin(phiK); + z = location[2] + radius * std::cos(thetaK); } void mParticle::InitialWithLargrangianPoints(const kernel& current_kernel){ - if (verbose) amrex::Print() << "mParticle::InitialWithLargrangianPoints " << std::endl; + if (verbose) amrex::Print() << "mParticle::InitialWithLargrangianPoints\n"; // Update the markers' locations for(mParIter pti(*this, euler_finest_level); pti.isValid(); ++pti){ @@ -252,7 +347,7 @@ void mParticle::InitialWithLargrangianPoints(const kernel& current_kernel){ current_kernel.thetaK.at(id), current_kernel.phiK.at(id), current_kernel.location, - current_kernel.radious, + current_kernel.radius, current_kernel.ml); } ); @@ -264,11 +359,11 @@ void mParticle::InitialWithLargrangianPoints(const kernel& current_kernel){ template > AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void VelocityInterpolation_cir(int p_iter, P const& p, Real& Up, Real& Vp, Real& Wp, - Array4 const& E, int EulerVIndex, - const int *lo, const int *hi, - GpuArray const& plo, - GpuArray const& dx, - DELTA_FUNCTION_TYPE type) + Array4 const& E, int EulerVIndex, + const int *lo, const int *hi, + GpuArray const& plo, + GpuArray const& dx, + DELTA_FUNCTION_TYPE type) { //std::cout << "lo " << lo[0] << " " << lo[1] << " "<< lo[2] << "\n"; @@ -317,7 +412,7 @@ void mParticle::VelocityInterpolation(MultiFab &EulerVel, DELTA_FUNCTION_TYPE type)// { - if (verbose) amrex::Print() << "mParticle::VelocityInterpolation " << std::endl; + if (verbose) amrex::Print() << "\tmParticle::VelocityInterpolation " << std::endl; //amrex::Print() << "euler_finest_level " << euler_finest_level << std::endl; const auto& gm = m_gdb->Geom(euler_finest_level); @@ -373,18 +468,18 @@ void mParticle::VelocityInterpolation(MultiFab &EulerVel, template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE 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, - GpuArray const& plo, - GpuArray const& dx, - DELTA_FUNCTION_TYPE type) + ParticleReal fxP, + ParticleReal fyP, + ParticleReal fzP, + Real & ib_fx, + Real & ib_fy, + Real & ib_fz, + Array4 const& E, + int EulerForceIndex, + Real dv, + GpuArray const& plo, + GpuArray const& dx, + DELTA_FUNCTION_TYPE type) { const Real d = AMREX_D_TERM(dx[0], *dx[1], *dx[2]); //plo to ii jj kk @@ -398,9 +493,9 @@ void ForceSpreading_cic (P const& p, Real Fx = fxP * dv; Real Fy = fyP * dv; Real Fz = fzP * dv; - ib_fx += Fx; - ib_fy += Fy; - ib_fz += Fz; + Gpu::Atomic::AddNoRet( &ib_fx, Fx); + Gpu::Atomic::AddNoRet(&ib_fy, Fy); + Gpu::Atomic::AddNoRet(&ib_fz, Fz); // calc_delta(i, j, k, dxi, rho); //lagrangian to Euler for(int ii = -2; ii < +3; ii++){ @@ -427,12 +522,12 @@ void mParticle::ForceSpreading(MultiFab & EulerForce, Real dv, DELTA_FUNCTION_TYPE type){ - if (verbose) amrex::Print() << "mParticle::ForceSpreading " << std::endl; + if (verbose) amrex::Print() << "\tmParticle::ForceSpreading\n"; ib_forece.scale(0); const auto& gm = m_gdb->Geom(euler_finest_level); auto plo = gm.ProbLoArray(); auto dxi = gm.CellSizeArray(); - const int EulerForceIndex = euler_force_index; + for(mParIter pti(*this, euler_finest_level); pti.isValid(); ++pti){ const Long np = pti.numParticles(); const auto& particles = pti.GetArrayOfStructs(); @@ -445,12 +540,12 @@ void mParticle::ForceSpreading(MultiFab & EulerForce, const auto& p_ptr = particles().data(); 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); + ForceSpreading_cic(p_ptr[i], fxP_ptr[i], fyP_ptr[i], fzP_ptr[i], ib_forece[0], ib_forece[1], ib_forece[2], Uarray, euler_force_index, dv, plo, dxi, type); }); } - EulerForce.SumBoundary(EulerForceIndex, 3, gm.periodicity()); + EulerForce.SumBoundary(euler_force_index, 3, gm.periodicity()); - if (verbose) { + if (false) { // Check the Multifab // Open a file stream for output std::ofstream outFile("EulerForce.txt"); @@ -483,76 +578,54 @@ void mParticle::ForceSpreading(MultiFab & EulerForce, } -void mParticle::UpdateParticles(const amrex::MultiFab& Euler, kernel& kernel, Real dt, Real alpha_k) +void mParticle::UpdateMarkers(kernel& current_kernel, Real dt) +{ +} + +void mParticle::UpdateParticles(const MultiFab& Euler_old, + const MultiFab& Euler, + const MultiFab& pvf, + kernel& kernel, Real dt) { - if (verbose) amrex::Print() << "mParticle::UpdateParticles " << std::endl; - - //update the kernel's infomation and cal body force - for(mParIter pti(*this, euler_finest_level); pti.isValid(); ++pti){ - auto &particles = pti.GetArrayOfStructs(); - auto *p_ptr = particles.data(); - auto &attri = pti.GetAttribs(); - auto *FxP = attri[P_ATTR::Fx_Marker].data(); - auto *FyP = attri[P_ATTR::Fy_Marker].data(); - auto *FzP = attri[P_ATTR::Fz_Marker].data(); - auto *UP = attri[P_ATTR::U_Marker].data(); - auto *VP = attri[P_ATTR::V_Marker].data(); - auto *WP = attri[P_ATTR::W_Marker].data(); - const Real Dv = kernel.dv; - const Long np = pti.numParticles(); - RealVect ForceDv{std::vector{0.0,0.0,0.0}}; - RealVect Moment{std::vector{0.0,0.0,0.0}}; - auto *ForceDv_ptr = &ForceDv; - auto *Moment_ptr = &Moment; - auto *location_ptr = &kernel.location; - auto *omega_ptr = &kernel.omega; - auto *velocity_ptr = &kernel.velocity; - auto *varphi_ptr = &kernel.varphi; - const Real rho_p = kernel.rho; - //sum - amrex::ParallelFor(np, [=] - AMREX_GPU_DEVICE (int i) noexcept{ - //calculate the force - //find current particle's lagrangian marker - *ForceDv_ptr += RealVect(AMREX_D_DECL(FxP[i],FyP[i],FzP[i])) * Dv; - *Moment_ptr += (RealVect(AMREX_D_DECL(p_ptr[i].pos(0),p_ptr[i].pos(1),p_ptr[i].pos(2))) - *location_ptr).crossProduct( - RealVect(AMREX_D_DECL(FxP[i],FyP[i],FzP[i]))) * Dv; - }); - RealVect oldVelocity = kernel.velocity; - RealVect oldOmega = kernel.omega; - kernel.velocity = kernel.velocity - - 2 * alpha_k * dt/( Math::pi() * 4 * Math::powi<3>(kernel.radious) / 3) / (kernel.rho - euler_fluid_rho) * (ForceDv);// + mVector(0.0, -9.8, 0.0)); - kernel.omega = kernel.omega - - 2 * alpha_k * dt * kernel.rho / cal_momentum(kernel.rho, kernel.radious) / (kernel.rho - euler_fluid_rho) * Moment; + //continue condition 6DOF + if((kernel.TLX + kernel.TLY + kernel.TLZ + kernel.RLX + kernel.RLY + kernel.RLZ) == 0) return; + //sum + CalculateSum_cir(kernel.sum_u, Euler_old, Euler, pvf, euler_velocity_index); + + {//reduce all kernel data + amrex::ParallelAllReduce::Sum(&kernel.sum_t[0], 3, ParallelDescriptor::Communicator()); + amrex::ParallelAllReduce::Sum(&kernel.sum_u[0], 3, ParallelDescriptor::Communicator()); + amrex::ParallelAllReduce::Sum(&kernel.ib_forece[0], 3, ParallelDescriptor::Communicator()); + } + //ioprocessor calculation + if(amrex::ParallelDescriptor::MyProc() == ParallelDescriptor::IOProcessorNumber()){ - auto deltaX = alpha_k * dt * (*velocity_ptr + oldVelocity); - *location_ptr = *location_ptr + deltaX; - *varphi_ptr = *varphi_ptr + alpha_k * dt * (*omega_ptr + oldOmega); - //sum - amrex::ParallelFor(np, [=] - AMREX_GPU_DEVICE (int i) noexcept{ - //calculate the force - //find current particle's lagrangian marker - RealVect tmp = (*omega_ptr).crossProduct(*location_ptr - RealVect(p_ptr[i].pos(0), - p_ptr[i].pos(1), - p_ptr[i].pos(2))); - FxP[i] = rho_p / dt *(UP[i] + tmp[0]); - FyP[i] = rho_p / dt *(VP[i] + tmp[1]); - FzP[i] = rho_p / dt *(WP[i] + tmp[2]); - p_ptr[i].pos(0) += deltaX[0]; - p_ptr[i].pos(1) += deltaX[1]; - p_ptr[i].pos(2) += deltaX[2]; - }); + Real Vp = Math::pi() * 4 / 3 * Math::powi<3>(kernel.radius); + //update the kernel's infomation and cal body force + //update kernel velocity + kernel.velocity = (kernel.velocity + + kernel.sum_u * euler_fluid_rho / dt + - euler_fluid_rho * kernel.ib_forece * kernel.dv + + m_gravity * (kernel.rho - euler_fluid_rho) * Vp + + kernel.Fcp) * dt / kernel.rho * Vp; + //kernel.omega = ; + kernel.velocity *= RealVect(kernel.TLX, kernel.TLY, kernel.TLZ); + kernel.omega *= RealVect(kernel.RLX, kernel.RLY, kernel.RLZ); + + kernel.location += kernel.velocity * dt; } - Redistribute(); + //Bcast velocity + ParallelDescriptor::Bcast(&kernel.location[0], 3, ParallelDescriptor::IOProcessorNumber()); + if (verbose) WriteAsciiFile(amrex::Concatenate("particle", 4)); } -void mParticle::ComputeLagrangianForce(Real dt, const kernel& kernel) +void mParticle::ComputeLagrangianForce(Real dt, + const kernel& kernel) { - if (verbose) amrex::Print() << "mParticle::ComputeLagrangianForce " << std::endl; + if (verbose) amrex::Print() << "\tmParticle::ComputeLagrangianForce\n"; Real Ub = kernel.velocity[0]; Real Vb = kernel.velocity[1]; @@ -579,6 +652,13 @@ void mParticle::ComputeLagrangianForce(Real dt, const kernel& kernel) if (verbose) WriteAsciiFile(amrex::Concatenate("particle", 3)); } + +void mParticle::VelocityCorrection(amrex::MultiFab &Euler, amrex::MultiFab &EulerForce, Real dt) +{ + if(verbose) amrex::Print() << "\tmParticle::VelocityCorrection\n"; + MultiFab::Saxpy(Euler, dt, EulerForce, ParticleProperties::euler_force_index, ParticleProperties::euler_velocity_index, 3, EulerForce.nGrow()); //VelocityCorrection +} + void mParticle::WriteParticleFile(int index) { WriteAsciiFile(amrex::Concatenate("particle", index)); @@ -587,12 +667,6 @@ void mParticle::WriteParticleFile(int index) void mParticle::WriteIBForce(int step, amrex::Real time, kernel& current_kernel) { {//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()); @@ -611,7 +685,7 @@ void mParticle::WriteIBForce(int step, amrex::Real time, kernel& current_kernel) 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; + amrex::Print() << "[Particle] 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] << "," @@ -626,10 +700,13 @@ void mParticle::WriteIBForce(int step, amrex::Real time, kernel& current_kernel) /* Particles member function */ /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ void Particles::create_particles(const Geometry &gm, - const DistributionMapping & dm, - const BoxArray & ba) + const DistributionMapping & dm, + const BoxArray & ba) { particle = new mParticle(gm, dm, ba); + ParticleProperties::plo = gm.ProbLoArray(); + ParticleProperties::phi = gm.ProbHiArray(); + ParticleProperties::dx = gm.CellSizeArray(); } mParticle* Particles::get_particles() @@ -637,34 +714,66 @@ mParticle* Particles::get_particles() return particle; } -void Particles::define_para(const Vector& x, - const Vector& y, - const Vector& z, - Real rho_s, - Real radious, - Real rho_f, - int force_index, - int velocity_index, - int finest_level) -{ - std::copy(x.begin(), x.end(), std::back_inserter(_x)); - std::copy(y.begin(), y.end(), std::back_inserter(_y)); - std::copy(z.begin(), z.end(), std::back_inserter(_z)); - - euler_fluid_rho = rho_f; - euler_solid_rho = rho_s; - - particle_radious = radious; - euler_finest_level = finest_level; - euler_velocity_index = velocity_index; - euler_force_index = force_index; +void Particles::init_particle(int level, Real gravity) +{ + if(particle != nullptr){ + particle->InitParticles( + ParticleProperties::_x, + ParticleProperties::_y, + ParticleProperties::_z, + ParticleProperties::_rho, + ParticleProperties::_Vx, + ParticleProperties::_Vy, + ParticleProperties::_Vz, + ParticleProperties::_TLX, + ParticleProperties::_TLY, + ParticleProperties::_TLZ, + ParticleProperties::_RLX, + ParticleProperties::_RLY, + ParticleProperties::_RLZ, + ParticleProperties::_radius, + ParticleProperties::euler_fluid_rho, + gravity, + ParticleProperties::euler_force_index, + ParticleProperties::euler_velocity_index, + level, + ParticleProperties::verbose, + ParticleProperties::loop_time); + } } -void Particles::init_particle() +void Particles::Initialize() { - if(particle != nullptr){ - particle->InitParticles(_x, _y, _z, euler_solid_rho, particle_radious, euler_fluid_rho, - euler_force_index, euler_velocity_index, euler_finest_level); + ParmParse pp("particle"); + + std::string particle_inputfile; + + pp.get("input",particle_inputfile); + + if(!particle_inputfile.empty()){ + amrex::Print() << "[Particle] : Reading partilces cfg file : " << particle_inputfile << "\n"; + ParmParse p_file(particle_inputfile); + p_file.getarr("x", ParticleProperties::_x); + p_file.getarr("y", ParticleProperties::_y); + p_file.getarr("z", ParticleProperties::_z); + p_file.getarr("rho",ParticleProperties::_rho); + p_file.get("radius", ParticleProperties::_radius); + p_file.getarr("velocity_x", ParticleProperties::_Vx); + p_file.getarr("velocity_y", ParticleProperties::_Vy); + p_file.getarr("velocity_z", ParticleProperties::_Vz); + p_file.getarr("TLX", ParticleProperties::_TLX); + p_file.getarr("TLY", ParticleProperties::_TLY); + p_file.getarr("TLZ", ParticleProperties::_TLZ); + p_file.getarr("RLX", ParticleProperties::_RLX); + p_file.getarr("RLY", ParticleProperties::_RLY); + p_file.getarr("RLZ", ParticleProperties::_RLZ); + p_file.get("LOOP", ParticleProperties::loop_time); + p_file.query("verbose", ParticleProperties::verbose); + p_file.get("euler_velocity_index", ParticleProperties::euler_velocity_index); + p_file.get("euler_force_index", ParticleProperties::euler_force_index); + p_file.get("euler_fluid_rho", ParticleProperties::euler_fluid_rho); + }else { + amrex::Abort("[Particle] : can't read particles settings, pls check your config file \"particle.input\""); } } \ No newline at end of file diff --git a/Source/NavierStokes.cpp b/Source/NavierStokes.cpp index 096e00b5..8b39da7d 100644 --- a/Source/NavierStokes.cpp +++ b/Source/NavierStokes.cpp @@ -10,6 +10,10 @@ #include #include +#ifdef AMREX_PARTICLES +#include +#endif + #ifdef BL_USE_VELOCITY #include #include @@ -2653,7 +2657,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(parent->levelSteps(0), time, S_new, EulerForce, 10, dt, 0.5); + Particles::get_particles()->InteractWithEuler(parent->levelSteps(0), time, S_new, EulerForce, dt); } //amrex::Abort("Stop here!"); #endif diff --git a/Source/NavierStokesBase.H b/Source/NavierStokesBase.H index 4cd71270..0c62139f 100644 --- a/Source/NavierStokesBase.H +++ b/Source/NavierStokesBase.H @@ -17,7 +17,6 @@ #include #ifdef AMREX_PARTICLES -#include #include #endif diff --git a/Source/NavierStokesBase.cpp b/Source/NavierStokesBase.cpp index 6b31e77d..e8cd22a4 100644 --- a/Source/NavierStokesBase.cpp +++ b/Source/NavierStokesBase.cpp @@ -16,6 +16,9 @@ #include #include +#ifdef AMREX_PARTICLES +#include +#endif #ifdef AMREX_USE_EB #include #include @@ -24,7 +27,6 @@ #include #include #endif - #ifdef AMREX_USE_TURBULENT_FORCING #include #endif @@ -352,12 +354,9 @@ void NavierStokesBase::define_workspace() #ifdef AMREX_PARTICLES if (level == parent->finestLevel()) { // amrex::Print() << "Check grids " << grids << " " << level << " " << parent->finestLevel() << std::endl; - Particles::create_particles(geom, dmap, grids); // Class constructor - Vector x {7.5}; // 7.5 - Vector y {3.0}; - Vector z {3.0}; - Particles::define_para(x, y, z, 1.0, 0.5, 1.0, 0, 0, level); - Particles::init_particle(); + Particles::Initialize(); + Particles::create_particles(geom, dmap, grids); // Class constructor + Particles::init_particle( level, gravity); } #endif } diff --git a/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere b/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere index 9a52901d..d0a8d972 100644 --- a/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere +++ b/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere @@ -16,6 +16,9 @@ ns.fixed_dt = 1.0 ns.cfl = 0.3 ns.init_iter = 0 +# Diffused IB input file +particle.input = particle_inputs + # Refinement criterion, use vorticity and presence of tracer amr.refinement_indicators = tracer @@ -155,3 +158,50 @@ nodal_proj.proj_abs_tol = 1.e-5 mac_proj.mac_tol = 1.e-4 mac_proj.mac_abs_tol = 1.e-5 + + +############################ +# # +# Diffused IB cfg file # +# # +############################ + +# particle's location +# x = p1x p2x p3x ... +particle_inputs.x = 7.5 +particle_inputs.y = 3.0 +particle_inputs.z = 3.0 + +# particle's density +# rho = p1r p2r p3r ... +particle_inputs.rho = 1.0 + +# particle's radius +# single +particle_inputs.radius = 0.5 + +# particle's velocity +# vx = p1vx p2vx p3vx ... +particle_inputs.velocity_x = 0.0 +particle_inputs.velocity_y = 0.0 +particle_inputs.velocity_z = 0.0 + +# particle's 6DOF +# TLX = p1tl p2tl ... +particle_inputs.TLX = 0 +particle_inputs.TLY = 0 +particle_inputs.TLZ = 0 +particle_inputs.RLX = 0 +particle_inputs.RLY = 0 +particle_inputs.RLZ = 0 + +# Ns loop time +particle_inputs.LOOP = 15 + +# msg print +particle_inputs.verbose = 1 + +# Euler cfg +particle_inputs.euler_velocity_index = 0 +particle_inputs.euler_force_index = 0 +particle_inputs.euler_fluid_rho = 1.0 From 13fcabe8c892f4f21788cc4f89f534a31b6351c2 Mon Sep 17 00:00:00 2001 From: S-Explorer <1246206018@qq.com> Date: Thu, 18 Apr 2024 09:13:59 +0800 Subject: [PATCH 2/9] fixed code, format --- Source/DiffusedIB.H | 2 +- Source/DiffusedIB.cpp | 95 ++++++++++--------- Source/NS_bcfill.H | 2 +- .../FlowPastSphere/inputs.3d.flow_past_sphere | 16 ++-- 4 files changed, 59 insertions(+), 56 deletions(-) diff --git a/Source/DiffusedIB.H b/Source/DiffusedIB.H index 7f88fe2f..5583d963 100644 --- a/Source/DiffusedIB.H +++ b/Source/DiffusedIB.H @@ -153,7 +153,7 @@ private: void ForceSpreading(amrex::MultiFab &Euler, RealVect& ib_forece, Real dv, DELTA_FUNCTION_TYPE type); - void VelocityCorrection(amrex::MultiFab &Euler, amrex::MultiFab &EulerForce, Real dt); + void VelocityCorrection(amrex::MultiFab &Euler, amrex::MultiFab &EulerForce, Real dt) const; void UpdateParticles(const MultiFab& Euler_old, const MultiFab& Euler, const MultiFab& pvf, kernel& kernel, Real dt); diff --git a/Source/DiffusedIB.cpp b/Source/DiffusedIB.cpp index ed889c89..886263e5 100644 --- a/Source/DiffusedIB.cpp +++ b/Source/DiffusedIB.cpp @@ -12,6 +12,8 @@ #include namespace fs = std::filesystem; +#define GHOST_CELLS 2 + using namespace amrex; /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ @@ -20,9 +22,9 @@ using namespace amrex; namespace ParticleProperties{ Vector _x{}, _y{}, _z{}, _rho{}; - Vector _Vx{}, _Vy{}, _Vz{}; + Vector Vx{}, Vy{}, Vz{}; Real _radius; - Vector _TLX{},_TLY{},_TLZ{},_RLX{},_RLY{},_RLZ{}; + Vector TLX{}, TLY{},TLZ{},RLX{},RLY{},RLZ{}; int euler_finest_level{0}; int euler_velocity_index{0}; int euler_force_index{0}; @@ -40,7 +42,7 @@ namespace ParticleProperties{ void nodal_phi_to_pvf(MultiFab& pvf, const MultiFab& phi_nodal) { - amrex::Print() << "In the nodal_phi_to_pvf " << std::endl; + amrex::Print() << "In the nodal_phi_to_pvf\n"; #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) @@ -77,6 +79,8 @@ void nodal_phi_to_pvf(MultiFab& pvf, const MultiFab& phi_nodal) void calculate_phi_nodal(MultiFab& phi_nodal, kernel& current_kernel) { + phi_nodal.setVal(0.0); + Real Xp2 = Math::powi<2>(current_kernel.location[0]); Real Yp2 = Math::powi<2>(current_kernel.location[1]); Real Zp2 = Math::powi<2>(current_kernel.location[2]); @@ -89,9 +93,9 @@ void calculate_phi_nodal(MultiFab& phi_nodal, kernel& current_kernel) amrex::ParallelFor(bx, [=, &pnfab] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - Real Xn = i * ParticleProperties::dx[0] + 0.5 * ParticleProperties::dx[0] + ParticleProperties::plo[0]; - Real Yn = j * ParticleProperties::dx[1] + 0.5 * ParticleProperties::dx[1] + ParticleProperties::plo[1]; - Real Zn = k * ParticleProperties::dx[2] + 0.5 * ParticleProperties::dx[2] + ParticleProperties::plo[2]; + Real Xn = i * ParticleProperties::dx[0] + ParticleProperties::plo[0]; + Real Yn = j * ParticleProperties::dx[1] + ParticleProperties::plo[1]; + Real Zn = k * ParticleProperties::dx[2] + ParticleProperties::plo[2]; pnfab(i,j,k) = std::sqrt((Math::powi<2>(Xn) - Xp2)/Rp2 + (Math::powi<2>(Yn) - Yp2)/Rp2 @@ -175,7 +179,7 @@ void mParticle::InteractWithEuler(int iStep, Real dt, DELTA_FUNCTION_TYPE type){ - if (verbose) amrex::Print() << "[Particle] mParticle::InteractWithEuler " << std::endl; + if (verbose) amrex::Print() << "[Particle] mParticle::InteractWithEuler\n"; for(kernel& kernel : particle_kernels){ InitialWithLargrangianPoints(kernel); // Initialize markers for a specific particle @@ -185,7 +189,7 @@ void mParticle::InteractWithEuler(int iStep, int loop = loop_time; while(loop > 0){ if(verbose) amrex::Print() << "[Particle] Ns loop index : " << loop + 1 << "\n"; - EulerForce.setVal(0.0, euler_force_index, 3, EulerForce.nGrow()); //clear Euler force + EulerForce.setVal(0.0, euler_force_index, 3, GHOST_CELLS); //clear Euler force VelocityInterpolation(EulerVel, type); ComputeLagrangianForce(dt, kernel); ForceSpreading(EulerForce, kernel.ib_forece, kernel.dv, type); @@ -220,7 +224,7 @@ void mParticle::InitParticles(const Vector& x, { verbose = _verbose; loop_time = _loop_time; - if (verbose) amrex::Print() << "[Particle] mParticle::InitParticles " << std::endl; + if (verbose) amrex::Print() << "[Particle] mParticle::InitParticles\n"; euler_finest_level = finest_level; euler_force_index = force_index; @@ -412,7 +416,7 @@ void mParticle::VelocityInterpolation(MultiFab &EulerVel, DELTA_FUNCTION_TYPE type)// { - if (verbose) amrex::Print() << "\tmParticle::VelocityInterpolation " << std::endl; + if (verbose) amrex::Print() << "\tmParticle::VelocityInterpolation\n"; //amrex::Print() << "euler_finest_level " << euler_finest_level << std::endl; const auto& gm = m_gdb->Geom(euler_finest_level); @@ -538,7 +542,7 @@ 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, [=, &ib_forece] + amrex::ParallelFor(np, [&] 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, euler_force_index, dv, plo, dxi, type); }); @@ -566,7 +570,7 @@ void mParticle::ForceSpreading(MultiFab & EulerForce, for (int j = bx.smallEnd(1); j <= bx.bigEnd(1); ++j) { for (int i = bx.smallEnd(0); i <= bx.bigEnd(0); ++i) { // This print statement is for demonstration and should not be used in actual GPU code. - outFile << "Processing i: " << i << ", j: " << j << ", k: " << k << " " << a(i,j,k,0) << " " << a(i,j,k,1) << " " << a(i,j,k,2) << std::endl; + outFile << "Processing i: " << i << ", j: " << j << ", k: " << k << " " << a(i,j,k,0) << " " << a(i,j,k,1) << " " << a(i,j,k,2) << "\n"; } } } @@ -587,19 +591,21 @@ void mParticle::UpdateParticles(const MultiFab& Euler_old, const MultiFab& pvf, kernel& kernel, Real dt) { - if (verbose) amrex::Print() << "mParticle::UpdateParticles " << std::endl; - //continue condition 6DOF - if((kernel.TLX + kernel.TLY + kernel.TLZ + kernel.RLX + kernel.RLY + kernel.RLZ) == 0) return; - //sum - CalculateSum_cir(kernel.sum_u, Euler_old, Euler, pvf, euler_velocity_index); - + if (verbose) amrex::Print() << "mParticle::UpdateParticles\n"; + {//reduce all kernel data amrex::ParallelAllReduce::Sum(&kernel.sum_t[0], 3, ParallelDescriptor::Communicator()); amrex::ParallelAllReduce::Sum(&kernel.sum_u[0], 3, ParallelDescriptor::Communicator()); amrex::ParallelAllReduce::Sum(&kernel.ib_forece[0], 3, ParallelDescriptor::Communicator()); } + //continue condition 6DOF + if((kernel.TLX + kernel.TLY + kernel.TLZ + kernel.RLX + kernel.RLY + kernel.RLZ) == 0) return; + //sum + CalculateSum_cir(kernel.sum_u, Euler_old, Euler, pvf, euler_velocity_index); + + //ioprocessor calculation - if(amrex::ParallelDescriptor::MyProc() == ParallelDescriptor::IOProcessorNumber()){ + // if(amrex::ParallelDescriptor::MyProc() == ParallelDescriptor::IOProcessorNumber()){ Real Vp = Math::pi() * 4 / 3 * Math::powi<3>(kernel.radius); //update the kernel's infomation and cal body force @@ -614,9 +620,9 @@ void mParticle::UpdateParticles(const MultiFab& Euler_old, kernel.omega *= RealVect(kernel.RLX, kernel.RLY, kernel.RLZ); kernel.location += kernel.velocity * dt; - } + // } //Bcast velocity - ParallelDescriptor::Bcast(&kernel.location[0], 3, ParallelDescriptor::IOProcessorNumber()); + // ParallelDescriptor::Bcast(&kernel.location[0], 3, ParallelDescriptor::IOProcessorNumber()); if (verbose) WriteAsciiFile(amrex::Concatenate("particle", 4)); } @@ -653,10 +659,10 @@ void mParticle::ComputeLagrangianForce(Real dt, } -void mParticle::VelocityCorrection(amrex::MultiFab &Euler, amrex::MultiFab &EulerForce, Real dt) +void mParticle::VelocityCorrection(amrex::MultiFab &Euler, amrex::MultiFab &EulerForce, Real dt) const { if(verbose) amrex::Print() << "\tmParticle::VelocityCorrection\n"; - MultiFab::Saxpy(Euler, dt, EulerForce, ParticleProperties::euler_force_index, ParticleProperties::euler_velocity_index, 3, EulerForce.nGrow()); //VelocityCorrection + MultiFab::Saxpy(Euler, dt, EulerForce, ParticleProperties::euler_force_index, ParticleProperties::euler_velocity_index, 3, GHOST_CELLS); //VelocityCorrection } void mParticle::WriteParticleFile(int index) @@ -666,11 +672,8 @@ void mParticle::WriteParticleFile(int index) void mParticle::WriteIBForce(int step, amrex::Real time, kernel& current_kernel) { - {//reduce all kernel data - 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()); - } + amrex::ParallelAllReduce::Sum(¤t_kernel.ib_forece[0], 3, ParallelDescriptor::Communicator()); + if(amrex::ParallelDescriptor::MyProc() != ParallelDescriptor::IOProcessorNumber()) return; std::string file("IB_Force_Particle_" + std::to_string(current_kernel.id) + ".csv"); @@ -723,15 +726,15 @@ void Particles::init_particle(int level, Real gravity) ParticleProperties::_y, ParticleProperties::_z, ParticleProperties::_rho, - ParticleProperties::_Vx, - ParticleProperties::_Vy, - ParticleProperties::_Vz, - ParticleProperties::_TLX, - ParticleProperties::_TLY, - ParticleProperties::_TLZ, - ParticleProperties::_RLX, - ParticleProperties::_RLY, - ParticleProperties::_RLZ, + ParticleProperties::Vx, + ParticleProperties::Vy, + ParticleProperties::Vz, + ParticleProperties::TLX, + ParticleProperties::TLY, + ParticleProperties::TLZ, + ParticleProperties::RLX, + ParticleProperties::RLY, + ParticleProperties::RLZ, ParticleProperties::_radius, ParticleProperties::euler_fluid_rho, gravity, @@ -759,15 +762,15 @@ void Particles::Initialize() p_file.getarr("z", ParticleProperties::_z); p_file.getarr("rho",ParticleProperties::_rho); p_file.get("radius", ParticleProperties::_radius); - p_file.getarr("velocity_x", ParticleProperties::_Vx); - p_file.getarr("velocity_y", ParticleProperties::_Vy); - p_file.getarr("velocity_z", ParticleProperties::_Vz); - p_file.getarr("TLX", ParticleProperties::_TLX); - p_file.getarr("TLY", ParticleProperties::_TLY); - p_file.getarr("TLZ", ParticleProperties::_TLZ); - p_file.getarr("RLX", ParticleProperties::_RLX); - p_file.getarr("RLY", ParticleProperties::_RLY); - p_file.getarr("RLZ", ParticleProperties::_RLZ); + p_file.getarr("velocity_x", ParticleProperties::Vx); + p_file.getarr("velocity_y", ParticleProperties::Vy); + p_file.getarr("velocity_z", ParticleProperties::Vz); + p_file.getarr("TLX", ParticleProperties::TLX); + p_file.getarr("TLY", ParticleProperties::TLY); + p_file.getarr("TLZ", ParticleProperties::TLZ); + p_file.getarr("RLX", ParticleProperties::RLX); + p_file.getarr("RLY", ParticleProperties::RLY); + p_file.getarr("RLZ", ParticleProperties::RLZ); p_file.get("LOOP", ParticleProperties::loop_time); p_file.query("verbose", ParticleProperties::verbose); p_file.get("euler_velocity_index", ParticleProperties::euler_velocity_index); diff --git a/Source/NS_bcfill.H b/Source/NS_bcfill.H index 5e146fe1..538c30d9 100644 --- a/Source/NS_bcfill.H +++ b/Source/NS_bcfill.H @@ -158,7 +158,7 @@ struct velFill if (bc.lo(idir) == BCType::ext_dir && iv[idir] < domain_box.smallEnd(idir)) { dest(i,j,k,dcomp+nc) = bcv[idir][orig_comp+nc]; - if ( (probtype == 102 || probtype == 1 ) && orig_comp == 0) { // shear flow inlet in x dir + if ( (probtype == 102 || probtype == 103 ) && orig_comp == 0) { // shear flow inlet in x dir dest(i,j,k,dcomp+nc) = ub + shearrate * y; // amrex::Print() << i << " " << j << " " << k << " " << ub << " " << shearrate << " " << y << " " << idir << std::endl; } diff --git a/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere b/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere index d0a8d972..1f749bfc 100644 --- a/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere +++ b/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere @@ -7,12 +7,12 @@ # Maximum number of coarse grid timesteps to be taken, if stop_time is # not reached first. -max_step = 1 +max_step = 10 # Time at which calculation stops, if max_step is not reached first. -stop_time = 4.0 +stop_time = 10.0 -ns.fixed_dt = 1.0 +ns.fixed_dt = 0.01 ns.cfl = 0.3 ns.init_iter = 0 @@ -35,7 +35,7 @@ amr.blocking_factor = 8 # Number of cells in each coordinate direction at the coarsest level # amr.n_cell = 16 8 8 -amr.n_cell = 32 16 16 +amr.n_cell = 192 96 96 # amr.n_cell = 288 128 128 amr.max_grid_size = 16 # amr.max_grid_size = 32 @@ -68,13 +68,13 @@ amr.v = 1 #******************************************************************************* # Interval (in number of coarse timesteps) between checkpoint(restart) files -amr.check_int = 1 +amr.check_int = 100 #amr.restart = chk01400 #******************************************************************************* # Interval (in number of coarse timesteps) between plot files -amr.plot_int = 1 +amr.plot_int = 50 #******************************************************************************* @@ -168,7 +168,7 @@ mac_proj.mac_abs_tol = 1.e-5 # particle's location # x = p1x p2x p3x ... -particle_inputs.x = 7.5 +particle_inputs.x = 6.0 particle_inputs.y = 3.0 particle_inputs.z = 3.0 @@ -196,7 +196,7 @@ particle_inputs.RLY = 0 particle_inputs.RLZ = 0 # Ns loop time -particle_inputs.LOOP = 15 +particle_inputs.LOOP = 2 # msg print particle_inputs.verbose = 1 From b99d6e833cc2f9285dad9233ce62e42fc294a010 Mon Sep 17 00:00:00 2001 From: yzeng Date: Mon, 22 Apr 2024 00:19:28 -0400 Subject: [PATCH 3/9] play with xuzhu's development branch --- Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere b/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere index 1f749bfc..0cfbe526 100644 --- a/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere +++ b/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere @@ -12,7 +12,7 @@ max_step = 10 # Time at which calculation stops, if max_step is not reached first. stop_time = 10.0 -ns.fixed_dt = 0.01 +ns.fixed_dt = 0.005 ns.cfl = 0.3 ns.init_iter = 0 @@ -153,11 +153,11 @@ ns.do_diffused_ib = 1 particles.do_nspc_particles = 0 -nodal_proj.proj_tol = 1.e-4 -nodal_proj.proj_abs_tol = 1.e-5 +nodal_proj.proj_tol = 1.e-7 +nodal_proj.proj_abs_tol = 1.e-8 -mac_proj.mac_tol = 1.e-4 -mac_proj.mac_abs_tol = 1.e-5 +mac_proj.mac_tol = 1.e-7 +mac_proj.mac_abs_tol = 1.e-8 ############################ From 3cc3e798ca41a29a246eeb1e90ed3787c0ddb7f4 Mon Sep 17 00:00:00 2001 From: S-Explorer <1246206018@qq.com> Date: Thu, 18 Apr 2024 09:13:59 +0800 Subject: [PATCH 4/9] fixed code, format --- Source/DiffusedIB.H | 2 +- Source/DiffusedIB.cpp | 99 ++++++++++--------- Source/NS_bcfill.H | 2 +- .../FlowPastSphere/inputs.3d.flow_past_sphere | 16 +-- 4 files changed, 61 insertions(+), 58 deletions(-) diff --git a/Source/DiffusedIB.H b/Source/DiffusedIB.H index 7f88fe2f..5583d963 100644 --- a/Source/DiffusedIB.H +++ b/Source/DiffusedIB.H @@ -153,7 +153,7 @@ private: void ForceSpreading(amrex::MultiFab &Euler, RealVect& ib_forece, Real dv, DELTA_FUNCTION_TYPE type); - void VelocityCorrection(amrex::MultiFab &Euler, amrex::MultiFab &EulerForce, Real dt); + void VelocityCorrection(amrex::MultiFab &Euler, amrex::MultiFab &EulerForce, Real dt) const; void UpdateParticles(const MultiFab& Euler_old, const MultiFab& Euler, const MultiFab& pvf, kernel& kernel, Real dt); diff --git a/Source/DiffusedIB.cpp b/Source/DiffusedIB.cpp index ed889c89..e7cea0d8 100644 --- a/Source/DiffusedIB.cpp +++ b/Source/DiffusedIB.cpp @@ -12,6 +12,8 @@ #include namespace fs = std::filesystem; +#define GHOST_CELLS 2 + using namespace amrex; /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ @@ -20,9 +22,9 @@ using namespace amrex; namespace ParticleProperties{ Vector _x{}, _y{}, _z{}, _rho{}; - Vector _Vx{}, _Vy{}, _Vz{}; + Vector Vx{}, Vy{}, Vz{}; Real _radius; - Vector _TLX{},_TLY{},_TLZ{},_RLX{},_RLY{},_RLZ{}; + Vector TLX{}, TLY{},TLZ{},RLX{},RLY{},RLZ{}; int euler_finest_level{0}; int euler_velocity_index{0}; int euler_force_index{0}; @@ -40,7 +42,7 @@ namespace ParticleProperties{ void nodal_phi_to_pvf(MultiFab& pvf, const MultiFab& phi_nodal) { - amrex::Print() << "In the nodal_phi_to_pvf " << std::endl; + amrex::Print() << "In the nodal_phi_to_pvf\n"; #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) @@ -77,6 +79,8 @@ void nodal_phi_to_pvf(MultiFab& pvf, const MultiFab& phi_nodal) void calculate_phi_nodal(MultiFab& phi_nodal, kernel& current_kernel) { + phi_nodal.setVal(0.0); + Real Xp2 = Math::powi<2>(current_kernel.location[0]); Real Yp2 = Math::powi<2>(current_kernel.location[1]); Real Zp2 = Math::powi<2>(current_kernel.location[2]); @@ -89,9 +93,9 @@ void calculate_phi_nodal(MultiFab& phi_nodal, kernel& current_kernel) amrex::ParallelFor(bx, [=, &pnfab] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - Real Xn = i * ParticleProperties::dx[0] + 0.5 * ParticleProperties::dx[0] + ParticleProperties::plo[0]; - Real Yn = j * ParticleProperties::dx[1] + 0.5 * ParticleProperties::dx[1] + ParticleProperties::plo[1]; - Real Zn = k * ParticleProperties::dx[2] + 0.5 * ParticleProperties::dx[2] + ParticleProperties::plo[2]; + Real Xn = i * ParticleProperties::dx[0] + ParticleProperties::plo[0]; + Real Yn = j * ParticleProperties::dx[1] + ParticleProperties::plo[1]; + Real Zn = k * ParticleProperties::dx[2] + ParticleProperties::plo[2]; pnfab(i,j,k) = std::sqrt((Math::powi<2>(Xn) - Xp2)/Rp2 + (Math::powi<2>(Yn) - Yp2)/Rp2 @@ -175,7 +179,7 @@ void mParticle::InteractWithEuler(int iStep, Real dt, DELTA_FUNCTION_TYPE type){ - if (verbose) amrex::Print() << "[Particle] mParticle::InteractWithEuler " << std::endl; + if (verbose) amrex::Print() << "[Particle] mParticle::InteractWithEuler\n"; for(kernel& kernel : particle_kernels){ InitialWithLargrangianPoints(kernel); // Initialize markers for a specific particle @@ -184,8 +188,8 @@ void mParticle::InteractWithEuler(int iStep, //for 1 -> Ns int loop = loop_time; while(loop > 0){ - if(verbose) amrex::Print() << "[Particle] Ns loop index : " << loop + 1 << "\n"; - EulerForce.setVal(0.0, euler_force_index, 3, EulerForce.nGrow()); //clear Euler force + if(verbose) amrex::Print() << "[Particle] Ns loop index : " << loop << "\n"; + EulerForce.setVal(0.0, euler_force_index, 3, GHOST_CELLS); //clear Euler force VelocityInterpolation(EulerVel, type); ComputeLagrangianForce(dt, kernel); ForceSpreading(EulerForce, kernel.ib_forece, kernel.dv, type); @@ -220,7 +224,7 @@ void mParticle::InitParticles(const Vector& x, { verbose = _verbose; loop_time = _loop_time; - if (verbose) amrex::Print() << "[Particle] mParticle::InitParticles " << std::endl; + if (verbose) amrex::Print() << "[Particle] mParticle::InitParticles\n"; euler_finest_level = finest_level; euler_force_index = force_index; @@ -412,7 +416,7 @@ void mParticle::VelocityInterpolation(MultiFab &EulerVel, DELTA_FUNCTION_TYPE type)// { - if (verbose) amrex::Print() << "\tmParticle::VelocityInterpolation " << std::endl; + if (verbose) amrex::Print() << "\tmParticle::VelocityInterpolation\n"; //amrex::Print() << "euler_finest_level " << euler_finest_level << std::endl; const auto& gm = m_gdb->Geom(euler_finest_level); @@ -493,7 +497,7 @@ void ForceSpreading_cic (P const& p, Real Fx = fxP * dv; Real Fy = fyP * dv; Real Fz = fzP * dv; - Gpu::Atomic::AddNoRet( &ib_fx, Fx); + Gpu::Atomic::AddNoRet(&ib_fx, Fx); Gpu::Atomic::AddNoRet(&ib_fy, Fy); Gpu::Atomic::AddNoRet(&ib_fz, Fz); // calc_delta(i, j, k, dxi, rho); @@ -538,7 +542,7 @@ 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, [=, &ib_forece] + amrex::ParallelFor(np, [&] 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, euler_force_index, dv, plo, dxi, type); }); @@ -566,7 +570,7 @@ void mParticle::ForceSpreading(MultiFab & EulerForce, for (int j = bx.smallEnd(1); j <= bx.bigEnd(1); ++j) { for (int i = bx.smallEnd(0); i <= bx.bigEnd(0); ++i) { // This print statement is for demonstration and should not be used in actual GPU code. - outFile << "Processing i: " << i << ", j: " << j << ", k: " << k << " " << a(i,j,k,0) << " " << a(i,j,k,1) << " " << a(i,j,k,2) << std::endl; + outFile << "Processing i: " << i << ", j: " << j << ", k: " << k << " " << a(i,j,k,0) << " " << a(i,j,k,1) << " " << a(i,j,k,2) << "\n"; } } } @@ -587,19 +591,21 @@ void mParticle::UpdateParticles(const MultiFab& Euler_old, const MultiFab& pvf, kernel& kernel, Real dt) { - if (verbose) amrex::Print() << "mParticle::UpdateParticles " << std::endl; - //continue condition 6DOF - if((kernel.TLX + kernel.TLY + kernel.TLZ + kernel.RLX + kernel.RLY + kernel.RLZ) == 0) return; - //sum - CalculateSum_cir(kernel.sum_u, Euler_old, Euler, pvf, euler_velocity_index); - + if (verbose) amrex::Print() << "mParticle::UpdateParticles\n"; + {//reduce all kernel data amrex::ParallelAllReduce::Sum(&kernel.sum_t[0], 3, ParallelDescriptor::Communicator()); amrex::ParallelAllReduce::Sum(&kernel.sum_u[0], 3, ParallelDescriptor::Communicator()); amrex::ParallelAllReduce::Sum(&kernel.ib_forece[0], 3, ParallelDescriptor::Communicator()); } + //continue condition 6DOF + if((kernel.TLX + kernel.TLY + kernel.TLZ + kernel.RLX + kernel.RLY + kernel.RLZ) == 0) return; + //sum + CalculateSum_cir(kernel.sum_u, Euler_old, Euler, pvf, euler_velocity_index); + + //ioprocessor calculation - if(amrex::ParallelDescriptor::MyProc() == ParallelDescriptor::IOProcessorNumber()){ + // if(amrex::ParallelDescriptor::MyProc() == ParallelDescriptor::IOProcessorNumber()){ Real Vp = Math::pi() * 4 / 3 * Math::powi<3>(kernel.radius); //update the kernel's infomation and cal body force @@ -614,9 +620,9 @@ void mParticle::UpdateParticles(const MultiFab& Euler_old, kernel.omega *= RealVect(kernel.RLX, kernel.RLY, kernel.RLZ); kernel.location += kernel.velocity * dt; - } + // } //Bcast velocity - ParallelDescriptor::Bcast(&kernel.location[0], 3, ParallelDescriptor::IOProcessorNumber()); + // ParallelDescriptor::Bcast(&kernel.location[0], 3, ParallelDescriptor::IOProcessorNumber()); if (verbose) WriteAsciiFile(amrex::Concatenate("particle", 4)); } @@ -653,10 +659,10 @@ void mParticle::ComputeLagrangianForce(Real dt, } -void mParticle::VelocityCorrection(amrex::MultiFab &Euler, amrex::MultiFab &EulerForce, Real dt) +void mParticle::VelocityCorrection(amrex::MultiFab &Euler, amrex::MultiFab &EulerForce, Real dt) const { if(verbose) amrex::Print() << "\tmParticle::VelocityCorrection\n"; - MultiFab::Saxpy(Euler, dt, EulerForce, ParticleProperties::euler_force_index, ParticleProperties::euler_velocity_index, 3, EulerForce.nGrow()); //VelocityCorrection + MultiFab::Saxpy(Euler, dt, EulerForce, ParticleProperties::euler_force_index, ParticleProperties::euler_velocity_index, 3, GHOST_CELLS); //VelocityCorrection } void mParticle::WriteParticleFile(int index) @@ -666,11 +672,8 @@ void mParticle::WriteParticleFile(int index) void mParticle::WriteIBForce(int step, amrex::Real time, kernel& current_kernel) { - {//reduce all kernel data - 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()); - } + amrex::ParallelAllReduce::Sum(¤t_kernel.ib_forece[0], 3, ParallelDescriptor::Communicator()); + if(amrex::ParallelDescriptor::MyProc() != ParallelDescriptor::IOProcessorNumber()) return; std::string file("IB_Force_Particle_" + std::to_string(current_kernel.id) + ".csv"); @@ -723,15 +726,15 @@ void Particles::init_particle(int level, Real gravity) ParticleProperties::_y, ParticleProperties::_z, ParticleProperties::_rho, - ParticleProperties::_Vx, - ParticleProperties::_Vy, - ParticleProperties::_Vz, - ParticleProperties::_TLX, - ParticleProperties::_TLY, - ParticleProperties::_TLZ, - ParticleProperties::_RLX, - ParticleProperties::_RLY, - ParticleProperties::_RLZ, + ParticleProperties::Vx, + ParticleProperties::Vy, + ParticleProperties::Vz, + ParticleProperties::TLX, + ParticleProperties::TLY, + ParticleProperties::TLZ, + ParticleProperties::RLX, + ParticleProperties::RLY, + ParticleProperties::RLZ, ParticleProperties::_radius, ParticleProperties::euler_fluid_rho, gravity, @@ -759,15 +762,15 @@ void Particles::Initialize() p_file.getarr("z", ParticleProperties::_z); p_file.getarr("rho",ParticleProperties::_rho); p_file.get("radius", ParticleProperties::_radius); - p_file.getarr("velocity_x", ParticleProperties::_Vx); - p_file.getarr("velocity_y", ParticleProperties::_Vy); - p_file.getarr("velocity_z", ParticleProperties::_Vz); - p_file.getarr("TLX", ParticleProperties::_TLX); - p_file.getarr("TLY", ParticleProperties::_TLY); - p_file.getarr("TLZ", ParticleProperties::_TLZ); - p_file.getarr("RLX", ParticleProperties::_RLX); - p_file.getarr("RLY", ParticleProperties::_RLY); - p_file.getarr("RLZ", ParticleProperties::_RLZ); + p_file.getarr("velocity_x", ParticleProperties::Vx); + p_file.getarr("velocity_y", ParticleProperties::Vy); + p_file.getarr("velocity_z", ParticleProperties::Vz); + p_file.getarr("TLX", ParticleProperties::TLX); + p_file.getarr("TLY", ParticleProperties::TLY); + p_file.getarr("TLZ", ParticleProperties::TLZ); + p_file.getarr("RLX", ParticleProperties::RLX); + p_file.getarr("RLY", ParticleProperties::RLY); + p_file.getarr("RLZ", ParticleProperties::RLZ); p_file.get("LOOP", ParticleProperties::loop_time); p_file.query("verbose", ParticleProperties::verbose); p_file.get("euler_velocity_index", ParticleProperties::euler_velocity_index); diff --git a/Source/NS_bcfill.H b/Source/NS_bcfill.H index 5e146fe1..538c30d9 100644 --- a/Source/NS_bcfill.H +++ b/Source/NS_bcfill.H @@ -158,7 +158,7 @@ struct velFill if (bc.lo(idir) == BCType::ext_dir && iv[idir] < domain_box.smallEnd(idir)) { dest(i,j,k,dcomp+nc) = bcv[idir][orig_comp+nc]; - if ( (probtype == 102 || probtype == 1 ) && orig_comp == 0) { // shear flow inlet in x dir + if ( (probtype == 102 || probtype == 103 ) && orig_comp == 0) { // shear flow inlet in x dir dest(i,j,k,dcomp+nc) = ub + shearrate * y; // amrex::Print() << i << " " << j << " " << k << " " << ub << " " << shearrate << " " << y << " " << idir << std::endl; } diff --git a/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere b/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere index d0a8d972..1f749bfc 100644 --- a/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere +++ b/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere @@ -7,12 +7,12 @@ # Maximum number of coarse grid timesteps to be taken, if stop_time is # not reached first. -max_step = 1 +max_step = 10 # Time at which calculation stops, if max_step is not reached first. -stop_time = 4.0 +stop_time = 10.0 -ns.fixed_dt = 1.0 +ns.fixed_dt = 0.01 ns.cfl = 0.3 ns.init_iter = 0 @@ -35,7 +35,7 @@ amr.blocking_factor = 8 # Number of cells in each coordinate direction at the coarsest level # amr.n_cell = 16 8 8 -amr.n_cell = 32 16 16 +amr.n_cell = 192 96 96 # amr.n_cell = 288 128 128 amr.max_grid_size = 16 # amr.max_grid_size = 32 @@ -68,13 +68,13 @@ amr.v = 1 #******************************************************************************* # Interval (in number of coarse timesteps) between checkpoint(restart) files -amr.check_int = 1 +amr.check_int = 100 #amr.restart = chk01400 #******************************************************************************* # Interval (in number of coarse timesteps) between plot files -amr.plot_int = 1 +amr.plot_int = 50 #******************************************************************************* @@ -168,7 +168,7 @@ mac_proj.mac_abs_tol = 1.e-5 # particle's location # x = p1x p2x p3x ... -particle_inputs.x = 7.5 +particle_inputs.x = 6.0 particle_inputs.y = 3.0 particle_inputs.z = 3.0 @@ -196,7 +196,7 @@ particle_inputs.RLY = 0 particle_inputs.RLZ = 0 # Ns loop time -particle_inputs.LOOP = 15 +particle_inputs.LOOP = 2 # msg print particle_inputs.verbose = 1 From 39b217bf0eb18208cef58c5194607a3967205dae Mon Sep 17 00:00:00 2001 From: S-Explorer <1246206018@qq.com> Date: Tue, 23 Apr 2024 23:06:18 +0800 Subject: [PATCH 5/9] plot DiffusedIB marker's information to plt file --- Source/NavierStokes.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/Source/NavierStokes.cpp b/Source/NavierStokes.cpp index 8b39da7d..48ae5e36 100644 --- a/Source/NavierStokes.cpp +++ b/Source/NavierStokes.cpp @@ -1155,6 +1155,12 @@ NavierStokes::writePlotFilePost (const std::string& dir, } #endif +#ifdef AMREX_PARTICLES + if(level == parent->finestLevel()){ + Particles::get_particles()->Checkpoint(dir, "particles"); + } +#endif + #ifdef AMREX_USE_EB if ( set_plot_coveredCell_val ) { From 8665786ea4606cc01c79f78222ce4b9f5a9201e2 Mon Sep 17 00:00:00 2001 From: ruohai0925 Date: Tue, 23 Apr 2024 19:43:16 -0700 Subject: [PATCH 6/9] debug the diffusedib code --- .../FlowPastSphere/inputs.3d.flow_past_sphere | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere b/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere index 0cfbe526..05172ce8 100644 --- a/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere +++ b/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere @@ -7,7 +7,7 @@ # Maximum number of coarse grid timesteps to be taken, if stop_time is # not reached first. -max_step = 10 +max_step = 1 # Time at which calculation stops, if max_step is not reached first. stop_time = 10.0 @@ -74,7 +74,7 @@ amr.check_int = 100 #******************************************************************************* # Interval (in number of coarse timesteps) between plot files -amr.plot_int = 50 +amr.plot_int = 100 #******************************************************************************* @@ -85,7 +85,7 @@ ns.vel_visc_coef = 0.01 #******************************************************************************* # Diffusion coefficient for first scalar -ns.scal_diff_coefs = 0.001 +ns.scal_diff_coefs = 0.0 #******************************************************************************* @@ -153,11 +153,11 @@ ns.do_diffused_ib = 1 particles.do_nspc_particles = 0 -nodal_proj.proj_tol = 1.e-7 -nodal_proj.proj_abs_tol = 1.e-8 +nodal_proj.proj_tol = 1.e-8 +nodal_proj.proj_abs_tol = 1.e-9 -mac_proj.mac_tol = 1.e-7 -mac_proj.mac_abs_tol = 1.e-8 +mac_proj.mac_tol = 1.e-8 +mac_proj.mac_abs_tol = 1.e-9 ############################ From c026199fc97401d96117b888489be475cf2a9ba1 Mon Sep 17 00:00:00 2001 From: ruohai0925 Date: Tue, 23 Apr 2024 21:26:41 -0700 Subject: [PATCH 7/9] update the diffusedib codes --- Source/DiffusedIB.H | 6 +-- Source/DiffusedIB.cpp | 37 +++++++++++-------- .../FlowPastSphere/inputs.3d.flow_past_sphere | 6 +-- 3 files changed, 27 insertions(+), 22 deletions(-) diff --git a/Source/DiffusedIB.H b/Source/DiffusedIB.H index 5583d963..3cf4c7fd 100644 --- a/Source/DiffusedIB.H +++ b/Source/DiffusedIB.H @@ -52,7 +52,7 @@ struct kernel{ Real rho{0.0}; int ml{0}; Real dv{0.0}; - RealVect ib_forece{0.0,0.0,0.0}; + RealVect ib_force{0.0,0.0,0.0}; std::vector thetaK; std::vector phiK; @@ -151,7 +151,7 @@ private: void ComputeLagrangianForce(Real dt, const kernel& kernel); - void ForceSpreading(amrex::MultiFab &Euler, RealVect& ib_forece, Real dv, DELTA_FUNCTION_TYPE type); + void ForceSpreading(amrex::MultiFab &Euler, RealVect& ib_force, Real dv, DELTA_FUNCTION_TYPE type); void VelocityCorrection(amrex::MultiFab &Euler, amrex::MultiFab &EulerForce, Real dt) const; @@ -159,7 +159,7 @@ private: static void WriteIBForce(int step, amrex::Real time, kernel& current_kernel); - uint32_t ib_forece_file_index = 0; + uint32_t ib_force_file_index = 0; RealVect m_gravity{0.0,0.0,0.0}; diff --git a/Source/DiffusedIB.cpp b/Source/DiffusedIB.cpp index 886263e5..db3b3d13 100644 --- a/Source/DiffusedIB.cpp +++ b/Source/DiffusedIB.cpp @@ -183,20 +183,25 @@ void mParticle::InteractWithEuler(int iStep, for(kernel& kernel : particle_kernels){ InitialWithLargrangianPoints(kernel); // Initialize markers for a specific particle - Redistribute(); - UpdateMarkers(kernel, dt); + //UpdateMarkers(kernel, dt); //for 1 -> Ns int loop = loop_time; while(loop > 0){ - if(verbose) amrex::Print() << "[Particle] Ns loop index : " << loop + 1 << "\n"; - EulerForce.setVal(0.0, euler_force_index, 3, GHOST_CELLS); //clear Euler force + if(verbose) amrex::Print() << "[Particle] Ns loop index : " << loop << "\n"; + VelocityInterpolation(EulerVel, type); ComputeLagrangianForce(dt, kernel); - ForceSpreading(EulerForce, kernel.ib_forece, kernel.dv, type); + + EulerForce.setVal(0.0, euler_force_index, 3, GHOST_CELLS); // clear Euler force + kernel.ib_force.scale(0); // clear kernel ib_force + ForceSpreading(EulerForce, kernel.ib_force, kernel.dv, type); + + WriteIBForce(iStep, time, kernel); VelocityCorrection(EulerVel, EulerForce, dt); + loop--; } - WriteIBForce(iStep, time, kernel); + // WriteIBForce(iStep, time, kernel); } } @@ -357,6 +362,7 @@ void mParticle::InitialWithLargrangianPoints(const kernel& current_kernel){ ); } // Redistribute the markers after updating their locations + Redistribute(); if (verbose) WriteAsciiFile(amrex::Concatenate("particle", 1)); } @@ -485,7 +491,7 @@ void ForceSpreading_cic (P const& p, GpuArray const& dx, DELTA_FUNCTION_TYPE type) { - const Real d = AMREX_D_TERM(dx[0], *dx[1], *dx[2]); + //const Real d = AMREX_D_TERM(dx[0], *dx[1], *dx[2]); //plo to ii jj kk Real lx = (p.pos(0) - plo[0]) / dx[0]; Real ly = (p.pos(1) - plo[1]) / dx[1]; @@ -497,7 +503,7 @@ void ForceSpreading_cic (P const& p, Real Fx = fxP * dv; Real Fy = fyP * dv; Real Fz = fzP * dv; - Gpu::Atomic::AddNoRet( &ib_fx, Fx); + Gpu::Atomic::AddNoRet(&ib_fx, Fx); Gpu::Atomic::AddNoRet(&ib_fy, Fy); Gpu::Atomic::AddNoRet(&ib_fz, Fz); // calc_delta(i, j, k, dxi, rho); @@ -522,12 +528,11 @@ void ForceSpreading_cic (P const& p, } void mParticle::ForceSpreading(MultiFab & EulerForce, - RealVect& ib_forece, + RealVect& ib_force, Real dv, DELTA_FUNCTION_TYPE type){ if (verbose) amrex::Print() << "\tmParticle::ForceSpreading\n"; - ib_forece.scale(0); const auto& gm = m_gdb->Geom(euler_finest_level); auto plo = gm.ProbLoArray(); auto dxi = gm.CellSizeArray(); @@ -544,7 +549,7 @@ void mParticle::ForceSpreading(MultiFab & EulerForce, 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], ib_forece[0], ib_forece[1], ib_forece[2], Uarray, euler_force_index, dv, plo, dxi, type); + ForceSpreading_cic(p_ptr[i], fxP_ptr[i], fyP_ptr[i], fzP_ptr[i], ib_force[0], ib_force[1], ib_force[2], Uarray, euler_force_index, dv, plo, dxi, type); }); } EulerForce.SumBoundary(euler_force_index, 3, gm.periodicity()); @@ -596,7 +601,7 @@ void mParticle::UpdateParticles(const MultiFab& Euler_old, {//reduce all kernel data amrex::ParallelAllReduce::Sum(&kernel.sum_t[0], 3, ParallelDescriptor::Communicator()); amrex::ParallelAllReduce::Sum(&kernel.sum_u[0], 3, ParallelDescriptor::Communicator()); - amrex::ParallelAllReduce::Sum(&kernel.ib_forece[0], 3, ParallelDescriptor::Communicator()); + amrex::ParallelAllReduce::Sum(&kernel.ib_force[0], 3, ParallelDescriptor::Communicator()); } //continue condition 6DOF if((kernel.TLX + kernel.TLY + kernel.TLZ + kernel.RLX + kernel.RLY + kernel.RLZ) == 0) return; @@ -612,7 +617,7 @@ void mParticle::UpdateParticles(const MultiFab& Euler_old, //update kernel velocity kernel.velocity = (kernel.velocity + kernel.sum_u * euler_fluid_rho / dt - - euler_fluid_rho * kernel.ib_forece * kernel.dv + - euler_fluid_rho * kernel.ib_force * kernel.dv + m_gravity * (kernel.rho - euler_fluid_rho) * Vp + kernel.Fcp) * dt / kernel.rho * Vp; //kernel.omega = ; @@ -662,7 +667,7 @@ void mParticle::ComputeLagrangianForce(Real dt, void mParticle::VelocityCorrection(amrex::MultiFab &Euler, amrex::MultiFab &EulerForce, Real dt) const { if(verbose) amrex::Print() << "\tmParticle::VelocityCorrection\n"; - MultiFab::Saxpy(Euler, dt, EulerForce, ParticleProperties::euler_force_index, ParticleProperties::euler_velocity_index, 3, GHOST_CELLS); //VelocityCorrection + MultiFab::Saxpy(Euler, dt, EulerForce, ParticleProperties::euler_force_index, ParticleProperties::euler_velocity_index, 3, 0); //VelocityCorrection } void mParticle::WriteParticleFile(int index) @@ -672,7 +677,7 @@ void mParticle::WriteParticleFile(int index) void mParticle::WriteIBForce(int step, amrex::Real time, kernel& current_kernel) { - amrex::ParallelAllReduce::Sum(¤t_kernel.ib_forece[0], 3, ParallelDescriptor::Communicator()); + amrex::ParallelAllReduce::Sum(¤t_kernel.ib_force[0], 3, ParallelDescriptor::Communicator()); if(amrex::ParallelDescriptor::MyProc() != ParallelDescriptor::IOProcessorNumber()) return; @@ -694,7 +699,7 @@ void mParticle::WriteIBForce(int step, amrex::Real time, kernel& current_kernel) << 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"; + << current_kernel.ib_force[0] << "," << current_kernel.ib_force[1] << "," << current_kernel.ib_force[2] << "\n"; } out_ib_force.close(); } diff --git a/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere b/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere index 05172ce8..a35d66da 100644 --- a/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere +++ b/Tutorials/FlowPastSphere/inputs.3d.flow_past_sphere @@ -68,13 +68,13 @@ amr.v = 1 #******************************************************************************* # Interval (in number of coarse timesteps) between checkpoint(restart) files -amr.check_int = 100 +amr.check_int = 1 #amr.restart = chk01400 #******************************************************************************* # Interval (in number of coarse timesteps) between plot files -amr.plot_int = 100 +amr.plot_int = 1 #******************************************************************************* @@ -196,7 +196,7 @@ particle_inputs.RLY = 0 particle_inputs.RLZ = 0 # Ns loop time -particle_inputs.LOOP = 2 +particle_inputs.LOOP = 1 # msg print particle_inputs.verbose = 1 From 1c702310aee27c38e17d2b89af1b5f9802011e8c Mon Sep 17 00:00:00 2001 From: ruohai0925 Date: Tue, 23 Apr 2024 22:15:17 -0700 Subject: [PATCH 8/9] fix the bug of 32 MPIs, always need to use redistribute while setting markerP.pos --- Source/DiffusedIB.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/Source/DiffusedIB.cpp b/Source/DiffusedIB.cpp index db3b3d13..e869e115 100644 --- a/Source/DiffusedIB.cpp +++ b/Source/DiffusedIB.cpp @@ -186,6 +186,7 @@ void mParticle::InteractWithEuler(int iStep, //UpdateMarkers(kernel, dt); //for 1 -> Ns int loop = loop_time; + BL_ASSERT(loop > 0); while(loop > 0){ if(verbose) amrex::Print() << "[Particle] Ns loop index : " << loop << "\n"; @@ -196,7 +197,9 @@ void mParticle::InteractWithEuler(int iStep, kernel.ib_force.scale(0); // clear kernel ib_force ForceSpreading(EulerForce, kernel.ib_force, kernel.dv, type); - WriteIBForce(iStep, time, kernel); + if (loop == loop_time) { + WriteIBForce(iStep, time, kernel); + } VelocityCorrection(EulerVel, EulerForce, dt); loop--; @@ -321,7 +324,7 @@ void mParticle::InitParticles(const Vector& x, particleTileTmp.push_back_real(Marker_attr); } } - // Redistribute(); // No need to redistribute here! + Redistribute(); // Still needs to redistribute here! if (verbose) WriteAsciiFile(amrex::Concatenate("particle", 0)); } From b832939aee5aa2b6e3b80ff7ce8bb37824a90e9d Mon Sep 17 00:00:00 2001 From: AoooE <32592710+S-Explorer@users.noreply.github.com> Date: Wed, 24 Apr 2024 13:33:06 +0800 Subject: [PATCH 9/9] Update DiffusedIB.cpp DiffusedIB.cpp line 607 : forece -> force --- Source/DiffusedIB.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/DiffusedIB.cpp b/Source/DiffusedIB.cpp index 56ecfa40..5cfcec50 100644 --- a/Source/DiffusedIB.cpp +++ b/Source/DiffusedIB.cpp @@ -604,7 +604,7 @@ void mParticle::UpdateParticles(const MultiFab& Euler_old, {//reduce all kernel data amrex::ParallelAllReduce::Sum(&kernel.sum_t[0], 3, ParallelDescriptor::Communicator()); amrex::ParallelAllReduce::Sum(&kernel.sum_u[0], 3, ParallelDescriptor::Communicator()); - amrex::ParallelAllReduce::Sum(&kernel.ib_forece[0], 3, ParallelDescriptor::Communicator()); + amrex::ParallelAllReduce::Sum(&kernel.ib_force[0], 3, ParallelDescriptor::Communicator()); } //continue condition 6DOF if((kernel.TLX + kernel.TLY + kernel.TLZ + kernel.RLX + kernel.RLY + kernel.RLZ) == 0) return; @@ -787,4 +787,4 @@ void Particles::Initialize() }else { amrex::Abort("[Particle] : can't read particles settings, pls check your config file \"particle.input\""); } -} \ No newline at end of file +}