radiation.hpp

#ifndef ABLATELIBRARY_RADIATION_HPP
#define ABLATELIBRARY_RADIATION_HPP
 
#include <petsc/private/dmimpl.h>
#include <petscdm.h>
#include <petscdmswarm.h>
#include <petscsf.h>
#include <memory>
#include <set>
#include <utility>
#include "eos/radiationProperties/radiationProperties.hpp"
#include "finiteVolume/compressibleFlowFields.hpp"
#include "finiteVolume/finiteVolumeSolver.hpp"
#include "io/interval/interval.hpp"
#include "monitors/logs/log.hpp"
#include "solver/cellSolver.hpp"
#include "solver/timeStepper.hpp"
#include "utilities/constants.hpp"
#include "utilities/loggable.hpp"
#include "utilities/mpiUtilities.hpp"
#include "utilities/petscUtilities.hpp"
 
namespace ablate::radiation {
 
class Radiation : protected utilities::Loggable<Radiation> {  
 
   public:
    Radiation(const std::string& solverId, const std::shared_ptr<domain::Region>& region, const PetscInt raynumber, std::shared_ptr<eos::radiationProperties::RadiationModel> radiationModelIn,
              std::shared_ptr<ablate::monitors::logs::Log> = {});
 
    virtual ~Radiation();
 
    struct Identifier {
        PetscInt originRank;
        PetscInt originRayId;
        PetscInt remoteRank;
        PetscInt remoteRayId;
        PetscInt nSegment;
    };
 
    struct Carrier {
        PetscReal Ij = 0;    
        PetscReal Krad = 1;  
    };
 
    static inline PetscReal GetBlackBodyTotalIntensity(PetscReal temperature, const PetscReal refractiveIndex) {
        return refractiveIndex * refractiveIndex * ablate::utilities::Constants::sbc * temperature * temperature * temperature * temperature / ablate::utilities::Constants::pi;
    }
 
    static inline PetscReal GetBlackBodyWavelengthIntensity(const PetscReal temperature, const PetscReal wavelength, const PetscReal refractiveIndex) {
        PetscReal E_lambda = (2 * ablate::utilities::Constants::pi * ablate::utilities::Constants::h * ablate::utilities::Constants::c * ablate::utilities::Constants::c);
        E_lambda /= refractiveIndex * refractiveIndex * wavelength * wavelength * wavelength * wavelength * wavelength;
        E_lambda /= exp((ablate::utilities::Constants::h * ablate::utilities::Constants::c) / (refractiveIndex * wavelength * ablate::utilities::Constants::k * temperature)) - 1;
        E_lambda /= ablate::utilities::Constants::pi;
 
        return E_lambda;
    }
 
    static inline std::string GetClassType() { return "Radiation"; }
 
    virtual void Setup(const ablate::domain::Range& cellRange, ablate::domain::SubDomain& subDomain);
 
    virtual void Initialize(const ablate::domain::Range& cellRange, ablate::domain::SubDomain& subDomain);
 
    inline void GetIntensity(PetscReal* intensity, PetscInt index, const domain::Range& cellRange, PetscReal temperature, PetscReal kappa) {
        for (int i = 0; i < (int)absorptivityFunction.propertySize; ++i) {                                                  // Compute the losses
            PetscReal netIntensity = -4.0 * ablate::utilities::Constants::pi * GetBlackBodyTotalIntensity(temperature, 1);  
 
            // add in precomputed gains
            netIntensity += evaluatedGains[absorptivityFunction.propertySize * (index - cellRange.start) + i];
 
            // scale by kappa
            netIntensity *= kappa;
 
            intensity[i] = abs(netIntensity) > ablate::utilities::Constants::large ? ablate::utilities::Constants::large * PetscSignReal(netIntensity) : netIntensity;
        }
    }
 
    //  Add a class called "radiation surface properties" which is used to compute the amount of intensity absorbed by the material. This should be an optional input in the radiation base class.
    // GetIntensity should compute this, because otherwise the losses will not be accounted for.
    // Just embed a function within the GetIntensity function that "Computes surface properties" and decides how much of the radiation to do.
    // Also. the losses are wavelength dependent based on what portion of the black body temperature they are emitting.
    // Then, the camera properties will be able to be viewed as a radiation surface property.
    // The radiation solver can take a vector of radiation surface properties.
    // Each radiation surface property will be evaluated on the EvaluateGains call.
    // Each property can be output separately.
 
    inline std::string GetId() { return solverId; };
 
    inline eos::ThermodynamicTemperatureFunction GetAbsorptionFunction() { return absorptivityFunction; };
 
    void EvaluateGains(Vec solVec, ablate::domain::Field temperatureField, Vec auxVec);
 
    virtual void ParticleStep(ablate::domain::SubDomain& subDomain, DM faceDM, const PetscScalar* faceGeomArray, DM radReturn, PetscInt nlocalpoints,
                              PetscInt nglobalpoints);  
 
    virtual void IdentifyNewRaysOnRank(domain::SubDomain& subDomain, DM radReturn, PetscInt nlocalpoints);
 
    virtual PetscReal SurfaceComponent(const PetscReal normal[], PetscInt iCell, PetscInt nphi, PetscInt ntheta);
 
    inline std::shared_ptr<eos::radiationProperties::RadiationModel> GetRadiationModel() { return radiationModel; }
 
   protected:
    DM radSearch = nullptr;
 
    Vec faceGeomVec = nullptr;
    Vec cellGeomVec = nullptr;
 
    MPI_Datatype carrierMpiType;
 
    struct CellSegment {
        PetscInt cell;
        PetscReal pathLength;
    };
 
    struct Virtualcoord {
        PetscReal x;
        PetscReal y;
        PetscReal z;
        PetscReal xdir;
        PetscReal ydir;
        PetscReal zdir;
        PetscReal hhere;
    };
 
 
    PetscReal FaceIntersect(PetscInt ip, Virtualcoord* virtualcoord, PetscFVFaceGeom* face) const;  
 
    void UpdateCoordinates(PetscInt ipart, Virtualcoord* virtualcoord, PetscReal* coord, PetscReal adv) const;
 
    void DeleteOutOfBounds(ablate::domain::SubDomain& subDomain);
 
    virtual void SetBoundary(CellSegment& raySegment, PetscInt index, Identifier identifier) {
        raySegment.cell = index;
        raySegment.pathLength = -1;
    }
 
    PetscInt dim = 0;  
    PetscInt nTheta;   
    PetscInt nPhi;     
    PetscReal minCellRadius{};
 
    std::vector<std::vector<CellSegment>> raySegments;
 
    std::vector<Carrier> raySegmentsCalculations;
 
    PetscInt numberOriginRays;
 
    PetscInt numberOriginCells;
 
    PetscInt raysPerCell;
 
    std::vector<unsigned short int> raySegmentsPerOriginRay;
 
    std::vector<Carrier> raySegmentSummary;
 
    std::vector<PetscReal> gainsFactor;
 
    std::vector<PetscScalar> evaluatedGains;
 
    PetscSF remoteAccess = nullptr;
 
    std::string solverId;
 
    const std::shared_ptr<domain::Region> region;
 
    const std::shared_ptr<eos::radiationProperties::RadiationModel> radiationModel;
 
    eos::ThermodynamicTemperatureFunction absorptivityFunction;
    eos::ThermodynamicTemperatureFunction emissivityFunction;
 
    // !Store a log used to output the required information
    const std::shared_ptr<ablate::monitors::logs::Log> log = nullptr;
    static inline constexpr char IdentifierField[] = "identifier";
    static inline constexpr char VirtualCoordField[] = "virtual coord";
 
   public:
    inline std::shared_ptr<domain::Region> GetRegion() const { return region; }
};
std::ostream& operator<<(std::ostream& os, const Radiation::Identifier& id);
 
}  // namespace ablate::radiation
#endif