Surface absorptivity addition (#509)

* Added the option for a surface absorptivity which defaults to the specified emissivity if not specified. A critical impact of this change is that the surface radiation emissivities now have to be specified as an input parameter and will need to be altered on all pre-existing input files using surface emission, or else it will default to 1!

* Need to change the regression test yaml inputs to account for this

* PETSC API and formatting changes.

* Bumped Version. Also included in this commit is the Fix of the variable change file as the variable and valid changes names were swapped.

CMakeLists.txt

@@ -4,7 +4,7 @@ cmake_minimum_required(VERSION 3.18.4)
 include(config/petscCompilers.cmake)
 
 # Set the project details
-project(ablateLibrary VERSION 0.12.28)
+project(ablateLibrary VERSION 0.12.29)
 
 # Load the Required 3rd Party Libaries
 pkg_check_modules(PETSc REQUIRED IMPORTED_TARGET GLOBAL PETSc)

src/boundarySolver/physics/sublimation.cpp

@@ -13,7 +13,7 @@ ablate::boundarySolver::physics::Sublimation::Sublimation(std::shared_ptr<subMod
                                                           std::shared_ptr<mathFunctions::MathFunction> additionalHeatFlux,
                                                           std::shared_ptr<finiteVolume::processes::PressureGradientScaling> pressureGradientScaling, bool diffusionFlame,
                                                           std::shared_ptr<ablate::radiation::SurfaceRadiation> radiationIn, const std::shared_ptr<io::interval::Interval> &intervalIn,
-                                                          double emissivityIn)
+                                                          const std::shared_ptr<ablate::parameters::Parameters> &parametersIn)
     : transportModel(std::move(transportModel)),
       eos(std::move(eos)),
       additionalHeatFlux(std::move(additionalHeatFlux)),
@@ -24,7 +24,8 @@ ablate::boundarySolver::physics::Sublimation::Sublimation(std::shared_ptr<subMod
       pressureGradientScaling(std::move(pressureGradientScaling)),
       radiation(std::move(radiationIn)),
       radiationInterval((intervalIn ? intervalIn : std::make_shared<io::interval::FixedInterval>())),
-      emissivity(emissivityIn == 0 ? 1.0 : emissivityIn),
+      emissivity(parametersIn ? parametersIn->Get<PetscReal>("emissivity", 1.) : 1.),
+      absorptivity((parametersIn && parametersIn->GetString("absorptivity")) ? parametersIn->Get<PetscReal>("absorptivity", -1.) : emissivity),
       sublimationModel(std::move(sublimationModel)) {}
 
 void ablate::boundarySolver::physics::Sublimation::Setup(ablate::boundarySolver::BoundarySolver &bSolver) {
@@ -165,7 +166,7 @@ PetscErrorCode ablate::boundarySolver::physics::Sublimation::SublimationFunction
     // Use the solution from the radiation solve.
     PetscReal radIntensity = 0;
     if (sublimation->radiation) {
-        sublimation->radiation->GetSurfaceIntensity(&radIntensity, fg->faceId, boundaryTemperature, sublimation->emissivity);
+        sublimation->radiation->GetSurfaceIntensity(&radIntensity, fg->faceId, boundaryTemperature, sublimation->emissivity, sublimation->absorptivity);
     }
 
     // compute the heat flux. Add the radiation heat flux for this face intensity if the radiation solver exists
@@ -313,7 +314,7 @@ PetscErrorCode ablate::boundarySolver::physics::Sublimation::SublimationOutputFu
 
     PetscReal radIntensity = 0;
     if (sublimation->radiation) {
-        sublimation->radiation->GetSurfaceIntensity(&radIntensity, fg->faceId, boundaryTemperature, sublimation->emissivity);
+        sublimation->radiation->GetSurfaceIntensity(&radIntensity, fg->faceId, boundaryTemperature, sublimation->emissivity, sublimation->absorptivity);
     }
 
     // compute the heat flux
@@ -455,7 +456,7 @@ PetscErrorCode ablate::boundarySolver::physics::Sublimation::UpdateBoundaryHeatT
     // Use the solution from the radiation solve.
     PetscReal radIntensity = 0;
     if (sublimation->radiation) {
-        sublimation->radiation->GetSurfaceIntensity(&radIntensity, fg->faceId, boundaryTemperature, sublimation->emissivity);
+        sublimation->radiation->GetSurfaceIntensity(&radIntensity, fg->faceId, boundaryTemperature, sublimation->emissivity, sublimation->absorptivity);
     }
 
     // compute the heat flux. Add the radiation heat flux for this face intensity if the radiation solver exists
@@ -487,4 +488,5 @@ REGISTER(ablate::boundarySolver::BoundaryProcess, ablate::boundarySolver::physic
          OPT(ablate::finiteVolume::processes::PressureGradientScaling, "pgs", "Pressure gradient scaling is used to scale the acoustic propagation speed and increase time step for low speed flows"),
          OPT(bool, "diffusionFlame", "disables contribution to the momentum equation. Should be true when advection is not solved. (Default is false)"),
          OPT(ablate::radiation::SurfaceRadiation, "radiation", "radiation instance for the sublimation solver to calculate heat flux"),
-         OPT(ablate::io::interval::Interval, "radiationInterval", "number of time steps between the radiation solves"), OPT(double, "emissivity", "radiation property of the fuel surface"));
+         OPT(ablate::io::interval::Interval, "radiationInterval", "number of time steps between the radiation solves"),
+         OPT(ablate::parameters::Parameters, "parameters", "emissivity and absorptivity parameters for the sublimation energy balance"));

src/boundarySolver/physics/sublimation.hpp

@@ -78,10 +78,15 @@ class Sublimation : public BoundaryProcess, public io::Serializable {
     const std::shared_ptr<io::interval::Interval> radiationInterval;
 
     /**
-     * Emissivity of the fuel surface. For the gray assumption, this indicates how much of the black body radiation is absorbed and emitted from the fuel surface.
+     * Emissivity of the fuel surface. For the gray assumption, this indicates how much of the black body radiation is emitted from the fuel surface.
      * */
     const double emissivity;
 
+    /**
+     * Absorptivity of the fuel surface. For the gray assumption, this indicates how much of the incoming radiative surface heat flux is absorbed at the fuel surface.
+     * */
+    const double absorptivity;
+
     /**
      * Set the species densityYi based upon the blowing rate.  Update the energy if needed to maintain temperature
      */
@@ -99,7 +104,8 @@ class Sublimation : public BoundaryProcess, public io::Serializable {
     explicit Sublimation(std::shared_ptr<subModels::SublimationModel> sublimationModel, std::shared_ptr<ablate::eos::transport::TransportModel> transportModel, std::shared_ptr<ablate::eos::EOS> eos,
                          const std::shared_ptr<ablate::mathFunctions::FieldFunction> & = {}, std::shared_ptr<mathFunctions::MathFunction> additionalHeatFlux = {},
                          std::shared_ptr<finiteVolume::processes::PressureGradientScaling> pressureGradientScaling = {}, bool diffusionFlame = false,
-                         std::shared_ptr<ablate::radiation::SurfaceRadiation> radiationIn = {}, const std::shared_ptr<io::interval::Interval> &intervalIn = {}, double emissivityIn = 1);
+                         std::shared_ptr<ablate::radiation::SurfaceRadiation> radiationIn = {}, const std::shared_ptr<io::interval::Interval> &intervalIn = {},
+                         const std::shared_ptr<ablate::parameters::Parameters> & = {});
 
     void Setup(ablate::boundarySolver::BoundarySolver &bSolver) override;
     void Initialize(ablate::boundarySolver::BoundarySolver &bSolver) override;

src/domain/subDomain.cpp

@@ -226,7 +226,7 @@ DM ablate::domain::SubDomain::GetSubDM() {
     // If the subDM has not been created, create one
     if (!subDM) {
         // filter by label
-        DMPlexFilter(GetDM(), label, labelValue, &subDM) >> utilities::PetscUtilities::checkError;
+        DMPlexFilter(GetDM(), label, labelValue, PETSC_FALSE, PETSC_FALSE, NULL, &subDM) >> utilities::PetscUtilities::checkError;
 
         // copy over all fields that were in the main dm
         for (auto& fieldInfo : GetFields()) {
@@ -757,7 +757,7 @@ void ablate::domain::SubDomain::CreateEmptySubDM(DM* inDM, std::shared_ptr<domai
         subDmValue = labelValue;
     }
     if (subDmLabel) {
-        DMPlexFilter(GetDM(), subDmLabel, subDmValue, inDM);
+        DMPlexFilter(GetDM(), subDmLabel, subDmValue, PETSC_FALSE, PETSC_FALSE, NULL, inDM);
     } else {
         DMClone(GetDM(), inDM);
     }

src/monitors/boundarySolverMonitor.cpp

@@ -66,7 +66,7 @@ void ablate::monitors::BoundarySolverMonitor::Register(std::shared_ptr<solver::S
     DMPlexLabelComplete(boundaryDm, boundaryFaceLabel) >> utilities::PetscUtilities::checkError;
 
     // Now create a sub dm with only the faces
-    DMPlexFilter(boundaryDm, boundaryFaceLabel, 1, &faceDm) >> utilities::PetscUtilities::checkError;
+    DMPlexFilter(boundaryDm, boundaryFaceLabel, 1, PETSC_FALSE, PETSC_FALSE, NULL, &faceDm) >> utilities::PetscUtilities::checkError;
 
     // Add each of the output components on each face in the faceDm
     for (const auto& field : boundarySolver->GetOutputComponents()) {

src/monitors/radiationFlux.cpp

@@ -26,7 +26,7 @@ void ablate::monitors::RadiationFlux::Register(std::shared_ptr<solver::Solver> s
     domain::Region::GetLabel(radiationFluxRegion, solverIn->GetSubDomain().GetDM(), radiationFluxRegionLabel, regionValue);
 
     // Now create a sub dm with only the faces
-    DMPlexFilter(solverIn->GetSubDomain().GetDM(), radiationFluxRegionLabel, regionValue, &fluxDm) >> utilities::PetscUtilities::checkError;
+    DMPlexFilter(solverIn->GetSubDomain().GetDM(), radiationFluxRegionLabel, regionValue, PETSC_FALSE, PETSC_FALSE, NULL, &fluxDm) >> utilities::PetscUtilities::checkError;
 
     /** Add each of the output components on each face in the fluxDm
      * the number of components should be equal to the number of ray tracers plus any ratio outputs?
@@ -180,7 +180,7 @@ PetscErrorCode ablate::monitors::RadiationFlux::Save(PetscViewer viewer, PetscIn
                      * Get the intensity calculated out of the ray tracer. Write it to the appropriate location in the face DM.
                      */
                     PetscReal wavelengths[radiation[rayTracerIndex]->GetAbsorptionFunction().propertySize];
-                    radiation[rayTracerIndex]->GetSurfaceIntensity(wavelengths, boundaryPt, 0, 1);
+                    radiation[rayTracerIndex]->GetSurfaceIntensity(wavelengths, boundaryPt, 0, 1, 1);
                     if (log) log->Printf("%i:", c);
                     for (int wavelengthIndex = 0; wavelengthIndex < radiation[rayTracerIndex]->GetAbsorptionFunction().propertySize; ++wavelengthIndex) {
                         globalFaceData[radiation[rayTracerIndex]->GetAbsorptionFunction().propertySize * rayTracerIndex + wavelengthIndex] = wavelengths[wavelengthIndex];

src/radiation/orthogonalRadiation.hpp

@@ -21,7 +21,7 @@ class OrthogonalRadiation : public ablate::radiation::SurfaceRadiation {
      */
     static inline std::string GetClassType() { return "OrthogonalRadiation"; }
 
-    inline void GetSurfaceIntensity(PetscReal* intensityReturn, PetscInt faceId, PetscReal temperature, PetscReal emissivity = 1.0) override {
+    inline void GetSurfaceIntensity(PetscReal* intensityReturn, PetscInt faceId, PetscReal temperature, PetscReal emissivity = 1.0, PetscReal absorptivity = 1.0) override {
         for (int i = 0; i < (int)absorptivityFunction.propertySize; ++i) {  // Compute the losses
             PetscReal netIntensity = evaluatedGains[absorptivityFunction.propertySize * indexLookup.GetAbsoluteIndex(faceId) + i];
             intensityReturn[i] = abs(netIntensity) > ablate::utilities::Constants::large ? ablate::utilities::Constants::large * PetscSignReal(netIntensity) : netIntensity;

src/radiation/surfaceRadiation.hpp

@@ -42,16 +42,16 @@ class SurfaceRadiation : public ablate::radiation::Radiation {
      * @param emissivity the emissivity of the surface
      * @return
      */
-    virtual inline void GetSurfaceIntensity(PetscReal* intensity, PetscInt faceId, PetscReal temperature, PetscReal emissivity = 1.0) {
+    virtual inline void GetSurfaceIntensity(PetscReal* intensity, PetscInt faceId, PetscReal temperature, PetscReal emissivity = 1.0, PetscReal absorptivity = 1.0) {
         // add in precomputed gains
         for (int i = 0; i < (int)absorptivityFunction.propertySize; ++i) {  // Compute the losses
-            PetscReal netIntensity = -ablate::utilities::Constants::sbc * temperature * temperature * temperature * temperature;
+            // Emitted From Surface = \epsilon \sigma T^4
+            PetscReal netIntensity = -emissivity * ablate::utilities::Constants::sbc * temperature * temperature * temperature * temperature;
 
-            netIntensity += evaluatedGains[absorptivityFunction.propertySize * indexLookup.GetAbsoluteIndex(faceId) + i];
-
-            // scale by kappa
-            netIntensity *= emissivity;
+            // Absorbed at the surface = q''_Rad * \alpha
+            netIntensity += absorptivity * evaluatedGains[absorptivityFunction.propertySize * indexLookup.GetAbsoluteIndex(faceId) + i];
 
+            // This line just ensures that the netIntenisty isn't going to extreme values, If this clipping is really needed, it probably should be fixed somewhere else -klb
             intensity[i] = abs(netIntensity) > ablate::utilities::Constants::large ? ablate::utilities::Constants::large * PetscSignReal(netIntensity) : netIntensity;
         }
     }

src/solver/criteria/validRange.cpp

@@ -1,44 +1,13 @@
-#include "variableChange.hpp"
+#include "validRange.hpp"
+#include "convergenceException.hpp"
 
 #include <domain/subDomain.hpp>
 #include <utility>
 
-ablate::solver::criteria::VariableChange::VariableChange(std::string variableName, double convergenceTolerance, ablate::utilities::MathUtilities::Norm convergenceNorm,
-                                                         std::shared_ptr<ablate::domain::Region> region)
-    : variableName(std::move(variableName)), convergenceTolerance(convergenceTolerance), convergenceNorm(convergenceNorm), region(std::move(region)) {}
+ablate::solver::criteria::ValidRange::ValidRange(std::string variableName, double lowerBound, double upperBound, std::shared_ptr<ablate::domain::Region> region)
+    : variableName(std::move(variableName)), lowerBound(lowerBound), upperBound(upperBound), region(std::move(region)) {}
 
-ablate::solver::criteria::VariableChange::~VariableChange() {
-    if (previousValues) {
-        VecDestroy(&previousValues) >> ablate::utilities::PetscUtilities::checkError;
-    }
-}
-
-void ablate::solver::criteria::VariableChange::Initialize(const ablate::domain::Domain& domain) {
-    // Get the subdomain for this region
-    auto subDomain = domain.GetSubDomain(region);
-
-    // Look up the field
-    const auto& field = subDomain->GetField(variableName);
-
-    // Get the sub vector
-    IS subIs;
-    DM subDm;
-    Vec locVec;
-    subDomain->GetFieldLocalVector(field, 0.0, &subIs, &locVec, &subDm) >> ablate::utilities::PetscUtilities::checkError;
-
-    // Create a copy of the local vec
-    VecDuplicate(locVec, &previousValues) >> ablate::utilities::PetscUtilities::checkError;
-    // Set the block size as one so that it only produces one norm value
-    VecSetBlockSize(previousValues, 1) >> ablate::utilities::PetscUtilities::checkError;
-
-    // copy over the current values to use as the previous state
-    VecCopy(locVec, previousValues) >> ablate::utilities::PetscUtilities::checkError;
-
-    // restore
-    subDomain->RestoreFieldLocalVector(field, &subIs, &locVec, &subDm) >> ablate::utilities::PetscUtilities::checkError;
-}
-
-bool ablate::solver::criteria::VariableChange::CheckConvergence(const ablate::domain::Domain& domain, PetscReal time, PetscInt step, const std::shared_ptr<ablate::monitors::logs::Log>& log) {
+bool ablate::solver::criteria::ValidRange::CheckConvergence(const ablate::domain::Domain& domain, PetscReal time, PetscInt step, const std::shared_ptr<ablate::monitors::logs::Log>& log) {
     // Get the subdomain for this region
     auto subDomain = domain.GetSubDomain(region);
 
@@ -51,31 +20,28 @@ bool ablate::solver::criteria::VariableChange::CheckConvergence(const ablate::do
     Vec locVec;
     subDomain->GetFieldLocalVector(field, 0.0, &subIs, &locVec, &subDm) >> ablate::utilities::PetscUtilities::checkError;
 
-    // Compute the norm
-    PetscReal norm;
-    ablate::utilities::MathUtilities::ComputeNorm(convergenceNorm, locVec, previousValues, &norm) >> ablate::utilities::PetscUtilities::checkError;
-
-    // Create a copy of the local vec
-    VecCopy(locVec, previousValues) >> ablate::utilities::PetscUtilities::checkError;
+    // Get the min and max values
+    PetscScalar min, max;
+    VecMin(locVec, nullptr, &min) >> ablate::utilities::PetscUtilities::checkError;
+    VecMax(locVec, nullptr, &max) >> ablate::utilities::PetscUtilities::checkError;
 
     // restore
     subDomain->RestoreFieldLocalVector(field, &subIs, &locVec, &subDm) >> ablate::utilities::PetscUtilities::checkError;
 
-    // Check to see if converged
-    bool converged = norm < convergenceTolerance;
-    if (log) {
-        if (converged) {
-            log->Printf("\tVariableChange %s converged to %g\n", variableName.c_str(), norm);
-        } else {
-            log->Printf("\tVariableChange %s error: %g\n", variableName.c_str(), norm);
+    // check upper and lower bounds
+    if (max < lowerBound || min > upperBound) {
+        std::string message = max < lowerBound ? variableName + " all values fall below the the lower bound " + std::to_string(lowerBound) + ". Maximum value is " + std::to_string(max)
+                                               : variableName + " all values exceed the the upper bound " + std::to_string(upperBound) + ". Minimum value is " + std::to_string(min);
+        if (log) {
+            log->Printf("%s\n", message.c_str());
         }
+        throw ConvergenceException(message);
     }
 
-    return converged;
+    return true;
 }
 
 #include "registrar.hpp"
-REGISTER(ablate::solver::criteria::ConvergenceCriteria, ablate::solver::criteria::VariableChange, "This class checks for a relative change in the specified variable between checks",
-         ARG(std::string, "name", "the variable to check"), ARG(double, "tolerance", "the tolerance to reach to be considered converged"),
-         ENUM(ablate::utilities::MathUtilities::Norm, "norm", "norm type ('l1','l1_norm','l2', 'linf', 'l2_norm')"),
-         OPT(ablate::domain::Region, "region", "the region to check for a converged variable"));
+REGISTER(ablate::solver::criteria::ConvergenceCriteria, ablate::solver::criteria::ValidRange, "This class will stop the convergence iterations if the bounds are exceeded",
+         ARG(std::string, "name", "the variable to check"), ARG(double, "lowerBound", "values lower than this will result in an exception"),
+         ARG(double, "upperBound", "values higher than this will result in an exception"), OPT(ablate::domain::Region, "region", "the region to check for a converged variable"));