diff --git a/CMakeLists.txt b/CMakeLists.txt index e27377f5c..1588327ae 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -50,6 +50,11 @@ include(config/findSLEPc.cmake) include(config/printVersion.cmake) find_package(HDF5 REQUIRED) +add_library(PAPI INTERFACE IMPORTED GLOBAL) +target_include_directories(PAPI INTERFACE "/p/lustre2/kolosret/papi/src/install/include") +target_link_libraries(PAPI INTERFACE "/p/lustre2/kolosret/papi/src/install/lib/libpapi.so") + + # specify in the required libaries target_link_libraries(ablateLibrary PUBLIC @@ -65,6 +70,7 @@ target_link_libraries(ablateLibrary nlohmann_json::nlohmann_json ${HDF5_LIBRARIES} ZERORK::zerork_cfd_plugin +# PAPI PRIVATE chrestCompilerFlags) diff --git a/config/findZerork.cmake b/config/findZerork.cmake index 469564373..481c7a739 100644 --- a/config/findZerork.cmake +++ b/config/findZerork.cmake @@ -75,7 +75,9 @@ if (NOT (DEFINED ENV{ZERORK_DIR})) endif () elseif (DEFINED ENV{ZERORK_DIR}) - message(STATUS "Found ZERORK_DIR, using prebuilt zerork") + + if (NOT (DEFINED ENV{ABLATE_GPU})) + message(STATUS "Found ZERORK_DIR, using prebuilt CPU zerork") add_library(zerork_cfd_plugin INTERFACE IMPORTED GLOBAL) target_include_directories(zerork_cfd_plugin INTERFACE "$ENV{ZERORK_DIR}/include") @@ -83,5 +85,17 @@ elseif (DEFINED ENV{ZERORK_DIR}) add_library(ZERORK::zerork_cfd_plugin ALIAS zerork_cfd_plugin) + elseif (DEFINED ENV{ABLATE_GPU}) + message(STATUS "Found ZERORK_DIR, using prebuilt GPU zerork") + + add_library(zerork_cfd_plugin_gpu INTERFACE IMPORTED GLOBAL) + target_include_directories(zerork_cfd_plugin_gpu INTERFACE "$ENV{ZERORK_DIR}/include") + target_link_libraries(zerork_cfd_plugin_gpu INTERFACE "$ENV{ZERORK_DIR}/lib/libzerork_cfd_plugin_gpu.so") + + add_library(ZERORK::zerork_cfd_plugin ALIAS zerork_cfd_plugin_gpu) + + endif () + + endif () diff --git a/src/boundarySolver/boundarySolver.cpp b/src/boundarySolver/boundarySolver.cpp index 2e9d499eb..6c4492f78 100644 --- a/src/boundarySolver/boundarySolver.cpp +++ b/src/boundarySolver/boundarySolver.cpp @@ -847,6 +847,14 @@ PetscErrorCode ablate::boundarySolver::BoundarySolver::PreRHSFunction(TS ts, Pet PetscCall(rhsFunction.first(*this, ts, time, initialStage, locX, rhsFunction.second)); } EndEvent(); + + + try { + // update any aux fields, including ghost cells + UpdateAuxFields(time, locX, subDomain->GetAuxVector()); + } catch (std::exception& exception) { + SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in UpdateAuxFields: %s", exception.what()); + } PetscFunctionReturn(0); } diff --git a/src/boundarySolver/lodi/isothermalWall.cpp b/src/boundarySolver/lodi/isothermalWall.cpp index 83ec724ae..2065772cb 100644 --- a/src/boundarySolver/lodi/isothermalWall.cpp +++ b/src/boundarySolver/lodi/isothermalWall.cpp @@ -11,6 +11,8 @@ void ablate::boundarySolver::lodi::IsothermalWall::Setup(ablate::boundarySolver: ablate::boundarySolver::lodi::LODIBoundary::Setup(bSolver); bSolver.RegisterFunction(IsothermalWallFunction, this, fieldNames, fieldNames, {}); + bSolver.RegisterPreRHSFunction(CorrectBoundaryEnergy, this); + if (nSpecEqs) { bSolver.RegisterFunction( MirrorSpecies, this, {finiteVolume::CompressibleFlowFields::EULER_FIELD, finiteVolume::CompressibleFlowFields::DENSITY_YI_FIELD}, {finiteVolume::CompressibleFlowFields::YI_FIELD}); @@ -37,14 +39,18 @@ PetscErrorCode ablate::boundarySolver::lodi::IsothermalWall::IsothermalWallFunct PetscReal boundarySpeedOfSound; PetscReal boundaryPressure; +// boundaryTemperature=300; + // Get the velocity and pressure on the surface { boundaryDensity = boundaryValues[uOff[isothermalWall->eulerId] + finiteVolume::CompressibleFlowFields::RHO]; for (PetscInt d = 0; d < dim; d++) { boundaryVel[d] = boundaryValues[uOff[isothermalWall->eulerId] + finiteVolume::CompressibleFlowFields::RHOU + d] / boundaryDensity; boundaryNormalVelocity += boundaryVel[d] * fg->normal[d]; + boundaryNormalVelocity += boundaryVel[d] * fg->normal[d]; } PetscCall(isothermalWall->computeTemperature.function(boundaryValues, &boundaryTemperature, isothermalWall->computeTemperature.context.get())); +// boundaryTemperature=300; PetscCall(isothermalWall->computeSpeedOfSound.function(boundaryValues, boundaryTemperature, &boundarySpeedOfSound, isothermalWall->computeSpeedOfSound.context.get())); PetscCall(isothermalWall->computePressureFromTemperature.function(boundaryValues, boundaryTemperature, &boundaryPressure, isothermalWall->computePressureFromTemperature.context.get())); } @@ -99,6 +105,7 @@ PetscErrorCode ablate::boundarySolver::lodi::IsothermalWall::IsothermalWallFunct // Get scriptL std::vector scriptL(isothermalWall->nEqs); scriptL[1 + dim] = lambda[1 + dim] * (dPdNorm - boundaryDensity * dVeldNorm * alpha2 * (velNormPrim - boundaryNormalVelocity - speedOfSoundPrim)); // Outgoing + scriptL[1 + dim] = 0; //L3 // acoustic // wave scriptL[0] = scriptL[1 + dim]; // Incoming acoustic wave @@ -135,9 +142,115 @@ PetscErrorCode ablate::boundarySolver::lodi::IsothermalWall::IsothermalWallFunct scriptL.data(), transformationMatrix, source); +// for (PetscInt d = 0; d < (2+dim+isothermalWall->nSpecEqs); d++) { +// source[d]=0; +// } + + PetscFunctionReturn(0); +} +PetscErrorCode ablate::boundarySolver::lodi::IsothermalWall::CorrectBoundaryEnergy( + ablate::boundarySolver::BoundarySolver& solver, + TS ts, PetscReal time, bool initialStage, Vec locX, void* ctx) { + + PetscFunctionBeginUser; + + auto isothermalWall = reinterpret_cast(ctx); + auto& subDomain = solver.GetSubDomain(); + + // Get field information + const auto& eulerField = subDomain.GetField(finiteVolume::CompressibleFlowFields::EULER_FIELD); + const auto& densityYiField = subDomain.GetField(finiteVolume::CompressibleFlowFields::DENSITY_YI_FIELD); + + // Get the DM and dimension + DM dm = subDomain.GetDM(); + PetscInt dim = subDomain.GetDimensions(); + + // Get array access to the local vector + PetscScalar* locXArray; + PetscCall(VecGetArray(locX, &locXArray)); + + // March over each boundary cell in the solver region + ablate::domain::Range cellRange; + solver.GetCellRange(cellRange); + + for (PetscInt c = cellRange.start; c < cellRange.end; ++c) { + PetscInt boundaryCell = cellRange.points ? cellRange.points[c] : c; + + // Get pointer to the boundary cell data + PetscScalar* cellData; + PetscCall(DMPlexPointLocalRef(dm, boundaryCell, locXArray, &cellData)); + + if (cellData) { + // Extract conserved variables + PetscReal rho_old = cellData[eulerField.offset + finiteVolume::CompressibleFlowFields::RHO]; +// PetscReal rhoE = cellData[eulerField.offset + finiteVolume::CompressibleFlowFields::RHOE]; + + if (rho_old<0){ + std::cout << "The density is negative for cell: " << boundaryCell << " \n"; + } + + PetscReal rho = PetscMax(0.005, rho_old); + rho = PetscMin(100, rho); + // TODO something smarter maybe? ... + + //Reset the density + cellData[eulerField.offset + finiteVolume::CompressibleFlowFields::RHO] = rho; +// cellData[eulerField.offset + fp::RHO] = rho; + //Reset the momentum + PetscReal mom; + for (PetscInt d = 0; d < dim; d++) { + mom=cellData[eulerField.offset + finiteVolume::CompressibleFlowFields::RHOU+d]; + cellData[eulerField.offset + finiteVolume::CompressibleFlowFields::RHOU+d] = mom / rho_old * rho; + } + + // Reset the Species + PetscReal yi; + for (PetscInt sp = 0; sp < isothermalWall->nSpecEqs; sp++) { + yi=cellData[densityYiField.offset+sp]; + cellData[densityYiField.offset+sp] = yi / rho_old * rho; + + } + + + // Calculate kinetic energy + PetscReal KE = 0.0; + for (PetscInt d = 0; d < dim; ++d) { + PetscReal momentum_d = cellData[eulerField.offset + finiteVolume::CompressibleFlowFields::RHOU + d]; + KE += momentum_d * momentum_d; + } + KE = 0.5 * KE / rho; + + // Sensible energy +// PetscReal e_current = rhoE / rho_old - KE; + + + PetscReal e_tmp; + PetscCall(isothermalWall->computeInternalEnergyFromTemperature.function(cellData + eulerField.offset, + isothermalWall->wallTemperature, &e_tmp, isothermalWall->computeInternalEnergyFromTemperature.context.get())); + +// PetscReal e_corrected = PetscMax(e_tmp, e_current); +// PetscReal e_corrected_total = e_corrected+ KE; + + e_tmp = e_tmp + KE; + + cellData[eulerField.offset + finiteVolume::CompressibleFlowFields::RHOE] = rho * (e_tmp); + + + + + + } + } + + // Don't forget to restore the range + solver.RestoreRange(cellRange); + + PetscCall(VecRestoreArray(locX, &locXArray)); PetscFunctionReturn(0); } + + PetscErrorCode ablate::boundarySolver::lodi::IsothermalWall::MirrorSpecies(PetscInt dim, const ablate::boundarySolver::BoundarySolver::BoundaryFVFaceGeom *fg, const PetscFVCellGeom *boundaryCell, const PetscInt *uOff, PetscScalar *boundaryValues, const PetscScalar *stencilValues, const PetscInt *aOff, PetscScalar *auxValues, const PetscScalar *stencilAuxValues, void *ctx) { @@ -155,6 +268,11 @@ PetscErrorCode ablate::boundarySolver::lodi::IsothermalWall::MirrorSpecies(Petsc boundaryValues[uOff[DENSITY_YI] + sp] = yi * boundaryDensity; auxValues[aOff[YI] + sp] = yi; } +// boundaryValues[uOff[EULER_FIELD] + RHO] = stencilValues[uOff[EULER_FIELD] + RHO]; +// boundaryValues[uOff[EULER_FIELD] + RHO+1] = stencilValues[uOff[EULER_FIELD] + RHO+1]; +// boundaryValues[uOff[EULER_FIELD] + RHO+2] = stencilValues[uOff[EULER_FIELD] + RHO+2]; +// boundaryValues[uOff[EULER_FIELD] + RHO+3] = stencilValues[uOff[EULER_FIELD] + RHO+3]; +// boundaryValues[uOff[EULER_FIELD] + RHO+4] = stencilValues[uOff[EULER_FIELD] + RHO+4]; PetscFunctionReturn(0); } diff --git a/src/boundarySolver/lodi/isothermalWall.hpp b/src/boundarySolver/lodi/isothermalWall.hpp index f156a9251..f908cf1c4 100644 --- a/src/boundarySolver/lodi/isothermalWall.hpp +++ b/src/boundarySolver/lodi/isothermalWall.hpp @@ -9,6 +9,9 @@ class IsothermalWall : public LODIBoundary { explicit IsothermalWall(std::shared_ptr eos, std::shared_ptr pressureGradientScaling = {}); void Setup(ablate::boundarySolver::BoundarySolver& bSolver) override; + static PetscErrorCode CorrectBoundaryEnergy(ablate::boundarySolver::BoundarySolver& solver, + TS ts, PetscReal time, bool initialStage, + Vec locX, void* ctx); static PetscErrorCode IsothermalWallFunction(PetscInt dim, const boundarySolver::BoundarySolver::BoundaryFVFaceGeom* fg, const PetscFVCellGeom* boundaryCell, const PetscInt uOff[], const PetscScalar* boundaryValues, const PetscScalar* stencilValues[], const PetscInt aOff[], const PetscScalar* auxValues, @@ -16,6 +19,7 @@ class IsothermalWall : public LODIBoundary { PetscScalar source[], void* ctx); private: + double wallTemperature = 300; static PetscErrorCode MirrorSpecies(PetscInt dim, const BoundarySolver::BoundaryFVFaceGeom* fg, const PetscFVCellGeom* boundaryCell, const PetscInt uOff[], PetscScalar* boundaryValues, const PetscScalar* stencilValues, const PetscInt aOff[], PetscScalar* auxValues, const PetscScalar* stencilAuxValues, void* ctx); }; diff --git a/src/boundarySolver/lodi/lodiBoundary.cpp b/src/boundarySolver/lodi/lodiBoundary.cpp index 06598b8fe..c0e135a4c 100644 --- a/src/boundarySolver/lodi/lodiBoundary.cpp +++ b/src/boundarySolver/lodi/lodiBoundary.cpp @@ -231,5 +231,6 @@ void ablate::boundarySolver::lodi::LODIBoundary::Setup(PetscInt dimsIn, PetscInt computeSpecificHeatConstantPressure = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::SpecificHeatConstantPressure, fields); computeSpecificHeatConstantVolume = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::SpecificHeatConstantVolume, fields); computeSensibleEnthalpyFunction = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::SensibleEnthalpy, fields); + computeInternalEnergyFromTemperature = eos->GetThermodynamicTemperatureFunction(eos::ThermodynamicProperty::InternalSensibleEnergy, fields); computePressure = eos->GetThermodynamicFunction(eos::ThermodynamicProperty::Pressure, fields); } \ No newline at end of file diff --git a/src/boundarySolver/lodi/lodiBoundary.hpp b/src/boundarySolver/lodi/lodiBoundary.hpp index 8a450c180..8b6e7b468 100644 --- a/src/boundarySolver/lodi/lodiBoundary.hpp +++ b/src/boundarySolver/lodi/lodiBoundary.hpp @@ -50,6 +50,7 @@ class LODIBoundary : public BoundaryProcess { eos::ThermodynamicTemperatureFunction computeSpecificHeatConstantVolume; eos::ThermodynamicTemperatureFunction computeSensibleEnthalpyFunction; eos::ThermodynamicFunction computePressure; + eos::ThermodynamicTemperatureFunction computeInternalEnergyFromTemperature; public: explicit LODIBoundary(std::shared_ptr eos, std::shared_ptr pressureGradientScaling = {}); diff --git a/src/eos/zerork.cpp b/src/eos/zerork.cpp index 4b639dba6..a20d1d254 100644 --- a/src/eos/zerork.cpp +++ b/src/eos/zerork.cpp @@ -129,15 +129,14 @@ PetscErrorCode ablate::eos::zerorkEOS::TemperatureTemperatureFunction(const Pets // compute the internal energy needed to compute temperature from the sensible enthalpy double sensibleenergy = conserved[functionContext->eulerOffset + ablate::finiteVolume::CompressibleFlowFields::RHOE] / density - 0.5 * speedSquare; - double enthalpymix = functionContext->mech->getMassEnthalpyFromTY(298.15, &reactorMassFrac[0]); - sensibleenergy += enthalpymix; + double energymix = functionContext->mech->getMassIntEnergyFromTY(298.15, &reactorMassFrac[0]); + sensibleenergy += energymix; // set the temperature from zerork *temperature = functionContext->mech->getTemperatureFromEY(sensibleenergy, &reactorMassFrac[0], temperatureGuess); PetscFunctionReturn(0); - PetscFunctionReturn(0); } PetscErrorCode ablate::eos::zerorkEOS::TemperatureMassFractionFunction(const PetscReal *conserved, const PetscReal *yi, PetscReal *property, void *ctx) { return TemperatureTemperatureMassFractionFunction(conserved, yi, 300, property, ctx); @@ -160,8 +159,8 @@ PetscErrorCode ablate::eos::zerorkEOS::TemperatureTemperatureMassFractionFunctio // compute the internal energy needed to compute temperature from the sensible enthalpy double sensibleenergy = conserved[functionContext->eulerOffset + ablate::finiteVolume::CompressibleFlowFields::RHOE] / density - 0.5 * speedSquare; - double enthalpyMixFormation = functionContext->mech->getMassEnthalpyFromTY(298.15, &reactorMassFrac[0]); - sensibleenergy += enthalpyMixFormation; + double energymix = functionContext->mech->getMassIntEnergyFromTY(298.15, &reactorMassFrac[0]); + sensibleenergy += energymix; // set the temperature from zerork *temperature = functionContext->mech->getTemperatureFromEY(sensibleenergy, &reactorMassFrac[0], temperatureGuess); @@ -196,9 +195,9 @@ PetscErrorCode ablate::eos::zerorkEOS::InternalSensibleEnergyTemperatureFunction FillreactorMassFracVectorFromDensityMassFractions(numSpc, density, conserved + functionContext->densityYiOffset, reactorMassFrac); double energyMix = functionContext->mech->getMassIntEnergyFromTY(temperature, &reactorMassFrac[0]); - double enthalpyMixFormation = functionContext->mech->getMassEnthalpyFromTY(298.15, &reactorMassFrac[0]); + double uMixFormation = functionContext->mech->getMassIntEnergyFromTY(298.15, &reactorMassFrac[0]); - *sensibleEnergyTemperature = energyMix - enthalpyMixFormation; + *sensibleEnergyTemperature = energyMix - uMixFormation; PetscFunctionReturn(0); } @@ -220,9 +219,9 @@ PetscErrorCode ablate::eos::zerorkEOS::InternalSensibleEnergyTemperatureMassFrac FillreactorMassFracVectorFromMassFractions(numSpc, yi, reactorMassFrac); double energyMix = functionContext->mech->getMassIntEnergyFromTY(temperature, &reactorMassFrac[0]); - double enthalpyMixFormation = functionContext->mech->getMassEnthalpyFromTY(298.15, &reactorMassFrac[0]); + double uMixFormation = functionContext->mech->getMassIntEnergyFromTY(298.15, &reactorMassFrac[0]); - *sensibleEnergyTemperature = energyMix - enthalpyMixFormation; + *sensibleEnergyTemperature = energyMix - uMixFormation; PetscFunctionReturn(0); } @@ -612,9 +611,9 @@ ablate::eos::EOSFunction ablate::eos::zerorkEOS::GetFieldFunctionFunction(const PetscReal density = mech->getDensityFromTPY(temperature, pressure, &reactorMassFrac[0]); double energyMix = mech->getMassIntEnergyFromTY(temperature, &reactorMassFrac[0]); - double enthalpyMixFormation = mech->getMassEnthalpyFromTY(298.15, &reactorMassFrac[0]); + double energyMixFormation = mech->getMassIntEnergyFromTY(298.15, &reactorMassFrac[0]); - double sensibleInternalEnergy = energyMix - enthalpyMixFormation; + double sensibleInternalEnergy = energyMix - energyMixFormation; // convert to total sensibleEnergy PetscReal kineticEnergy = 0; @@ -643,9 +642,9 @@ ablate::eos::EOSFunction ablate::eos::zerorkEOS::GetFieldFunctionFunction(const std::vector reactorMassFrac(nSpc, 0.); FillreactorMassFracVectorFromMassFractions(nSpc, yi, reactorMassFrac); - double enthalpyMixFormation = mech->getMassEnthalpyFromTY(298.15, &reactorMassFrac[0]); + double energyMixFormation = mech->getMassIntEnergyFromTY(298.15, &reactorMassFrac[0]); - double internalEnergy = sensibleInternalEnergy + enthalpyMixFormation; + double internalEnergy = sensibleInternalEnergy + energyMixFormation; double temperature = mech->getTemperatureFromEY(internalEnergy, &reactorMassFrac[0], 300); PetscReal density = mech->getDensityFromTPY(temperature, pressure, &reactorMassFrac[0]); @@ -703,9 +702,9 @@ ablate::eos::EOSFunction ablate::eos::zerorkEOS::GetFieldFunctionFunction(const std::vector reactorMassFrac(nSpc, 0.); FillreactorMassFracVectorFromMassFractions(nSpc, yi, reactorMassFrac); - double enthalpyMixFormation = mech->getMassEnthalpyFromTY(298.15, &reactorMassFrac[0]); + double energyMixFormation = mech->getMassIntEnergyFromTY(298.15, &reactorMassFrac[0]); - double internalEnergy = sensibleInternalEnergy + enthalpyMixFormation; + double internalEnergy = sensibleInternalEnergy + energyMixFormation; double temperature = mech->getTemperatureFromEY(internalEnergy, &reactorMassFrac[0], 300); PetscReal density = mech->getDensityFromTPY(temperature, pressure, &reactorMassFrac[0]); diff --git a/src/eos/zerork/sourceCalculatorZeroRK.cpp b/src/eos/zerork/sourceCalculatorZeroRK.cpp index 98de51c63..a5d76b464 100644 --- a/src/eos/zerork/sourceCalculatorZeroRK.cpp +++ b/src/eos/zerork/sourceCalculatorZeroRK.cpp @@ -6,6 +6,12 @@ #include "utilities/mpiUtilities.hpp" #include "utilities/stringUtilities.hpp" +#include +#include "domain/range.hpp" +#include +#include + + void ablate::eos::zerorkeos::SourceCalculator::ChemistryConstraints::Set(const std::shared_ptr& options) { if (options) { verbose = options->Get("verbose", verbose); @@ -24,6 +30,8 @@ void ablate::eos::zerorkeos::SourceCalculator::ChemistryConstraints::Set(const s errorhandle = options->Get("errorhandle", errorhandle); dumpreactor = options->Get("dumpreactor", dumpreactor); dumpfailed = options->Get("dumpfailed", dumpfailed); + n_reactors_max = options->Get("n_reactors_max", n_reactors_max); + n_reactors_min = options->Get("n_reactors_min", n_reactors_min); cvode_num_retries = options->Get("cvode_num_retries", cvode_num_retries); cvode_retry_absolute_tolerance_adjustment = options->Get("cvode_retry_absolute_tolerance_adjustment", cvode_retry_absolute_tolerance_adjustment); cvode_retry_relative_tolerance_adjustment = options->Get("cvode_retry_relative_tolerance_adjustment", cvode_retry_relative_tolerance_adjustment); @@ -39,6 +47,51 @@ ablate::eos::zerorkeos::SourceCalculator::SourceCalculator(const std::vector(numberCells * (eosIn->mech->getNumSpecies() + 1)); + auto chemEos = std::dynamic_pointer_cast(eos); + + auto speciesElementInformation = chemEos->GetSpeciesElementalInformation(); + auto elementInformation = chemEos->GetElementInformation(); // element -> atomic mass + auto speciesMolecularMass = chemEos->GetSpeciesMolecularMass(); // species -> MW + const auto& species = chemEos->GetSpeciesVariables(); + + const std::vector trackingElements{"C", "H"}; + std::map massFractionsFuel; + std::map massFractionsOxidizer; + for (const auto& spName : species) { + massFractionsFuel[spName] = 0.0; + massFractionsOxidizer[spName] = 0.0; + } + + const std::string fuelName = "MMA"; + const std::string oxidizerName = "O2"; + + if (!massFractionsFuel.count(fuelName) || !massFractionsOxidizer.count(oxidizerName)) { + throw std::invalid_argument("SourceCalculator mixture fraction: could not find PMMA or O2 species in mechanism."); + } + + massFractionsFuel[fuelName] = 1.0; + massFractionsOxidizer[oxidizerName] = 1.0; + + zMixCoefficients.assign(species.size(), 0.0); + for (std::size_t s = 0; s < species.size(); ++s) { + const auto& spName = species[s]; + auto& speciesElement = speciesElementInformation[spName]; + for (const auto& elem : trackingElements) { + zMixCoefficients[s] += speciesElement[elem] * elementInformation[elem] / speciesMolecularMass[spName]; + } + } + + + zMixFuel = 0.0; + zMixOxidizer = 0.0; + for (std::size_t s = 0; s < species.size(); ++s) { + const auto& spName = species[s]; + zMixFuel += zMixCoefficients[s] * massFractionsFuel[spName]; + zMixOxidizer += zMixCoefficients[s] * massFractionsOxidizer[spName]; + } + + + // Look for the euler field auto eulerField = std::find_if(fields.begin(), fields.end(), [](const auto& field) { return field.name == ablate::finiteVolume::CompressibleFlowFields::EULER_FIELD; }); if (eulerField == fields.end()) { @@ -65,6 +118,12 @@ ablate::eos::zerorkeos::SourceCalculator::SourceCalculator(const std::vector> utilities::PetscUtilities::checkError; + std::size_t nCells = cellRange.end - cellRange.start; + reactorzMix.resize(nCells, 0.0); + reactorzMixGrad.resize(nCells * dim, 0.0); + // zerork state load up int nSpc = eos->mech->getNumSpecies(); // Number of Species int nState = nSpc + 1; @@ -168,21 +234,26 @@ void ablate::eos::zerorkeos::SourceCalculator::ComputeSource(const ablate::domai std::vector reactorP(nReactors); std::vector density2(nReactors); std::vector sensibleenergy(nReactors); - std::vector velmag2(nReactors); +// std::vector velmag2(nReactors); std::vector reactorMassFrac(nReactors * nSpc); std::vector enthalpyOfFormation(nSpc); std::vector reactorEval(nReactors, 0); + std::vector reactoriDs(nReactors); + std::vector reactorLocalIdx(nReactors, -1); + // Set up the vectors that are actually evaluated with the temperature threshold std::vector reactorTEval(nReactors); std::vector reactorPEval(nReactors); std::vector reactorMassFracEval(nReactors * nSpc); // get the current state from petsc - int p = 0; +// std::vector cellToLocal(nCells, -1); + int p = 0; //active reactor index for (int i = cellRange.start; i < cellRange.end; ++i) { const PetscInt cell = cellRange.points ? cellRange.points[i] : i; const std::size_t k = i - cellRange.start; +// cellToLocal[cell] = k; const PetscScalar* eulerField = nullptr; DMPlexPointLocalFieldRead(solutionDm, cell, eulerId, flowArray, &eulerField) >> utilities::PetscUtilities::checkError; @@ -225,24 +296,289 @@ void ablate::eos::zerorkeos::SourceCalculator::ComputeSource(const ablate::domai reactorT[k] = eos->mech->getTemperatureFromEY(sensibleenergy[k], &reactorMassFrac[k * nSpc], 2000); reactorP[k] = eos->mech->getPressureFromTVY(reactorT[k], 1 / density, &reactorMassFrac[k * nSpc]); + reactorzMix[k]=ComputeMixtureFraction(&reactorMassFrac[k * nSpc]); // Set up the vector that is actually being solved if (reactorT[k] > chemistryConstraints.thresholdTemperature) { reactorEval[k] = 1; reactorTEval[p] = reactorT[k]; reactorPEval[p] = reactorP[k]; + reactoriDs[p]=cell; + reactorLocalIdx[p] = static_cast(k); for (int s = 0; s < nSpc; s++) { reactorMassFracEval[p * nSpc + s] = reactorMassFrac[k * nSpc + s]; } +// reactorzMixEval[p]=ComputeMixtureFraction(&reactorMassFracEval[p * nSpc]); p += 1; } } + ComputeMixtureFractionGradients(solutionDm, cellRange, reactorzMix, reactorzMixGrad, dim); +// double maxelem = *std::max_element(reactorzMixGrad.begin(),reactorzMixGrad.end()); +// std::cout << "The max gradient is: "<< maxelem < reactorzMixGradMag(p, 0.0); + std::vector reactorChiEval(p, -1.0); + + int OHind=12; //Eventually this shouldnt be hardcoded + double OHcutoff=0.001; + double Xin_lo=-3.; + double Xin_high=1.; + double Xout_low=0.0; + double Xout_high=1.0; + + for (int i = 0; i= static_cast(nCells)) { + throw std::runtime_error("ComputeSource: invalid local cell index for reactor."); + } + diff = SutherlandDiff(reactorTEval[i], density2[k]); + //Calculate the square of the gradient +// diff= SutherlandDiff(reactorTEval[i],density2[reactoriDs[i]]); +// for (int s = 0; s < dim; s++) { +// reactorzMixGradMag[i] += reactorzMixGrad[i * dim + s]*reactorzMixGrad[i * dim + s]; +// } + for (int s = 0; s < dim; s++) { + double g = reactorzMixGrad[static_cast(k) * dim + s]; + reactorzMixGradMag[i] += g * g; + } + + if (weight && atoi(weight) == 1) { + reactorChiEval[i] = 2 * diff * reactorzMixGradMag[i]; + }else if (weight && atoi(weight) == 2) { + double chi = 2 * diff * reactorzMixGradMag[i]; + reactorChiEval[i] = chi*chi; + }else if (weight && atoi(weight) == 3) { + reactorChiEval[i] = 2 * diff * reactorzMixGradMag[i]; + }else if (weight && atoi(weight) == 4) { + double chi =2 * diff * reactorzMixGradMag[i]; + chi=std::max(chi,0.0001); + //Using the fitting constants + double u = std::log10(chi); + double v = p0 + p1 * u + p2 * u * u; + reactorChiEval[i] = std::pow(10.0, v); + }else if (weight && atoi(weight) == 5) { + double chi =2 * diff * reactorzMixGradMag[i]; + chi=std::max(chi,0.0001); + //Using the fitting constants + double u = std::log10(chi); + double v = pt0 + pt1 * u + pt2 * u * u; + reactorChiEval[i] = std::pow(10.0, v); + }else if (weight && atoi(weight) == 6) { + + double chi = 2 * diff * reactorzMixGradMag[i]; + //limit chi + chi=std::max(chi,0.0001); + + double Yi_OH = reactorMassFracEval[i * nSpc + OHind]; + //Using the fitting constants + double u = std::log10(chi); + double logchi = u; + if (u >= Xin_lo && u <= Xin_high && Yi_OH > OHcutoff){ + logchi= Xout_low + (u - Xin_lo)*(Xout_high - Xout_low)/(Xin_high - Xin_lo); +// std::cout<< "Chi was adjusted from: "<= Xin_lo && u <= Xin_high && Yi_OH > OHcutoff){ + logchi= Xout_low + (u - Xin_lo)*(Xout_high - Xout_low)/(Xin_high - Xin_lo); + // std::cout<< "Chi was adjusted from: "< 0) { +// double maxChi = *std::max_element(reactorChiEval.begin(),reactorChiEval.end()); +// double minChi = *std::min_element(reactorChiEval.begin(),reactorChiEval.end()); +// std::cout << "The max chi is: "<< maxChi < reactorWeights(N); + // save Yi for source terms std::vector ys = reactorMassFracEval; + if (weight && atoi(weight) == 1) { + std::cout <<"Setting Chi as the weight"< 0) { + zerork_reactor_set_aux_field_pointer(ZERORK_FIELD_COST, &reactorChiEval[0], zrm_handle); + } + } + else if(weight && atoi(weight) == 2){ + std::cout <<"Using option 2 for the weights"< 0) { + zerork_reactor_set_aux_field_pointer(ZERORK_FIELD_COST, &reactorChiEval[0], zrm_handle); + } + }else if(weight && atoi(weight) == 3){ + std::cout <<"Using option 3 for the weights"< idx(N); + std::iota(idx.begin(), idx.end(), 0); + + // Sort indices by descending Chi + std::sort(idx.begin(), idx.end(), [&](int a, int b) { + return reactorChiEval[a] > reactorChiEval[b]; + }); + + // Max Chi value + double chi_max = reactorChiEval[idx[0]]; + + // Assign MAX Chi to top-N hardest (batch size) + int n = std::min(chemistryConstraints.n_reactors_max, N); + for(int i = 0; i < n; ++i){ + int id = idx[i]; + reactorWeights[id] = chi_max; + } + if (N > 0) { + zerork_reactor_set_aux_field_pointer(ZERORK_FIELD_COST, &reactorWeights[0], zrm_handle); + } + + }else if(weight && atoi(weight) == 4){ + std::cout <<"Using option 4 for the weights"< 0) { + zerork_reactor_set_aux_field_pointer(ZERORK_FIELD_COST, &reactorChiEval[0], zrm_handle); + } + }else if(weight && atoi(weight) == 5){ + std::cout <<"Using option 5 for the weights"< 0) { + double minVal = *std::min_element(reactorChiEval.begin(), reactorChiEval.end()); + minVal = std::max(minVal,0.00001); + for (auto& v : reactorChiEval) { + v /= minVal; + } + zerork_reactor_set_aux_field_pointer(ZERORK_FIELD_COST, &reactorChiEval[0], zrm_handle); + } + }else if(weight && atoi(weight) == 6){ + std::cout <<"Using option 6 for the weights"< 0) { + zerork_reactor_set_aux_field_pointer(ZERORK_FIELD_COST, &reactorChiEval[0], zrm_handle); + } + }else if(weight && atoi(weight) == 7){ + std::cout <<"Using option 6 for the weights"< 0) { + zerork_reactor_set_aux_field_pointer(ZERORK_FIELD_COST, &reactorChiEval[0], zrm_handle); + } + }else if(weight && atoi(weight) == 8){ + std::cout <<"Using option 5 for the weights"< 0) { + double minVal = *std::min_element(reactorChiEval.begin(), reactorChiEval.end()); + minVal = std::max(minVal,0.00001); + for (auto& v : reactorChiEval) { + v /= minVal; + } + zerork_reactor_set_aux_field_pointer(ZERORK_FIELD_COST, &reactorChiEval[0], zrm_handle); + } + }else if(weight && atoi(weight) == 9){ + std::cout <<"Using option 5 for the weights"< 0) { + double minVal = *std::min_element(reactorChiEval.begin(), reactorChiEval.end()); + minVal = std::max(minVal,0.00001); + for (auto& v : reactorChiEval) { + v /= minVal; + } + zerork_reactor_set_aux_field_pointer(ZERORK_FIELD_COST, &reactorChiEval[0], zrm_handle); + } + } + else{ + std::cout <<"No weigths are set"<> utilities::PetscUtilities::checkError; EndEvent(); } void ablate::eos::zerorkeos::SourceCalculator::AddSource(const ablate::domain::Range& cellRange, Vec, Vec locFVec) { @@ -372,6 +709,16 @@ void ablate::eos::zerorkeos::SourceCalculator::AddSource(const ablate::domain::R EndEvent(); } +double ablate::eos::zerorkeos::SourceCalculator::ComputeMixtureFraction(const double* yi) const { + double zMix = 0.0; + for (std::size_t s = 0; s < zMixCoefficients.size(); ++s) { + zMix += zMixCoefficients[s] * yi[s]; + } + if (zMix < 0.0) zMix = 0.0; + if (zMix > 1.0) zMix = 1.0; + return (zMix - zMixOxidizer) / (zMixFuel - zMixOxidizer); +} + std::ostream& ablate::eos::zerorkeos::operator<<(std::ostream& os, const ablate::eos::zerorkeos::SourceCalculator::ReactorType& v) { switch (v) { case ablate::eos::zerorkeos::SourceCalculator::ReactorType::ConstantPressure: @@ -402,4 +749,161 @@ std::istream& ablate::eos::zerorkeos::operator>>(std::istream& is, ablate::eos:: " Default is Contstant volume"); } return is; -} \ No newline at end of file +} + + +double ablate::eos::zerorkeos::SourceCalculator::SutherlandDiff(double T,double rho){ + double Diff=0; + double mu0=1.716e-5; // Pa·s at T0 (air) + double T0=273.15; // K + double S=110.4; // Sutherland constant for air [K] + double Sc=0.7; // Schmidt number (approx for air) + + double mu = mu0 * std::pow((T / T0), 1.5) * (T0 + S) / (T + S); + return Diff = mu / (rho * Sc); +} + + + + + + + + + +void ablate::eos::zerorkeos::SourceCalculator::ComputeMixtureFractionGradients( + DM dm, + const ablate::domain::Range& cellRange, + const std::vector& zMixValues, + std::vector& zMixGrad, + PetscInt dim) { + + // number of local cells in this range + const PetscInt nCells = cellRange.end - cellRange.start; + if ((PetscInt)zMixValues.size() != nCells) { + throw std::runtime_error("ComputeMixtureFractionGradients: zMixValues size must equal number of cells in cellRange."); + } + + // get spatial dimension if dim <= 0 + PetscInt dmDim; + DMGetDimension(dm, &dmDim) >> utilities::PetscUtilities::checkError; + if (dim <= 0) dim = dmDim; + + // geometry: cell and face + Vec cellGeomVec = nullptr, faceGeomVec = nullptr; + DMPlexComputeGeometryFVM(dm, &cellGeomVec, &faceGeomVec) >> utilities::PetscUtilities::checkError; + + DM dmFace = nullptr, dmCell = nullptr; + const PetscScalar* faceGeomArray = nullptr; + const PetscScalar* cellGeomArray = nullptr; + VecGetDM(faceGeomVec, &dmFace) >> utilities::PetscUtilities::checkError; + VecGetDM(cellGeomVec, &dmCell) >> utilities::PetscUtilities::checkError; + VecGetArrayRead(faceGeomVec, &faceGeomArray) >> utilities::PetscUtilities::checkError; + VecGetArrayRead(cellGeomVec, &cellGeomArray) >> utilities::PetscUtilities::checkError; + + // ghost label (for faces) + DMLabel ghostLabel = nullptr; + DMGetLabel(dm, "ghost", &ghostLabel) >> utilities::PetscUtilities::checkError; + + // build map from DMPLEX cell point -> local index [0..nCells) + PetscInt cStart, cEnd; + DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd) >> utilities::PetscUtilities::checkError; + std::vector cellToLocal(cEnd, -1); + for (PetscInt i = cellRange.start; i < cellRange.end; ++i) { + const PetscInt cell = cellRange.points ? cellRange.points[i] : i; + const PetscInt localIdx = i - cellRange.start; + cellToLocal[cell] = localIdx; + } + + // zero output gradients + zMixGrad.assign((std::size_t)nCells * dim, 0.0); + + // loop over cells in the range + for (PetscInt i = cellRange.start; i < cellRange.end; ++i) { + const PetscInt cell = cellRange.points ? cellRange.points[i] : i; + const PetscInt localIdx = i - cellRange.start; + + // cell geometry + const PetscFVCellGeom* cg = nullptr; + DMPlexPointLocalRead(dmCell, cell, cellGeomArray, &cg) >> utilities::PetscUtilities::checkError; + if (!cg || cg->volume <= 0.0) continue; + + const double zCell = zMixValues[(std::size_t)localIdx]; + + // accumulator for ∑ z_f * n_f A_f + double gradAccum[3] = {0.0, 0.0, 0.0}; + + // faces of this cell + PetscInt numFaces = 0; + const PetscInt* faces = nullptr; + DMPlexGetConeSize(dm, cell, &numFaces) >> utilities::PetscUtilities::checkError; + DMPlexGetCone(dm, cell, &faces) >> utilities::PetscUtilities::checkError; + + for (PetscInt fi = 0; fi < numFaces; ++fi) { + const PetscInt face = faces[fi]; + + // skip ghost / boundary / children + PetscInt ghost = -1; + PetscBool boundary; + PetscInt numChildren; + if (ghostLabel) { + DMLabelGetValue(ghostLabel, face, &ghost) >> utilities::PetscUtilities::checkError; + } + DMIsBoundaryPoint(dm, face, &boundary) >> utilities::PetscUtilities::checkError; + DMPlexGetTreeChildren(dm, face, &numChildren, nullptr) >> utilities::PetscUtilities::checkError; + if (ghost >= 0 || boundary || numChildren > 0) continue; + + // get neighbors of this face + const PetscInt* cells = nullptr; + PetscInt supportSize = 0; + DMPlexGetSupportSize(dm, face, &supportSize) >> utilities::PetscUtilities::checkError; + DMPlexGetSupport(dm, face, &cells) >> utilities::PetscUtilities::checkError; + if (supportSize == 0) continue; // should not happen + + const PetscInt c0 = cells[0]; + const PetscInt c1 = (supportSize > 1) ? cells[1] : -1; + + // orientation of normal: use sign to orient fg->normal wrt this cell + double sign = 0.0; + if (cell == c0) { + sign = 1.0; + } else if (cell == c1) { + sign = -1.0; + } else { + continue; // this face is not attached to this cell (paranoia) + } + + // face geometry + const PetscFVFaceGeom* fg = nullptr; + DMPlexPointLocalRead(dmFace, face, faceGeomArray, &fg) >> utilities::PetscUtilities::checkError; + if (!fg) continue; + + // mixture fraction at face: central average if neighbor exists in our range + double zFace = zCell; + if (c1 >= 0) { + const PetscInt otherCell = (cell == c0) ? c1 : c0; + const PetscInt otherLocal = (otherCell >= 0 && otherCell < cEnd) ? cellToLocal[otherCell] : -1; + if (otherLocal >= 0) { + const double zOther = zMixValues[(std::size_t)otherLocal]; + zFace = 0.5 * (zCell + zOther); + } + } + + // accumulate: z_f * (outward normal * area) + for (PetscInt d = 0; d < dim; ++d) { + gradAccum[d] += zFace * sign * static_cast(fg->normal[d]); + } + } + + const double invVol = 1.0 / static_cast(cg->volume); + for (PetscInt d = 0; d < dim; ++d) { + zMixGrad[(std::size_t)localIdx * dim + d] = gradAccum[d] * invVol; + } + } + + // cleanup geometry vecs + VecRestoreArrayRead(faceGeomVec, &faceGeomArray) >> utilities::PetscUtilities::checkError; + VecRestoreArrayRead(cellGeomVec, &cellGeomArray) >> utilities::PetscUtilities::checkError; + VecDestroy(&faceGeomVec) >> utilities::PetscUtilities::checkError; + VecDestroy(&cellGeomVec) >> utilities::PetscUtilities::checkError; +} diff --git a/src/eos/zerork/sourceCalculatorZeroRK.hpp b/src/eos/zerork/sourceCalculatorZeroRK.hpp index 1dafa5eec..59bee4462 100644 --- a/src/eos/zerork/sourceCalculatorZeroRK.hpp +++ b/src/eos/zerork/sourceCalculatorZeroRK.hpp @@ -47,6 +47,10 @@ class SourceCalculator : public ChemistryModel::SourceCalculator, private utilit // max iterations for cvode int dumpfailed = 0; + int n_reactors_max = 10000; + + int n_reactors_min = 100; + // max iterations for cvode int maxiteration = 10000; @@ -99,7 +103,32 @@ class SourceCalculator : public ChemistryModel::SourceCalculator, private utilit */ void AddSource(const ablate::domain::Range& cellRange, Vec localXVec, Vec localFVec) override; + //MMixture fraction calculator for ZeroRK load balancing + double ComputeMixtureFraction(const double* yi) const; + + double SutherlandDiff(double T,double rho); + + void ComputeMixtureFractionGradients( + DM dm, + const ablate::domain::Range& cellRange, + const std::vector& zMixValues, + std::vector& zMixGrad, + PetscInt dim); + private: + + //member variables holding required stuff for zmix calculator + std::vector zMixCoefficients; + double zMixFuel = 0.0; + double zMixOxidizer = 0.0; + + +// std::vector mixtureFraction; + std::vector reactorzMix; + std::vector reactorzMixGrad; + + DM dmZMixGrad = nullptr; // DM for mixture fraction gradient + std::vector sourceZeroRKAtI; zerork_handle zrm_handle; //! copy of constraints diff --git a/src/finiteVolume/cellInterpolant.cpp b/src/finiteVolume/cellInterpolant.cpp index 4f69fd6c2..24671afcb 100644 --- a/src/finiteVolume/cellInterpolant.cpp +++ b/src/finiteVolume/cellInterpolant.cpp @@ -3,8 +3,9 @@ #include ablate::finiteVolume::CellInterpolant::CellInterpolant(std::shared_ptr subDomainIn, const std::shared_ptr& solverRegion, Vec faceGeomVec, Vec cellGeomVec, - double maxGradIn) - : subDomain(std::move(std::move(subDomainIn))), maxLimGrad(maxGradIn) { + double maxGradIn,const ablate::domain::Range& faceRange) + : subDomain(std::move(std::move(subDomainIn))), maxLimGrad(maxGradIn),flowLabelVec(numLabel * (faceRange.end - faceRange.start), 0) { +// StartEvent("FiniteVolumeSolver::CellInterpolant::CompluteRHS::constructor"); auto getGradientDm = [this, solverRegion, faceGeomVec, cellGeomVec](const domain::Field& fieldInfo, std::vector& gradDMs) { auto petscField = subDomain->GetPetscFieldObject(fieldInfo); auto petscFieldFV = (PetscFV)petscField; @@ -26,10 +27,84 @@ ablate::finiteVolume::CellInterpolant::CellInterpolant(std::shared_ptrGetDM(), regionLabel, regionValue); + + auto dm = subDomain->GetDM(); + + DMLabel ghostLabel; + DMGetLabel(dm, "ghost", &ghostLabel) >> utilities::PetscUtilities::checkError; + + DM faceDM, cellDM; + VecGetDM(faceGeomVec, &faceDM) >> utilities::PetscUtilities::checkError; + VecGetDM(cellGeomVec, &cellDM) >> utilities::PetscUtilities::checkError; + + const PetscScalar* cellGeomArray = nullptr; + const PetscScalar* faceGeomArray = nullptr; + VecGetArrayRead(cellGeomVec, &cellGeomArray) >> utilities::PetscUtilities::checkError; + VecGetArrayRead(faceGeomVec, &faceGeomArray) >> utilities::PetscUtilities::checkError; + + + +// for (PetscInt f = faceRange.start; f < faceRange.end; ++f) { +// const PetscInt face = faceRange.points ? faceRange.points[f] : f; +// +// //! The indicies are in order per face: +// // ghost, nsupp, nchild, leftFlowLabelValue, rightFlowLabelValue +// +// +// // make sure that this is a valid face +// PetscInt ghost, nsupp, nchild; +// DMLabelGetValue(ghostLabel, face, &ghost) >> utilities::PetscUtilities::checkError; +// DMPlexGetSupportSize(dm, face, &nsupp) >> utilities::PetscUtilities::checkError; +// DMPlexGetTreeChildren(dm, face, &nchild, nullptr) >> utilities::PetscUtilities::checkError; +// +// flowLabelVec[numLabel*f]=ghost; +// flowLabelVec[numLabel*f+1]=nsupp; +// flowLabelVec[numLabel*f+2]=nchild; +// +// +// // Get the face geometry +// const PetscInt* faceCells; +// PetscFVFaceGeom* fg; +// PetscFVCellGeom *cgL, *cgR; +// DMPlexPointLocalRead(faceDM, face, faceGeomArray, &fg) >> utilities::PetscUtilities::checkError; +// DMPlexGetSupport(dm, face, &faceCells) >> utilities::PetscUtilities::checkError; +// DMPlexPointLocalRead(cellDM, faceCells[0], cellGeomArray, &cgL) >> utilities::PetscUtilities::checkError; +// DMPlexPointLocalRead(cellDM, faceCells[1], cellGeomArray, &cgR) >> utilities::PetscUtilities::checkError; +// +// PetscInt leftFlowLabelValue = regionValue; +// PetscInt rightFlowLabelValue = regionValue; +// // start = MPI_Wtime(); +// if (regionLabel) { +// DMLabelGetValue(regionLabel, faceCells[0], &leftFlowLabelValue); +// DMLabelGetValue(regionLabel, faceCells[1], &rightFlowLabelValue); +// } +// flowLabelVec[numLabel*f+3]=leftFlowLabelValue; +// flowLabelVec[numLabel*f+4]=rightFlowLabelValue; +// +// //Check if the cells are ghost cells +// if (ghostLabel) { +// DMLabelGetValue(ghostLabel, faceCells[0], &leftFlowLabelValue); +// DMLabelGetValue(ghostLabel, faceCells[1], &rightFlowLabelValue); +// } +// flowLabelVec[numLabel*f+5]=leftFlowLabelValue; +// flowLabelVec[numLabel*f+6]=rightFlowLabelValue; +// +// } + + + // Compute the gradient dm for each field that supports it for (const auto& fieldInfo : subDomain->GetFields()) { getGradientDm(fieldInfo, gradientCellDms); } +// EndEvent(); } ablate::finiteVolume::CellInterpolant::~CellInterpolant() { @@ -43,9 +118,11 @@ ablate::finiteVolume::CellInterpolant::~CellInterpolant() { void ablate::finiteVolume::CellInterpolant::ComputeRHS(PetscReal time, Vec locXVec, Vec locAuxVec, Vec locFVec, const std::shared_ptr& solverRegion, std::vector& rhsFunctions, const ablate::domain::Range& faceRange, const ablate::domain::Range& cellRange, Vec cellGeomVec, Vec faceGeomVec) { + StartEvent("FiniteVolumeSolver::CellInterpolant::ComputeRHS::getDMsize"); auto dm = subDomain->GetDM(); auto dmAux = subDomain->GetAuxDM(); + /* 1: Get sizes from dm and dmAux */ PetscSection section = nullptr; DMGetLocalSection(dm, §ion) >> utilities::PetscUtilities::checkError; @@ -56,6 +133,7 @@ void ablate::finiteVolume::CellInterpolant::ComputeRHS(PetscReal time, Vec locXV PetscDSGetNumFields(ds, &nf) >> utilities::PetscUtilities::checkError; PetscDSGetTotalDimension(ds, &totDim) >> utilities::PetscUtilities::checkError; + // Check to see if the dm has an auxVec/auxDM associated with it. If it does, extract it PetscDS dsAux = subDomain->GetAuxDiscreteSystem(); PetscInt naf = 0, totDimAux = 0; @@ -63,7 +141,9 @@ void ablate::finiteVolume::CellInterpolant::ComputeRHS(PetscReal time, Vec locXV PetscDSGetTotalDimension(dsAux, &totDimAux) >> utilities::PetscUtilities::checkError; PetscDSGetNumFields(dsAux, &naf) >> utilities::PetscUtilities::checkError; } + EndEvent(); + StartEvent("FiniteVolumeSolver::CellInterpolant::ComputeRHS::getgeoData"); /* 2: Get geometric data */ // We can use a single call for the geometry data because it does not depend on the fv object const PetscScalar* cellGeomArray = nullptr; @@ -87,19 +167,28 @@ void ablate::finiteVolume::CellInterpolant::ComputeRHS(PetscReal time, Vec locXV // there must be a separate gradient vector/dm for field because they can be different sizes std::vector locGradVecs(nf, nullptr); + EndEvent(); + /* Reconstruct and limit cell gradients */ // for each field compute the gradient in the localGrads vector for (const auto& field : subDomain->GetFields()) { + // This is logged in the function ComputeFieldGradients(field, locXVec, locGradVecs[field.subId], gradientCellDms[field.subId], cellGeomVec, faceGeomVec, faceRange, cellRange); } +// StartEvent("FiniteVolumeSolver::CellInterpolant::ComputeRHS::setlocalvec"); std::vector locGradArrays(nf, nullptr); for (const auto& field : subDomain->GetFields()) { if (locGradVecs[field.subId]) { VecGetArrayRead(locGradVecs[field.subId], &locGradArrays[field.subId]) >> utilities::PetscUtilities::checkError; } } +// EndEvent(); + if (time <= 1.1E-8) std::cout << "The mesh has " << faceRange.end << " faces \n" ; + if (time <= 1.1E-8) std::cout << "The mesh has " << cellRange.end << " cells \n" ; + + //This is measured inside the function ComputeFluxSourceTerms(dm, ds, totDim, @@ -119,7 +208,7 @@ void ablate::finiteVolume::CellInterpolant::ComputeRHS(PetscReal time, Vec locXV rhsFunctions, faceRange, cellRange); - + StartEvent("FiniteVolumeSolver::CellInterpolant::ComputeRHS::cleanup"); // clean up cell grads for (const auto& field : subDomain->GetFields()) { if (locGradVecs[field.subId]) { @@ -137,6 +226,7 @@ void ablate::finiteVolume::CellInterpolant::ComputeRHS(PetscReal time, Vec locXV VecRestoreArray(locFVec, &locFArray) >> utilities::PetscUtilities::checkError; VecRestoreArrayRead(faceGeomVec, (const PetscScalar**)&faceGeomArray) >> utilities::PetscUtilities::checkError; VecRestoreArrayRead(cellGeomVec, (const PetscScalar**)&cellGeomArray) >> utilities::PetscUtilities::checkError; + EndEvent(); } void ablate::finiteVolume::CellInterpolant::ComputeRHS(PetscReal time, Vec locXVec, Vec locAuxVec, Vec locFVec, const std::shared_ptr& solverRegion, @@ -329,13 +419,14 @@ static PetscErrorCode DMPlexApplyLimiter_Internal(DM dm, DM dmCell, PetscLimiter void ablate::finiteVolume::CellInterpolant::ComputeFieldGradients(const domain::Field& field, Vec xLocalVec, Vec& gradLocVec, DM& dmGrad, Vec cellGeomVec, Vec faceGeomVec, const ablate::domain::Range& faceRange, const ablate::domain::Range& cellRange) { - // get the FVM petsc field associated with this field - auto fvm = (PetscFV)subDomain->GetPetscFieldObject(field); - auto dm = subDomain->GetFieldDM(field); +// StartEvent("FiniteVolumeSolver::CellInterpolant::ComputeRHS::ComputeFieldGradients"); +// get the FVM petsc field associated with this field + auto fvm = (PetscFV)subDomain->GetPetscFieldObject(field); // 1 read + auto dm = subDomain->GetFieldDM(field); // 1 read // Get the dm for this grad field // If there is no grad, return - if (!dmGrad) { + if (!dmGrad) { // 1 write return; } @@ -343,13 +434,17 @@ void ablate::finiteVolume::CellInterpolant::ComputeFieldGradients(const domain:: DMGetLocalVector(dmGrad, &gradLocVec) >> utilities::PetscUtilities::checkError; // Get the correct sized vec (gradient for this field) - Vec gradGlobVec; - DMGetGlobalVector(dmGrad, &gradGlobVec) >> utilities::PetscUtilities::checkError; - VecZeroEntries(gradGlobVec) >> utilities::PetscUtilities::checkError; + Vec gradGlobVec; // 1 write + DMGetGlobalVector(dmGrad, &gradGlobVec) >> utilities::PetscUtilities::checkError; // N read + VecZeroEntries(gradGlobVec) >> utilities::PetscUtilities::checkError; // vector length writes + +// PetscInt size; +// VecGetSize(gradGlobVec, &size)>> utilities::PetscUtilities::checkError; +// std::cout << "Global size of gradGlobVec: " << size << std::endl; // check to see if there is a ghost label DMLabel ghostLabel; - DMGetLabel(dm, "ghost", &ghostLabel) >> utilities::PetscUtilities::checkError; + DMGetLabel(dm, "ghost", &ghostLabel) >> utilities::PetscUtilities::checkError; // 1 read, 1 write // Get the face geometry DM dmFace; @@ -363,7 +458,7 @@ void ablate::finiteVolume::CellInterpolant::ComputeFieldGradients(const domain:: // extract the global grad array PetscScalar* gradGlobArray; - VecGetArray(gradGlobVec, &gradGlobArray); + VecGetArray(gradGlobVec, &gradGlobArray); //This returns just the pointer, no copies. // Get the dof and dim PetscInt dim = subDomain->GetDimensions(); @@ -376,7 +471,7 @@ void ablate::finiteVolume::CellInterpolant::ComputeFieldGradients(const domain:: PetscBool boundary; PetscInt ghost = -1; if (ghostLabel) { - DMLabelGetValue(ghostLabel, face, &ghost); + DMLabelGetValue(ghostLabel, face, &ghost); //1 read } DMIsBoundaryPoint(dm, face, &boundary); PetscInt numChildren; @@ -408,14 +503,16 @@ void ablate::finiteVolume::CellInterpolant::ComputeFieldGradients(const domain:: for (PetscInt d = 0; d < dim; ++d) { if (cgrad[0]) cgrad[0][pd * dim + d] += fg->grad[0][d] * delta; if (cgrad[1]) cgrad[1][pd * dim + d] -= fg->grad[1][d] * delta; +// std::cout << "The index is: " << pd * dim + d << std::endl; + } } } - +// EndEvent(); // Check for a limiter the limiter PetscLimiter lim; PetscFVGetLimiter(fvm, &lim) >> utilities::PetscUtilities::checkError; - if (lim) { + if (false) { /* Limit interior gradients (using cell-based loop because it generalizes better to vector limiters) */ // Get the cell geometry DM dmCell; @@ -434,7 +531,7 @@ void ablate::finiteVolume::CellInterpolant::ComputeFieldGradients(const domain:: PetscScalar* cx; // cell solution/field values PetscFVCellGeom* cg; // cell geometry PetscScalar* cgrad; // cell gradient - PetscInt coneSize; // cell conectivity + PetscInt coneSize; // cell connectivity DMPlexGetConeSize(dm, cell, &coneSize) >> utilities::PetscUtilities::checkError; DMPlexGetCone(dm, cell, &cellFaces) >> utilities::PetscUtilities::checkError; @@ -485,6 +582,8 @@ void ablate::finiteVolume::CellInterpolant::ComputeFieldGradients(const domain:: DMRestoreWorkArray(dm, dof, MPIU_REAL, &cellPhi) >> utilities::PetscUtilities::checkError; VecRestoreArrayRead(cellGeomVec, &cellGeometryArray); } + + StartEvent("FiniteVolumeSolver::CellInterpolant::ComputeRHS::GradientComm"); // Communicate gradient values VecRestoreArray(gradGlobVec, &gradGlobArray) >> utilities::PetscUtilities::checkError; DMGlobalToLocalBegin(dmGrad, gradGlobVec, INSERT_VALUES, gradLocVec) >> utilities::PetscUtilities::checkError; @@ -494,6 +593,7 @@ void ablate::finiteVolume::CellInterpolant::ComputeFieldGradients(const domain:: VecRestoreArrayRead(xLocalVec, &xLocalArray) >> utilities::PetscUtilities::checkError; VecRestoreArrayRead(faceGeomVec, &faceGeometryArray) >> utilities::PetscUtilities::checkError; DMRestoreGlobalVector(dmGrad, &gradGlobVec) >> utilities::PetscUtilities::checkError; + EndEvent(); } void ablate::finiteVolume::CellInterpolant::ComputeFluxSourceTerms(DM dm, PetscDS ds, PetscInt totDim, const PetscScalar* xArray, DM dmAux, PetscDS dsAux, PetscInt totDimAux, @@ -502,8 +602,10 @@ void ablate::finiteVolume::CellInterpolant::ComputeFluxSourceTerms(DM dm, PetscD const std::shared_ptr& solverRegion, std::vector& rhsFunctions, const ablate::domain::Range& faceRange, const ablate::domain::Range& cellRange) { +// StartEvent("FiniteVolumeSolver::CellInterpolant::ComputeRHS::ComputeSource::Setup"); PetscInt dim = subDomain->GetDimensions(); + // Size up the work arrays (uL, uR, gradL, gradR, auxL, auxR, gradAuxL, gradAuxR), these are only sized for one face at a time PetscScalar* flux; DMGetWorkArray(dm, totDim, MPIU_SCALAR, &flux) >> utilities::PetscUtilities::checkError; @@ -557,85 +659,272 @@ void ablate::finiteVolume::CellInterpolant::ComputeFluxSourceTerms(DM dm, PetscD DMLabel regionLabel = nullptr; PetscInt regionValue = 0; domain::Region::GetLabel(solverRegion, subDomain->GetDM(), regionLabel, regionValue); - // March over each face in this region - for (PetscInt f = faceRange.start; f < faceRange.end; ++f) { - const PetscInt face = faceRange.points ? faceRange.points[f] : f; - // make sure that this is a valid face - PetscInt ghost, nsupp, nchild; - DMLabelGetValue(ghostLabel, face, &ghost) >> utilities::PetscUtilities::checkError; - DMPlexGetSupportSize(dm, face, &nsupp) >> utilities::PetscUtilities::checkError; - DMPlexGetTreeChildren(dm, face, &nchild, nullptr) >> utilities::PetscUtilities::checkError; - if (ghost >= 0 || nsupp > 2 || nchild > 0) continue; + if (regionValue!=1){ + throw std::runtime_error("The regionValue in cellinterpollant is " + std::to_string(regionValue) +", the main is assumed to have a regionValue of 1"); - // Get the face geometry - const PetscInt* faceCells; - PetscFVFaceGeom* fg; - PetscFVCellGeom *cgL, *cgR; - DMPlexPointLocalRead(faceDM, face, faceGeomArray, &fg) >> utilities::PetscUtilities::checkError; - DMPlexGetSupport(dm, face, &faceCells) >> utilities::PetscUtilities::checkError; - DMPlexPointLocalRead(cellDM, faceCells[0], cellGeomArray, &cgL) >> utilities::PetscUtilities::checkError; - DMPlexPointLocalRead(cellDM, faceCells[1], cellGeomArray, &cgR) >> utilities::PetscUtilities::checkError; - - PetscInt leftFlowLabelValue = regionValue; - PetscInt rightFlowLabelValue = regionValue; - if (regionLabel) { - DMLabelGetValue(regionLabel, faceCells[0], &leftFlowLabelValue); - DMLabelGetValue(regionLabel, faceCells[1], &rightFlowLabelValue); - } - // compute the left/right face values - ProjectToFace(subDomain->GetFields(), ds, *fg, faceCells[0], *cgL, dm, xArray, dmGrads, locGradArrays, uL, gradL, leftFlowLabelValue == regionValue); - ProjectToFace(subDomain->GetFields(), ds, *fg, faceCells[1], *cgR, dm, xArray, dmGrads, locGradArrays, uR, gradR, rightFlowLabelValue == regionValue); + } - // determine the left/right cells - if (auxArray) { - // Get the field values at this cell - DMPlexPointLocalRead(dmAux, faceCells[0], auxArray, &auxL) >> utilities::PetscUtilities::checkError; - DMPlexPointLocalRead(dmAux, faceCells[1], auxArray, &auxR) >> utilities::PetscUtilities::checkError; - } +// double totalTime1=0; +// int callCount1=0; + // double totalTime2=0; + // int callCount2=0; + // double totalTime3=0; + // int callCount3=0; +// int output_it=1000; +// +// double start = MPI_Wtime(); + // start = MPI_Wtime(); + + + + +// StartEvent("FiniteVolumeSolver::CellInterpolant::ComputeRHS::ComputeSource::Project"); +// for (PetscInt f = faceRange.start; f < faceRange.end; ++f) { +// const PetscInt face = faceRange.points ? faceRange.points[f] : f; +// +// // make sure that this is a valid facel +// PetscInt ghost, nsupp, nchild; +// DMLabelGetValue(ghostLabel, face, &ghost) >> utilities::PetscUtilities::checkError; +// DMPlexGetSupportSize(dm, face, &nsupp) >> utilities::PetscUtilities::checkError; +// DMPlexGetTreeChildren(dm, face, &nchild, nullptr) >> utilities::PetscUtilities::checkError; +// +// if (ghost >= 0 || nsupp > 2 || nchild > 0) continue; +// // Get the face geometry +// const PetscInt* faceCells; +// PetscFVFaceGeom* fg; +// PetscFVCellGeom *cgL, *cgR; +// DMPlexPointLocalRead(faceDM, face, faceGeomArray, &fg) >> utilities::PetscUtilities::checkError; +// DMPlexGetSupport(dm, face, &faceCells) >> utilities::PetscUtilities::checkError; +// DMPlexPointLocalRead(cellDM, faceCells[0], cellGeomArray, &cgL) >> utilities::PetscUtilities::checkError; +// DMPlexPointLocalRead(cellDM, faceCells[1], cellGeomArray, &cgR) >> utilities::PetscUtilities::checkError; +// +// // compute the left/right face values +//// ProjectToFace(subDomain->GetFields(), ds, *fg, faceCells[0], *cgL, dm, xArray, dmGrads, locGradArrays, uL, gradL, 1); +//// ProjectToFace(subDomain->GetFields(), ds, *fg, faceCells[1], *cgR, dm, xArray, dmGrads, locGradArrays, uR, gradR, 1); +// +// ProjectToFace(subDomain->GetFields(), ds, *fg, faceCells[0], *cgL, dm, xArray, dmGrads, locGradArrays, uL, gradL,faceCells[1], *cgR, 1); +// ProjectToFace(subDomain->GetFields(), ds, *fg, faceCells[1], *cgR, dm, xArray, dmGrads, locGradArrays, uR, gradR, faceCells[0], *cgL, 1); +// } +// EndEvent(); + + + + +// StartEvent("FiniteVolumeSolver::CellInterpolant::ComputeRHS::ComputeSource::SourceTotal"); + + bool oldcode=true; + if (oldcode) { + for (PetscInt f = faceRange.start; f < faceRange.end; ++f) { + const PetscInt face = faceRange.points ? faceRange.points[f] : f; + + // make sure that this is a valid facel + PetscInt ghost, nsupp, nchild; + DMLabelGetValue(ghostLabel, face, &ghost) >> utilities::PetscUtilities::checkError; + DMPlexGetSupportSize(dm, face, &nsupp) >> utilities::PetscUtilities::checkError; + DMPlexGetTreeChildren(dm, face, &nchild, nullptr) >> utilities::PetscUtilities::checkError; + + if (ghost >= 0 || nsupp > 2 || nchild > 0) continue; + // Get the face geometry + const PetscInt* faceCells; + PetscFVFaceGeom* fg; + PetscFVCellGeom *cgL, *cgR; + DMPlexPointLocalRead(faceDM, face, faceGeomArray, &fg) >> utilities::PetscUtilities::checkError; + DMPlexGetSupport(dm, face, &faceCells) >> utilities::PetscUtilities::checkError; - // March over each source function - for (std::size_t fun = 0; fun < rhsFunctions.size(); fun++) { - PetscInt fluxOffset = 0; // Flux offset for the function ( Currently calculated by just adding the number of components of the previous fields) - PetscArrayzero(flux, totDim) >> utilities::PetscUtilities::checkError; - const auto& rhsFluxFunctionDescription = rhsFunctions[fun]; - rhsFluxFunctionDescription.function(dim, fg, uOff[fun].data(), uL, uR, aOff[fun].data(), auxL, auxR, flux, rhsFluxFunctionDescription.context) >> utilities::PetscUtilities::checkError; - // add the fluxes back to the cell - for (std::size_t updateFieldIdx = 0; updateFieldIdx < rhsFunctions[fun].updateFields.size(); updateFieldIdx++) { - PetscInt cellLabelValue = regionValue; - PetscScalar *fL = nullptr, *fR = nullptr; - DMLabelGetValue(ghostLabel, faceCells[0], &ghost) >> utilities::PetscUtilities::checkError; - if (regionLabel) { - DMLabelGetValue(regionLabel, faceCells[0], &cellLabelValue) >> utilities::PetscUtilities::checkError; - } - if (ghost <= 0 && regionValue == cellLabelValue) { - DMPlexPointLocalFieldRef(dm, faceCells[0], fluxId[fun][updateFieldIdx], locFArray, &fL) >> utilities::PetscUtilities::checkError; - } + DMPlexPointLocalRead(cellDM, faceCells[0], cellGeomArray, &cgL) >> utilities::PetscUtilities::checkError; + DMPlexPointLocalRead(cellDM, faceCells[1], cellGeomArray, &cgR) >> utilities::PetscUtilities::checkError; - cellLabelValue = regionValue; - DMLabelGetValue(ghostLabel, faceCells[1], &ghost) >> utilities::PetscUtilities::checkError; - if (regionLabel) { - DMLabelGetValue(regionLabel, faceCells[1], &cellLabelValue) >> utilities::PetscUtilities::checkError; - } - if (ghost <= 0 && regionValue == cellLabelValue) { - DMPlexPointLocalFieldRef(dm, faceCells[1], fluxId[fun][updateFieldIdx], locFArray, &fR) >> utilities::PetscUtilities::checkError; + PetscInt leftFlowLabelValue = regionValue; + PetscInt rightFlowLabelValue = regionValue; + + if (regionLabel) { + DMLabelGetValue(regionLabel, faceCells[0], &leftFlowLabelValue); + DMLabelGetValue(regionLabel, faceCells[1], &rightFlowLabelValue); + } + +// ProjectToFace(subDomain->GetFields(), ds, *fg, faceCells[0], *cgL, dm, xArray, dmGrads, locGradArrays, uL, gradL, leftFlowLabelValue == regionValue); +// ProjectToFace(subDomain->GetFields(), ds, *fg, faceCells[1], *cgR, dm, xArray, dmGrads, locGradArrays, uR, gradR, rightFlowLabelValue == regionValue); + + ProjectToFace(subDomain->GetFields(), ds, *fg, faceCells[0], *cgL, dm, xArray, dmGrads, locGradArrays, uL, gradL,faceCells[1], *cgR, leftFlowLabelValue == regionValue); + ProjectToFace(subDomain->GetFields(), ds, *fg, faceCells[1], *cgR, dm, xArray, dmGrads, locGradArrays, uR, gradR, faceCells[0], *cgL, rightFlowLabelValue == regionValue); + + // determine the left/right cells + if (auxArray) { + // Get the field values at this cell + DMPlexPointLocalRead(dmAux, faceCells[0], auxArray, &auxL) >> utilities::PetscUtilities::checkError; + DMPlexPointLocalRead(dmAux, faceCells[1], auxArray, &auxR) >> utilities::PetscUtilities::checkError; + } + + for (std::size_t fun = 0; fun < rhsFunctions.size(); fun++) { + PetscInt fluxOffset = 0; // Flux offset for the function ( Currently calculated by just adding the number of components of the previous fields) + PetscArrayzero(flux, totDim) >> utilities::PetscUtilities::checkError; + const auto& rhsFluxFunctionDescription = rhsFunctions[fun]; + rhsFluxFunctionDescription.function(dim, fg, uOff[fun].data(), uL, uR, aOff[fun].data(), auxL, auxR, flux, rhsFluxFunctionDescription.context) >> utilities::PetscUtilities::checkError; + // add the fluxes back to the cell + for (std::size_t updateFieldIdx = 0; updateFieldIdx < rhsFunctions[fun].updateFields.size(); updateFieldIdx++) { + PetscInt cellLabelValue = regionValue; + PetscScalar *fL = nullptr, *fR = nullptr; + DMLabelGetValue(ghostLabel, faceCells[0], &ghost) >> utilities::PetscUtilities::checkError; + if (regionLabel) { + DMLabelGetValue(regionLabel, faceCells[0], &cellLabelValue) >> utilities::PetscUtilities::checkError; + } + if (ghost <= 0 && regionValue == cellLabelValue) { + DMPlexPointLocalFieldRef(dm, faceCells[0], fluxId[fun][updateFieldIdx], locFArray, &fL) >> utilities::PetscUtilities::checkError; + } + + cellLabelValue = regionValue; + DMLabelGetValue(ghostLabel, faceCells[1], &ghost) >> utilities::PetscUtilities::checkError; + if (regionLabel) { + DMLabelGetValue(regionLabel, faceCells[1], &cellLabelValue) >> utilities::PetscUtilities::checkError; + } + if (ghost <= 0 && regionValue == cellLabelValue) { + DMPlexPointLocalFieldRef(dm, faceCells[1], fluxId[fun][updateFieldIdx], locFArray, &fR) >> utilities::PetscUtilities::checkError; + } + + for (PetscInt d = 0; d < (fluxComponentSize[fun][updateFieldIdx]); ++d) { + if (fL) fL[d] -= flux[fluxOffset + d] / cgL->volume; + if (fR) fR[d] += flux[fluxOffset + d] / cgR->volume; + } + fluxOffset += fluxComponentSize[fun][updateFieldIdx]; } + // EndEvent(); + } + // ++callCount3; + // if (callCount3==output_it) { + // PetscPrintf(PETSC_COMM_WORLD, "Section3 total %d time: %f s\n", callCount3, totalTime3); } + + // EndEvent(); + } + } else { + for (PetscInt f = faceRange.start; f < faceRange.end; ++f) { + const PetscInt face = faceRange.points ? faceRange.points[f] : f; + + // make sure that this is a valid facel +// PetscInt ghost2, nsupp2, nchild2; +// DMLabelGetValue(ghostLabel, face, &ghost2) >> utilities::PetscUtilities::checkError; +// DMPlexGetSupportSize(dm, face, &nsupp2) >> utilities::PetscUtilities::checkError; +// DMPlexGetTreeChildren(dm, face, &nchild2, nullptr) >> utilities::PetscUtilities::checkError; +// + PetscInt ghost = flowLabelVec[numLabel*f]; + PetscInt nsupp = flowLabelVec[numLabel*f+1]; + PetscInt nchild = flowLabelVec[numLabel*f+2]; + +// if (ghost2!=ghost || nsupp2!=nsupp || nchild2!=nchild){ +// throw std::runtime_error("either ghost or nsupp or child in not matching for face:" + std::to_string(face) ); +// } + + if (ghost >= 0 || nsupp > 2 || nchild > 0) continue; + // Get the face geometry + const PetscInt* faceCells; + PetscFVFaceGeom* fg; + PetscFVCellGeom *cgL, *cgR; + DMPlexPointLocalRead(faceDM, face, faceGeomArray, &fg) >> utilities::PetscUtilities::checkError; + DMPlexGetSupport(dm, face, &faceCells) >> utilities::PetscUtilities::checkError; + + DMPlexPointLocalRead(cellDM, faceCells[0], cellGeomArray, &cgL) >> utilities::PetscUtilities::checkError; + DMPlexPointLocalRead(cellDM, faceCells[1], cellGeomArray, &cgR) >> utilities::PetscUtilities::checkError; + + PetscInt leftFlowLabelValue = flowLabelVec[numLabel*f+3]; + PetscInt rightFlowLabelValue = flowLabelVec[numLabel*f+4]; + +// PetscInt leftFlowLabelValue2 = regionValue; +// PetscInt rightFlowLabelValue2 = regionValue; +// if (regionLabel) { +// DMLabelGetValue(regionLabel, faceCells[0], &leftFlowLabelValue2); +// DMLabelGetValue(regionLabel, faceCells[1], &rightFlowLabelValue2); +// } +// if (leftFlowLabelValue!=leftFlowLabelValue2 || rightFlowLabelValue2!=rightFlowLabelValue ){ +// throw std::runtime_error("one of the face labels dont match matching for face:" + std::to_string(face) ); +// } + +// ProjectToFace(subDomain->GetFields(), ds, *fg, faceCells[0], *cgL, dm, xArray, dmGrads, locGradArrays, uL, gradL, leftFlowLabelValue == regionValue); +// ProjectToFace(subDomain->GetFields(), ds, *fg, faceCells[1], *cgR, dm, xArray, dmGrads, locGradArrays, uR, gradR, rightFlowLabelValue == regionValue); + + ProjectToFace(subDomain->GetFields(), ds, *fg, faceCells[0], *cgL, dm, xArray, dmGrads, locGradArrays, uL, gradL,faceCells[1], *cgR, leftFlowLabelValue == regionValue); + ProjectToFace(subDomain->GetFields(), ds, *fg, faceCells[1], *cgR, dm, xArray, dmGrads, locGradArrays, uR, gradR, faceCells[0], *cgL, rightFlowLabelValue == regionValue); + + + // determine the left/right cells + if (auxArray) { + // Get the field values at this cell + DMPlexPointLocalRead(dmAux, faceCells[0], auxArray, &auxL) >> utilities::PetscUtilities::checkError; + DMPlexPointLocalRead(dmAux, faceCells[1], auxArray, &auxR) >> utilities::PetscUtilities::checkError; + } + + for (std::size_t fun = 0; fun < rhsFunctions.size(); fun++) { + PetscInt fluxOffset = 0; // Flux offset for the function ( Currently calculated by just adding the number of components of the previous fields) + PetscArrayzero(flux, totDim) >> utilities::PetscUtilities::checkError; + const auto& rhsFluxFunctionDescription = rhsFunctions[fun]; + rhsFluxFunctionDescription.function(dim, fg, uOff[fun].data(), uL, uR, aOff[fun].data(), auxL, auxR, flux, rhsFluxFunctionDescription.context) >> utilities::PetscUtilities::checkError; + // add the fluxes back to the cell + for (std::size_t updateFieldIdx = 0; updateFieldIdx < rhsFunctions[fun].updateFields.size(); updateFieldIdx++) { + PetscInt cellLabelValue = regionValue; +// + PetscScalar *fL = nullptr, *fR = nullptr; + +// PetscInt cellLabelValue2 = regionValue; +// DMLabelGetValue(ghostLabel, faceCells[0], &ghost2) >> utilities::PetscUtilities::checkError; +// if (regionLabel) { +// DMLabelGetValue(regionLabel, faceCells[0], &cellLabelValue2) >> utilities::PetscUtilities::checkError; +// } + + ghost = flowLabelVec[numLabel*f+5]; + cellLabelValue = flowLabelVec[numLabel*f+3]; + +// if (ghost!=ghost2 || cellLabelValue2!=cellLabelValue ){ +// throw std::runtime_error("either ghost or cellLabelValue2 (left) is not good for source terms, for face:" + std::to_string(face) ); +// } + + if (ghost <= 0 && regionValue == cellLabelValue) { + DMPlexPointLocalFieldRef(dm, faceCells[0], fluxId[fun][updateFieldIdx], locFArray, &fL) >> utilities::PetscUtilities::checkError; + } + + cellLabelValue = regionValue; +// +// DMLabelGetValue(ghostLabel, faceCells[1], &ghost2) >> utilities::PetscUtilities::checkError; +// if (regionLabel) { +// DMLabelGetValue(regionLabel, faceCells[1], &cellLabelValue2) >> utilities::PetscUtilities::checkError; +// } + + + ghost = flowLabelVec[numLabel*f+6]; + cellLabelValue = flowLabelVec[numLabel*f+4]; - for (PetscInt d = 0; d < (fluxComponentSize[fun][updateFieldIdx]); ++d) { - if (fL) fL[d] -= flux[fluxOffset + d] / cgL->volume; - if (fR) fR[d] += flux[fluxOffset + d] / cgR->volume; +// if (ghost!=ghost2 || cellLabelValue2!=cellLabelValue ){ +// throw std::runtime_error("either ghost or cellLabelValue2 (right) is not good for source terms, for face:" + std::to_string(face) ); +// } + + if (ghost <= 0 && regionValue == cellLabelValue) { + DMPlexPointLocalFieldRef(dm, faceCells[1], fluxId[fun][updateFieldIdx], locFArray, &fR) >> utilities::PetscUtilities::checkError; + } + + for (PetscInt d = 0; d < (fluxComponentSize[fun][updateFieldIdx]); ++d) { + if (fL) fL[d] -= flux[fluxOffset + d] / cgL->volume; + if (fR) fR[d] += flux[fluxOffset + d] / cgR->volume; + } + fluxOffset += fluxComponentSize[fun][updateFieldIdx]; } - fluxOffset += fluxComponentSize[fun][updateFieldIdx]; + // EndEvent(); } + // ++callCount3; + // if (callCount3==output_it) { + // PetscPrintf(PETSC_COMM_WORLD, "Section3 total %d time: %f s\n", callCount3, totalTime3); } + + // EndEvent(); } } +// EndEvent(); // cleanup +// StartEvent("FiniteVolumeSolver::CellInterpolant::ComputeRHS::ComputeSource::Cleanup"); DMRestoreWorkArray(dm, totDim, MPIU_SCALAR, &flux) >> utilities::PetscUtilities::checkError; DMRestoreWorkArray(dm, totDim, MPIU_SCALAR, &uL) >> utilities::PetscUtilities::checkError; DMRestoreWorkArray(dm, totDim, MPIU_SCALAR, &uR) >> utilities::PetscUtilities::checkError; DMRestoreWorkArray(dm, dim * totDim, MPIU_SCALAR, &gradL) >> utilities::PetscUtilities::checkError; DMRestoreWorkArray(dm, dim * totDim, MPIU_SCALAR, &gradR) >> utilities::PetscUtilities::checkError; +// EndEvent(); + } static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, DMLabel regionLabel, PetscInt regionValue, PetscFV fvm, DM dmFace, PetscScalar* fgeom, DM dmCell, PetscScalar* cgeom) { @@ -887,10 +1176,31 @@ PetscErrorCode ablate::finiteVolume::CellInterpolant::ComputeGradientFVM(DM dm, PetscCall(PetscSectionDestroy(§ionGrad)); PetscFunctionReturn(0); } -void ablate::finiteVolume::CellInterpolant::ProjectToFace(const std::vector& fields, PetscDS ds, const PetscFVFaceGeom& faceGeom, PetscInt cellId, const PetscFVCellGeom& cellGeom, + + +void ablate::finiteVolume::CellInterpolant::ProjectToFace(const std::vector& fields, PetscDS ds, const PetscFVFaceGeom& faceGeom, PetscInt cellId, const PetscFVCellGeom& cellGeom, DM dm, const PetscScalar* xArray, const std::vector& dmGrads, const std::vector& gradArrays, PetscScalar* u, PetscScalar* grad, bool projectField) { +// StartEvent("FiniteVolumeSolver::CellInterpolant::ComputeRHS::ProjectToFace"); + + //Timing +// double start = MPI_Wtime(); + + + //Papi low level +// int EventSet = PAPI_NULL; +// long long values[1]; +// PAPI_create_eventset(&EventSet); +// PAPI_add_event(EventSet, PAPI_DP_OPS); +// PAPI_start(EventSet); + + + //PAPI high level +// int retval; +// retval = PAPI_hl_region_begin("project"); + const auto dim = subDomain->GetDimensions(); + // [R: 1] — Read subDomain // Keep track of derivative offset PetscInt* offsets; @@ -898,8 +1208,10 @@ void ablate::finiteVolume::CellInterpolant::ProjectToFace(const std::vector> utilities::PetscUtilities::checkError; PetscDSGetComponentDerivativeOffsets(ds, &dirOffsets) >> utilities::PetscUtilities::checkError; + // March over each field for (const auto& field : fields) { + // [R: 1], [W: 1] per loop variable PetscReal dx[3]; PetscScalar* xCell; PetscScalar* gradCell; @@ -918,7 +1230,9 @@ void ablate::finiteVolume::CellInterpolant::ProjectToFace(const std::vector& fields, PetscDS ds, const PetscFVFaceGeom& faceGeom, PetscInt cellId, const PetscFVCellGeom& cellGeom, + DM dm, const PetscScalar* xArray, const std::vector& dmGrads, const std::vector& gradArrays, PetscScalar* u, + PetscScalar* grad, PetscInt neighborCellId, const PetscFVCellGeom& neighborCellGeom, bool projectField) { + + const auto dim = subDomain->GetDimensions(); + + // Keep track of derivative offset + PetscInt* offsets; + PetscInt* dirOffsets; + PetscDSGetComponentOffsets(ds, &offsets) >> utilities::PetscUtilities::checkError; + PetscDSGetComponentDerivativeOffsets(ds, &dirOffsets) >> utilities::PetscUtilities::checkError; + + // March over each field + for (const auto& field : fields) { + PetscReal dx[3]; + PetscScalar* xCell; + PetscScalar* gradCell; + + // Get the field values at this cell + DMPlexPointLocalFieldRead(dm, cellId, field.subId, xArray, &xCell) >> utilities::PetscUtilities::checkError; + + // If we need to project the field + if (projectField && dmGrads[field.subId]) { + DMPlexPointLocalRead(dmGrads[field.subId], cellId, gradArrays[field.subId], &gradCell) >> utilities::PetscUtilities::checkError; + DMPlex_WaxpyD_Internal(dim, -1, cellGeom.centroid, faceGeom.centroid, dx); + + // Apply limiter to the gradient + auto fvm = (PetscFV)subDomain->GetPetscFieldObject(field); + PetscLimiter lim; + PetscFVGetLimiter(fvm, &lim) >> utilities::PetscUtilities::checkError; + + PetscReal limitedGrad[dim * field.numberComponents]; + if (lim && neighborCellId >= 0 ) { + // Get neighbor cell values + PetscScalar* neighborXCell; + DMPlexPointLocalFieldRead(dm, neighborCellId, field.subId, xArray, &neighborXCell) >> utilities::PetscUtilities::checkError; + + // Calculate distance vector between cell centers + PetscReal cellDistance[3]; + // v_i = NeighborCentroid_i - ThisCentroid_i = dx_i + DMPlex_WaxpyD_Internal(dim, -1, cellGeom.centroid, neighborCellGeom.centroid, cellDistance); + + // Apply limiting per component based on this specific face + for (PetscInt c = 0; c < field.numberComponents; ++c) { + PetscReal phi = 1.0; // Default to no limiting + + // Calculate the denom(gradient dot distance) + PetscReal denom = 0.0; + for (PetscInt d = 0; d < dim; ++d) { + denom += gradCell[c * dim + d] * cellDistance[d]; + } + + if (PetscAbsReal(denom) > PETSC_SMALL) { + // Use symmetric slope limiter form (Berger, Aftosmis, and Murman 2005) + PetscReal flim = 0.5 * PetscRealPart(neighborXCell[c] - xCell[c]) / denom; + PetscLimiterLimit(lim, flim, &phi) >> utilities::PetscUtilities::checkError; + } + +// // Apply additional maxGradient limiting if needed +// PetscBool cancelGrad = PETSC_FALSE; +// if (phi == 0.0) { +// for (PetscInt d = 0; d < dim; d++) { +// if (PetscAbsReal(gradCell[c * dim + d]) > maxLimGrad) { +// cancelGrad = PETSC_TRUE; +// break; +// } +// } +// } +// // Apply limiting to gradient components +// for (PetscInt d = 0; d < dim; ++d) { +// if (cancelGrad) { +// limitedGrad[c * dim + d] = 0.0; +// } else { +// limitedGrad[c * dim + d] = phi * gradCell[c * dim + d]; +// } +// } +// } + + // Apply limiting to gradient components + for (PetscInt d = 0; d < dim; ++d) { + limitedGrad[c * dim + d] = phi * gradCell[c * dim + d]; + } + } + + for (PetscInt c = 0; c < field.numberComponents; ++c) { + PetscReal projection = 0.0; + for (PetscInt d = 0; d < dim; ++d) { + projection += limitedGrad[c * dim + d] * dx[d]; + } + u[offsets[field.subId] + c] = xCell[c] + projection; + + // Copy the limited gradient into the grad vector + for (PetscInt d = 0; d < dim; d++) { + grad[dirOffsets[field.subId] + c * dim + d] = limitedGrad[c * dim + d]; + } + + + double xlim_min=0.30841; + double xlim_max=1; //0.351656//0.341656 is the end plane of the nozzle + //Overwrite with cell centered values. This is just so the rocket doesnt break in the nozzle... 11/1/25 + if(xlim_min> utilities::PetscUtilities::checkError; + for (PetscInt c = 0; c < field.numberComponents; ++c) { + u[offsets[field.subId] + c] = xCell[c]; + for (PetscInt d = 0; d < dim; d++) { + grad[dirOffsets[field.subId] + c * dim + d] = gradCell[c * dim + d]; + } + } + + } else { + // Just copy the cell centered value on to the face + for (PetscInt c = 0; c < field.numberComponents; ++c) { + u[offsets[field.subId] + c] = xCell[c]; + + // fill the grad with NAN to prevent use + for (PetscInt d = 0; d < dim; d++) { + grad[dirOffsets[field.subId] + c * dim + d] = NAN; + } + } + } + } +} \ No newline at end of file diff --git a/src/finiteVolume/cellInterpolant.hpp b/src/finiteVolume/cellInterpolant.hpp index 15c8025ac..a18316e5e 100644 --- a/src/finiteVolume/cellInterpolant.hpp +++ b/src/finiteVolume/cellInterpolant.hpp @@ -6,9 +6,13 @@ #include "domain/range.hpp" #include "domain/region.hpp" #include "domain/subDomain.hpp" -namespace ablate::finiteVolume { +//#include "/p/lustre2/kolosret/papi/src/install/include/papi.h" +#include +#include -class CellInterpolant { +namespace ablate::finiteVolume { +constexpr int numLabel = 7; +class CellInterpolant : private utilities::Loggable { public: /** * Function assumes that the left/right solution and aux variables are discontinuous across the interface @@ -43,6 +47,14 @@ class CellInterpolant { std::vector auxFields; }; +// void handle_error (int retval) +// { +// printf("PAPI error %d: %s\n", retval, PAPI_strerror(retval)); +// exit(1); +// } + + double totalTime=0.0; + int callCount=0; private: //! use the subDomain to setup the problem std::shared_ptr subDomain; @@ -53,6 +65,11 @@ class CellInterpolant { // Maximum value for gradients for the multi-direction flux limiter const double maxLimGrad; + + //! Vector to hold all the labels, I dont know the size right now + std::vector flowLabelVec; + + /** * Function to compute the flux source terms */ @@ -67,6 +84,10 @@ class CellInterpolant { void ProjectToFace(const std::vector& fields, PetscDS ds, const PetscFVFaceGeom& faceGeom, PetscInt cellId, const PetscFVCellGeom& cellGeom, DM dm, const PetscScalar* xArray, const std::vector& dmGrads, const std::vector& gradArrays, PetscScalar* u, PetscScalar* grad, bool projectField = true); + void ProjectToFace(const std::vector& fields, PetscDS ds, const PetscFVFaceGeom& faceGeom, PetscInt cellId, const PetscFVCellGeom& cellGeom, DM dm, const PetscScalar* xArray, + const std::vector& dmGrads, const std::vector& gradArrays, PetscScalar* u, PetscScalar* grad, PetscInt neighborCellId , const PetscFVCellGeom& neighborCellGeom, bool projectField = true); + + /** * computes the cell gradients * @param field @@ -102,7 +123,7 @@ class CellInterpolant { * @param faceGeomVec * @param cellGeomVec */ - CellInterpolant(std::shared_ptr subDomain, const std::shared_ptr& solverRegion, Vec faceGeomVec, Vec cellGeomVec, double maxGradIn); + CellInterpolant(std::shared_ptr subDomain, const std::shared_ptr& solverRegion, Vec faceGeomVec, Vec cellGeomVec, double maxGradIn,const ablate::domain::Range& faceRange); ~CellInterpolant(); /** @@ -125,6 +146,35 @@ class CellInterpolant { const ablate::domain::Range& cellRange, Vec cellGeomVec); }; + + +//void get_memory_usage() { +// struct rusage usage; +// if (getrusage(RUSAGE_SELF, &usage) == 0) { +// printf("Max resident set size: %ld KB\n", usage.ru_maxrss); +// } else { +// perror("getrusage failed"); +// } +//} + } // namespace ablate::finiteVolume + +// void get_memory_usage() { +// FILE *file = fopen("/proc/self/statm", "r"); +// if (file == NULL) { +// perror("Error opening /proc/self/statm"); +// return; +// } +// +// long pages; +// if (fscanf(file, "%ld", &pages) == 1) { +// long page_size = sysconf(_SC_PAGESIZE); // Get page size in bytes +// printf("Memory usage: %ld KB\n", (pages * page_size) / 1024); +// } +// +// fclose(file); +// } + + #endif // ABLATELIBRARY_CELLINTERPOLANT_HPP diff --git a/src/finiteVolume/faceInterpolant.cpp b/src/finiteVolume/faceInterpolant.cpp index b9680e25c..5d23a0ffc 100644 --- a/src/finiteVolume/faceInterpolant.cpp +++ b/src/finiteVolume/faceInterpolant.cpp @@ -262,12 +262,27 @@ void ablate::finiteVolume::FaceInterpolant::RestoreInterpolatedFaceVectors(Vec, void ablate::finiteVolume::FaceInterpolant::ComputeRHS(PetscReal time, Vec locXVec, Vec locAuxVec, Vec locFVec, const std::shared_ptr& solverRegion, std::vector& rhsFunctions, const ablate::domain::Range& faceRange, Vec cellGeomVec, Vec faceGeomVec) { +// PetscInt globalSizeX, globalSizeAux; +// VecGetSize(locXVec, &globalSizeX) >> utilities::PetscUtilities::checkError; + if (locAuxVec) { + PetscInt globalSizeAux; + VecGetSize(locAuxVec, &globalSizeAux) >> utilities::PetscUtilities::checkError; + PetscInt globalSizeX; + VecGetSize(locAuxVec, &globalSizeX) >> utilities::PetscUtilities::checkError; + + } + +// std::printf("The length of the locXVec is: %d \n",globalSizeX); +// StartEvent("FiniteVolumeSolver::FaceInterpolant::ComputeRHS::Gradient"); // get the dm auto dm = subDomain->GetDM(); // interpolate to the faces Vec faceSolutionVec, faceAuxVec, faceSolutionGradVec, faceAuxGradVec; GetInterpolatedFaceVectors(locXVec, locAuxVec, faceSolutionVec, faceAuxVec, faceSolutionGradVec, faceAuxGradVec); +// EndEvent(); +// StartEvent("FiniteVolumeSolver::FaceInterpolant::ComputeRHS::Setup"); + // check for ghost cells DMLabel ghostLabel; @@ -351,7 +366,8 @@ void ablate::finiteVolume::FaceInterpolant::ComputeRHS(PetscReal time, Vec locXV } } } - +// EndEvent(); +// StartEvent("FiniteVolumeSolver::FaceInterpolant::ComputeRHS::Sourcecalc"); // march over each face for (PetscInt f = faceRange.start; f < faceRange.end; f++) { PetscInt face = faceRange.points ? faceRange.points[f] : f; @@ -385,9 +401,11 @@ void ablate::finiteVolume::FaceInterpolant::ComputeRHS(PetscReal time, Vec locXV PetscFVFaceGeom* fg; DMPlexPointLocalRead(faceDM, face, faceGeomArray, &fg); +// EndEvent(); // March over each source function for (std::size_t fun = 0; fun < rhsFunctions.size(); fun++) { +// StartEvent("FiniteVolumeSolver::FaceInterpolant::ComputeRHS::diffusion"); PetscArrayzero(flux.data(), totDim) >> utilities::PetscUtilities::checkError; PetscInt fluxOffset = 0; // Flux offset for the function ( Currently calculated by just adding the number of components of the previous fields) const auto& rhsFluxFunctionDescription = rhsFunctions[fun]; @@ -404,6 +422,8 @@ void ablate::finiteVolume::FaceInterpolant::ComputeRHS(PetscReal time, Vec locXV flux.data(), rhsFluxFunctionDescription.context) >> utilities::PetscUtilities::checkError; +// EndEvent(); +// StartEvent("FiniteVolumeSolver::FaceInterpolant::ComputeRHS::Addflux"); for (std::size_t updateFieldIdx = 0; updateFieldIdx < rhsFunctions[fun].updateFields.size(); updateFieldIdx++) { // add the flux back to the cell PetscScalar *fL = nullptr, *fR = nullptr; @@ -431,9 +451,11 @@ void ablate::finiteVolume::FaceInterpolant::ComputeRHS(PetscReal time, Vec locXV } fluxOffset += fluxComponentSize[fun][updateFieldIdx]; } +// EndEvent(); } } - +// EndEvent(); +// StartEvent("FiniteVolumeSolver::FaceInterpolant::ComputeRHS::Cleanup"); VecRestoreArrayRead(faceSolutionVec, &faceSolutionArray); VecRestoreArrayRead(faceSolutionGradVec, &faceSolutionGradArray); @@ -445,4 +467,5 @@ void ablate::finiteVolume::FaceInterpolant::ComputeRHS(PetscReal time, Vec locXV VecRestoreArrayRead(cellGeomVec, &cellGeomArray) >> utilities::PetscUtilities::checkError; VecRestoreArrayRead(faceGeomVec, &faceGeomArray) >> utilities::PetscUtilities::checkError; RestoreInterpolatedFaceVectors(locXVec, locAuxVec, faceSolutionVec, faceAuxVec, faceSolutionGradVec, faceAuxGradVec); +// EndEvent(); } diff --git a/src/finiteVolume/faceInterpolant.hpp b/src/finiteVolume/faceInterpolant.hpp index 0a0cabd9c..28bc5e4c2 100644 --- a/src/finiteVolume/faceInterpolant.hpp +++ b/src/finiteVolume/faceInterpolant.hpp @@ -8,7 +8,7 @@ namespace ablate::finiteVolume { -class FaceInterpolant { +class FaceInterpolant : private utilities::Loggable{ private: //! use the subDomain to setup the problem std::shared_ptr subDomain; diff --git a/src/finiteVolume/finiteVolumeSolver.cpp b/src/finiteVolume/finiteVolumeSolver.cpp index 6c21ed1be..502f7213d 100644 --- a/src/finiteVolume/finiteVolumeSolver.cpp +++ b/src/finiteVolume/finiteVolumeSolver.cpp @@ -235,13 +235,16 @@ PetscErrorCode ablate::finiteVolume::FiniteVolumeSolver::ComputeRHSFunction(Pets ablate::domain::Range faceRange, cellRange; GetFaceRange(faceRange); GetCellRange(cellRange); + +// PAPI_library_init(PAPI_VER_CURRENT); + try { +// MPI_Barrier(PETSC_COMM_WORLD); StartEvent("FiniteVolumeSolver::ComputeRHSFunction::discontinuousFluxFunction"); if (!discontinuousFluxFunctionDescriptions.empty()) { if (cellInterpolant == nullptr) { - cellInterpolant = std::make_unique(subDomain, GetRegion(), faceGeomVec, cellGeomVec, maxlimit); + cellInterpolant = std::make_unique(subDomain, GetRegion(), faceGeomVec, cellGeomVec, maxlimit,faceRange); } - cellInterpolant->ComputeRHS(time, locXVec, subDomain->GetAuxVector(), locFVec, GetRegion(), discontinuousFluxFunctionDescriptions, faceRange, cellRange, cellGeomVec, faceGeomVec); } EndEvent(); @@ -249,14 +252,15 @@ PetscErrorCode ablate::finiteVolume::FiniteVolumeSolver::ComputeRHSFunction(Pets SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in CellInterpolant discontinuousFluxFunction: %s", exception.what()); } + try { StartEvent("FiniteVolumeSolver::ComputeRHSFunction::pointFunction"); if (!pointFunctionDescriptions.empty()) { if (cellInterpolant == nullptr) { - cellInterpolant = std::make_unique(subDomain, GetRegion(), faceGeomVec, cellGeomVec, maxlimit); +// cellInterpolant = std::make_unique(subDomain, GetRegion(), faceGeomVec, cellGeomVec, maxlimit); } - cellInterpolant->ComputeRHS(time, locXVec, subDomain->GetAuxVector(), locFVec, GetRegion(), pointFunctionDescriptions, cellRange, cellGeomVec); +// cellInterpolant->ComputeRHS(time, locXVec, subDomain->GetAuxVector(), locFVec, GetRegion(), pointFunctionDescriptions, cellRange, cellGeomVec); } EndEvent(); } catch (std::exception& exception) { @@ -482,9 +486,28 @@ PetscErrorCode ablate::finiteVolume::FiniteVolumeSolver::PreRHSFunction(TS ts, P PetscCall(rhsFunction.first(*this, ts, time, initialStage, locX, rhsFunction.second)); } EndEvent(); + + StartEvent("FiniteVolumeSolver::PreRHSFunction::FixEnergy"); + try { + // update any aux fields, including ghost cells +// FixEnergy(subDomain->GetAuxVector(),subDomain->GetSolutionVector()); + } catch (std::exception& exception) { + SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in Energy fix: %s", exception.what()); + } + EndEvent(); + PetscFunctionReturn(0); } +//PetscErrorCode FixEnergy(Vec locAuxField,Vec locSolField){ +// +// double a=1; +// +// PetscFunctionReturn(0); +//} + + + #include "registrar.hpp" REGISTER(ablate::solver::Solver, ablate::finiteVolume::FiniteVolumeSolver, "finite volume solver", ARG(std::string, "id", "the name of the flow field"), OPT(ablate::domain::Region, "region", "the region to apply this solver. Default is entire domain"), diff --git a/src/finiteVolume/finiteVolumeSolver.hpp b/src/finiteVolume/finiteVolumeSolver.hpp index 88352e50a..69fabc7bb 100644 --- a/src/finiteVolume/finiteVolumeSolver.hpp +++ b/src/finiteVolume/finiteVolumeSolver.hpp @@ -13,7 +13,11 @@ #include "solver/timeStepper.hpp" #include "utilities/constants.hpp" #include "utilities/vectorUtilities.hpp" - +//#include "/p/lustre2/kolosret/papi/src/install/include/papi.h" +#include +//#include +//#include +#include namespace ablate::finiteVolume { // forward declare the FlowProcesses @@ -35,7 +39,15 @@ class FiniteVolumeSolver : public solver::CellSolver, //! store an enum for the fields in the meshCharacteristicsDm enum MeshCharacteristics { MIN_CELL_RADIUS = 0, MAX_CELL_RADIUS }; +// //Additional post step to fix temperature +// void PostStep(TS ts); + private: + + //Function to limit the energy + PetscErrorCode FixEnergy(Vec locXVec,Vec locSolVec); + + /** * struct to describe the compute timestamp functions */ diff --git a/src/finiteVolume/fluxCalculator/ausmpUp.cpp b/src/finiteVolume/fluxCalculator/ausmpUp.cpp index c372227fa..1c300797a 100644 --- a/src/finiteVolume/fluxCalculator/ausmpUp.cpp +++ b/src/finiteVolume/fluxCalculator/ausmpUp.cpp @@ -1,100 +1,261 @@ #include "ausmpUp.hpp" +double ablate::finiteVolume::fluxCalculator::AusmpUp::totalTime = 0.0; +int ablate::finiteVolume::fluxCalculator::AusmpUp::callCount = 0; + ablate::finiteVolume::fluxCalculator::AusmpUp::AusmpUp(double mInf, std::shared_ptr pgs) : pgs(pgs), mInf(mInf) {} +//ablate::finiteVolume::fluxCalculator::Direction ablate::finiteVolume::fluxCalculator::AusmpUp::AusmpUpFunction(void* ctx, PetscReal uL, PetscReal aL, PetscReal rhoL, PetscReal pL, PetscReal uR, +// PetscReal aR, PetscReal rhoR, PetscReal pR, PetscReal* massFlux, PetscReal* p12) { +// +// // extract pgs/minf if provided +// auto ausmUp = (ablate::finiteVolume::fluxCalculator::AusmpUp*)ctx; +// PetscReal pgsAlpha = ausmUp->pgs ? ausmUp->pgs->GetAlpha() : 1.0; +// PetscReal mInf = ausmUp->mInf; +// +// // Compute the density at the interface +// PetscReal rho12 = (0.5) * (rhoL + rhoR); +// +// // compute the speed of sound at a12 +// PetscReal a12 = 0.5 * (aL + aR) / pgsAlpha; // Simple average of aL and aR. This can be replaced with eq. 30; +// +// // Compute the left and right mach numbers +// PetscReal mL = uL / a12; +// PetscReal mR = uR / a12; +// +// // Compute mBar2 (eq 70) +// PetscReal mBar2 = (PetscSqr(uL) + PetscSqr(uR)) / (2.0 * a12 * a12); +// +// // compute mInf2 or set fa to unity +// PetscReal fa = 1.0; +// if (mInf > 0) { +// PetscReal mInf2 = PetscSqr(mInf); +// +// PetscReal mO2 = PetscMin(1.0, PetscMax(mBar2, mInf2)); +// PetscReal mO = PetscSqrtReal(mO2); +// fa = mO * (2.0 - mO); +// } +// // compute the mach number on the interface +// PetscReal m12 = M4Plus(mL) + M4Minus(mR) - (Kp / fa) * PetscMax(1.0 - (sigma * mBar2), 0) * (pR - pL) / (rho12 * a12 * a12 * pgsAlpha * pgsAlpha); +// +// // store the mass flux; +// Direction direction; +// if (m12 > 0) { +// direction = LEFT; +// *massFlux = a12 * m12 * rhoL; +// } else { +// direction = RIGHT; +// *massFlux = a12 * m12 * rhoR; +// } +// +// // Pressure +// if (p12) { +// double p5Plus = P5Plus(mL, fa); +// double p5Minus = P5Minus(mR, fa); +// +// *p12 = p5Plus * pL + p5Minus * pR - Ku * p5Plus * p5Minus * rho12 * fa * a12 * a12 * pgsAlpha * pgsAlpha * (mR - mL); +// *p12 /= PetscSqr(pgsAlpha); +// } +// return direction; +//} ablate::finiteVolume::fluxCalculator::Direction ablate::finiteVolume::fluxCalculator::AusmpUp::AusmpUpFunction(void* ctx, PetscReal uL, PetscReal aL, PetscReal rhoL, PetscReal pL, PetscReal uR, + PetscReal aR, PetscReal rhoR, PetscReal pR, PetscReal* massFlux, PetscReal* p12) { - // extract pgs/minf if provided - auto ausmUp = (ablate::finiteVolume::fluxCalculator::AusmpUp*)ctx; - PetscReal pgsAlpha = ausmUp->pgs ? ausmUp->pgs->GetAlpha() : 1.0; - PetscReal mInf = ausmUp->mInf; + + + //PAPI high level +// int retval; +// retval = PAPI_hl_region_begin("ausmup"); + +// //PAPI low level +// int EventSet = PAPI_NULL; +// long long values[1]; +// // Create event set +// PAPI_create_eventset(&EventSet); +// PAPI_add_event(EventSet, PAPI_DP_OPS); +// PAPI_start(EventSet); + + + //Timing +// double start = MPI_Wtime(); + + + auto ausmUp = (ablate::finiteVolume::fluxCalculator::AusmpUp*)ctx; // 1 read (ctx), 1 write (ausmUp pointer assignment) + PetscReal pgsAlpha = ausmUp->pgs ? ausmUp->pgs->GetAlpha() : 1.0; // 1 read (ausmUp->pgs), 1 conditional read (GetAlpha), 1 write + PetscReal mInf = ausmUp->mInf; // 1 read, 1 write + // Compute the density at the interface - PetscReal rho12 = (0.5) * (rhoL + rhoR); + PetscReal rho12 = 0.5 * (rhoL + rhoR); // 2 reads (rhoL, rhoR), 1 write // compute the speed of sound at a12 - PetscReal a12 = 0.5 * (aL + aR) / pgsAlpha; // Simple average of aL and aR. This can be replaced with eq. 30; + PetscReal a12 = 0.5 * (aL + aR) / pgsAlpha; // 3 reads (aL, aR, pgsAlpha), 1 write // Compute the left and right mach numbers - PetscReal mL = uL / a12; - PetscReal mR = uR / a12; + PetscReal mL = uL / a12; // 2 reads (uL, a12), 1 write + PetscReal mR = uR / a12; // 2 reads (uR, a12), 1 write + //R 12, W 7 so far // Compute mBar2 (eq 70) - PetscReal mBar2 = (PetscSqr(uL) + PetscSqr(uR)) / (2.0 * a12 * a12); + PetscReal mBar2 = (PetscSqr(uL) + PetscSqr(uR)) / (2.0 * a12 * a12); // 3 reads (uL, uR, a12), 1 write // compute mInf2 or set fa to unity - PetscReal fa = 1.0; - if (mInf > 0) { - PetscReal mInf2 = PetscSqr(mInf); + PetscReal fa = 1.0; // 1 write + if (mInf > 0) { // 1 read (mInf) + PetscReal mInf2 = PetscSqr(mInf); // 1 read, 1 write - PetscReal mO2 = PetscMin(1.0, PetscMax(mBar2, mInf2)); - PetscReal mO = PetscSqrtReal(mO2); - fa = mO * (2.0 - mO); + PetscReal mO2 = PetscMin(1.0, PetscMax(mBar2, mInf2)); // 2 reads (mBar2, mInf2), 1 write + PetscReal mO = PetscSqrtReal(mO2); // 1 read, 1 write + fa = mO * (2.0 - mO); // 1 read, 1 write } + // compute the mach number on the interface PetscReal m12 = M4Plus(mL) + M4Minus(mR) - (Kp / fa) * PetscMax(1.0 - (sigma * mBar2), 0) * (pR - pL) / (rho12 * a12 * a12 * pgsAlpha * pgsAlpha); + // ~12 reads (mL, mR, fa, mBar2, pR, pL, rho12, a12 x2, pgsAlpha x2, Kp, sigma), 1 write // store the mass flux; - Direction direction; - if (m12 > 0) { - direction = LEFT; - *massFlux = a12 * m12 * rhoL; + // doesnt matter which branch for memory + Direction direction; // 1 write + if (m12 > 0) { // 1 read + direction = LEFT; // 1 write + *massFlux = a12 * m12 * rhoL; // 3 reads, 1 write } else { - direction = RIGHT; - *massFlux = a12 * m12 * rhoR; + direction = RIGHT; // 1 write + *massFlux = a12 * m12 * rhoR; // 3 reads, 1 write } + //32 R + 2*M4_R, 13 W+2*M4_R so far + 1 if loop + // Pressure - if (p12) { - double p5Plus = P5Plus(mL, fa); - double p5Minus = P5Minus(mR, fa); + if (p12) { // 1 read + double p5Plus = P5Plus(mL, fa); // 2 reads, 1 write + double p5Minus = P5Minus(mR, fa); // 2 reads, 1 write *p12 = p5Plus * pL + p5Minus * pR - Ku * p5Plus * p5Minus * rho12 * fa * a12 * a12 * pgsAlpha * pgsAlpha * (mR - mL); - *p12 /= PetscSqr(pgsAlpha); + // ~11 reads (p5Plus, p5Minus, pL, pR, Ku, rho12, fa, a12 x2, pgsAlpha x2, mR, mL), 1 write + *p12 /= PetscSqr(pgsAlpha); // 2 reads (p12, pgsAlpha), 1 write } - return direction; + + +// PAPI low level +// PAPI_stop(EventSet, values); +// printf("FLOPs counted: %lld\n", values[0]); + + + //PAPI high level +// retval = PAPI_hl_region_end("ausmup"); +// if ( retval != PAPI_OK ){ +// std::cout << "abc" << std::endl; +// } + + + + //Timing function +// totalTime += MPI_Wtime() - start; +// ++callCount; +// if (callCount == 100000) { +// PetscPrintf( +// PETSC_COMM_WORLD, "StaticFunction called %d times, total time: %f s\n", ablate::finiteVolume::fluxCalculator::AusmpUp::callCount, ablate::finiteVolume::fluxCalculator::AusmpUp::totalTime); +// } + return direction; // 1 read } -PetscReal ablate::finiteVolume::fluxCalculator::AusmpUp::M1Plus(PetscReal m) { return 0.5 * (m + PetscAbs(m)); } +PetscReal ablate::finiteVolume::fluxCalculator::AusmpUp::M1Plus(PetscReal m) { + return 0.5 * (m + PetscAbs(m)); // 2 reads (m), 1 write (return) +} + +PetscReal ablate::finiteVolume::fluxCalculator::AusmpUp::M2Plus(PetscReal m) { + return 0.25 * PetscSqr(m + 1); // 1 read (m), 1 write (return) +} -PetscReal ablate::finiteVolume::fluxCalculator::AusmpUp::M2Plus(PetscReal m) { return 0.25 * PetscSqr(m + 1); } +PetscReal ablate::finiteVolume::fluxCalculator::AusmpUp::M1Minus(PetscReal m) { + return 0.5 * (m - PetscAbs(m)); // 2 reads (m), 1 write (return) +} -PetscReal ablate::finiteVolume::fluxCalculator::AusmpUp::M1Minus(PetscReal m) { return 0.5 * (m - PetscAbs(m)); } -PetscReal ablate::finiteVolume::fluxCalculator::AusmpUp::M2Minus(PetscReal m) { return -0.25 * PetscSqr(m - 1); } +PetscReal ablate::finiteVolume::fluxCalculator::AusmpUp::M2Minus(PetscReal m) { + return -0.25 * PetscSqr(m - 1); // 1 read (m), 1 write (return) +} PetscReal ablate::finiteVolume::fluxCalculator::AusmpUp::M4Plus(PetscReal m) { - if (PetscAbs(m) >= 1.0) { - return M1Plus(m); + if (PetscAbs(m) >= 1.0) { // 1 read (m) + return M1Plus(m); // 1 read (m), 1 write (return) } else { return M2Plus(m) * (1.0 - 16.0 * beta * M2Minus(m)); + // 2 reads (m), 1 read (beta), 1 write (return) } } + PetscReal ablate::finiteVolume::fluxCalculator::AusmpUp::M4Minus(PetscReal m) { - if (PetscAbs(m) >= 1.0) { - return M1Minus(m); + if (PetscAbs(m) >= 1.0) { // 1 read (m) + return M1Minus(m); // 1 read (m), 1 write (return) } else { return M2Minus(m) * (1.0 + 16.0 * beta * M2Plus(m)); + // 2 reads (m), 1 read (beta), 1 write (return) } } + PetscReal ablate::finiteVolume::fluxCalculator::AusmpUp::P5Plus(PetscReal m, double fa) { - if (PetscAbs(m) >= 1.0) { - return (M1Plus(m) / (m + 1E-30)); + if (PetscAbs(m) >= 1.0) { // 1 read (m) + return (M1Plus(m) / (m + 1E-30)); // 2 reads (m), 1 write (return) } else { - // compute alpha - double alpha = 3.0 / 16.0 * (-4.0 + 5 * fa * fa); + double alpha = 3.0 / 16.0 * (-4.0 + 5 * fa * fa); // 2 reads (fa), 1 write (alpha) return (M2Plus(m) * ((2.0 - m) - 16. * alpha * m * M2Minus(m))); + // 3 reads (m, alpha, fa), 1 write (return) } } + PetscReal ablate::finiteVolume::fluxCalculator::AusmpUp::P5Minus(PetscReal m, double fa) { - if (PetscAbs(m) >= 1.0) { - return (M1Minus(m) / (m + 1E-30)); + if (PetscAbs(m) >= 1.0) { // 1 read (m) + return (M1Minus(m) / (m + 1E-30)); // 2 reads (m), 1 write (return) } else { - double alpha = 3.0 / 16.0 * (-4.0 + 5 * fa * fa); + double alpha = 3.0 / 16.0 * (-4.0 + 5 * fa * fa); // 2 reads (fa), 1 write (alpha) return (M2Minus(m) * ((-2.0 - m) + 16. * alpha * m * M2Plus(m))); + // 3 reads (m, alpha, fa), 1 write (return) } } +//PetscReal ablate::finiteVolume::fluxCalculator::AusmpUp::M1Plus(PetscReal m) { return 0.5 * (m + PetscAbs(m)); } +// +//PetscReal ablate::finiteVolume::fluxCalculator::AusmpUp::M2Plus(PetscReal m) { return 0.25 * PetscSqr(m + 1); } +// +//PetscReal ablate::finiteVolume::fluxCalculator::AusmpUp::M1Minus(PetscReal m) { return 0.5 * (m - PetscAbs(m)); } +//PetscReal ablate::finiteVolume::fluxCalculator::AusmpUp::M2Minus(PetscReal m) { return -0.25 * PetscSqr(m - 1); } +// +//PetscReal ablate::finiteVolume::fluxCalculator::AusmpUp::M4Plus(PetscReal m) { +// if (PetscAbs(m) >= 1.0) { +// return M1Plus(m); +// } else { +// return M2Plus(m) * (1.0 - 16.0 * beta * M2Minus(m)); +// } +//} +//PetscReal ablate::finiteVolume::fluxCalculator::AusmpUp::M4Minus(PetscReal m) { +// if (PetscAbs(m) >= 1.0) { +// return M1Minus(m); +// } else { +// return M2Minus(m) * (1.0 + 16.0 * beta * M2Plus(m)); +// } +//} +//PetscReal ablate::finiteVolume::fluxCalculator::AusmpUp::P5Plus(PetscReal m, double fa) { +// if (PetscAbs(m) >= 1.0) { +// return (M1Plus(m) / (m + 1E-30)); +// } else { +// // compute alpha +// double alpha = 3.0 / 16.0 * (-4.0 + 5 * fa * fa); +// +// return (M2Plus(m) * ((2.0 - m) - 16. * alpha * m * M2Minus(m))); +// } +//} +//PetscReal ablate::finiteVolume::fluxCalculator::AusmpUp::P5Minus(PetscReal m, double fa) { +// if (PetscAbs(m) >= 1.0) { +// return (M1Minus(m) / (m + 1E-30)); +// } else { +// double alpha = 3.0 / 16.0 * (-4.0 + 5 * fa * fa); +// return (M2Minus(m) * ((-2.0 - m) + 16. * alpha * m * M2Plus(m))); +// } +//} + + #include "registrar.hpp" REGISTER(ablate::finiteVolume::fluxCalculator::FluxCalculator, ablate::finiteVolume::fluxCalculator::AusmpUp, "A sequel to AUSM, Part II: AUSM+-up for all speeds, Meng-Sing Liou, Pages 137-170, 2006", OPT(double, "mInf", "the reference mach number"), diff --git a/src/finiteVolume/fluxCalculator/ausmpUp.hpp b/src/finiteVolume/fluxCalculator/ausmpUp.hpp index 47b5d4bda..5994a4b0c 100644 --- a/src/finiteVolume/fluxCalculator/ausmpUp.hpp +++ b/src/finiteVolume/fluxCalculator/ausmpUp.hpp @@ -8,11 +8,12 @@ namespace ablate::finiteVolume::fluxCalculator { /** * A sequel to AUSM, Part II: AUSM+-up for all speeds */ -class AusmpUp : public fluxCalculator::FluxCalculator { +class AusmpUp : public fluxCalculator::FluxCalculator,private utilities::Loggable { private: // AusmUp uses a pgs if provided const std::shared_ptr pgs; + static Direction AusmpUpFunction(void*, PetscReal uL, PetscReal aL, PetscReal rhoL, PetscReal pL, PetscReal uR, PetscReal aR, PetscReal rhoR, PetscReal pR, PetscReal* massFlux, PetscReal* p12); static PetscReal M1Plus(PetscReal m); @@ -24,10 +25,17 @@ class AusmpUp : public fluxCalculator::FluxCalculator { const inline static PetscReal Ku = 0.75; const inline static PetscReal sigma = 0.25; + + + // The reference infinity mach number const double mInf; public: + static double totalTime; + static int callCount; + + explicit AusmpUp(double mInf, std::shared_ptr = {}); AusmpUp(AusmpUp const&) = delete; AusmpUp& operator=(AusmpUp const&) = delete; diff --git a/src/finiteVolume/processes/compactCompressibleNSSpeciesTransport.cpp b/src/finiteVolume/processes/compactCompressibleNSSpeciesTransport.cpp index bff9a76a2..48fa929b7 100644 --- a/src/finiteVolume/processes/compactCompressibleNSSpeciesTransport.cpp +++ b/src/finiteVolume/processes/compactCompressibleNSSpeciesTransport.cpp @@ -8,6 +8,17 @@ #include "utilities/constants.hpp" #include "utilities/mathUtilities.hpp" #include "utilities/petscUtilities.hpp" +double ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::totalTimeDiff = 0.0; +int ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::callCountDiff = 0; +double ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::totalTimeDiffEner = 0.0; +int ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::callCountDiffEner = 0; +double ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::totalTimeDiffSpec = 0.0; +int ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::callCountDiffSpec = 0; + +double ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::totalTimeState = 0.0; +int ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::callCountState = 0; +double ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::totalTimeAUSM = 0.0; +int ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::callCountAUSM = 0; ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::CompactCompressibleNSSpeciesTransport(const std::shared_ptr& parametersIn, std::shared_ptr eosIn, @@ -163,16 +174,15 @@ PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTran const PetscScalar* fieldR, const PetscInt* aOff, const PetscScalar* auxL, const PetscScalar* auxR, PetscScalar* flux, void* ctx) { PetscFunctionBeginUser; + auto advectionData = (AdvectionData*)ctx; // 1 read, 1 write - auto advectionData = (AdvectionData*)ctx; - - const int EULER_FIELD = 0; - const int RHOYI_FIELD = 1; + const int EULER_FIELD = 0; // 1 write + const int RHOYI_FIELD = 1; // 1 write // Compute the norm PetscReal norm[3]; - utilities::MathUtilities::NormVector(dim, fg->normal, norm); - const PetscReal areaMag = utilities::MathUtilities::MagVector(dim, fg->normal); + utilities::MathUtilities::NormVector(dim, fg->normal, norm); // D reads (fg->normal), D writes (norm) + const PetscReal areaMag = utilities::MathUtilities::MagVector(dim, fg->normal); // 3*D + 1 reads , D writes // Decode the left and right states PetscReal densityL; @@ -181,23 +191,42 @@ PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTran PetscReal internalEnergyL; PetscReal aL; PetscReal pL; - +// double start = MPI_Wtime(); // decode the left side { - densityL = fieldL[uOff[EULER_FIELD] + CompressibleFlowFields::RHO]; + densityL = fieldL[uOff[EULER_FIELD] + CompressibleFlowFields::RHO]; // 2 reads, 1 write PetscReal temperatureL; + PetscCall(advectionData->computeTemperature.function(fieldL, auxL[aOff[0]] * .67 + .33 * auxR[aOff[0]], &temperatureL, advectionData->computeTemperature.context.get())); + // State evaluation: Temperature, Energy, a, P + if(std::isnan(temperatureL)) { + std::cout << "temperatureL is NaN." << std::endl; + temperatureL=273; + } // Get the velocity in this direction - normalVelocityL = 0.0; - for (PetscInt d = 0; d < dim; d++) { - velocityL[d] = fieldL[uOff[EULER_FIELD] + CompressibleFlowFields::RHOU + d] / densityL; - normalVelocityL += velocityL[d] * norm[d]; + normalVelocityL = 0.0; // 0 reads, 1 write + for (PetscInt d = 0; d < dim; d++) { // 1 reads, 1 write + velocityL[d] = fieldL[uOff[EULER_FIELD] + CompressibleFlowFields::RHOU + d] / densityL; // 3 reads, 1 write + normalVelocityL += velocityL[d] * norm[d]; // 2 reads, 1 write } PetscCall(advectionData->computeInternalEnergy.function(fieldL, temperatureL, &internalEnergyL, advectionData->computeInternalEnergy.context.get())); PetscCall(advectionData->computeSpeedOfSound.function(fieldL, temperatureL, &aL, advectionData->computeSpeedOfSound.context.get())); PetscCall(advectionData->computePressure.function(fieldL, temperatureL, &pL, advectionData->computePressure.context.get())); + +// double yiSumL=0; +// for (PetscInt sp = 0; sp < advectionData->numberSpecies; sp++) { +// // Limit the bounds +// PetscScalar yi = fieldL[uOff[RHOYI_FIELD] + sp] / densityL; +//// yi = PetscMax(0.0, yi); +//// yi = PetscMin(1.0, yi); +// yiSumL += yi; +// +// } +// printf("The sum of mass fractions from the LEFT state is %f \n", yiSumL); + + } PetscReal densityR; @@ -207,41 +236,65 @@ PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTran PetscReal aR; PetscReal pR; + //In terms of memory same as left { // decode right state densityR = fieldR[uOff[EULER_FIELD] + CompressibleFlowFields::RHO]; PetscReal temperatureR; PetscCall(advectionData->computeTemperature.function(fieldR, auxR[aOff[0]] * .67 + .33 * auxL[aOff[0]], &temperatureR, advectionData->computeTemperature.context.get())); - + if(std::isnan(temperatureR)) { + std::cout << "temperatureR is NaN." << std::endl; + temperatureR=273; + } // Get the velocity in this direction normalVelocityR = 0.0; for (PetscInt d = 0; d < dim; d++) { velocityR[d] = fieldR[uOff[EULER_FIELD] + CompressibleFlowFields::RHOU + d] / densityR; normalVelocityR += velocityR[d] * norm[d]; } +// double yiSumR=0; +// for (PetscInt sp = 0; sp < advectionData->numberSpecies; sp++) { +// // Limit the bounds +// PetscScalar yi = fieldR[uOff[RHOYI_FIELD] + sp] / densityR; +//// yi = PetscMax(0.0, yi); +//// yi = PetscMin(1.0, yi); +// yiSumR += yi; +// +// } +// printf("The sum of mass fractions from the RIGHT state is %f \n", yiSumR); PetscCall(advectionData->computeInternalEnergy.function(fieldR, temperatureR, &internalEnergyR, advectionData->computeInternalEnergy.context.get())); PetscCall(advectionData->computeSpeedOfSound.function(fieldR, temperatureR, &aR, advectionData->computeSpeedOfSound.context.get())); PetscCall(advectionData->computePressure.function(fieldR, temperatureR, &pR, advectionData->computePressure.context.get())); } +// totalTimeState += MPI_Wtime() - start; +// ++callCountState; +// if (callCountState==1000000) { +// PetscPrintf(PETSC_COMM_WORLD, "State total %d time: %f s\n", callCountState, totalTimeState); } + // get the face values PetscReal massFlux; PetscReal p12; +// start = MPI_Wtime(); + fluxCalculator::Direction direction = advectionData->fluxCalculatorFunction(advectionData->fluxCalculatorCtx, normalVelocityL, aL, densityL, pL, normalVelocityR, aR, densityR, pR, &massFlux, &p12); if (direction == fluxCalculator::LEFT) { - flux[CompressibleFlowFields::RHO] = massFlux * areaMag; - PetscReal velMagL = utilities::MathUtilities::MagVector(dim, velocityL); - PetscReal HL = internalEnergyL + velMagL * velMagL / 2.0 + pL / densityL; - flux[CompressibleFlowFields::RHOE] = HL * massFlux * areaMag; - for (PetscInt n = 0; n < dim; n++) { - flux[CompressibleFlowFields::RHOU + n] = velocityL[n] * massFlux * areaMag + p12 * fg->normal[n]; + flux[CompressibleFlowFields::RHO] = massFlux * areaMag; // 3 reads, 1 write + PetscReal velMagL = utilities::MathUtilities::MagVector(dim, velocityL); // 3*D + 1 reads , D writes + PetscReal HL = internalEnergyL + velMagL * velMagL / 2.0 + pL / densityL; // 5 reads, 1 write + flux[CompressibleFlowFields::RHOE] = HL * massFlux * areaMag; // 4 reads, 1 write + for (PetscInt n = 0; n < dim; n++) { // 1 reads, 1 write + flux[CompressibleFlowFields::RHOU + n] = velocityL[n] * massFlux * areaMag + p12 * fg->normal[n]; // 7 reads, 1 write + } + for (PetscInt ns = 0; ns < advectionData->numberSpecies; ns++) { // 1 reads, 1 write + flux[uOff[RHOYI_FIELD] + ns] = massFlux * fieldL[uOff[RHOYI_FIELD] + ns] / densityL * areaMag; // 7 reads, 1 write } - for (PetscInt ns = 0; ns < advectionData->numberSpecies; ns++) flux[uOff[RHOYI_FIELD] + ns] = massFlux * fieldL[uOff[RHOYI_FIELD] + ns] / densityL * areaMag; } else if (direction == fluxCalculator::RIGHT) { + flux[CompressibleFlowFields::RHO] = massFlux * areaMag; PetscReal velMagR = utilities::MathUtilities::MagVector(dim, velocityR); PetscReal HR = internalEnergyR + velMagR * velMagR / 2.0 + pR / densityR; @@ -249,7 +302,10 @@ PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTran for (PetscInt n = 0; n < dim; n++) { flux[CompressibleFlowFields::RHOU + n] = velocityR[n] * massFlux * areaMag + p12 * fg->normal[n]; } - for (PetscInt ns = 0; ns < advectionData->numberSpecies; ns++) flux[uOff[RHOYI_FIELD] + ns] = massFlux * fieldR[uOff[RHOYI_FIELD] + ns] / densityR * areaMag; + for (PetscInt ns = 0; ns < advectionData->numberSpecies; ns++){ + flux[uOff[RHOYI_FIELD] + ns] = massFlux * fieldR[uOff[RHOYI_FIELD] + ns] / densityR * areaMag; + + } } else { flux[CompressibleFlowFields::RHO] = massFlux * areaMag; @@ -267,9 +323,127 @@ PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTran flux[uOff[RHOYI_FIELD] + ns] = massFlux * 0.5 * (fieldR[uOff[RHOYI_FIELD] + ns] + fieldL[uOff[RHOYI_FIELD] + ns]) / (0.5 * (densityL + densityR)) * areaMag; } +// totalTimeAUSM += MPI_Wtime() - start; +// ++callCountAUSM; +// if (callCountAUSM==1000000) { +// PetscPrintf(PETSC_COMM_WORLD, "AUSM total %d time: %f s\n", callCountAUSM, totalTimeAUSM); } + PetscFunctionReturn(0); } + + + +//PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::AdvectionFlux(PetscInt dim, const PetscFVFaceGeom* fg, const PetscInt* uOff, const PetscScalar* fieldL, +// const PetscScalar* fieldR, const PetscInt* aOff, const PetscScalar* auxL, const PetscScalar* auxR, +// PetscScalar* flux, void* ctx) { +// PetscFunctionBeginUser; +// auto advectionData = (AdvectionData*)ctx; +// +// const int EULER_FIELD = 0; +// const int RHOYI_FIELD = 1; +// +// // Compute the norm +// PetscReal norm[3]; +// utilities::MathUtilities::NormVector(dim, fg->normal, norm); +// const PetscReal areaMag = utilities::MathUtilities::MagVector(dim, fg->normal); +// +// // Decode the left and right states +// PetscReal densityL; +// PetscReal normalVelocityL; +// PetscReal velocityL[3]; +// PetscReal internalEnergyL; +// PetscReal aL; +// PetscReal pL; +// +// // decode the left side +// { +// densityL = fieldL[uOff[EULER_FIELD] + CompressibleFlowFields::RHO]; +// PetscReal temperatureL; +// +// PetscCall(advectionData->computeTemperature.function(fieldL, auxL[aOff[0]] * .67 + .33 * auxR[aOff[0]], &temperatureL, advectionData->computeTemperature.context.get())); +// +// // Get the velocity in this direction +// normalVelocityL = 0.0; +// for (PetscInt d = 0; d < dim; d++) { +// velocityL[d] = fieldL[uOff[EULER_FIELD] + CompressibleFlowFields::RHOU + d] / densityL; +// normalVelocityL += velocityL[d] * norm[d]; +// } +// PetscCall(advectionData->computeInternalEnergy.function(fieldL, temperatureL, &internalEnergyL, advectionData->computeInternalEnergy.context.get())); +// PetscCall(advectionData->computeSpeedOfSound.function(fieldL, temperatureL, &aL, advectionData->computeSpeedOfSound.context.get())); +// PetscCall(advectionData->computePressure.function(fieldL, temperatureL, &pL, advectionData->computePressure.context.get())); +// } +// +// PetscReal densityR; +// PetscReal normalVelocityR; +// PetscReal velocityR[3]; +// PetscReal internalEnergyR; +// PetscReal aR; +// PetscReal pR; +// +// { // decode right state +// densityR = fieldR[uOff[EULER_FIELD] + CompressibleFlowFields::RHO]; +// PetscReal temperatureR; +// +// PetscCall(advectionData->computeTemperature.function(fieldR, auxR[aOff[0]] * .67 + .33 * auxL[aOff[0]], &temperatureR, advectionData->computeTemperature.context.get())); +// +// // Get the velocity in this direction +// normalVelocityR = 0.0; +// for (PetscInt d = 0; d < dim; d++) { +// velocityR[d] = fieldR[uOff[EULER_FIELD] + CompressibleFlowFields::RHOU + d] / densityR; +// normalVelocityR += velocityR[d] * norm[d]; +// } +// +// PetscCall(advectionData->computeInternalEnergy.function(fieldR, temperatureR, &internalEnergyR, advectionData->computeInternalEnergy.context.get())); +// PetscCall(advectionData->computeSpeedOfSound.function(fieldR, temperatureR, &aR, advectionData->computeSpeedOfSound.context.get())); +// PetscCall(advectionData->computePressure.function(fieldR, temperatureR, &pR, advectionData->computePressure.context.get())); +// } +// +// // get the face values +// PetscReal massFlux; +// PetscReal p12; +// +// fluxCalculator::Direction direction = +// advectionData->fluxCalculatorFunction(advectionData->fluxCalculatorCtx, normalVelocityL, aL, densityL, pL, normalVelocityR, aR, densityR, pR, &massFlux, &p12); +// +// if (direction == fluxCalculator::LEFT) { +// flux[CompressibleFlowFields::RHO] = massFlux * areaMag; +// PetscReal velMagL = utilities::MathUtilities::MagVector(dim, velocityL); +// PetscReal HL = internalEnergyL + velMagL * velMagL / 2.0 + pL / densityL; +// flux[CompressibleFlowFields::RHOE] = HL * massFlux * areaMag; +// for (PetscInt n = 0; n < dim; n++) { +// flux[CompressibleFlowFields::RHOU + n] = velocityL[n] * massFlux * areaMag + p12 * fg->normal[n]; +// } +// for (PetscInt ns = 0; ns < advectionData->numberSpecies; ns++) flux[uOff[RHOYI_FIELD] + ns] = massFlux * fieldL[uOff[RHOYI_FIELD] + ns] / densityL * areaMag; +// } else if (direction == fluxCalculator::RIGHT) { +// flux[CompressibleFlowFields::RHO] = massFlux * areaMag; +// PetscReal velMagR = utilities::MathUtilities::MagVector(dim, velocityR); +// PetscReal HR = internalEnergyR + velMagR * velMagR / 2.0 + pR / densityR; +// flux[CompressibleFlowFields::RHOE] = HR * massFlux * areaMag; +// for (PetscInt n = 0; n < dim; n++) { +// flux[CompressibleFlowFields::RHOU + n] = velocityR[n] * massFlux * areaMag + p12 * fg->normal[n]; +// } +// for (PetscInt ns = 0; ns < advectionData->numberSpecies; ns++) flux[uOff[RHOYI_FIELD] + ns] = massFlux * fieldR[uOff[RHOYI_FIELD] + ns] / densityR * areaMag; +// } else { +// flux[CompressibleFlowFields::RHO] = massFlux * areaMag; +// +// PetscReal velMagL = utilities::MathUtilities::MagVector(dim, velocityL); +// PetscReal HL = internalEnergyL + velMagL * velMagL / 2.0 + pL / densityL; +// +// PetscReal velMagR = utilities::MathUtilities::MagVector(dim, velocityR); +// PetscReal HR = internalEnergyR + velMagR * velMagR / 2.0 + pR / densityR; +// +// flux[CompressibleFlowFields::RHOE] = 0.5 * (HL + HR) * massFlux * areaMag; +// for (PetscInt n = 0; n < dim; n++) { +// flux[CompressibleFlowFields::RHOU + n] = 0.5 * (velocityL[n] + velocityR[n]) * massFlux * areaMag + p12 * fg->normal[n]; +// } +// for (PetscInt ns = 0; ns < advectionData->numberSpecies; ns++) +// flux[uOff[RHOYI_FIELD] + ns] = massFlux * 0.5 * (fieldR[uOff[RHOYI_FIELD] + ns] + fieldL[uOff[RHOYI_FIELD] + ns]) / (0.5 * (densityL + densityR)) * areaMag; +// } +// +// PetscFunctionReturn(0); +//} + double ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport::ComputeCflTimeStep(TS ts, ablate::finiteVolume::FiniteVolumeSolver& flow, void* ctx) { // Get the dm and current solution vector DM dm; @@ -587,6 +761,7 @@ PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTran const PetscInt aOff_x[], const PetscScalar aux[], const PetscScalar gradAux[], PetscScalar flux[], void* ctx) { PetscFunctionBeginUser; +// double start = MPI_Wtime(); // this order is based upon the order that they are passed into RegisterRHSFunction const int T = 0; const int VEL = 1; @@ -636,6 +811,11 @@ PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTran // zero out the density flux flux[CompressibleFlowFields::RHO] = 0.0; +// totalTimeDiff += MPI_Wtime() - start; +// ++callCountDiff; +// if (callCountDiff==1000000) { +// PetscPrintf(PETSC_COMM_WORLD, "Diff total %d time: %f s\n", callCountDiff, totalTimeDiff); } + PetscFunctionReturn(0); } @@ -644,6 +824,7 @@ PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTran const PetscInt aOff_x[], const PetscScalar aux[], const PetscScalar gradAux[], PetscScalar flux[], void* ctx) { PetscFunctionBeginUser; +// double start = MPI_Wtime(); // this order is based upon the order that they are passed into RegisterRHSFunction const int yi = 0; const int euler = 0; @@ -678,6 +859,10 @@ PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTran flux[CompressibleFlowFields::RHOE] += speciesFlux; } } +// totalTimeDiffEner += MPI_Wtime() - start; +// ++callCountDiffEner; +// if (callCountDiffEner==1000000) { +// PetscPrintf(PETSC_COMM_WORLD, "DiffEner total %d time: %f s\n", callCountDiffEner, totalTimeDiffEner); } PetscFunctionReturn(0); } @@ -730,6 +915,7 @@ PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTran const PetscInt aOff_x[], const PetscScalar aux[], const PetscScalar gradAux[], PetscScalar flux[], void* ctx) { PetscFunctionBeginUser; +// double start = MPI_Wtime(); // this order is based upon the order that they are passed into RegisterRHSFunction const int yi = 0; const int euler = 0; @@ -756,7 +942,12 @@ PetscErrorCode ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTran flux[sp] += speciesFlux; } } - +// totalTimeDiffSpec += MPI_Wtime() - start; +// ++callCountDiffSpec; +// int rank; +// MPI_Comm_rank(PETSC_COMM_WORLD, &rank); +// if (callCountDiffSpec==1000000 && rank == 0) { +// PetscPrintf(PETSC_COMM_WORLD, "Diffspecies total %d time: %f s on rank: %d\n", callCountDiffSpec, totalTimeDiffSpec,rank); } PetscFunctionReturn(0); } diff --git a/src/finiteVolume/processes/compactCompressibleNSSpeciesTransport.hpp b/src/finiteVolume/processes/compactCompressibleNSSpeciesTransport.hpp index 91cdb47ac..504dddebc 100644 --- a/src/finiteVolume/processes/compactCompressibleNSSpeciesTransport.hpp +++ b/src/finiteVolume/processes/compactCompressibleNSSpeciesTransport.hpp @@ -11,6 +11,17 @@ namespace ablate::finiteVolume::processes { class CompactCompressibleNSSpeciesTransport : public FlowProcess { private: + static double totalTimeDiff; + static int callCountDiff; + static double totalTimeDiffEner; + static int callCountDiffEner; + static double totalTimeDiffSpec; + static int callCountDiffSpec; + static double totalTimeState; + static int callCountState; + static double totalTimeAUSM; + static int callCountAUSM; + // store ctx needed for function advection function that is passed into Petsc struct AdvectionData { // flow CFL