diff --git a/DEMSystems/sphereDEMSystem/thermalSphereDEMSystem.cpp b/DEMSystems/sphereDEMSystem/thermalSphereDEMSystem.cpp new file mode 100644 index 000000000..b0dae6355 --- /dev/null +++ b/DEMSystems/sphereDEMSystem/thermalSphereDEMSystem.cpp @@ -0,0 +1,252 @@ +/*------------------------------- phasicFlow --------------------------------- + O C enter of + O O E ngineering and + O O M ultiscale modeling of + OOOOOOO F luid flow +------------------------------------------------------------------------------ + Copyright (C): www.cemf.ir + email: hamid.r.norouzi AT gmail.com +------------------------------------------------------------------------------ +Licence: + This file is part of phasicFlow code. It is a free software for simulating + granular and multiphase flows. You can redistribute it and/or modify it under + the terms of GNU General Public License v3 or any other later versions. + + phasicFlow is distributed to help others in their research in the field of + granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the + implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + +-----------------------------------------------------------------------------*/ + +#include "thermalSphereDEMSystem.hpp" +#include "vocabs.hpp" + +namespace pFlow +{ + +//----------------------------- protected methods ----------------------------- + +// ========================================================================= // +// Section 3: Core Physics Loop +// ========================================================================= // + +bool thermalSphereDEMSystem::loop() +{ + do + { + // 3.1 Handle particle injection triggers + if (!insertion_().insertParticles( + Control().time().currentIter(), + Control().time().currentTime(), + Control().time().dt())) + { + fatalError << "Particle insertion failed " + << "in thermalSphereDEMSystem.\n"; + return false; + } + + // 3.2 Initialize physics accumulators + geometry_->beforeIteration(); + interaction_->beforeIteration(); + particles_->beforeIteration(); + + // 3.3 Mechanical collision evaluation + interaction_->iterate(); + + // 3.4 Thermodynamic evaluation (Q_pp, Q_pfp, Q_rad) + if (thermalInteraction_) + { + thermalInteraction_->iterate(); + } + + // 3.5 Equations of motion and explicit energy integration + particles_->iterate(); + + // 3.6 Clean up and state finalization + geometry_->iterate(); + particles_->afterIteration(); + geometry_->afterIteration(); + + } while(Control()++); + + return true; +} + +//----------------------------- constructors ---------------------------------- + +// ========================================================================= // +// Section 1: Constructors +// ========================================================================= // + +thermalSphereDEMSystem::thermalSphereDEMSystem( + word demSystemName, + const std::vector& domains, + int argc, + char* argv[], + bool requireRVel) +: + sphereDEMSystem(demSystemName, domains, argc, argv, requireRVel) +{ + REPORT(0) << "\nInitializing thermal DEM components..." << END_REPORT; + + // Reset base instances for thermal override + interaction_.reset(); + insertion_.reset(); + particles_.reset(); + spheres_.reset(); + + // 1.1 Load thermal properties + auto thermalProps = thermalProperty( + propertyFile__, + Control().caseSetup().path()); + + // 1.2 Initialize thermal shapes + auto* combinedShape = new thermalSphereShape( + shapeFile__, + &Control().caseSetup(), + thermalProps); + + spheres_ = uniquePtr(combinedShape); + + // 1.3 Initialize thermal particles on GPU + auto* tp = new thermalSphereParticles( + Control(), + *combinedShape, + *combinedShape); + + particles_ = uniquePtr(tp); + thermalParticles_ = tp; + + // 1.4 Reconstruct insertion mechanism + insertion_ = makeUnique( + particles_(), + particles_().spheres()); + + if (!thermalParticles_->initializeThermalParticles()) + { + fatalError << "Failed to initialize thermal properties " + << "for particles.\n"; + } + + // 1.5 Reconstruct mechanical interactions + interaction_ = interaction::create( + Control(), + Particles(), + Geometry()); + + // 1.6 Initialize Unified Thermal Interaction Model + REPORT(0) << "Creating thermal interactions " + << "(Conduction, Radiation, PFP)..." << END_REPORT; + + box localDomain = domains.empty() ? box() : domains[0]; + + thermalInteraction_ = makeUnique( + Control(), + *thermalParticles_, + localDomain); + + // 1.7 Update distribution boundaries + real minD, maxD; + particles_->boundingSphereMinMax(minD, maxD); + particleDistribution_ = makeUnique(domains, maxD); +} + +//---------------------------- public methods --------------------------------- + +// ========================================================================= // +// Section 2: Time Integration Constraints +// ========================================================================= // + +bool thermalSphereDEMSystem::iterate( + real upToTime, + real timeToWrite, + word timeName) +{ + Control().time().setStopAt(upToTime); + Control().time().setOutputToFile(timeToWrite, timeName); + + return loop(); +} + +bool thermalSphereDEMSystem::iterate(real upToTime) +{ + Control().time().setStopAt(upToTime); + + return loop(); +} + +// ========================================================================= // +// Section 4: Data Exchange Interfaces (CFD-DEM Coupling) +// ========================================================================= // + +span thermalSphereDEMSystem::temperature() +{ + auto& hVec = thermalParticles_->temperatureHost(); + return span(hVec.data(), hVec.size()); +} + +span thermalSphereDEMSystem::emissivity() +{ + auto& hVec = thermalParticles_->emissivityHost(); + return span(hVec.data(), hVec.size()); +} + +span thermalSphereDEMSystem::radSumTemp() +{ + auto& hVec = thermalParticles_->radSumTempHost(); + return span(hVec.data(), hVec.size()); +} + +span thermalSphereDEMSystem::radNumPrt() +{ + auto& hVec = thermalParticles_->radNumPrtHost(); + return span(hVec.data(), hVec.size()); +} + +span thermalSphereDEMSystem::parFluidHeatSourceConv() +{ + auto& hVec = thermalParticles_->heatSourceConvHost(); + return span(hVec.data(), hVec.size()); +} + +span thermalSphereDEMSystem::parFluidHeatSourceRad() +{ + auto& hVec = thermalParticles_->heatSourceRadHost(); + return span(hVec.data(), hVec.size()); +} + +bool thermalSphereDEMSystem::sendFluidHeatSourcesToDEM() +{ + thermalParticles_->heatSourcesHostUpdatedSync(); + return true; +} + +// ========================================================================= // +// Section 5: PFP Pipeline Exchange +// ========================================================================= // + +span thermalSphereDEMSystem::parFluidKappa() +{ + auto& hVec = thermalParticles_->fluidKappaHost(); + return span(hVec.data(), hVec.size()); +} + +span thermalSphereDEMSystem::parFluidAlpha() +{ + auto& hVec = thermalParticles_->fluidAlphaHost(); + return span(hVec.data(), hVec.size()); +} + +bool thermalSphereDEMSystem::sendFluidPropertiesToDEM() +{ + thermalParticles_->fluidPropertiesHostUpdatedSync(); + return true; +} + +//+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +} // pFlow + + + + diff --git a/DEMSystems/sphereDEMSystem/thermalSphereDEMSystem.hpp b/DEMSystems/sphereDEMSystem/thermalSphereDEMSystem.hpp new file mode 100644 index 000000000..6b9d31005 --- /dev/null +++ b/DEMSystems/sphereDEMSystem/thermalSphereDEMSystem.hpp @@ -0,0 +1,129 @@ +/*------------------------------- phasicFlow --------------------------------- + O C enter of + O O E ngineering and + O O M ultiscale modeling of + OOOOOOO F luid flow +------------------------------------------------------------------------------ + Copyright (C): www.cemf.ir + email: hamid.r.norouzi AT gmail.com +------------------------------------------------------------------------------ +Licence: + This file is part of phasicFlow code. It is a free software for simulating + granular and multiphase flows. You can redistribute it and/or modify it under + the terms of GNU General Public License v3 or any other later versions. + + phasicFlow is distributed to help others in their research in the field of + granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the + implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + +-----------------------------------------------------------------------------*/ + +#ifndef pFlow_thermalSphereDEMSystem_hpp +#define pFlow_thermalSphereDEMSystem_hpp + +#include "sphereDEMSystem.hpp" +#include "thermalSphereParticles.hpp" +#include "thermalInteraction.hpp" +#include "thermalProperty.hpp" +#include "thermalSphereShape.hpp" + +namespace pFlow +{ + +/** + * @brief Extends the base mechanical DEM solver to handle thermodynamic + * physics (Conduction, PFP, and Radiation). + */ +class thermalSphereDEMSystem +: + public sphereDEMSystem +{ +protected: + + //- protected members + + /// @brief Direct host-view access for thermal particles. + thermalSphereParticles* thermalParticles_ = nullptr; + + /// @brief Manages inter-particle conduction, fluid-bridge heat, radiation. + uniquePtr thermalInteraction_ = nullptr; + + //- protected methods + + /** + * @brief Core integration loop encompassing mechanical & thermal updates. + * @return True upon successful execution. + */ + bool loop(); + +public: + + //- Type info + + TypeInfo("thermalSphereDEMSystem"); + + //- constructors + + thermalSphereDEMSystem( + word demSystemName, + const std::vector& domains, + int argc, + char* argv[], + bool requireRVel = false); + + ~thermalSphereDEMSystem() override = default; + + //- public methods + + add_vCtor( + DEMSystem, + thermalSphereDEMSystem, + word + ); + + bool iterate( + real upToTime, + real timeToWrite, + word timeName) override; + + bool iterate(real upToTime) override; + + span temperature() override; + + span emissivity() override; + + span radSumTemp() override; + + span radNumPrt() override; + + span parFluidHeatSourceConv() override; + + span parFluidHeatSourceRad() override; + + bool sendFluidHeatSourcesToDEM() override; + + span parFluidKappa() override; + + span parFluidAlpha() override; + + bool sendFluidPropertiesToDEM() override; + + /** + * @brief Evaluates radiation module availability. + * @return True if the module is active and enabled by user dictionary. + */ + inline + bool hasRadiation() const override + { + return thermalInteraction_ != nullptr && + thermalInteraction_->isRadiationEnabled(); + } + +}; // thermalSphereDEMSystem + +} // pFlow + +#endif // pFlow_thermalSphereDEMSystem_hpp + + + diff --git a/solvers/heatSphereGranFlow/CMakeLists.txt b/solvers/heatSphereGranFlow/CMakeLists.txt new file mode 100644 index 000000000..43e85c51c --- /dev/null +++ b/solvers/heatSphereGranFlow/CMakeLists.txt @@ -0,0 +1,7 @@ + +set(source_files +heatSphereGranFlow.cpp +) +set(link_lib Kokkos::kokkos phasicFlow Particles Geometry Property Interaction Interaction Utilities) + +pFlow_make_executable_install(heatSphereGranFlow source_files link_lib) diff --git a/solvers/heatSphereGranFlow/createDEMComponents.hpp b/solvers/heatSphereGranFlow/createDEMComponents.hpp new file mode 100644 index 000000000..647d28d04 --- /dev/null +++ b/solvers/heatSphereGranFlow/createDEMComponents.hpp @@ -0,0 +1,121 @@ +/*------------------------------- phasicFlow --------------------------------- + O C enter of + O O E ngineering and + O O M ultiscale modeling of + OOOOOOO F luid flow +------------------------------------------------------------------------------ + Copyright (C): www.cemf.ir + email: hamid.r.norouzi AT gmail.com +------------------------------------------------------------------------------ +Licence: + This file is part of phasicFlow code. It is a free software for simulating + granular and multiphase flows. You can redistribute it and/or modify it under + the terms of GNU General Public License v3 or any other later versions. + + phasicFlow is distributed to help others in their research in the field of + granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the + implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + +-----------------------------------------------------------------------------*/ + +/** + * @file createDEMComponents.hpp + * @brief Initialization sequence for the thermal DEM simulation objects. + * + * @details + * Instantiates the thermal shape, particle container, insertion + * mechanism, and mechanical + thermal interaction models required + * before the main time loop begins. + */ + +// ========================================================================= // +// Section 1: Shape & Material Initialization +// ========================================================================= // + +REPORT(0) << "Reading thermal shapes dictionary..." << END_REPORT; + +/** + * @brief Geometric + thermal shape dictionary. + * Binds per-material thermal data from thermalProperty + * to per-shape geometry properties. + */ +pFlow::thermalSphereShape spheres +( + pFlow::shapeFile__, + &Control.caseSetup(), + proprties // thermalProperty instance from the main solver +); + +// ========================================================================= // +// Section 2: Particle Container Initialization +// ========================================================================= // + +REPORT(0) << "\nReading thermal sphere particles . . ." << END_REPORT; + +/** + * @brief Main GPU-backed thermal particle container. + */ +pFlow::thermalSphereParticles sphParticles +( + Control, + spheres, + spheres +); + +// ========================================================================= // +// Section 3: Particle Insertion Mechanism +// ========================================================================= // + +REPORT(0) << "\nCreating particle insertion object . . ." << END_REPORT; + +/** + * @brief Time-triggered particle injector. + */ +auto sphInsertion = pFlow::sphereInsertion +( + sphParticles, + sphParticles.spheres() +); + +// ========================================================================= // +// Section 4: Mechanical Interaction Model +// ========================================================================= // + +REPORT(0) << "\nCreating interaction model for sphere-sphere contact . . ." + << END_REPORT; + +/** + * @brief Factory-instantiated contact-force model (e.g., Hertz-Mindlin). + * Handles particle-particle and particle-wall collisions. + */ +auto interactionPtr = pFlow::interaction::create +( + Control, + sphParticles, + surfGeometry +); + +auto& sphInteraction = interactionPtr(); + +// ========================================================================= // +// Section 5: Thermal Interaction Model (Fixed for Standalone Mode) +// ========================================================================= // + +REPORT(0) << "\nCreating unified thermal interaction model " + << "(Conduction, PFP, Radiation) . . ." << END_REPORT; + +/** + * @brief Thermal physics dispatcher. + * Computes Q_pp (Batchelor-O'Brien) and particle-particle radiation. + */ +auto thermalIntPtr = pFlow::makeUnique +( + Control, + sphParticles, + pFlow::box() +); + +auto& thermalInt = thermalIntPtr(); + + + diff --git a/solvers/heatSphereGranFlow/heatSphereGranFlow.cpp b/solvers/heatSphereGranFlow/heatSphereGranFlow.cpp new file mode 100644 index 000000000..0696fb706 --- /dev/null +++ b/solvers/heatSphereGranFlow/heatSphereGranFlow.cpp @@ -0,0 +1,163 @@ +/*------------------------------- phasicFlow --------------------------------- + O C enter of + O O E ngineering and + O O M ultiscale modeling of + OOOOOOO F luid flow +------------------------------------------------------------------------------ + Copyright (C): www.cemf.ir + email: hamid.r.norouzi AT gmail.com +------------------------------------------------------------------------------ +Licence: + This file is part of phasicFlow code. It is a free software for simulating + granular and multiphase flows. You can redistribute it and/or modify it under + the terms of GNU General Public License v3 or any other later versions. + + phasicFlow is distributed to help others in their research in the field of + granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the + implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + +-----------------------------------------------------------------------------*/ + +/** + * @file heatSphereGranFlow.cpp + * @brief Standalone DEM solver for granular flow with heat transfer. + * + * @details + * This solver simulates the granular flow of cohesion-less, spherical + * particles while additionally solving the explicit particle energy + * equation (contact conduction, radiation, and particle-fluid-particle + * sub-grid heat transfer) via Kokkos kernels. It carries no chemical + * reaction capability; for reacting flows, use multiSpeciesGranFlow, + * which extends this same thermal layer. + * + * Note: Q_conv (convection) and Q_pfp (fluid bridge) remain 0.0 since + * the Eulerian fluid mesh does not exist in standalone mode. This mode + * is ideal for unit-testing conduction (Q_pp) and radiation. + */ + +#include "vocabs.hpp" +#include "phasicFlowKokkos.hpp" +#include "systemControl.hpp" +#include "commandLine.hpp" +#include "property.hpp" +#include "geometry.hpp" +#include "sphereParticles.hpp" +#include "interaction.hpp" +#include "Insertions.hpp" + +// --- Thermal Additions --- +#include "thermalProperty.hpp" +#include "thermalSphereShape.hpp" +#include "thermalSphereParticles.hpp" +#include "thermalInteraction.hpp" + +/** + * @brief Main execution entry point for the standalone thermal DEM + * solver. + */ +int main(int argc, char* argv[]) +{ + // ===================================================================== // + // Section 1: Initialization & CLI Parsing + // ===================================================================== // + pFlow::commandLine cmds + ( + "heatSphereGranFlow", + "DEM solver for non-cohesive spherical particles with heat " + "transfer, particle insertion mechanism, and moving geometry." + ); + + bool isCoupling = false; + + if (!cmds.parse(argc, argv)) return 0; + + // this should be palced in each main + pFlow::processors::initProcessors(argc, argv); + pFlow::initialize_pFlowProcessors(); + + #include "initialize_Control.hpp" + + // ===================================================================== // + // Section 2: Material & Geometry Setup + // ===================================================================== // + + /// Read global thermal properties from the case directory. + auto proprties = pFlow::thermalProperty + ( + pFlow::propertyFile__, + Control.caseSetup().path() + ); + + #include "setSurfaceGeometry.hpp" + + #include "createDEMComponents.hpp" + + // ===================================================================== // + // Section 3: Solver Capabilities Notice + // ===================================================================== // + REPORT(0) + << "\n[INFO] Standalone Thermal Mode Active.\n" + << " Q_pp (contact conduction) : Computed via Kokkos kernel\n" + << " Q_rad (radiation) : Computed via Kokkos kernel\n" + << " Q_conv / Q_pfp = 0 (No Eulerian fluid mesh " + << "present)\n" + << " No chemical reaction capability in this solver.\n" + << " Use multiSpeciesGranFlow for reacting flows, or\n" + << " unresolvedHeatSpherePFPlus for coupled CFD-DEM heat " + << "transfer.\n" + << END_REPORT; + + // ===================================================================== // + // Section 4: Main Transient Time Loop + // ===================================================================== // + REPORT(0) << "\nStart of time loop . . .\n" << END_REPORT; + + do + { + // 4.1 Particle insertion phase + if (!sphInsertion.insertParticles( + Control.time().currentIter(), + Control.time().currentTime(), + Control.time().dt())) + { + fatalError + << "particle insertion failed in heatSphereGranFlow " + << "solver.\n"; + return 1; + } + + // 4.2 Pre-processing updates (reset forces, predict, etc.) + surfGeometry.beforeIteration(); + sphParticles.beforeIteration(); + sphInteraction.beforeIteration(); + + // 4.3 Evaluate contact interactions (particle-particle, wall) + sphInteraction.iterate(); + + // 4.4 Thermal physics (conduction, radiation, PFP) + thermalInt.iterate(); + + // 4.5 Update boundary kinematics + surfGeometry.iterate(); + + // 4.6 Update particle kinematics and integrate temperature + sphParticles.iterate(); + + // 4.7 Post-processing cleanups + sphInteraction.afterIteration(); + surfGeometry.afterIteration(); + sphParticles.afterIteration(); + + } while (Control++); + + REPORT(0) << "\nEnd of time loop.\n" << END_REPORT; + + // this should be palced in each main + #include "finalize.hpp" + pFlow::processors::finalizeProcessors(); + + return 0; +} + + + diff --git a/src/Interaction/thermalInteraction/thermalInteraction.cpp b/src/Interaction/thermalInteraction/thermalInteraction.cpp new file mode 100644 index 000000000..b965f75cc --- /dev/null +++ b/src/Interaction/thermalInteraction/thermalInteraction.cpp @@ -0,0 +1,298 @@ +/*------------------------------- phasicFlow --------------------------------- + O C enter of + O O E ngineering and + O O M ultiscale modeling of + OOOOOOO F luid flow +------------------------------------------------------------------------------ + Copyright (C): www.cemf.ir + email: hamid.r.norouzi AT gmail.com +------------------------------------------------------------------------------ +Licence: + This file is part of phasicFlow code. It is a free software for simulating + granular and multiphase flows. You can redistribute it and/or modify it under + the terms of GNU General Public License v3 or any other later versions. + + phasicFlow is distributed to help others in their research in the field of + granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the + implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + +-----------------------------------------------------------------------------*/ + +#include "thermalInteraction.hpp" +#include "thermalInteractionKernels.hpp" + +namespace pFlow +{ + +//----------------------------- constructors ---------------------------------- + +// ========================================================================= // +// Constructor +// ========================================================================= // + +thermalInteraction::thermalInteraction( + systemControl& control, + const thermalSphereParticles& particles, + const box& domainBox) +: + control_(control), + particles_(particles), + thermalTimer_("thermalInteraction", &control.timers()) +{ + dictionary thermoDict( + "thermoPhysicalInteraction", + control_.caseSetup().path() + "thermoPhysicalInteraction"); + + // ---------------------------------------------------------------------- // + // Radiation + // ---------------------------------------------------------------------- // + if (!thermoDict.containsDataEntry("enableRadiation")) + { + fatalErrorInFunction + << "Missing MANDATORY entry 'enableRadiation' " + << "in thermoPhysicalInteraction dictionary." << endl; + fatalExit; + } + + Logical enableRad = thermoDict.getVal("enableRadiation"); + + if (enableRad) + { + enableRadiation_ = true; + + // radCut determines which particles are treated as radiating + // neighbours of one another. A silently-defaulted value would + // switch radiation off in effect while still reporting it as + // enabled, so it must be supplied explicitly. + if (!thermoDict.containsDataEntry("radCut")) + { + fatalErrorInFunction + << "Parameter 'radCut' is mandatory when enableRadiation " + << "is true.\nPlease add it to the thermoPhysicalInteraction " + << "dictionary." << endl; + fatalExit; + } + radCut_ = thermoDict.getVal("radCut"); + + radUpdateInterval_ = thermoDict.getValOrSet( + "radUpdateInterval", + 1); + + if (radUpdateInterval_ == 0) + { + fatalErrorInFunction + << "'radUpdateInterval' must be a positive integer, got 0." + << endl; + fatalExit; + } + + REPORT(0) << "Creating Radiation interaction model . . ." << END_REPORT; + } + else + { + enableRadiation_ = false; + REPORT(0) << Yellow_Text(" -> Radiation is disabled by user.") + << END_REPORT; + } + + // ---------------------------------------------------------------------- // + // Collisional Heat Conduction (Q_pp) + // ---------------------------------------------------------------------- // + if (!thermoDict.containsDataEntry("enableConduction")) + { + fatalErrorInFunction + << "Missing MANDATORY entry 'enableConduction' " + << "in thermoPhysicalInteraction dictionary." << endl; + fatalExit; + } + + Logical enableCond = thermoDict.getVal("enableConduction"); + enableConduction_ = enableCond ? true : false; + + if (enableConduction_) + { + REPORT(0) << "Creating Collisional Heat Transfer (Q_p-p) model . . ." + << END_REPORT; + } + else + { + REPORT(0) << Yellow_Text(" -> Collisional Heat Transfer is disabled.") + << END_REPORT; + } + + // ---------------------------------------------------------------------- // + // Particle-Fluid-Particle (PFP) Sub-grid Heat Transfer + // ---------------------------------------------------------------------- // + if (!thermoDict.containsDataEntry("enablePFP")) + { + fatalErrorInFunction + << "Missing MANDATORY entry 'enablePFP' " + << "in thermoPhysicalInteraction dictionary." << endl; + fatalExit; + } + + Logical enablePfpFlag = thermoDict.getVal("enablePFP"); + enablePFP_ = enablePfpFlag ? true : false; + + if (enablePFP_) + { + REPORT(0) << "Creating Particle-Fluid-Particle (PFP) " + << "sub-grid Heat Transfer model . . ." << END_REPORT; + } + else + { + REPORT(0) << Yellow_Text(" -> PFP Heat Transfer is disabled.") + << END_REPORT; + } + + // ---------------------------------------------------------------------- // + // Hertzian simulation-scale Young's modulus. + // + // Used to compute the mechanical contact radius whenever two particles + // touch. That contact radius feeds both the collisional conduction + // rate (Q_pp) and the PFP contact-limit radius r_sij, regardless of + // which of the two mechanisms triggered the calculation, so it must + // be supplied whenever either is enabled. + // ---------------------------------------------------------------------- // + if (enableConduction_ || enablePFP_) + { + if (thermoDict.containsDataEntry("simYoungsModulus")) + { + simYoungsModulus_ = thermoDict.getVal("simYoungsModulus"); + } + else + { + fatalErrorInFunction + << "Parameter 'simYoungsModulus' is mandatory when " + << "enableConduction or enablePFP is true.\n" + << "Please add it to the thermoPhysicalInteraction dictionary." + << endl; + fatalExit; + } + } + + // ---------------------------------------------------------------------- // + // Determine the neighbor search cell size. + // ---------------------------------------------------------------------- // + real pfpCut = 3.0 * particles_.getShapes().maxBoundingSphere(); + real searchCut = radCut_; + + if (enablePFP_ && pfpCut > searchCut) + { + searchCut = pfpCut; + } + + real cellSize = (searchCut > 1e-12) + ? searchCut + : 3.0 * particles_.getShapes().maxBoundingSphere(); + + REPORT(1) << " Thermal interaction search cell size: " << cellSize << " m" + << END_REPORT; + + mapper_ = makeUnique( + domainBox, + cellSize, + particles_.pointPosition().deviceViewAll(), + particles_.dynPointStruct().activePointsMaskDevice(), + false, + true); +} + +//---------------------------- public methods --------------------------------- + +// ========================================================================= // +// Core Thermal Iteration +// ========================================================================= // + +void thermalInteraction::iterate() +{ + if (!enableRadiation_ && !enableConduction_ && !enablePFP_) + { + return; + } + + thermalTimer_.start(); + bool boxChanged = false; + + bool mapperBuiltOk = mapper_->build( + particles_.pointPosition().deviceViewAll(), + particles_.dynPointStruct().activePointsMaskDevice(), + boxChanged); + + if (!mapperBuiltOk) + { + output + << "\n" + << Yellow_Text("[thermalInteraction] WARNING — mapperNBS failed " + "to build") + << " at step " << stepCounter_ << ".\n" + << " Likely cause: a burned-out ghost particle has an extreme " + << "position\n" + << " that forces the search box beyond allocatable limits.\n" + << " Thermal interactions (Q_pp, Q_pfp, radiation) are SKIPPED " + << "this step.\n" + << endl; + + thermalTimer_.end(); + stepCounter_++; + return; + } + + auto searchBox = mapper_->getSearchCells(); + auto domainMin = searchBox.domainBox().minPoint(); + auto cellSize = searchBox.cellSize(); + int32x3 numCells(searchBox.nx(), searchBox.ny(), searchBox.nz()); + + bool calcRad = enableRadiation_ && (stepCounter_ % radUpdateInterval_ == 0); + bool calcCond = enableConduction_; + bool calcPFP = enablePFP_; + + if (calcCond) + { + Kokkos::deep_copy(particles_.heatSourceCondPP().deviceViewAll(), 0.0); + } + + if (calcPFP) + { + Kokkos::deep_copy(particles_.heatSourcePFP().deviceViewAll(), 0.0); + } + + thermalInteractionKernels::calcThermalInteractions( + particles_.dynPointStruct().activePointsMaskDevice(), + particles_.pointPosition().deviceViewAll(), + particles_.velocity().deviceViewAll(), + particles_.rVelocity().deviceViewAll(), + particles_.diameter().deviceViewAll(), + particles_.mass().deviceViewAll(), + particles_.temperature().deviceViewAll(), + particles_.Cp().deviceViewAll(), + particles_.conductivity().deviceViewAll(), + particles_.E0().deviceViewAll(), + particles_.nu().deviceViewAll(), + particles_.fluidKappa().deviceViewAll(), + particles_.fluidAlpha().deviceViewAll(), + mapper_->getCellIterator(), + domainMin, + cellSize, + numCells, + radCut_, + simYoungsModulus_, + calcRad, + calcCond, + calcPFP, + particles_.heatSourceCondPP().deviceViewAll(), + particles_.heatSourcePFP().deviceViewAll(), + particles_.radSumTemp().deviceViewAll(), + particles_.radNumPrt().deviceViewAll()); + + thermalTimer_.end(); + stepCounter_++; +} + +//+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +} // pFlow + + + + diff --git a/src/Interaction/thermalInteraction/thermalInteraction.hpp b/src/Interaction/thermalInteraction/thermalInteraction.hpp new file mode 100644 index 000000000..2d11f8ca3 --- /dev/null +++ b/src/Interaction/thermalInteraction/thermalInteraction.hpp @@ -0,0 +1,120 @@ +/*------------------------------- phasicFlow --------------------------------- + O C enter of + O O E ngineering and + O O M ultiscale modeling of + OOOOOOO F luid flow +------------------------------------------------------------------------------ + Copyright (C): www.cemf.ir + email: hamid.r.norouzi AT gmail.com +------------------------------------------------------------------------------ +Licence: + This file is part of phasicFlow code. It is a free software for simulating + granular and multiphase flows. You can redistribute it and/or modify it under + the terms of GNU General Public License v3 or any other later versions. + + phasicFlow is distributed to help others in their research in the field of + granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the + implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + +-----------------------------------------------------------------------------*/ + +#ifndef pFlow_thermalInteraction_hpp +#define pFlow_thermalInteraction_hpp + +#include "systemControl.hpp" +#include "thermalSphereParticles.hpp" +#include "mapperNBS.hpp" +#include "Timer.hpp" + +namespace pFlow +{ + +/** + * @brief Dispatcher for intra-phase thermodynamic interactions. + * * Manages the calculation of particle-particle conduction (Q_pp), + * sub-grid Particle-Fluid-Particle heat transfer (Q_pfp), and + * local radiation neighbourhood sums. + */ +class thermalInteraction +{ +public: + + //- Type info + + TypeInfo("thermalInteraction"); + +private: + + //- private members + + // --- Section 1: System References --- + + systemControl& control_; + + const thermalSphereParticles& particles_; + + uniquePtr mapper_ = nullptr; + + // --- Section 2: Physics Control Flags --- + + /// @brief Toggles radiation neighbourhood calculations. + bool enableRadiation_ = false; + + uint32 radUpdateInterval_ = 1; + + real radCut_ = 0.0; + + /// @brief Toggles direct particle-particle contact conduction. + bool enableConduction_ = false; + + real simYoungsModulus_ = 1e7; + + /// @brief Toggles sub-grid fluid bridge heat transfer (PFP). + bool enablePFP_ = false; + + // --- Section 3: Performance & Tracking --- + + uint32 stepCounter_ = 0; + + Timer thermalTimer_; + +public: + + //- constructors + + /** + * @brief Constructs the thermal interaction manager. + */ + thermalInteraction( + systemControl& control, + const thermalSphereParticles& prtcl, + const box& domainBox); + + ~thermalInteraction() = default; + + //- public methods + + // --- Section 4: Public Interface --- + + /** + * @brief Checks if radiation physics is actively executing. + * @return True if radiation is globally enabled by the user. + */ + inline + bool isRadiationEnabled() const + { + return enableRadiation_; + } + + /** + * @brief Executes the neighbor-search and thermodynamic kernels. + */ + void iterate(); + +}; // thermalInteraction + +} // pFlow + +#endif // pFlow_thermalInteraction_hpp + + diff --git a/src/Interaction/thermalInteraction/thermalInteractionKernels.cpp b/src/Interaction/thermalInteraction/thermalInteractionKernels.cpp new file mode 100644 index 000000000..1b6777bfb --- /dev/null +++ b/src/Interaction/thermalInteraction/thermalInteractionKernels.cpp @@ -0,0 +1,404 @@ +/*------------------------------- phasicFlow --------------------------------- + O C enter of + O O E ngineering and + O O M ultiscale modeling of + OOOOOOO F luid flow +------------------------------------------------------------------------------ + Copyright (C): www.cemf.ir + email: hamid.r.norouzi AT gmail.com +------------------------------------------------------------------------------ +Licence: + This file is part of phasicFlow code. It is a free software for simulating + granular and multiphase flows. You can redistribute it and/or modify it under + the terms of GNU General Public License v3 or any other later versions. + + phasicFlow is distributed to help others in their research in the field of + granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the + implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + +-----------------------------------------------------------------------------*/ + +#include "thermalInteractionKernels.hpp" +#include + +namespace pFlow +{ +namespace thermalInteractionKernels +{ + +// Gauss-Legendre quadrature abscissae/weights for the PFP flux integral. +// Kept at file scope (rather than re-built per thread) to reduce +// per-thread register pressure on the GPU. +KOKKOS_INLINE_FUNCTION +constexpr real t_GL[5] = { + -0.9061798459, -0.5384693101, 0.0, 0.5384693101, 0.9061798459 +}; + +KOKKOS_INLINE_FUNCTION +constexpr real w_GL[5] = { + 0.2369268850, 0.4786286705, 0.5688888889, 0.4786286705, 0.2369268850 +}; + +// Dynamic scheduling: neighbour counts vary strongly between dense and +// dilute regions of the particle bed, so static chunking would leave some +// threads idle while others are still sweeping crowded cells. +using policy = Kokkos::RangePolicy< + pFlow::DefaultExecutionSpace, + Kokkos::Schedule, + Kokkos::IndexType>; + +void calcThermalInteractions( + const pFlagTypeDevice& m, + const deviceViewType1D& pos, + const deviceViewType1D& tvel, + const deviceViewType1D& rvel, + const deviceViewType1D& diameter, + const deviceViewType1D& mass, + const deviceViewType1D& temperature, + const deviceViewType1D& Cp, + const deviceViewType1D& K, + const deviceViewType1D& E0, + const deviceViewType1D& nu, + const deviceViewType1D& fluidKappa, + const deviceViewType1D& fluidAlpha, + const mapperNBS::CellIterator& cellIter, + const realx3& domainMin, + const real& cellSize, + const int32x3& numCells, + const real radCut, + const real simYoungsModulus, + const bool calcRad, + const bool calcCond, + const bool calcPFP, + deviceViewType1D Q_pp, + deviceViewType1D Q_pfp, + deviceViewType1D radSumTemp, + deviceViewType1D radNumPrt) +{ + auto r = m.activeRange(); + + Kokkos::parallel_for( + "calcThermalInteractions", + policy(r.start(), r.end()), + KOKKOS_LAMBDA(uint32 i) + { + if (m(i)) + { + realx3 p_i = pos[i]; + real R_i = 0.5 * diameter[i]; + real T_i = temperature[i]; + + real pi = Kokkos::numbers::pi_v; + + real radCutSq = radCut * radCut; + real sumT = 0.0; + uint32 count = 0; + + int32 c_x = static_cast( + (p_i.x() - domainMin.x()) / cellSize); + int32 c_y = static_cast( + (p_i.y() - domainMin.y()) / cellSize); + int32 c_z = static_cast( + (p_i.z() - domainMin.z()) / cellSize); + + // Sweep the immediate 27-cell neighborhood + for (int32 cx = c_x - 1; cx <= c_x + 1; ++cx) + { + for (int32 cy = c_y - 1; cy <= c_y + 1; ++cy) + { + for (int32 cz = c_z - 1; cz <= c_z + 1; ++cz) + { + if (cx >= 0 && cx < numCells.x() && + cy >= 0 && cy < numCells.y() && + cz >= 0 && cz < numCells.z()) + { + uint32 j = cellIter.start(cx, cy, cz); + + while (j != mapperNBS::CellIterator::NoPos) + { + if (i != j && m(j)) + { + real dx = p_i.x() - pos[j].x(); + real dy = p_i.y() - pos[j].y(); + real dz = p_i.z() - pos[j].z(); + + real distSq = dx*dx + dy*dy + dz*dz; + + // ===================================== + // 1. Radiation Check (Asymmetric + // execution, all i-j pairs) + // ===================================== + if (calcRad && distSq <= radCutSq) + { + sumT += temperature[j]; + count++; + } + + // ===================================== + // 2. Conduction & PFP Checks (Symmetric + // execution: i < j only). Computing + // both forces here halves the + // computational load. + // ===================================== + if ((calcCond || calcPFP) && i < j) + { + real R_j = 0.5 * diameter[j]; + real sumRadiiSq = + (R_i + R_j) * (R_i + R_j); + + bool isContact = + (distSq < sumRadiiSq); + real dist = sqrt(distSq); + real rc_real = 0.0; + + // --- 2.A: Static Contact + // Conduction + // (Eq. 6.159) --- + if (isContact && dist > 1e-12) + { + real R_eff = + (R_i * R_j) / (R_i + R_j); + + // Inverse-modulus quantities + // for the simulation's Young's + // modulus and the material's + // real Young's modulus. + real term_E_sim = + (1.0 - nu[i]*nu[i]) / + simYoungsModulus + + (1.0 - nu[j]*nu[j]) / + simYoungsModulus; + + real term_E_real = + (1.0 - nu[i]*nu[i]) / E0[i]+ + (1.0 - nu[j]*nu[j]) / E0[j]; + + // Geometric contact radius from + // the normal overlap of the two + // spheres: a^2 = R_eff*overlap, + // the standard Hertzian + // relationship between + // interpenetration and + // contact-patch radius. Purely + // geometric - no relative + // velocity or contact-time + // model involved, matching the + // static/sustained contact + // regime of Eq. 6.159. + real overlap = + (R_i + R_j) - dist; + real rc_geom = + sqrt(R_eff * overlap); + + // Contact radius correction + // (Eq. 6.166 & 6.167): the + // overlap above was produced + // using the simulation's + // Young's modulus, so rc_geom + // overstates the contact radius + // a real, stiffer material + // would give for the same + // geometric interpenetration. + // c = (E_sim/E_real)^0.2; in + // terms of the inverse-modulus + // quantities above (term_E=1/E) + // this is: + // c=(term_E_real/term_E_sim)^.2 + real c_corr = + pow(term_E_real / + term_E_sim, + 0.2); + + rc_real = c_corr * rc_geom; + + if (calcCond) + { + real tempDiff = + temperature[j] - T_i; + real num = + 4.0 * rc_real * + tempDiff; + real den = + (1.0 / K[i]) + + (1.0 / K[j]); + real Q_rate = num / den; + + // Atomic accumulation for + // thread safety + Kokkos::atomic_add( + &Q_pp[i], Q_rate); + Kokkos::atomic_add( + &Q_pp[j], -Q_rate); + } + } + + // --- 2.B: Particle-Fluid-Particle + // Sub-grid Heat Transfer + // (Eq. 6.160) --- + if (calcPFP && dist > 1e-12) + { + real R_star = 0.5 * (R_i + R_j); + real H = 0.5 * + (dist - R_i - R_j); + real k_f = 0.5 * + (fluidKappa[i] + + fluidKappa[j]); + + // Cut-off rule: Ignored if + // H/R* > 0.5 + if (H / R_star <= 0.5 && + k_f > 1e-12) + { + // Eq. 6.164: r_ij + // evaluation based on + // local porosity + real r_sij = + isContact ? + rc_real : 0.0; + + real eps_avg = 0.5 * + (fluidAlpha[i] + + fluidAlpha[j]); + + real solid_frac = + 1.0 - eps_avg; + + if (solid_frac < 0.01) + { + // Clamp to avoid inf + solid_frac = 0.01; + } + + real r_ij = 0.56 * R_star * + pow(solid_frac, + -1.0/3.0); + + // Eq. 6.162: r_sf (Upper + // limit of integration) + real R_H = R_star + + (H > 0.0 ? H : 0.0); + + real r_sf = (R_star * r_ij)/ + sqrt(r_ij*r_ij + + R_H*R_H); + + if (r_sf > r_sij) + { + // Loop-free 5-point + // Gauss-Legendre + // Quadrature + real A = r_sij; + real B = r_sf; + real c1 = 0.5 * (B - A); + real c2 = 0.5 * (A + B); + + real integral = 0.0; + + // Explicit unrolled + // loop for GPU + // registers + // (Thread-safe) + for (int k=0; k<5; ++k) + { + real r_pt = + c1 * t_GL[k] + + c2; + real r2 = r_pt*r_pt; + + real Ri2 = R_i*R_i; + real root_i = 0.0; + if (Ri2 > r2) + { + root_i = + sqrt(Ri2 - + r2); + } + + real Rj2 = R_j*R_j; + real root_j = 0.0; + if (Rj2 > r2) + { + root_j = + sqrt(Rj2 - + r2); + } + + // Lens gap physical + // thickness + real gap = dist - + root_i - root_j; + + if (gap < 0.0) + { + gap = 0.0; + } + + real term_i = + (R_i - root_i) / + K[i]; + real term_j = + (R_j - root_j) / + K[j]; + real term_f = + gap / k_f; + + real R_th = term_i + + term_j + term_f; + + // Avoid division by + // zero at perfect + // rigid contact + // centers + if (R_th > 1e-12) + { + real F = + (2.0 * pi * + r_pt)/R_th; + integral += + w_GL[k] * F; + } + } + integral *= c1; + + // Apply integrated PFP + // flux symmetrically + real Q_pfp_val = + integral * + (temperature[j] - + T_i); + + Kokkos::atomic_add( + &Q_pfp[i], + Q_pfp_val); + Kokkos::atomic_add( + &Q_pfp[j], + -Q_pfp_val); + } + } + } + } + } + j = cellIter.next(j); + } + } + } + } + } + + // Finalize Radiation (Saves sum to memory) + if (calcRad) + { + radSumTemp[i] = sumT; + radNumPrt[i] = count; + } + } + }); + + Kokkos::fence(); +} + +} // thermalInteractionKernels +} // pFlow + + + diff --git a/src/Interaction/thermalInteraction/thermalInteractionKernels.hpp b/src/Interaction/thermalInteraction/thermalInteractionKernels.hpp new file mode 100644 index 000000000..46cbfb8b0 --- /dev/null +++ b/src/Interaction/thermalInteraction/thermalInteractionKernels.hpp @@ -0,0 +1,140 @@ +/*------------------------------- phasicFlow --------------------------------- + O C enter of + O O E ngineering and + O O M ultiscale modeling of + OOOOOOO F luid flow +------------------------------------------------------------------------------ + Copyright (C): www.cemf.ir + email: hamid.r.norouzi AT gmail.com +------------------------------------------------------------------------------ +Licence: + This file is part of phasicFlow code. It is a free software for simulating + granular and multiphase flows. You can redistribute it and/or modify it under + the terms of GNU General Public License v3 or any other later versions. + + phasicFlow is distributed to help others in their research in the field of + granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the + implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + +-----------------------------------------------------------------------------*/ + +#ifndef pFlow_thermalInteractionKernels_hpp +#define pFlow_thermalInteractionKernels_hpp + +#include "types.hpp" +#include "pointFlag.hpp" +#include "mapperNBS.hpp" + +namespace pFlow +{ +namespace thermalInteractionKernels +{ + +/** + * @brief GPU kernel computing collisional (Q_pp), Particle-Fluid-Particle + * (Q_pfp), and radiation heat transfer for all active particle pairs. + * + * @details + * Dispatched by thermalInteraction::iterate() once per DEM sub-step. + * The kernel performs three independent calculations in a single neighbour- + * search sweep: + * + * ### 1. Static-contact conduction (Q_pp) + * Heat conduction through solid–solid contact between particles in + * sustained/static contact. The contact radius is computed purely from + * geometry: rc = sqrt(R_eff * overlap), the standard Hertzian relationship + * between the normal interpenetration of two spheres and their circular + * contact-patch radius. No relative velocity or contact-time model is + * involved. Because the simulation typically uses a softened Young's + * modulus to allow a larger DEM timestep, the resulting overlap - and + * hence rc - overstates what a real, stiffer material would produce; a + * correction factor derived from the ratio of simulated to real Young's + * modulus scales rc back down before it is used in the heat transfer + * rate. Enabled when `calcCond = true`. + * + * ### 2. Particle-Fluid-Particle sub-grid heat transfer (Q_pfp) + * When two particles are within a dimensionless gap H/R* <= 0.5, the fluid + * between them forms a thin thermal bridge. The effective conductance is: + * + * Q_pfp = 2pi k_f R* \int_0^{H/R*} r^2/(h + r^2/2R*) dr + * + * evaluated using 5-point Gauss-Legendre quadrature (Rong & Horio, 1999). + * Enabled when `calcPFP = true`. + * + * #### Sub-grid energy approximation + * The fluid energy equation (TEqn on the CFD side) does **not** contain + * a corresponding sink term for the energy consumed by the PFP bridge. + * This is an intentional sub-grid approximation based on the following: + * + * 1. PFP operates below CFD mesh resolution in the unresolved regime. + * 2. The heat redistributes between solid particles; the fluid acts only + * as a passive conduit — energy given to particle i is received from + * particle j and vice versa. Net fluid energy change approx 0 at the + * cell scale. + * 3. For dilute-to-moderate packing (alpha_s < 0.45), PFP is typically one + * order of magnitude smaller than convective heat transfer, introducing + * < 5 % error in the fluid energy budget. + * + * **Validity conditions:** + * - Cell volume >> particle volume (unresolved regime, V_cell/V_p > 10) + * - Solid conductivity >> fluid conductivity (k_s/k_f >> 1) + * - Solid volume fraction alpha_s < 0.45 + * + * For dense beds (alpha_s > 0.5) or particles with low k_s/k_f, the fluid-side + * energy sink from PFP should be evaluated and added to TEqn explicitly. + * + * #### Note on fluid property sampling + * `fluidKappa` and `fluidAlpha` are sampled at the particle's Eulerian cell + * centre (not at the bridge midpoint). The resulting discretisation error + * is O(d_p / Delta x), negligible for d_p << Delta x. + * + * ### 3. Radiation (neighbourhood-based) + * Accumulates radiating-neighbour temperature sums and counts for the + * linearised radiation model in sphereHeatTransfer. Radiative exchange + * happens only between solid particles; the carrier fluid is radiatively + * transparent and never appears in this accumulation. Enabled when + * `calcRad = true`. + * + * @param calcRad Enable radiation neighbourhood accumulation. + * @param calcCond Enable particle-particle static-contact conduction (Q_pp). + * @param calcPFP Enable particle-fluid-particle sub-grid transfer (Q_pfp). + */ +void calcThermalInteractions( + const pFlagTypeDevice& m, + const deviceViewType1D& pos, + const deviceViewType1D& tvel, + const deviceViewType1D& rvel, + const deviceViewType1D& diameter, + const deviceViewType1D& mass, + const deviceViewType1D& temperature, + const deviceViewType1D& Cp, + const deviceViewType1D& K, + const deviceViewType1D& E0, + const deviceViewType1D& nu, + // local fluid kappa at particle cell [W/(m.K)] + const deviceViewType1D& fluidKappa, + // local fluid porosity alpha at particle cell [-] + const deviceViewType1D& fluidAlpha, + const mapperNBS::CellIterator& cellIter, + const realx3& domainMin, + const real& cellSize, + const int32x3& numCells, + const real radCut, + const real simYoungsModulus, + const bool calcRad, + const bool calcCond, + const bool calcPFP, + // static-contact conduction output [W] + deviceViewType1D Q_pp, + // PFP sub-grid output [W] + deviceViewType1D Q_pfp, + deviceViewType1D radSumTemp, + deviceViewType1D radNumPrt); + +} // thermalInteractionKernels +} // pFlow + +#endif // pFlow_thermalInteractionKernels_hpp + + + diff --git a/src/Particles/SphereParticles/thermalSphereParticles/thermalSphereParticles.cpp b/src/Particles/SphereParticles/thermalSphereParticles/thermalSphereParticles.cpp new file mode 100644 index 000000000..e93f899e8 --- /dev/null +++ b/src/Particles/SphereParticles/thermalSphereParticles/thermalSphereParticles.cpp @@ -0,0 +1,464 @@ +/*------------------------------- phasicFlow --------------------------------- + O C enter of + O O E ngineering and + O O M ultiscale modeling of + OOOOOOO F luid flow +------------------------------------------------------------------------------ + Copyright (C): www.cemf.ir + email: hamid.r.norouzi AT gmail.com +------------------------------------------------------------------------------ +Licence: + This file is part of phasicFlow code. It is a free software for simulating + granular and multiphase flows. You can redistribute it and/or modify it under + the terms of GNU General Public License v3 or any other later versions. + + phasicFlow is distributed to help others in their research in the field of + granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the + implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + +-----------------------------------------------------------------------------*/ + +#include "thermalSphereParticles.hpp" +#include "thermalSphereParticlesKernels.hpp" +#include + +namespace pFlow +{ + +//----------------------------- protected methods ----------------------------- + +// ========================================================================= // +// Section 1: Memory Management +// ========================================================================= // + +void thermalSphereParticles::checkHostMemory() +{ + sphereFluidParticles::checkHostMemory(); + + if (temperature_.size() != temperatureHost_.size()) + { + size_t oldSize = temperatureHost_.size(); + size_t newSize = temperature_.size(); + + resizeNoInit(temperatureHost_, newSize); + resizeNoInit(heatSourceConvHost_, newSize); + resizeNoInit(heatSourceRadHost_, newSize); + resizeNoInit(heatSourceCondPPHost_, newSize); + resizeNoInit(emissivityHost_, newSize); + resizeNoInit(radSumTempHost_, newSize); + resizeNoInit(radNumPrtHost_, newSize); + resizeNoInit(fluidKappaHost_, newSize); + resizeNoInit(fluidAlphaHost_, newSize); + + for (size_t i = oldSize; i < newSize; ++i) + { + temperatureHost_[i] = temperature_.field()[i]; + heatSourceConvHost_[i] = heatSourceConv_.field()[i]; + heatSourceRadHost_[i] = heatSourceRad_.field()[i]; + heatSourceCondPPHost_[i] = heatSourceCondPP_.field()[i]; + emissivityHost_[i] = emissivity_.field()[i]; + radSumTempHost_[i] = radSumTemp_.field()[i]; + radNumPrtHost_[i] = radNumPrt_.field()[i]; + fluidKappaHost_[i] = fluidKappa_.field()[i]; + fluidAlphaHost_[i] = fluidAlpha_.field()[i]; + } + } +} + +//----------------------------- constructors ---------------------------------- + +// ========================================================================= // +// Section 2: Constructor & Initialization +// ========================================================================= // + +thermalSphereParticles::thermalSphereParticles( + systemControl& control, + const sphereShape& shpShape, + const thermalSphereShape& thShpShape) +: + sphereFluidParticles(control, shpShape), + thSpheres_(thShpShape), + temperature_( + objectFile( + "temperature", + "", + objectFile::READ_ALWAYS, + objectFile::WRITE_ALWAYS), + dynPointStruct(), + 0.0), + Cp_( + objectFile( + "Cp", + "", + objectFile::READ_NEVER, + objectFile::WRITE_NEVER), + dynPointStruct(), + 1.0), + conductivity_( + objectFile( + "conductivity", + "", + objectFile::READ_NEVER, + objectFile::WRITE_NEVER), + dynPointStruct(), + 1.0), + heatSourceConv_( + objectFile( + "heatSourceConv", + "", + objectFile::READ_ALWAYS, + objectFile::WRITE_ALWAYS), + dynPointStruct(), + 0.0), + heatSourceRad_( + objectFile( + "heatSourceRad", + "", + objectFile::READ_ALWAYS, + objectFile::WRITE_ALWAYS), + dynPointStruct(), + 0.0), + heatSourceCondPP_( + objectFile( + "heatSourceCondPP", + "", + objectFile::READ_ALWAYS, + objectFile::WRITE_ALWAYS), + dynPointStruct(), + 0.0), + heatSourcePFP_( + objectFile( + "heatSourcePFP", + "", + objectFile::READ_ALWAYS, + objectFile::WRITE_ALWAYS), + dynPointStruct(), + 0.0), + emissivity_( + objectFile( + "emissivity", + "", + objectFile::READ_ALWAYS, + objectFile::WRITE_ALWAYS), + dynPointStruct(), + 0.0), + radSumTemp_( + objectFile( + "radSumTemp", + "", + objectFile::READ_ALWAYS, + objectFile::WRITE_ALWAYS), + dynPointStruct(), + 0.0), + radNumPrt_( + objectFile( + "radNumPrt", + "", + objectFile::READ_ALWAYS, + objectFile::WRITE_ALWAYS), + dynPointStruct(), + 0u), + E0_( + objectFile( + "E0", + "", + objectFile::READ_NEVER, + objectFile::WRITE_NEVER), + dynPointStruct(), + 1e9), + nu_( + objectFile( + "nu", + "", + objectFile::READ_NEVER, + objectFile::WRITE_NEVER), + dynPointStruct(), + 0.3), + temperatureRate_( + objectFile( + "temperatureRate", + "", + objectFile::READ_NEVER, + objectFile::WRITE_NEVER), + dynPointStruct(), + 0.0), + fluidKappa_( + objectFile( + "fluidKappa", + "", + objectFile::READ_ALWAYS, + objectFile::WRITE_ALWAYS), + dynPointStruct(), + 0.0), + fluidAlpha_( + objectFile( + "fluidAlpha", + "", + objectFile::READ_ALWAYS, + objectFile::WRITE_ALWAYS), + dynPointStruct(), + 0.0), + heatTransferTimer_("heatTransfer", &this->timers()), + temperatureIntegrationTimer_("tempInt", &this->timers()) +{ + initializeThermalParticles(); + checkHostMemory(); + + temperatureHostUpdatedSync(); + radiationDataHostUpdatedSync(); +} + +//---------------------------- public methods --------------------------------- + +bool thermalSphereParticles::initializeThermalParticles() +{ + auto activeMask = this->dynPointStruct().activePointsMaskDevice(); + + realVector h_Cp = thSpheres_.heatCapacities(); + realVector h_K = thSpheres_.heatConductivities(); + realVector h_Eps = thSpheres_.emissivities(); + realVector h_E0 = thSpheres_.realYoungsModuli(); + realVector h_Nu = thSpheres_.poissonRatios(); + + deviceViewType1D d_Cp ("dCp", h_Cp.size()); + deviceViewType1D d_K ("dK", h_K.size()); + deviceViewType1D d_Eps("dEps", h_Eps.size()); + deviceViewType1D d_E0 ("dE0", h_E0.size()); + deviceViewType1D d_Nu ("dNu", h_Nu.size()); + + auto m_Cp = Kokkos::create_mirror_view(d_Cp); + auto m_K = Kokkos::create_mirror_view(d_K); + auto m_Eps = Kokkos::create_mirror_view(d_Eps); + auto m_E0 = Kokkos::create_mirror_view(d_E0); + auto m_Nu = Kokkos::create_mirror_view(d_Nu); + + for (size_t i = 0; i < h_Cp.size(); ++i) + { + m_Cp (i) = h_Cp [i]; + m_K (i) = h_K [i]; + m_Eps(i) = h_Eps[i]; + m_E0 (i) = h_E0 [i]; + m_Nu (i) = h_Nu [i]; + } + + Kokkos::deep_copy(d_Cp, m_Cp); + Kokkos::deep_copy(d_K, m_K); + Kokkos::deep_copy(d_Eps, m_Eps); + Kokkos::deep_copy(d_E0, m_E0); + Kokkos::deep_copy(d_Nu, m_Nu); + + thermalSphereParticlesKernels::initThermalProperties( + activeMask, + shapeIndex().deviceViewAll(), + Cp_.deviceViewAll(), + conductivity_.deviceViewAll(), + emissivity_.deviceViewAll(), + E0_.deviceViewAll(), + nu_.deviceViewAll(), + d_Cp, + d_K, + d_Eps, + d_E0, + d_Nu); + + return true; +} + +// ========================================================================= // +// Section 3: Core Iteration Logic +// ========================================================================= // + +bool thermalSphereParticles::beforeIteration() +{ + sphereFluidParticles::beforeIteration(); + checkHostMemory(); + + if (heatSourceConvHost_.size() == heatSourceConv_.deviceView().size()) + { + Kokkos::deep_copy(heatSourceConv_.deviceView(), heatSourceConvHost_); + } + + if (heatSourceRadHost_.size() == heatSourceRad_.deviceView().size()) + { + Kokkos::deep_copy(heatSourceRad_.deviceView(), heatSourceRadHost_); + } + + temperatureRate_.field().fill(0.0); + + temperatureHostUpdatedSync(); + radiationDataHostUpdatedSync(); + + return true; +} + +bool thermalSphereParticles::iterate() +{ + if (!sphereFluidParticles::iterate()) + { + return false; + } + + auto mask = dynPointStruct().activePointsMaskDevice(); + + heatTransferTimer_.start(); + + thermalSphereParticlesKernels::calcFluidParticleHeatTransfer( + mask, + diameter().deviceViewAll(), + mass().deviceViewAll(), + Cp().deviceViewAll(), + temperature().deviceViewAll(), + heatSourceConv_.deviceViewAll(), + heatSourceRad_.deviceViewAll(), + heatSourceCondPP_.deviceViewAll(), + heatSourcePFP_.deviceViewAll(), + temperatureRate_.deviceViewAll()); + + heatTransferTimer_.end(); + + temperatureIntegrationTimer_.start(); + + thermalSphereParticlesKernels::integrateTemperature( + mask, + control().time().dt(), + temperature().deviceViewAll(), + temperatureRate_.deviceViewAll()); + + temperatureIntegrationTimer_.end(); + + return true; +} + +// ========================================================================= // +// Section 4: Particle Insertion +// ========================================================================= // + +bool thermalSphereParticles::insertParticles( + const realx3Vector& pos, + const wordVector& names, + const anyList& vars) +{ + anyList nv(vars); + + realVector cpV ("Cp"); + realVector kV ("k"); + realVector epsV ("emissivity"); + realVector e0V ("E0"); + realVector nuV ("nu"); + realVector tV ("T"); + realVector kappaV ("fluidKappa"); + realVector alphaV ("fluidAlpha"); + realVector pfpV ("heatSourcePFP"); + + for (const auto& name : names) + { + uint32 i; + if (thSpheres_.shapeNameToIndex(name, i)) + { + cpV .push_back(thSpheres_.heatCapacity(i)); + kV .push_back(thSpheres_.heatConductivity(i)); + epsV.push_back(thSpheres_.emissivity(i)); + e0V .push_back(thSpheres_.realYoungsModulus(i)); + nuV .push_back(thSpheres_.poissonRatio(i)); + + kappaV.push_back(0.0); + alphaV.push_back(0.0); + pfpV .push_back(0.0); + + // Fetch insertion temperature directly from properties + tV.push_back(thSpheres_.insertionTemperature()); + } + } + + nv.emplaceBack(Cp_.name() + "Vector", std::move(cpV)); + nv.emplaceBack(conductivity_.name() + "Vector", std::move(kV)); + nv.emplaceBack(emissivity_.name() + "Vector", std::move(epsV)); + nv.emplaceBack(E0_.name() + "Vector", std::move(e0V)); + nv.emplaceBack(nu_.name() + "Vector", std::move(nuV)); + nv.emplaceBack(temperature_.name() + "Vector", std::move(tV)); + nv.emplaceBack(fluidKappa_.name() + "Vector", std::move(kappaV)); + nv.emplaceBack(fluidAlpha_.name() + "Vector", std::move(alphaV)); + nv.emplaceBack(heatSourcePFP_.name() + "Vector", std::move(pfpV)); + + return sphereFluidParticles::insertParticles(pos, names, nv); +} + +// ========================================================================= // +// Section 5: MPI Synchronisation Routines (Host <-> Device) +// ========================================================================= // + +void thermalSphereParticles::heatSourcesHostUpdatedSync() +{ + checkHostMemory(); + + bool sizeConv = + (heatSourceConvHost_.size() == heatSourceConv_.deviceView().size()); + bool sizeRad = + (heatSourceRadHost_.size() == heatSourceRad_.deviceView().size()); + bool sizeCond = + (heatSourceCondPPHost_.size() == heatSourceCondPP_.deviceView().size()); + + if (sizeConv && sizeRad && sizeCond) + { + Kokkos::deep_copy( + heatSourceConv_.deviceView(), + heatSourceConvHost_); + + Kokkos::deep_copy( + heatSourceRad_.deviceView(), + heatSourceRadHost_); + + Kokkos::deep_copy( + heatSourceCondPP_.deviceView(), + heatSourceCondPPHost_); + } +} + +void thermalSphereParticles::fluidPropertiesHostUpdatedSync() +{ + checkHostMemory(); + + bool sizeKappa = + (fluidKappaHost_.size() == fluidKappa_.deviceView().size()); + bool sizeAlpha = + (fluidAlphaHost_.size() == fluidAlpha_.deviceView().size()); + + if (sizeKappa && sizeAlpha) + { + Kokkos::deep_copy(fluidKappa_.deviceView(), fluidKappaHost_); + Kokkos::deep_copy(fluidAlpha_.deviceView(), fluidAlphaHost_); + } +} + +void thermalSphereParticles::temperatureHostUpdatedSync() +{ + checkHostMemory(); + + if (temperatureHost_.size() == temperature_.deviceView().size()) + { + Kokkos::deep_copy(temperatureHost_, temperature_.deviceView()); + } +} + +void thermalSphereParticles::radiationDataHostUpdatedSync() +{ + checkHostMemory(); + + bool sizeEps = + (emissivityHost_.size() == emissivity_.deviceView().size()); + bool sizeSum = + (radSumTempHost_.size() == radSumTemp_.deviceView().size()); + bool sizeNum = + (radNumPrtHost_.size() == radNumPrt_.deviceView().size()); + + if (sizeEps && sizeSum && sizeNum) + { + Kokkos::deep_copy(emissivityHost_, emissivity_.deviceView()); + Kokkos::deep_copy(radSumTempHost_, radSumTemp_.deviceView()); + Kokkos::deep_copy(radNumPrtHost_, radNumPrt_.deviceView()); + } +} + +//+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +} // pFlow diff --git a/src/Particles/SphereParticles/thermalSphereParticles/thermalSphereParticles.hpp b/src/Particles/SphereParticles/thermalSphereParticles/thermalSphereParticles.hpp new file mode 100644 index 000000000..24a01a87f --- /dev/null +++ b/src/Particles/SphereParticles/thermalSphereParticles/thermalSphereParticles.hpp @@ -0,0 +1,367 @@ +/*------------------------------- phasicFlow --------------------------------- + O C enter of + O O E ngineering and + O O M ultiscale modeling of + OOOOOOO F luid flow +------------------------------------------------------------------------------ + Copyright (C): www.cemf.ir + email: hamid.r.norouzi AT gmail.com +------------------------------------------------------------------------------ +Licence: + This file is part of phasicFlow code. It is a free software for simulating + granular and multiphase flows. You can redistribute it and/or modify it under + the terms of GNU General Public License v3 or any other later versions. + + phasicFlow is distributed to help others in their research in the field of + granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the + implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + +-----------------------------------------------------------------------------*/ + +#ifndef pFlow_thermalSphereParticles_hpp +#define pFlow_thermalSphereParticles_hpp + +#include "sphereFluidParticles.hpp" +#include "thermalSphereShape.hpp" + +namespace pFlow +{ + +/** + * @brief Manages the thermal state and thermodynamic properties of + * spherical particles on the GPU. + * + * @details + * Extends sphereFluidParticles by introducing device (Kokkos) memory for: + * - Temperatures and integration rates (Explicit Euler). + * - Thermodynamic properties (heat capacities, conductivities, emissivities). + * - Multi-mode heat sources (Convection, Radiation, Conduction, PFP). + * - Host (CPU) mirror arrays used for MPI and OpenFOAM coupling synchronization. + */ +class thermalSphereParticles +: + public sphereFluidParticles +{ +public: + + //- Type info + + TypeInfo("thermalSphereParticles"); + +private: + + //- private members + + // --- Section 1: Shape Reference --- + + const thermalSphereShape& thSpheres_; + + // --- Section 2: Device Fields (GPU/Kokkos) --- + + realPointField_D temperature_; + realPointField_D Cp_; + realPointField_D conductivity_; + realPointField_D temperatureRate_; + + realPointField_D heatSourceConv_; + realPointField_D heatSourceRad_; + realPointField_D heatSourceCondPP_; + realPointField_D heatSourcePFP_; + + realPointField_D emissivity_; + realPointField_D radSumTemp_; + uint32PointField_D radNumPrt_; + + realPointField_D E0_; + realPointField_D nu_; + + realPointField_D fluidKappa_; + realPointField_D fluidAlpha_; + + // --- Section 3: Performance Timers --- + + Timer heatTransferTimer_; + Timer temperatureIntegrationTimer_; + + // --- Section 4: Host Mirror Fields (CPU RAM) --- + + hostViewType1D temperatureHost_; + hostViewType1D heatSourceConvHost_; + hostViewType1D heatSourceRadHost_; + hostViewType1D heatSourceCondPPHost_; + hostViewType1D emissivityHost_; + hostViewType1D radSumTempHost_; + hostViewType1D radNumPrtHost_; + hostViewType1D fluidKappaHost_; + hostViewType1D fluidAlphaHost_; + +protected: + + //- protected methods + + // --- Section 5: Memory Management --- + + /** + * @brief Ensures host arrays are sized to match their corresponding + * device arrays and initializes newly allocated memory slots. + */ + void checkHostMemory(); + +public: + + //- constructors + + // --- Section 6: Constructor and Initialization --- + + thermalSphereParticles( + systemControl& control, + const sphereShape& shpShape, + const thermalSphereShape& thShpShape); + + ~thermalSphereParticles() override = default; + + //- public methods + + /** + * @brief Scatters per-material thermal properties to individual + * particle slots on the GPU. + * @return True upon successful mapping. + */ + bool initializeThermalParticles(); + + // --- Section 7: Core Iteration Hooks --- + + bool beforeIteration() override; + + bool iterate() override; + + bool insertParticles( + const realx3Vector& pos, + const wordVector& names, + const anyList& vars) override; + + // --- Section 8: Device Accessors --- + + inline + const realPointField_D& temperature() const + { + return temperature_; + } + + inline + realPointField_D& temperature() + { + return temperature_; + } + + inline + const realPointField_D& Cp() const + { + return Cp_; + } + + inline + const realPointField_D& conductivity() const + { + return conductivity_; + } + + inline + realPointField_D& conductivity() + { + return conductivity_; + } + + inline + const realPointField_D& heatSourceConv() const + { + return heatSourceConv_; + } + + inline + realPointField_D& heatSourceConv() + { + return heatSourceConv_; + } + + inline + const realPointField_D& heatSourceRad() const + { + return heatSourceRad_; + } + + inline + realPointField_D& heatSourceRad() + { + return heatSourceRad_; + } + + inline + const realPointField_D& heatSourceCondPP() const + { + return heatSourceCondPP_; + } + + inline + realPointField_D& heatSourceCondPP() + { + return heatSourceCondPP_; + } + + inline + const realPointField_D& heatSourcePFP() const + { + return heatSourcePFP_; + } + + inline + realPointField_D& heatSourcePFP() + { + return heatSourcePFP_; + } + + inline + const realPointField_D& emissivity() const + { + return emissivity_; + } + + inline + realPointField_D& emissivity() + { + return emissivity_; + } + + inline + const realPointField_D& radSumTemp() const + { + return radSumTemp_; + } + + inline + realPointField_D& radSumTemp() + { + return radSumTemp_; + } + + inline + const uint32PointField_D& radNumPrt() const + { + return radNumPrt_; + } + + inline + uint32PointField_D& radNumPrt() + { + return radNumPrt_; + } + + inline + const realPointField_D& E0() const + { + return E0_; + } + + inline + const realPointField_D& nu() const + { + return nu_; + } + + inline + const realPointField_D& fluidKappa() const + { + return fluidKappa_; + } + + inline + realPointField_D& fluidKappa() + { + return fluidKappa_; + } + + inline + const realPointField_D& fluidAlpha() const + { + return fluidAlpha_; + } + + inline + realPointField_D& fluidAlpha() + { + return fluidAlpha_; + } + + // --- Section 9: Host Accessors --- + + inline + auto& temperatureHost() + { + return temperatureHost_; + } + + inline + auto& heatSourceConvHost() + { + return heatSourceConvHost_; + } + + inline + auto& heatSourceRadHost() + { + return heatSourceRadHost_; + } + + inline + auto& heatSourceCondPPHost() + { + return heatSourceCondPPHost_; + } + + inline + auto& emissivityHost() + { + return emissivityHost_; + } + + inline + auto& radSumTempHost() + { + return radSumTempHost_; + } + + inline + auto& radNumPrtHost() + { + return radNumPrtHost_; + } + + inline + auto& fluidKappaHost() + { + return fluidKappaHost_; + } + + inline + auto& fluidAlphaHost() + { + return fluidAlphaHost_; + } + + // --- Section 10: Synchronisation Routines (Host <-> Device) --- + + void heatSourcesHostUpdatedSync(); + + void fluidPropertiesHostUpdatedSync(); + + void temperatureHostUpdatedSync(); + + void radiationDataHostUpdatedSync(); + +}; // thermalSphereParticles + +} // pFlow + +#endif // pFlow_thermalSphereParticles_hpp + + diff --git a/src/Particles/SphereParticles/thermalSphereParticles/thermalSphereParticlesKernels.cpp b/src/Particles/SphereParticles/thermalSphereParticles/thermalSphereParticlesKernels.cpp new file mode 100644 index 000000000..dc3b20bdd --- /dev/null +++ b/src/Particles/SphereParticles/thermalSphereParticles/thermalSphereParticlesKernels.cpp @@ -0,0 +1,170 @@ +/*------------------------------- phasicFlow --------------------------------- + O C enter of + O O E ngineering and + O O M ultiscale modeling of + OOOOOOO F luid flow +------------------------------------------------------------------------------ + Copyright (C): www.cemf.ir + email: hamid.r.norouzi AT gmail.com +------------------------------------------------------------------------------ +Licence: + This file is part of phasicFlow code. It is a free software for simulating + granular and multiphase flows. You can redistribute it and/or modify it under + the terms of GNU General Public License v3 or any other later versions. + + phasicFlow is distributed to help others in their research in the field of + granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the + implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + +-----------------------------------------------------------------------------*/ + +#include "thermalSphereParticlesKernels.hpp" + +namespace pFlow +{ +namespace thermalSphereParticlesKernels +{ + +// ========================================================================= // +// Section 1: Execution Policy Definition +// Schedule is used for optimal load balancing on uniform particle +// arrays +// ========================================================================= // + +using policy = Kokkos::RangePolicy< + pFlow::DefaultExecutionSpace, + Kokkos::Schedule, + Kokkos::IndexType>; + +// ========================================================================= // +// Section 2: Property Initialization Kernel +// ========================================================================= // + +void initThermalProperties( + const pFlagTypeDevice& m, + const deviceViewType1D& idx, + deviceViewType1D Cp, + deviceViewType1D K, + deviceViewType1D emissivity, + deviceViewType1D E0, + deviceViewType1D nu, + const deviceViewType1D& sCp, + const deviceViewType1D& sK, + const deviceViewType1D& sEps, + const deviceViewType1D& sE0, + const deviceViewType1D& sNu) +{ + auto r = m.activeRange(); + + Kokkos::parallel_for( + "initThermalProps", + policy(r.start(), r.end()), + KOKKOS_LAMBDA(uint32 i) + { + if (m(i)) + { + // Map the particle to its material type + uint32 j = idx[i]; + + // Scatter properties to the main particle arrays + Cp[i] = sCp[j]; + K[i] = sK[j]; + emissivity[i] = sEps[j]; + E0[i] = sE0[j]; + nu[i] = sNu[j]; + } + }); + + Kokkos::fence(); +} + +// ========================================================================= // +// Section 3: Energy Equation Kernel +// ========================================================================= // + +void calcFluidParticleHeatTransfer( + const pFlagTypeDevice& m, + const deviceViewType1D& d, + const deviceViewType1D& ms, + const deviceViewType1D& Cp, + const deviceViewType1D& T, + const deviceViewType1D& Q_conv, + const deviceViewType1D& Q_rad, + const deviceViewType1D& Q_pp, + const deviceViewType1D& Q_pfp, + deviceViewType1D TR) +{ + auto r = m.activeRange(); + + Kokkos::parallel_for( + "calcHeatTransferRate", + policy(r.start(), r.end()), + KOKKOS_LAMBDA(uint32 i) + { + if (m(i)) + { + // --------------------------------------------------------- // + // Particle Energy Equation (Lumped Capacitance Model) + // m * Cp * dT/dt = Q_conv + Q_rad + Q_pp + Q_pfp + // --------------------------------------------------------- // + + // Thermal Inertia = Mass [kg] * Specific Heat Capacity + // [J/(kg.K)] = [J/K] + real thermalInertia = ms[i] * Cp[i]; + + // Protection against division by zero (Thermal Inertia Guard) + if (thermalInertia > 1e-12) + { + TR[i] = (Q_conv[i] + Q_rad[i] + Q_pp[i] + Q_pfp[i]) / + thermalInertia; + } + else + { + TR[i] = 0.0; + } + } + }); + + Kokkos::fence(); +} + +// ========================================================================= // +// Section 4: Time Integration Kernel +// ========================================================================= // + +void integrateTemperature( + const pFlagTypeDevice& m, + const real dt, + deviceViewType1D T, + const deviceViewType1D& TR) +{ + auto r = m.activeRange(); + + Kokkos::parallel_for( + "integrateTemp", + policy(r.start(), r.end()), + KOKKOS_LAMBDA(uint32 i) + { + if (m(i)) + { + // --------------------------------------------------------- // + // Explicit Euler Integration for particle temperature. + // Stability is naturally enhanced because Q_rad is + // analytically linearized in the CFD solver, preventing + // the stiff T^4 equation from causing numerical explosions. + // --------------------------------------------------------- // + + T[i] += TR[i] * dt; + } + }); + + Kokkos::fence(); +} + +//+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +} // thermalSphereParticlesKernels +} // pFlow + + + diff --git a/src/Particles/SphereParticles/thermalSphereParticles/thermalSphereParticlesKernels.hpp b/src/Particles/SphereParticles/thermalSphereParticles/thermalSphereParticlesKernels.hpp new file mode 100644 index 000000000..7ace1ca8b --- /dev/null +++ b/src/Particles/SphereParticles/thermalSphereParticles/thermalSphereParticlesKernels.hpp @@ -0,0 +1,130 @@ +/*------------------------------- phasicFlow --------------------------------- + O C enter of + O O E ngineering and + O O M ultiscale modeling of + OOOOOOO F luid flow +------------------------------------------------------------------------------ + Copyright (C): www.cemf.ir + email: hamid.r.norouzi AT gmail.com +------------------------------------------------------------------------------ +Licence: + This file is part of phasicFlow code. It is a free software for simulating + granular and multiphase flows. You can redistribute it and/or modify it under + the terms of GNU General Public License v3 or any other later versions. + + phasicFlow is distributed to help others in their research in the field of + granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the + implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + +-----------------------------------------------------------------------------*/ + +#ifndef pFlow_thermalSphereParticlesKernels_hpp +#define pFlow_thermalSphereParticlesKernels_hpp + +#include "types.hpp" +#include "pointFlag.hpp" + +namespace pFlow +{ + +/** + * @namespace pFlow::thermalSphereParticlesKernels + * @brief High-performance GPU/CPU kernels for Lagrangian particle + * thermodynamics. + * + * @details + * This namespace isolates the raw numerical integration and physical + * calculations for particle heat transfer from the memory management classes. + * By using Kokkos `deviceViewType1D`, these functions map directly to + * massively parallel execution spaces (like NVIDIA CUDA or AMD HIP), ensuring + * that thermodynamic updates for millions of particles occur efficiently. + */ +namespace thermalSphereParticlesKernels +{ + + /** + * @brief Maps macroscopic shape properties to individual particle arrays. + * + * @param mask Active particle flag. + * @param shapeIndex The index linking a particle to its material type. + * @param Cp [OUT] Specific heat capacity array of particles. + * @param K [OUT] Thermal conductivity array of particles. + * @param emissivity [OUT] Surface emissivity array of particles. + * @param E0 [OUT] Real Young's Modulus array of particles. + * @param nu [OUT] Poisson's Ratio array of particles. + * @param shapeCp Dictionary-loaded heat capacities per material type. + * @param shapeK Dictionary-loaded conductivities per material type. + * @param shapeEps Dictionary-loaded emissivities per material type. + * @param shapeE0 Dictionary-loaded Real Young's Moduli per material. + * @param shapeNu Dictionary-loaded Poisson's Ratios per material type. + */ + void initThermalProperties( + const pFlagTypeDevice& mask, + const deviceViewType1D& shapeIndex, + deviceViewType1D Cp, + deviceViewType1D K, + deviceViewType1D emissivity, + deviceViewType1D E0, + deviceViewType1D nu, + const deviceViewType1D& shapeCp, + const deviceViewType1D& shapeK, + const deviceViewType1D& shapeEps, + const deviceViewType1D& shapeE0, + const deviceViewType1D& shapeNu); + + /** + * @brief Evaluates the First Law of Thermodynamics for each particle. + * + * @details + * Uses the Lumped Capacitance Model to calculate the temporal + * temperature derivative: + * dT/dt = (Q_conv + Q_rad + Q_pp + Q_pfp) / (m * Cp) + * + * @param mask Active particle flag. + * @param diameter Particle diameter array. + * @param mass Particle mass array. + * @param Cp Particle specific heat capacity array. + * @param temperature Current particle temperature array. + * @param Q_conv Convective heat source array [W]. + * @param Q_rad Radiative heat source array [W]. + * @param Q_pp Collisional heat transfer array [W]. + * @param Q_pfp Particle-Fluid-Particle sub-grid heat [W]. + * @param temperatureRate [OUT] The resulting rate of temperature change. + */ + void calcFluidParticleHeatTransfer( + const pFlagTypeDevice& mask, + const deviceViewType1D& diameter, + const deviceViewType1D& mass, + const deviceViewType1D& Cp, + const deviceViewType1D& temperature, + const deviceViewType1D& Q_conv, + const deviceViewType1D& Q_rad, + const deviceViewType1D& Q_pp, + const deviceViewType1D& Q_pfp, + deviceViewType1D temperatureRate); + + /** + * @brief Marches the particle temperatures forward in time. + * + * @details + * Applies Explicit Euler integration: + * T(t+dt) = T(t) + (dT/dt) * dt + * + * @param mask Active particle flag. + * @param dt The physical time step size [s]. + * @param temperature [IN/OUT] Particle temperature array to be updated. + * @param temperatureRate Computed rate of temperature change (dT/dt). + */ + void integrateTemperature( + const pFlagTypeDevice& mask, + const real dt, + deviceViewType1D temperature, + const deviceViewType1D& temperatureRate); + +} // thermalSphereParticlesKernels +} // pFlow + +#endif // pFlow_thermalSphereParticlesKernels_hpp + + + diff --git a/src/Particles/SphereParticles/thermalSphereShape/thermalSphereShape.cpp b/src/Particles/SphereParticles/thermalSphereShape/thermalSphereShape.cpp new file mode 100644 index 000000000..55b5ccd6e --- /dev/null +++ b/src/Particles/SphereParticles/thermalSphereShape/thermalSphereShape.cpp @@ -0,0 +1,141 @@ +/*------------------------------- phasicFlow --------------------------------- + O C enter of + O O E ngineering and + O O M ultiscale modeling of + OOOOOOO F luid flow +------------------------------------------------------------------------------ + Copyright (C): www.cemf.ir + email: hamid.r.norouzi AT gmail.com +------------------------------------------------------------------------------ +Licence: + This file is part of phasicFlow code. It is a free software for simulating + granular and multiphase flows. You can redistribute it and/or modify it under + the terms of GNU General Public License v3 or any other later versions. + + phasicFlow is distributed to help others in their research in the field of + granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the + implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + +-----------------------------------------------------------------------------*/ + +#include "thermalSphereShape.hpp" +#include "thermalProperty.hpp" + +namespace pFlow +{ + +//----------------------------- private methods ------------------------------- + +bool thermalSphereShape::readThermalProperties() +{ + // Retrieve the array that links each shape index to a specific material ID + auto pids = shapePropertyIds(); + + // Allocate memory for shape-specific thermal property vectors + cp_ = realVector("Cp", numShapes()); + k_ = realVector("k", numShapes()); + emissivity_ = realVector("emissivity", numShapes()); + E0_ = realVector("realYoungsModuli", numShapes()); + nu_ = realVector("poissonRatios", numShapes()); + + // ---------------------------------------------------------------------- // + // Direct access to the already-constructed base property object, which + // holds all dictionary values parsed from the case file. No separate + // fileDictionary instantiation or hardcoded path is needed here. + // ---------------------------------------------------------------------- // + const thermalProperty* tProps = + dynamic_cast(&properties()); + + if (!tProps) + { + fatalErrorInFunction + << "Provided property object is not a thermalProperty!" + << endl; + fatalExit; + } + + const realVector& allCp = tProps->heatCapacities(); + const realVector& allK = tProps->heatConductivities(); + const realVector& allEps = tProps->emissivities(); + const realVector& allE0 = tProps->realYoungsModuli(); + const realVector& allNu = tProps->poissonRatios(); + + // Map the global material properties to the specific local shapes + for (uint32 i = 0; i < numShapes(); ++i) + { + cp_[i] = allCp [pids[i]]; + k_[i] = allK [pids[i]]; + emissivity_[i] = allEps[pids[i]]; + E0_[i] = allE0 [pids[i]]; + nu_[i] = allNu [pids[i]]; + } + + // ---------------------------------------------------------------------- // + // Initial temperature assigned to newly inserted particles. + // + // A silently-defaulted value here (e.g. a stale 300 K) would mean any + // particle inserted through dynamic insertion (rather than read from + // an initial positions file) enters the domain at an unintended + // temperature with no warning at all, so this must be supplied + // explicitly rather than falling back to a built-in default. + // ---------------------------------------------------------------------- // + if (!properties().containsDataEntry("insertionTemperature")) + { + fatalErrorInFunction + << "Missing MANDATORY entry 'insertionTemperature' in the " + << "interaction dictionary." << endl; + fatalExit; + } + + insertionTemperature_ = properties().getVal("insertionTemperature"); + + return true; +} + +//---------------------------- protected methods ------------------------------ + +bool thermalSphereShape::writeToDict(dictionary& dict) const +{ + bool isWritten = sphereShape::writeToDict(dict) + && dict.add("heatCapacities", cp_) + && dict.add("heatConductivities", k_) + && dict.add("emissivities", emissivity_) + && dict.add("realYoungsModuli", E0_) + && dict.add("poissonRatios", nu_); + + return isWritten; +} + +//----------------------------- constructors ---------------------------------- + +// ========================================================================= // +// Constructors +// ========================================================================= // + +thermalSphereShape::thermalSphereShape( + const word& fileName, + repository* owner, + const property& prop) +: + sphereShape(fileName, owner, prop) +{ + readThermalProperties(); +} + +thermalSphereShape::thermalSphereShape( + const word& shapeType, + const word& fileName, + repository* owner, + const property& prop) +: + thermalSphereShape(fileName, owner, prop) +{ + // Body intentionally empty — delegates to primary constructor +} + +//+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +} // pFlow + + + diff --git a/src/Particles/SphereParticles/thermalSphereShape/thermalSphereShape.hpp b/src/Particles/SphereParticles/thermalSphereShape/thermalSphereShape.hpp new file mode 100644 index 000000000..749bd8ba7 --- /dev/null +++ b/src/Particles/SphereParticles/thermalSphereShape/thermalSphereShape.hpp @@ -0,0 +1,222 @@ +/*------------------------------- phasicFlow --------------------------------- + O C enter of + O O E ngineering and + O O M ultiscale modeling of + OOOOOOO F luid flow +------------------------------------------------------------------------------ + Copyright (C): www.cemf.ir + email: hamid.r.norouzi AT gmail.com +------------------------------------------------------------------------------ +Licence: + This file is part of phasicFlow code. It is a free software for simulating + granular and multiphase flows. You can redistribute it and/or modify it under + the terms of GNU General Public License v3 or any other later versions. + + phasicFlow is distributed to help others in their research in the field of + granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the + implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + +-----------------------------------------------------------------------------*/ + +#ifndef pFlow_thermalSphereShape_hpp +#define pFlow_thermalSphereShape_hpp + +#include "sphereShape.hpp" + +namespace pFlow +{ + +/** + * @class thermalSphereShape + * @brief Maps global material thermal properties to specific discrete + * particle shapes. + * + * @details + * While `thermalProperty` acts as a global database of material properties + * (e.g., "Steel", "Glass"), this class assigns those macroscopic properties + * to specific geometric entities (e.g., "Small_Steel", "Large_Glass"). + * It inherits from the mechanical `sphereShape` and adds thermodynamic data + * arrays (Cp, k, emissivity, E0, nu) tailored to the number of defined shapes + * in the simulation. + */ +class thermalSphereShape +: + public sphereShape +{ +public: + + //- Type info + + TypeInfo("shape"); + +private: + + //- private members + + /// @brief Heat capacity mapped to each specific particle shape + /// [J/(kg.K)]. + realVector cp_; + + /// @brief Thermal conductivity mapped to each specific particle shape + /// [W/(m.K)]. + realVector k_; + + /// @brief Surface emissivity mapped to each specific particle shape + /// (dimensionless). + realVector emissivity_; + + /// @brief Real Young's Modulus mapped to each specific particle shape + /// [Pa]. + realVector E0_; + + /// @brief Poisson's ratio mapped to each specific particle shape + /// (dimensionless). + realVector nu_; + + /** + * @brief Initial temperature [K] assigned to newly inserted particles. + * + * When a batch of particles is inserted and no existing particles are + * present to sample from, this value is used as the initial + * temperature. + * + * Read from the property dictionary key 'insertionTemperature', which + * is mandatory: the case file must state this value explicitly rather + * than relying on a built-in default. + */ + real insertionTemperature_ = real(300); + + //- private methods + + /** + * @brief Populates the shape-specific thermal arrays. + * @details Reads the master material properties from the simulation + * dictionary and maps them to the local shape arrays using the + * shape-to-material ID index. + * @return True if mapping is successful. + */ + bool readThermalProperties(); + +protected: + + //- protected methods + + /** + * @brief Serializes both mechanical and thermal shape data to a + * dictionary. + * @param dict The target phasicFlow dictionary object. + * @return True if all data is successfully written. + */ + bool writeToDict(dictionary& dict) const override; + +public: + + //- constructors + + thermalSphereShape( + const word& fileName, + repository* owner, + const property& prop); + + thermalSphereShape( + const word& shapeType, + const word& fileName, + repository* owner, + const property& prop); + + ~thermalSphereShape() override = default; + + //- public methods + + // ================================================================= // + // Accessor Methods (Vector Level) + // ================================================================= // + + inline + realVector heatCapacities() const + { + return cp_; + } + + inline + realVector heatConductivities() const + { + return k_; + } + + inline + realVector emissivities() const + { + return emissivity_; + } + + inline + realVector realYoungsModuli() const + { + return E0_; + } + + inline + realVector poissonRatios() const + { + return nu_; + } + + // ================================================================= // + // Accessor Methods (Scalar Level for specific shape indices) + // ================================================================= // + + inline + real heatCapacity(uint32 i) const + { + return cp_[i]; + } + + inline + real heatConductivity(uint32 i) const + { + return k_[i]; + } + + inline + real emissivity(uint32 i) const + { + return emissivity_[i]; + } + + inline + real realYoungsModulus(uint32 i) const + { + return E0_[i]; + } + + inline + real poissonRatio(uint32 i) const + { + return nu_[i]; + } + + /** + * @brief Initial temperature for newly inserted particles [K]. + * + * Used by thermalSphereParticles::insertParticles() when the + * temperature field is empty (first insertion event) and no existing + * particle can be sampled from. Configured via the mandatory + * 'insertionTemperature' key in the property dictionary. + */ + inline + real insertionTemperature() const + { + return insertionTemperature_; + } + + add_vCtor(shape, thermalSphereShape, word); + +}; // thermalSphereShape + +} // pFlow + +#endif // pFlow_thermalSphereShape_hpp + + + diff --git a/src/Property/thermalProperty/thermalProperty.cpp b/src/Property/thermalProperty/thermalProperty.cpp new file mode 100644 index 000000000..e7db25b82 --- /dev/null +++ b/src/Property/thermalProperty/thermalProperty.cpp @@ -0,0 +1,171 @@ +/*------------------------------- phasicFlow --------------------------------- + O C enter of + O O E ngineering and + O O M ultiscale modeling of + OOOOOOO F luid flow +------------------------------------------------------------------------------ + Copyright (C): www.cemf.ir + email: hamid.r.norouzi AT gmail.com +------------------------------------------------------------------------------ +Licence: + This file is part of phasicFlow code. It is a free software for simulating + granular and multiphase flows. You can redistribute it and/or modify it under + the terms of GNU General Public License v3 or any other later versions. + + phasicFlow is distributed to help others in their research in the field of + granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the + implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + +-----------------------------------------------------------------------------*/ + +#include "thermalProperty.hpp" +#include "fileDictionary.hpp" + +namespace pFlow +{ + +//----------------------------- private methods ------------------------------- + +// ========================================================================= // +// Section 1: Dictionary I/O +// ========================================================================= // + +bool thermalProperty::readDictionary() +{ + uniquePtr thermoDictPtr = nullptr; + + // Dynamic path resolution + if (p_dir_ != nullptr) + { + thermoDictPtr = makeUnique( + "thermoPhysicalInteraction", + *p_dir_); + } + else + { + // Safe fallback for legacy code calling the default constructor + thermoDictPtr = makeUnique( + "thermoPhysicalInteraction", + fileSystem("caseSetup")); + } + + auto& thermoDict = thermoDictPtr(); + + // Read thermal properties + heatCapacities_ = + thermoDict.getVal("heatCapacities"); + + heatConductivities_ = + thermoDict.getVal("heatConductivities"); + + emissivities_ = + thermoDict.getVal("emissivities"); + + // Read mechanical properties from the base interaction dictionary + realYoungsModuli_ = + this->getVal("realYoungsModuli"); + + poissonRatios_ = + this->getVal("poissonRatios"); + + // Validate sizes against the materials list + bool isValid = + (materials().size() == heatCapacities_.size()) && + (materials().size() == heatConductivities_.size()) && + (materials().size() == emissivities_.size()) && + (materials().size() == realYoungsModuli_.size()) && + (materials().size() == poissonRatios_.size()); + + if (!isValid) + { + fatalErrorInFunction + << " Mismatch in the number of material properties between " + << "'interaction' and 'thermoPhysicalInteraction' dictionaries." + << endl; + fatalExit; + } + + return isValid; +} + +bool thermalProperty::writeDictionary() +{ + bool isWritten = + add("heatCapacities", heatCapacities_) && + add("heatConductivities", heatConductivities_) && + add("emissivities", emissivities_) && + add("realYoungsModuli", realYoungsModuli_) && + add("poissonRatios", poissonRatios_); + + if (!isWritten) + { + fatalErrorInFunction + << " Error in writing thermal properties to dictionary " + << globalName() << endl; + } + + return isWritten; +} + +//----------------------------- constructors ---------------------------------- + +// ========================================================================= // +// Section 2: Constructors +// ========================================================================= // + +thermalProperty::thermalProperty( + const word& fileName, + repository* owner) +: + property(fileName, owner) +{ + if (!readDictionary()) + { + fatalExit; + } +} + +thermalProperty::thermalProperty( + const word& fileName, + const fileSystem& dir) +: + property(fileName, dir), + p_dir_(&dir) +{ + if (!readDictionary()) + { + fatalExit; + } +} + +thermalProperty::thermalProperty( + const word& fileName, + const wordVector& materials, + const realVector& densities, + const realVector& heatCapacities, + const realVector& heatConductivities, + const realVector& emissivities, + const realVector& realYoungsModuli, + const realVector& poissonRatios, + repository* owner) +: + property(fileName, materials, densities, owner), + heatCapacities_(heatCapacities), + heatConductivities_(heatConductivities), + emissivities_(emissivities), + realYoungsModuli_(realYoungsModuli), + poissonRatios_(poissonRatios) +{ + if (!writeDictionary()) + { + fatalExit; + } +} + +//+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +} // pFlow + + + + diff --git a/src/Property/thermalProperty/thermalProperty.hpp b/src/Property/thermalProperty/thermalProperty.hpp new file mode 100644 index 000000000..6e150c64b --- /dev/null +++ b/src/Property/thermalProperty/thermalProperty.hpp @@ -0,0 +1,183 @@ +/*------------------------------- phasicFlow --------------------------------- + O C enter of + O O E ngineering and + O O M ultiscale modeling of + OOOOOOO F luid flow +------------------------------------------------------------------------------ + Copyright (C): www.cemf.ir + email: hamid.r.norouzi AT gmail.com +------------------------------------------------------------------------------ +Licence: + This file is part of phasicFlow code. It is a free software for simulating + granular and multiphase flows. You can redistribute it and/or modify it under + the terms of GNU General Public License v3 or any other later versions. + + phasicFlow is distributed to help others in their research in the field of + granular and multiphase flows, but WITHOUT ANY WARRANTY; without even the + implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + +-----------------------------------------------------------------------------*/ + +#ifndef pFlow_thermalProperty_hpp +#define pFlow_thermalProperty_hpp + +#include "property.hpp" + +namespace pFlow +{ + +/** + * @class thermalProperty + * @brief Manages macroscopic thermal properties of materials in the simulation. + * + * @details + * Acts as the centralized database for thermodynamic material properties + * (e.g., heat capacity, conductivity, emissivity) loaded from the + * case dictionaries. It ensures that the mechanical and thermal + * definitions of materials remain strictly synchronized. + */ +class thermalProperty +: + public property +{ +public: + + //- Type info + + TypeInfo("thermalProperty"); + +private: + + //- private members + + // --- Section 2: Material Property Arrays --- + + realVector heatCapacities_; + + realVector heatConductivities_; + + realVector emissivities_; + + realVector realYoungsModuli_; + + realVector poissonRatios_; + + //- private methods + + // --- Section 3: File I/O --- + + bool readDictionary(); + + bool writeDictionary(); + +protected: + + //- protected members + + // --- Section 1: Internal Path Resolution --- + + /// @brief Safely caches the dictionary directory path. + const fileSystem* p_dir_ = nullptr; + +public: + + //- constructors + + // --- Section 4: Constructors --- + + explicit thermalProperty( + const word& fileName, + repository* owner = nullptr); + + thermalProperty( + const word& fileName, + const fileSystem& dir); + + thermalProperty( + const word& fileName, + const wordVector& materials, + const realVector& densities, + const realVector& heatCapacities, + const realVector& heatConductivities, + const realVector& emissivities, + const realVector& realYoungsModuli, + const realVector& poissonRatios, + repository* owner = nullptr); + + ~thermalProperty() override = default; + + //- public methods + + // --- Section 5: Vector Accessor Methods --- + + inline + const auto& heatCapacities() const + { + return heatCapacities_; + } + + inline + const auto& heatConductivities() const + { + return heatConductivities_; + } + + inline + const auto& emissivities() const + { + return emissivities_; + } + + inline + const auto& realYoungsModuli() const + { + return realYoungsModuli_; + } + + inline + const auto& poissonRatios() const + { + return poissonRatios_; + } + + // --- Section 6: Scalar Accessor Methods --- + + inline + real heatCapacity(uint32 i) const + { + return heatCapacities_[i]; + } + + inline + real heatConductivity(uint32 i) const + { + return heatConductivities_[i]; + } + + inline + real emissivity(uint32 i) const + { + return emissivities_[i]; + } + + inline + real realYoungsModulus(uint32 i) const + { + return realYoungsModuli_[i]; + } + + inline + real poissonRatio(uint32 i) const + { + return poissonRatios_[i]; + } + +}; // thermalProperty + +} // pFlow + +#endif // pFlow_thermalProperty_hpp + + + +