diff --git a/apps/benchmarks/MaterialTransport/thermal_PSM/thermal_PSM_particles/benchmark.prm b/apps/benchmarks/MaterialTransport/thermal_PSM/thermal_PSM_particles/benchmark.prm
index eccd61943f476346a8fc1b436916d8fce7a39070..2a64635fb6d22d7275674894d43c5edb1c99fa10 100644
--- a/apps/benchmarks/MaterialTransport/thermal_PSM/thermal_PSM_particles/benchmark.prm
+++ b/apps/benchmarks/MaterialTransport/thermal_PSM/thermal_PSM_particles/benchmark.prm
@@ -2,11 +2,11 @@
NumericalSetup
{
// product of number of blocks should be equal to number of used processes
- numXBlocks 2;
+ numXBlocks 1;
numYBlocks 1;
- numZBlocks 2;
+ numZBlocks 1;
- cellsPerBlockPerDirection 100;
+ cellsPerBlockPerDirection 32;
periodicInY false;
periodicInZ false;
@@ -19,22 +19,22 @@ NumericalSetup
Thot 1; // ThotSI = ThotLB (concentration conversion from SI to LB)
Tcold 0;
PrandtlNumber 0.71;
- RayleighNumber 1000000;
- Mach 0.01;
+ RayleighNumber 100000;
+ Mach 0.1;
relaxationRate 1.5;
- timeSteps 50000;
+ timeSteps 1000;
// particle distribution in LBM units
useParticles true; // if true, PSM/particle mapping/velocity computation/hydrodynamic force reduction is used, else LBM is used
particleDiameter 10.0;
- particleGenerationSpacing 11.0;
+ particleGenerationSpacing 5.0;
generationDomainFraction <0.8, 1.0, 1.0>; // fraction of the domain where particles are generated
generationPointOfReferenceOffset <0, 0, 0>; // offset of point of reference from domain center, see SCIterator.h
useParticleOffset true; // offset every second particle layer in flow direction
- particleSubBlockSize <10, 10, 10>;
+ particleSubBlockSize <4, 4, 4>;
// fluid quantities in LBM units
uInflow 0.00008;
@@ -42,7 +42,7 @@ NumericalSetup
Output
{
- vtkSpacing 0;
+ vtkSpacing 1;
vtkFolder vtk_out;
performanceLogFrequency 10;
diff --git a/apps/benchmarks/MaterialTransport/thermal_PSM/thermal_PSM_particles/thermalPSMCodegenParticles.py b/apps/benchmarks/MaterialTransport/thermal_PSM/thermal_PSM_particles/thermalPSMCodegenParticles.py
index 6ef8d01434abb9cec639f9a6e9d9c67e812e5186..0add80ff7edfe269b1ed7f311847934b53edeec3 100644
--- a/apps/benchmarks/MaterialTransport/thermal_PSM/thermal_PSM_particles/thermalPSMCodegenParticles.py
+++ b/apps/benchmarks/MaterialTransport/thermal_PSM/thermal_PSM_particles/thermalPSMCodegenParticles.py
@@ -46,7 +46,7 @@ const bool infoCsePdfs = {cse_pdfs};
with CodeGeneration() as ctx:
data_type = "float64" if ctx.double_accuracy else "float32"
stencil_fluid = LBStencil(Stencil.D3Q19)
- stencil_concentration = LBStencil(Stencil.D3Q7)
+ stencil_concentration = LBStencil(Stencil.D3Q27)
omega = sp.Symbol("omega") # for now same for both the sweeps
init_density_fluid = sp.Symbol("init_density_fluid")
init_density_concentration = sp.Symbol("init_density_concentration")
@@ -101,13 +101,14 @@ with CodeGeneration() as ctx:
f"particle_v({MaxParticlesPerCell * stencil_fluid.D}), particle_f({MaxParticlesPerCell * stencil_fluid.D}), Bs({MaxParticlesPerCell}): {data_type}[3D]",
layout=layout,
)
+ particle_temperatures = ps.fields(f"particle_t({MaxParticlesPerCell}) :{data_type}[3D]",layout=layout)
# Solid fraction field
B = ps.fields(f"b({1}): {data_type}[3D]", layout=layout)
- #force_concentration_on_fluid = sp.Matrix([0, (rho_0)*alpha*(concentration_field.center - T0)*gravity_LBM,0])
- force_concentration_on_fluid = sp.Matrix([0, 0,(rho_0)*alpha*(concentration_field.center - T0)*gravity_LBM])
+ force_concentration_on_fluid = sp.Matrix([0, (rho_0)*alpha*(concentration_field.center - T0)*gravity_LBM,0])
+ #force_concentration_on_fluid = sp.Matrix([0, 0,(rho_0)*alpha*(concentration_field.center - T0)*gravity_LBM])
# Fluid LBM optimisation
lbm_fluid_opt = LBMOptimisation(
@@ -132,8 +133,8 @@ with CodeGeneration() as ctx:
relaxation_rate=omega_f,
#relaxation_rates=[0,1,1,Sv,Sv,Sv,Sq,Sq,Sv],
output={"velocity": velocity_field},
- #force= force_concentration_on_fluid,
- #force_model=ForceModel.GUO,
+ force= force_concentration_on_fluid,
+ force_model=ForceModel.GUO,
compressible=True,
)
@@ -163,7 +164,7 @@ with CodeGeneration() as ctx:
stencil=stencil_concentration,
method=Method.MRT,
relaxation_rates=[1,qk,qk,qk,qe,qe,qe], # MRT 2
- #relaxation_rate=omega,
+ #relaxation_rate=omega_c,
velocity_input=velocity_field,
output={"density": concentration_field},
compressible=True,
@@ -178,6 +179,7 @@ with CodeGeneration() as ctx:
MaxParticlesPerCell=MaxParticlesPerCell,
individual_fraction_field=Bs,
particle_force_field=None,
+ particle_temperature_field=particle_temperatures,
)
psm_concentration_config = LBMConfig(
@@ -186,7 +188,6 @@ with CodeGeneration() as ctx:
relaxation_rate=omega_c,
velocity_input=velocity_field,
output={"density": concentration_field},
- #force_model=ForceModel.LUO,
compressible=True,
psm_config=psm_config_C,
)
@@ -199,7 +200,7 @@ with CodeGeneration() as ctx:
# =====================
method_fluid = create_lb_method(lbm_config=psm_fluid_config)
- method_concentration = create_lb_method(lbm_config=lbm_concentration_config)
+ method_concentration = create_lb_method(lbm_config=psm_concentration_config)
init_velocity = sp.symbols("init_velocity_:3")
pdfs_fluid_setter = macroscopic_values_setter(
method_fluid, density=init_density_fluid, velocity=init_velocity, pdfs=pdfs_fluid.center_vector
@@ -210,16 +211,29 @@ with CodeGeneration() as ctx:
)
# Use average velocity of all intersecting particles when setting PDFs (mandatory for SC=3)
+ rhs = []
for i, sub_exp in enumerate(pdfs_fluid_setter.subexpressions[-3:]):
- rhs = []
- for summand in sub_exp.rhs.args:
- rhs.append(summand * (1.0 - B.center))
+ if(len(sub_exp.rhs.args) > 0):
+ for summand in (sub_exp.rhs.args):
+ rhs.append(summand * (1.0 - B.center))
+ else:
+ rhs.append(sub_exp.rhs * (1.0 - B.center))
for p in range(2):
rhs.append(particle_velocities(p * stencil_fluid.D + i) * Bs.center(p))
pdfs_fluid_setter.subexpressions.remove(sub_exp)
pdfs_fluid_setter.subexpressions.append(Assignment(sub_exp.lhs, Add(*rhs)))
+ rhs = []
+ # Use average temperature of all intersecting particles when setting PDFs (mandatory for SC=3)
+ sub_exp_con = pdfs_concentration_setter.subexpressions[0]
+ rhs_con = []
+ rhs_con.append(sub_exp_con.rhs * (1.0 - B.center))
+ for p in range(2):
+ rhs_con.append(particle_temperatures(p) * Bs.center(p))
+ pdfs_concentration_setter.subexpressions.remove(sub_exp_con)
+ pdfs_concentration_setter.subexpressions.append(Assignment(sub_exp_con.lhs, Add(*rhs_con)))
+ print("con setter after manip ", pdfs_concentration_setter.subexpressions)
# specify the target
if ctx.gpu:
diff --git a/apps/benchmarks/MaterialTransport/thermal_PSM/thermal_PSM_particles/thermalPSMParticles.cpp b/apps/benchmarks/MaterialTransport/thermal_PSM/thermal_PSM_particles/thermalPSMParticles.cpp
index 1fbb2002ba41312543d7a75399330a17d59174b2..9cb01d6330a41d8a1df69c5503b83c496d6de6e7 100644
--- a/apps/benchmarks/MaterialTransport/thermal_PSM/thermal_PSM_particles/thermalPSMParticles.cpp
+++ b/apps/benchmarks/MaterialTransport/thermal_PSM/thermal_PSM_particles/thermalPSMParticles.cpp
@@ -30,6 +30,7 @@
#include "core/grid_generator/SCIterator.h"
#include "core/logging/all.h"
#include "core/timing/RemainingTimeLogger.h"
+#include "core/math/all.h"
#include "field/AddToStorage.h"
#include "field/vtk/all.h"
@@ -189,6 +190,7 @@ int main(int argc, char** argv)
const real_t rho_0 = 1.0;
const real_t omega_f = lbm::collision_model::omegaFromViscosity(kinematicViscosityLB);
const real_t omega_t = lbm::collision_model::omegaFromViscosity(thermalDiffusivityLB);
+ const real_t omega_c = (omega_t);
const real_t qk = 1/(4*thermalDiffusivityLB + 0.5);
const real_t qe = 1.0;
// calculations for verification and correctness
@@ -256,17 +258,11 @@ int main(int argc, char** argv)
simulationDomain.yMax() * (real_t(1) + generationDomainFraction[1]) / real_t(2),
simulationDomain.zMax() * (real_t(1) + generationDomainFraction[2]) / real_t(2)));
real_t particleOffset = particleGenerationSpacing / real_t(2);
- for (auto pt :
- grid_generator::SCGrid(generationDomain, generationDomain.center() + generationPointOfReferenceOffset,
- particleGenerationSpacing))
- {
- // Offset every second particle layer in flow direction to avoid channels in flow direction
- if (useParticleOffset &&
- uint_t(round(math::abs(generationDomain.center()[0] - pt[0]) / (particleGenerationSpacing))) % uint_t(2) !=
- uint_t(0))
- {
- pt = pt + Vector3(real_t(0), particleOffset, particleOffset);
- }
+ /*for (auto pt :
+ grid_generator::SCGrid(generationDomain, generationDomain.center(),
+ particleGenerationSpacing))*/
+ //{
+ auto pt = Vector3< real_t >(simulationDomain.xMax()/2,simulationDomain.yMax()/2,simulationDomain.zMax()/2);
if (rpdDomain->isContainedInProcessSubdomain(uint_c(mpi::MPIManager::instance()->rank()), pt))
{
mesa_pd::data::Particle&& p = *ps->create();
@@ -275,9 +271,14 @@ int main(int argc, char** argv)
p.setOwner(mpi::MPIManager::instance()->rank());
p.setShapeID(sphereShape);
p.setType(0);
- p.setTemperature(0.4);
+ p.setTemperature(math::realRandom(real_t(0.2),real_t(1)));
+
+ //p.setTemperature(0.308382);
+ //p.setLinearVelocity(Vector3<real_t>(0.1,0.1,0.1));
+
+
}
- }
+ //}
}
////////////////////////
@@ -367,8 +368,8 @@ int main(int argc, char** argv)
boundariesBlockString += "\n BoundariesConcentration";
boundariesBlockString += "{"
- "Border { direction W; walldistance -1; flag Neumann_Concentration; }"
- "Border { direction E; walldistance -1; flag Neumann_Concentration; }";
+ "Border { direction W; walldistance -1; flag Density_Concentration_east; }"
+ "Border { direction E; walldistance -1; flag Density_Concentration_east; }";
if (!periodicInY)
{
@@ -378,8 +379,8 @@ int main(int argc, char** argv)
if (!periodicInZ)
{
- boundariesBlockString += "Border { direction T; walldistance -1; flag Density_Concentration_west; }"
- "Border { direction B; walldistance -1; flag Density_Concentration_east; }";
+ boundariesBlockString += "Border { direction T; walldistance -1; flag Neumann_Concentration; }"
+ "Border { direction B; walldistance -1; flag Neumann_Concentration; }";
}
boundariesBlockString += "}";
WALBERLA_ROOT_SECTION()
@@ -451,9 +452,19 @@ int main(int argc, char** argv)
particleAndVolumeFractionSoA.BsFieldID, particleAndVolumeFractionSoA.BFieldID,densityConcentrationFieldCPUGPUID,
particleAndVolumeFractionSoA.particleVelocitiesFieldID, pdfFieldFluidCPUGPUID,velFieldFluidCPUGPUID,T0,alphaLB,gravityLB,real_t(1.0),rho_0);*/
- pystencils::InitializeFluidDomain pdfSetterFluid(particleAndVolumeFractionSoA.BsFieldID,particleAndVolumeFractionSoA.BFieldID,densityConcentrationFieldCPUGPUID,particleAndVolumeFractionSoA.particleVelocitiesFieldID,pdfFieldFluidCPUGPUID,T0,alphaLB,gravityLB,real_t(1),real_t(0),rho_0);
- pystencils::InitializeConcentrationDomain pdfSetterConcentration(
- densityConcentrationFieldCPUGPUID, pdfFieldConcentrationCPUGPUID, velFieldFluidCPUGPUID);
+ pystencils::InitializeFluidDomain pdfSetterFluid(particleAndVolumeFractionSoA.BsFieldID,particleAndVolumeFractionSoA.BFieldID,densityConcentrationFieldCPUGPUID,particleAndVolumeFractionSoA.particleVelocitiesFieldID,pdfFieldFluidCPUGPUID,T0,alphaLB,gravityLB,real_t(1),real_t(0),real_t(0),real_t(0),rho_0);
+ //pystencils::InitializeFluidDomain pdfSetterFluid(particleAndVolumeFractionSoA.BFieldID,densityConcentrationFieldCPUGPUID,pdfFieldFluidCPUGPUID,T0,alphaLB,gravityLB,real_t(1),real_t(0),real_t(0),real_t(0),rho_0);
+ //pystencils::InitializeFluidDomain pdfSetterFluid(particleAndVolumeFractionSoA.BFieldID,densityConcentrationFieldCPUGPUID,pdfFieldFluidCPUGPUID,T0,alphaLB,gravityLB,real_t(1),real_t(0),rho_0);
+
+ WALBERLA_LOG_INFO_ON_ROOT("reached here fluid initi done");
+ pystencils::InitializeConcentrationDomain pdfSetterConcentration(particleAndVolumeFractionSoA.BsFieldID,particleAndVolumeFractionSoA.BFieldID,
+ densityConcentrationFieldCPUGPUID,particleAndVolumeFractionSoA.particleTemperaturesFieldID, pdfFieldConcentrationCPUGPUID, velFieldFluidCPUGPUID);
+
+ /*pystencils::InitializeConcentrationDomain pdfSetterConcentration(particleAndVolumeFractionSoA.BsFieldID,
+ particleAndVolumeFractionSoA.particleTemperaturesFieldID, pdfFieldConcentrationCPUGPUID, velFieldFluidCPUGPUID);*/
+ //pystencils::InitializeConcentrationDomain pdfSetterConcentration(pdfFieldConcentrationCPUGPUID, velFieldFluidCPUGPUID);
+
+ WALBERLA_LOG_INFO_ON_ROOT("reached here concentr initi done");
#endif
@@ -471,11 +482,14 @@ int main(int argc, char** argv)
//initializeFluidDomain(&(*blockIt));
if (useParticles) {
psmSweepCollection.setParticleVelocitiesSweep(&(*blockIt));
+ psmSweepCollection.setParticleTemperaturesSweep(&(*blockIt));
+ WALBERLA_LOG_INFO_ON_ROOT("set the particle temps");
}
pdfSetterFluid(&(*blockIt));
- WALBERLA_LOG_INFO_ON_ROOT("reached till here 4");
+ WALBERLA_LOG_INFO_ON_ROOT("set the fluid pdfs");
pdfSetterConcentration(&(*blockIt));
+ WALBERLA_LOG_INFO_ON_ROOT("set the conce pdfs");
}
@@ -517,8 +531,9 @@ int main(int argc, char** argv)
pdfFieldFluidID, velFieldFluidID, T0, alphaLB, gravityLB,
rho_0);
#else
- pystencils::FluidMacroGetter getterSweep_fluid(particleAndVolumeFractionSoA.BFieldID,densityConcentrationFieldID, densityFluidFieldID,
- pdfFieldFluidCPUGPUID, velFieldFluidID, T0, alphaLB, gravityLB,
+ WALBERLA_LOG_INFO_ON_ROOT("reachd till hereeeeeee in ceod");
+ pystencils::FluidMacroGetter getterSweep_fluid(particleAndVolumeFractionSoA.BFieldID,densityConcentrationFieldCPUGPUID, densityFluidFieldID,
+ pdfFieldFluidCPUGPUID, velFieldFluidCPUGPUID, T0, alphaLB, gravityLB,
rho_0);
#endif
@@ -526,7 +541,7 @@ int main(int argc, char** argv)
#ifdef WALBERLA_BUILD_WITH_GPU_SUPPORT
pystencils::ConcentrationMacroGetter getterSweep_concentration(densityConcentrationFieldID, pdfFieldConcentrationID);
#else
- pystencils::ConcentrationMacroGetter getterSweep_concentration(densityConcentrationFieldID,
+ pystencils::ConcentrationMacroGetter getterSweep_concentration(densityConcentrationFieldCPUGPUID,
pdfFieldConcentrationCPUGPUID);
#endif
@@ -587,24 +602,29 @@ int main(int argc, char** argv)
make_shared< field::VTKWriter< VelocityField_fluid_T > >(velFieldFluidID, "Fluid Velocity"));
#else
vtkOutput_Fluid->addCellDataWriter(
- make_shared< field::VTKWriter< VelocityField_fluid_T > >(velFieldFluidID, "Fluid Velocity"));
+ make_shared< field::VTKWriter< VelocityField_fluid_T > >(velFieldFluidCPUGPUID, "Fluid Velocity"));
#endif
vtkOutput_Fluid->addCellDataWriter(
make_shared< field::VTKWriter< DensityField_fluid_T > >(densityFluidFieldID, "Fluid Density"));
vtkOutput_Fluid->addCellDataWriter(
make_shared< field::VTKWriter< FlagField_T > >(flagFieldFluidID, "FluidFlagField"));
+ vtkOutput_Fluid->addCellDataWriter(
+ make_shared< field::VTKWriter< BField_T > >(particleAndVolumeFractionSoA.BFieldID, "OverlapFraction"));
+ vtkOutput_Fluid->addCellDataWriter(
+ make_shared< field::VTKWriter< particleTemperaturesField_T > >(particleAndVolumeFractionSoA.particleTemperaturesFieldID, "particle temp filed"));
#ifdef WALBERLA_BUILD_WITH_GPU_SUPPORT
vtkOutput_Concentration->addCellDataWriter(
make_shared< field::VTKWriter< DensityField_concentration_T > >(densityConcentrationFieldID, "Concentration"));
#else
vtkOutput_Concentration->addCellDataWriter(
- make_shared< field::VTKWriter< DensityField_concentration_T > >(densityConcentrationFieldID, "Concentration"));
+ make_shared< field::VTKWriter< DensityField_concentration_T > >(densityConcentrationFieldCPUGPUID, "Concentration"));
#endif
vtkOutput_Concentration->addCellDataWriter(
make_shared< field::VTKWriter< FlagField_T > >(flagFieldConcentrationID, "ConcentrationFlagField"));
+ //WALBERLA_LOG_INFO("reache here in vtk ");
timeloop.addFuncBeforeTimeStep(vtk::writeFiles(vtkOutput_Fluid), "VTK output Fluid");
timeloop.addFuncBeforeTimeStep(vtk::writeFiles(vtkOutput_Concentration), "VTK output Concentration");
}
@@ -624,14 +644,23 @@ int main(int argc, char** argv)
particleAndVolumeFractionSoA.particleVelocitiesFieldID,
pdfFieldFluidCPUGPUID,velFieldFluidCPUGPUID,T0,alphaLB,gravityLB,omega_f,rho_0,frameWidth);
- pystencils::LBMFluidSweep lbmFluidSweep(pdfFieldFluidCPUGPUID,velFieldFluidCPUGPUID,omega_f);
+ pystencils::PSMConcentrationSweep psmConcentrationSweep(
+ particleAndVolumeFractionSoA.BsFieldID, particleAndVolumeFractionSoA.BFieldID, densityConcentrationFieldCPUGPUID,
+ particleAndVolumeFractionSoA.particleTemperaturesFieldID,particleAndVolumeFractionSoA.particleVelocitiesFieldID,pdfFieldConcentrationCPUGPUID, velFieldFluidCPUGPUID,omega_c);
+
+ /*pystencils::PSMConcentrationSweep psmConcentrationSweep(
+ particleAndVolumeFractionSoA.BsFieldID, particleAndVolumeFractionSoA.BFieldID, densityConcentrationFieldCPUGPUID,
+ particleAndVolumeFractionSoA.particleVelocitiesFieldID,pdfFieldConcentrationCPUGPUID, velFieldFluidCPUGPUID,omega_c);
+*/
- pystencils::LBMFluidSplitSweep lbmFluidSplitSweep(pdfFieldFluidCPUGPUID,velFieldFluidCPUGPUID,omega_f,frameWidth);
+ pystencils::LBMFluidSweep lbmFluidSweep(densityConcentrationFieldCPUGPUID,pdfFieldFluidCPUGPUID,velFieldFluidCPUGPUID,T0,alphaLB,gravityLB,omega_f,rho_0);
+
+ pystencils::LBMFluidSplitSweep lbmFluidSplitSweep(densityConcentrationFieldCPUGPUID,pdfFieldFluidCPUGPUID,velFieldFluidCPUGPUID,T0,alphaLB,gravityLB,omega_f,rho_0,frameWidth);
pystencils::LBMConcentrationSweep lbmConcentrationSweep(densityConcentrationFieldCPUGPUID, pdfFieldConcentrationCPUGPUID,
velFieldFluidCPUGPUID,qe,qk);
- pystencils::LBMConcentrationSplitSweep lbmConcentrationSplitSweep(densityConcentrationFieldCPUGPUID,pdfFieldConcentrationCPUGPUID, velFieldFluidCPUGPUID,qe,qk, frameWidth);
+
// Add LBM (fluid and concentration) communication function and boundary handling sweep
@@ -660,56 +689,37 @@ int main(int argc, char** argv)
if (useParticles)
{
- if(useCommunicationHiding){
- commTimeloop.add() << BeforeFunction([&]() { com_fluid.startCommunication(); })
- << Sweep(deviceSyncWrapper(psmSweepCollection.particleMappingSweep), "Particle mapping");
- commTimeloop.add() << Sweep(deviceSyncWrapper(psmSweepCollection.setParticleVelocitiesSweep),
- "Set particle velocities");
- commTimeloop.add() << Sweep(deviceSyncWrapper(psmFluidSplitSweep.getInnerSweep()), "PSM inner sweep")
- << AfterFunction([&]() { com_fluid.wait(); }, "LBM Communication (wait)");
- timeloop.add() << Sweep(deviceSyncWrapper(psmFluidSplitSweep.getOuterSweep()), "PSM outer sweep");
-
- commTimeloop.add() << BeforeFunction([&]() { com_concentration.startCommunication(); }, "LBM concentration Communication (start)")
- << Sweep(deviceSyncWrapper(lbmConcentrationSplitSweep.getInnerSweep()), "LBM concentration inner sweep")
- << AfterFunction([&]() { com_concentration.wait(); }, "LBM concentration Communication (wait)");
- timeloop.add() << Sweep(deviceSyncWrapper(lbmConcentrationSplitSweep.getOuterSweep()), "LBM concentration outer sweep");
- // after both the sweeps, reduce the particle forces.
- timeloop.add() << Sweep(deviceSyncWrapper(psmSweepCollection.reduceParticleForcesSweep),
- "Reduce particle forces");
- }
- else
- {
timeloop.add() << Sweep(deviceSyncWrapper(psmSweepCollection.particleMappingSweep), "Particle mapping");
timeloop.add() << Sweep(deviceSyncWrapper(psmSweepCollection.setParticleVelocitiesSweep),
"Set particle velocities");
timeloop.add() << Sweep(deviceSyncWrapper(psmFluidSweep), "PSM Fluid sweep");
- timeloop.add() << Sweep(deviceSyncWrapper(lbmConcentrationSweep), "LBM Concentration sweep");
+ //timeloop.add() << Sweep(deviceSyncWrapper(lbmConcentrationSweep), "LBM Concentration sweep");
+ timeloop.add() << Sweep(deviceSyncWrapper(psmSweepCollection.setParticleTemperaturesSweep),
+ "Set particle temperatures");
+ timeloop.add() << Sweep(deviceSyncWrapper(psmConcentrationSweep), "PSM Concentration sweep");
// after both the sweeps, reduce the particle forces.
timeloop.add() << Sweep(deviceSyncWrapper(psmSweepCollection.reduceParticleForcesSweep),
"Reduce particle forces");
- }
+
}
else{
- if (useCommunicationHiding)
- {
- commTimeloop.add() << BeforeFunction([&]() { com_fluid.startCommunication(); })
- << Sweep(deviceSyncWrapper(lbmFluidSplitSweep.getInnerSweep()), "PSM inner sweep")
- << AfterFunction([&]() { com_fluid.wait(); }, "LBM Communication (wait)");
- timeloop.add() << Sweep(deviceSyncWrapper(lbmFluidSplitSweep.getOuterSweep()), "PSM outer sweep");
-
- commTimeloop.add() << BeforeFunction([&]() { com_concentration.startCommunication(); }, "LBM concentration Communication (start)")
- << Sweep(deviceSyncWrapper(lbmConcentrationSplitSweep.getInnerSweep()), "LBM concentration inner sweep")
- << AfterFunction([&]() { com_concentration.wait(); }, "LBM concentration Communication (wait)");
- timeloop.add() << Sweep(deviceSyncWrapper(lbmConcentrationSplitSweep.getOuterSweep()), "LBM concentration outer sweep");
- }
- else
- {
+ /*WALBERLA_LOG_INFO_ON_ROOT("No particles used so executing normal fluid and concentration sweeps");
timeloop.add() << Sweep(deviceSyncWrapper(lbmFluidSweep), "LBM Fluid sweep");
- timeloop.add() << Sweep(deviceSyncWrapper(lbmConcentrationSweep), "LBM Concentration sweep");
- }
+ timeloop.add() << Sweep(deviceSyncWrapper(lbmConcentrationSweep), "LBM Concentration sweep");*/
+ timeloop.add() << Sweep(deviceSyncWrapper(psmSweepCollection.particleMappingSweep), "Particle mapping");
+ timeloop.add() << Sweep(deviceSyncWrapper(psmSweepCollection.setParticleVelocitiesSweep),
+ "Set particle velocities");
+ timeloop.add() << Sweep(deviceSyncWrapper(psmFluidSweep), "PSM Fluid sweep");
+ //timeloop.add() << Sweep(deviceSyncWrapper(lbmConcentrationSweep), "LBM Concentration sweep");
+ timeloop.add() << Sweep(deviceSyncWrapper(psmConcentrationSweep), "PSM Concentration sweep");
+
+ // after both the sweeps, reduce the particle forces.
+ timeloop.add() << Sweep(deviceSyncWrapper(psmSweepCollection.reduceParticleForcesSweep),
+ "Reduce particle forces");
+
}
diff --git a/src/lbm_mesapd_coupling/DataTypesCodegen.h b/src/lbm_mesapd_coupling/DataTypesCodegen.h
index a49c9ea1ade8275f22935b9725b9e1e820304c3c..6d5edb31bb67e11781566f5aba808ef54697ff8c 100644
--- a/src/lbm_mesapd_coupling/DataTypesCodegen.h
+++ b/src/lbm_mesapd_coupling/DataTypesCodegen.h
@@ -56,6 +56,7 @@ using idxField_T = GhostLayerField< size_t, MaxParticlesPerCel
using BField_T = GhostLayerField< real_t, 1 >;
using particleVelocitiesField_T = GhostLayerField< real_t, MaxParticlesPerCell * 3 >;
using particleForcesField_T = GhostLayerField< real_t, MaxParticlesPerCell * 3 >;
+using particleTemperaturesField_T = GhostLayerField< real_t, MaxParticlesPerCell*1 >;
#ifdef WALBERLA_BUILD_WITH_GPU_SUPPORT
using nOverlappingParticlesFieldGPU_T = walberla::gpu::GPUField< uint_t >;
using BsFieldGPU_T = walberla::gpu::GPUField< real_t >;
@@ -75,6 +76,7 @@ struct ParticleAndVolumeFractionSoA_T
BlockDataID BFieldID;
BlockDataID particleVelocitiesFieldID;
BlockDataID particleForcesFieldID;
+ BlockDataID particleTemperaturesFieldID;
// relaxation rate omega is used for Weighting_T != 1
real_t omega_;
// UIDs of the particles are stored during mapping, and it is checked that they are the same during the PSM kernel.
@@ -109,6 +111,8 @@ struct ParticleAndVolumeFractionSoA_T
bs, "particle velocities field CPU", real_t(0), field::fzyx, uint_t(1), true);
particleForcesFieldID = field::addToStorage< particleForcesField_T >(bs, "particle forces field CPU", real_t(0),
field::fzyx, uint_t(1), true);
+ particleTemperaturesFieldID = field::addToStorage< particleTemperaturesField_T >(bs, "particle temperatures field CPU", real_t(0),
+ field::fzyx, uint_t(1), true);
#endif
omega_ = omega;
}
diff --git a/src/lbm_mesapd_coupling/partially_saturated_cells_method/codegen/PSMSweepCollection.h b/src/lbm_mesapd_coupling/partially_saturated_cells_method/codegen/PSMSweepCollection.h
index f4d8400fae806d82f28f0332646e2d7aa5069c82..07feccf1c693f4fb98eb7d6455d0a0f114220440 100644
--- a/src/lbm_mesapd_coupling/partially_saturated_cells_method/codegen/PSMSweepCollection.h
+++ b/src/lbm_mesapd_coupling/partially_saturated_cells_method/codegen/PSMSweepCollection.h
@@ -61,11 +61,15 @@ class PSMSweepCollection
setParticleVelocitiesSweep(SetParticleVelocitiesSweep< ParticleAccessor_T, ParticleSelector_T, Weighting_T >(
bs, ac, ps, particleAndVolumeFractionSoA)),
reduceParticleForcesSweep(ReduceParticleForcesSweep< ParticleAccessor_T, ParticleSelector_T, Weighting_T >(
+ bs, ac, ps, particleAndVolumeFractionSoA)),
+ setParticleTemperaturesSweep(SetParticleTemperaturesSweep< ParticleAccessor_T, ParticleSelector_T, Weighting_T >(
bs, ac, ps, particleAndVolumeFractionSoA))
+
{}
SphereFractionMappingSweep< ParticleAccessor_T, ParticleSelector_T, Weighting_T > particleMappingSweep;
SetParticleVelocitiesSweep< ParticleAccessor_T, ParticleSelector_T, Weighting_T > setParticleVelocitiesSweep;
ReduceParticleForcesSweep< ParticleAccessor_T, ParticleSelector_T, Weighting_T > reduceParticleForcesSweep;
+ SetParticleTemperaturesSweep< ParticleAccessor_T, ParticleSelector_T, Weighting_T >setParticleTemperaturesSweep;
};
template< typename SweepCollection, typename PSMSweep >
diff --git a/src/lbm_mesapd_coupling/partially_saturated_cells_method/codegen/PSMWrapperSweepsCPU.h b/src/lbm_mesapd_coupling/partially_saturated_cells_method/codegen/PSMWrapperSweepsCPU.h
index 3eb79994f76de9183005e2016eef54fd42da196f..56b26f6ae993f9ca72e8bcd4cee4c77fdc13cc9d 100644
--- a/src/lbm_mesapd_coupling/partially_saturated_cells_method/codegen/PSMWrapperSweepsCPU.h
+++ b/src/lbm_mesapd_coupling/partially_saturated_cells_method/codegen/PSMWrapperSweepsCPU.h
@@ -46,6 +46,73 @@ namespace psm
namespace gpu
{
+template< typename ParticleAccessor_T, typename ParticleSelector_T, int Weighting_T >
+class SetParticleTemperaturesSweep
+{
+ public:
+ SetParticleTemperaturesSweep(const shared_ptr< StructuredBlockStorage >& bs,
+ const shared_ptr< ParticleAccessor_T >& ac,
+ const ParticleSelector_T& mappingParticleSelector,
+ ParticleAndVolumeFractionSoA_T< Weighting_T >& particleAndVolumeFractionSoA)
+ : bs_(bs), ac_(ac), mappingParticleSelector_(mappingParticleSelector),
+ particleAndVolumeFractionSoA_(particleAndVolumeFractionSoA)
+ {}
+ void operator()(IBlock* block)
+ {
+ // Check that uids of the particles have not changed since the last mapping to avoid incorrect indices
+ std::vector< walberla::id_t > currentUIDs;
+ for (size_t idx = 0; idx < ac_->size(); ++idx)
+ {
+ if (mappingParticleSelector_(idx, *ac_)) { currentUIDs.push_back(ac_->getUid(idx)); }
+ }
+ WALBERLA_ASSERT(particleAndVolumeFractionSoA_.mappingUIDs == currentUIDs);
+
+ size_t numMappedParticles = 0;
+ for (size_t idx = 0; idx < ac_->size(); ++idx)
+ {
+ if (mappingParticleSelector_(idx, *ac_)) { numMappedParticles++; }
+ }
+
+ if (numMappedParticles == uint_t(0)) return;
+
+ size_t arraySizes = numMappedParticles * sizeof(real_t) * 1;
+
+ // Allocate unified memory for the particle information required for computing the temperature (used in
+ // the solid collision operator for temperature field)
+ real_t* temperatures = (real_t*) malloc(arraySizes);
+ memset(temperatures, 0, arraySizes);
+
+ // Store particle information inside memory
+ size_t idxMapped = 0;
+ for (size_t idx = 0; idx < ac_->size(); ++idx)
+ {
+ if (mappingParticleSelector_(idx, *ac_))
+ {
+ temperatures[idxMapped] = ac_->getTemperature(idx);
+ WALBERLA_LOG_INFO("temperature from accessor is " << ac_->getTemperature(idx));
+ idxMapped++;
+ }
+ }
+
+ auto nOverlappingParticlesField =
+ block->getData< nOverlappingParticlesField_T >(particleAndVolumeFractionSoA_.nOverlappingParticlesFieldID);
+ auto idxField = block->getData< idxField_T >(particleAndVolumeFractionSoA_.idxFieldID);
+ auto particleTemperaturesField =
+ block->getData< particleTemperaturesField_T >(particleAndVolumeFractionSoA_.particleTemperaturesFieldID);
+ WALBERLA_FOR_ALL_CELLS_XYZ(
+ particleTemperaturesField, for (uint_t p = 0; p < nOverlappingParticlesField->get(x, y, z); p++) {
+ particleTemperaturesField->get(x, y, z, p) = temperatures[idxField->get(x, y, z, p)];
+ })
+ free(temperatures);
+ }
+
+ private:
+ shared_ptr< StructuredBlockStorage > bs_;
+ const shared_ptr< ParticleAccessor_T > ac_;
+ const ParticleSelector_T& mappingParticleSelector_;
+ ParticleAndVolumeFractionSoA_T< Weighting_T >& particleAndVolumeFractionSoA_;
+};
+
template< typename ParticleAccessor_T, typename ParticleSelector_T, int Weighting_T >
class SetParticleVelocitiesSweep
{