Petsc profile

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)
 

config/findZerork.cmake

@@ -75,13 +75,27 @@ 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")
     target_link_libraries(zerork_cfd_plugin INTERFACE "$ENV{ZERORK_DIR}/lib/libzerork_cfd_plugin.so")
 
     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 ()
 

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);
 }
 

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<PetscReal> 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<IsothermalWall*>(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);
 }

src/boundarySolver/lodi/isothermalWall.hpp

@@ -9,13 +9,17 @@ class IsothermalWall : public LODIBoundary {
     explicit IsothermalWall(std::shared_ptr<eos::EOS> eos, std::shared_ptr<finiteVolume::processes::PressureGradientScaling> 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,
                                                  const PetscScalar* stencilAuxValues[], PetscInt stencilSize, const PetscInt stencil[], const PetscScalar stencilWeights[], const PetscInt sOff[],
                                                  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);
 };

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);
 }

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::EOS> eos, std::shared_ptr<finiteVolume::processes::PressureGradientScaling> pressureGradientScaling = {});

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<double> 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<double> 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]);