Adding ability to update multiple Fields at once with cell interpolants and face interpolants,

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.13.01)
+project(ablateLibrary VERSION 0.13.02)
 
 # Load the Required 3rd Party Libaries
 pkg_check_modules(PETSc REQUIRED IMPORTED_TARGET GLOBAL PETSc)

CMakePresets.json

@@ -34,36 +34,6 @@
         "PKG_CONFIG_PATH": "$env{PETSC_DIR}/arch-ablate-opt/lib/pkgconfig:$penv{PKG_CONFIG_PATH}"
       }
     },
-    {
-      "name": "OptMainPetsc",
-      "displayName": "CLion Opt Config",
-      "description": "Default build for ABLATE in CLion",
-      "binaryDir": "${sourceDir}/cmake-build-MainPetsc",
-      "cacheVariables": {
-        "CMAKE_BUILD_TYPE": "Release",
-        "COMPILE_MPI_COMMAND": "$env{PETSC_DIR}/arch-ablate-opt/bin/mpirun"
-      },
-      "environment": {
-        "PETSC_DIR": "/home/klbud/Software/PetscUpdating/petsc",
-        "PETSC_ARCH": "arch-ablate-opt",
-        "PKG_CONFIG_PATH": "$env{PETSC_DIR}/arch-ablate-opt/lib/pkgconfig:$penv{PKG_CONFIG_PATH}"
-      }
-    },
-    {
-      "name": "OptFluentPetsc",
-      "displayName": "CLion Opt Config",
-      "description": "Default build for ABLATE in CLion",
-      "binaryDir": "${sourceDir}/cmake-build-FluentPetsc",
-      "cacheVariables": {
-        "CMAKE_BUILD_TYPE": "Release",
-        "COMPILE_MPI_COMMAND": "$env{PETSC_DIR}/arch-ablate-opt/bin/mpirun"
-      },
-      "environment": {
-        "PETSC_DIR": "/home/klbud/Software/PetscUpdating/petscW",
-        "PETSC_ARCH": "arch-ablate-opt",
-        "PKG_CONFIG_PATH": "$env{PETSC_DIR}/arch-ablate-opt/lib/pkgconfig:$penv{PKG_CONFIG_PATH}"
-      }
-    },
     {
       "name": "local-ablate-opt-info",
       "displayName": "CLion RelWithDebugInfo Config",

src/finiteVolume/cellInterpolant.cpp

@@ -99,7 +99,6 @@ void ablate::finiteVolume::CellInterpolant::ComputeRHS(PetscReal time, Vec locXV
             VecGetArrayRead(locGradVecs[field.subId], &locGradArrays[field.subId]) >> utilities::PetscUtilities::checkError;
         }
     }
-
     ComputeFluxSourceTerms(dm,
                            ds,
                            totDim,
@@ -499,8 +498,8 @@ void ablate::finiteVolume::CellInterpolant::ComputeFluxSourceTerms(DM dm, PetscD
     PetscScalar *auxL = nullptr, *auxR = nullptr;
 
     // Precompute the offsets to pass into the rhsFluxFunctionDescriptions
-    std::vector<PetscInt> fluxComponentSize(rhsFunctions.size());
-    std::vector<PetscInt> fluxId(rhsFunctions.size());
+    std::vector<std::vector<PetscInt>> fluxComponentSize(rhsFunctions.size());
+    std::vector<std::vector<PetscInt>> fluxId(rhsFunctions.size());
     std::vector<std::vector<PetscInt>> uOff(rhsFunctions.size());
     std::vector<std::vector<PetscInt>> aOff(rhsFunctions.size());
 
@@ -509,9 +508,11 @@ void ablate::finiteVolume::CellInterpolant::ComputeFluxSourceTerms(DM dm, PetscD
     PetscDSGetComponentOffsets(ds, &uOffTotal) >> utilities::PetscUtilities::checkError;
 
     for (std::size_t fun = 0; fun < rhsFunctions.size(); fun++) {
-        const auto& field = subDomain->GetField(rhsFunctions[fun].field);
-        fluxComponentSize[fun] = field.numberComponents;
-        fluxId[fun] = field.id;
+        for (std::size_t f = 0; f < rhsFunctions[fun].updateFields.size(); f++) {
+            const auto& field = subDomain->GetField(rhsFunctions[fun].updateFields[f]);
+            fluxComponentSize[fun].push_back(field.numberComponents);
+            fluxId[fun].push_back(field.id);
+        }
         for (std::size_t f = 0; f < rhsFunctions[fun].inputFields.size(); f++) {
             uOff[fun].push_back(uOffTotal[rhsFunctions[fun].inputFields[f]]);
         }
@@ -534,7 +535,6 @@ 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;
@@ -574,33 +574,36 @@ void ablate::finiteVolume::CellInterpolant::ComputeFluxSourceTerms(DM dm, PetscD
 
         // 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;
+                }
 
-            // add the flux back to the cell
-            PetscScalar *fL = nullptr, *fR = nullptr;
-            PetscInt cellLabelValue = regionValue;
-            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], 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], locFArray, &fR) >> 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]; ++d) {
-                if (fL) fL[d] -= flux[d] / cgL->volume;
-                if (fR) fR[d] += flux[d] / cgR->volume;
+                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];
             }
         }
     }

src/finiteVolume/cellInterpolant.hpp

@@ -26,7 +26,7 @@ class CellInterpolant {
         DiscontinuousFluxFunction function;
         void* context;
 
-        PetscInt field;
+        std::vector<PetscInt> updateFields;
         std::vector<PetscInt> inputFields;
         std::vector<PetscInt> auxFields;
     };

src/finiteVolume/compressibleFlowSolver.cpp

@@ -1,37 +1,55 @@
 #include "compressibleFlowSolver.hpp"
 #include <utility>
 #include "compressibleFlowFields.hpp"
+#include "finiteVolume/processes/compactCompressibleNSSpeciesSingleProgressTransport.hpp"
+#include "finiteVolume/processes/compactCompressibleNSSpeciesTransport.hpp"
 #include "finiteVolume/processes/evTransport.hpp"
 #include "finiteVolume/processes/navierStokesTransport.hpp"
 #include "finiteVolume/processes/speciesTransport.hpp"
 #include "utilities/vectorUtilities.hpp"
 
+/**
+ * The compact argument in this constructor is used to define whether we use seperate processes for the Euler (NS) transport, species transport, and EV transport.
+ * a value of 0, is all seperate
+ * a value of 1, is combined Euler and species transport, with seperate EV transport
+ * a value of 2, is combined Euler and species and The single field Progress EV transport (soot number density in this case)
+ */
 ablate::finiteVolume::CompressibleFlowSolver::CompressibleFlowSolver(std::string solverId, std::shared_ptr<domain::Region> region, std::shared_ptr<parameters::Parameters> options,
                                                                      const std::shared_ptr<eos::EOS>& eosIn, const std::shared_ptr<parameters::Parameters>& parameters,
                                                                      const std::shared_ptr<eos::transport::TransportModel>& transport,
                                                                      const std::shared_ptr<fluxCalculator::FluxCalculator>& fluxCalculatorIn,
                                                                      std::vector<std::shared_ptr<processes::Process>> additionalProcesses,
                                                                      std::vector<std::shared_ptr<boundaryConditions::BoundaryCondition>> boundaryConditions,
-                                                                     const std::shared_ptr<eos::transport::TransportModel>& evTransport)
-    : FiniteVolumeSolver(std::move(solverId), std::move(region), std::move(options),
-                         utilities::VectorUtilities::Merge(
-                             {
-                                 // create assumed processes for compressible flow
-                                 std::make_shared<ablate::finiteVolume::processes::NavierStokesTransport>(
-                                     parameters, eosIn, fluxCalculatorIn, transport, utilities::VectorUtilities::Find<ablate::finiteVolume::processes::PressureGradientScaling>(additionalProcesses)),
-                                 std::make_shared<ablate::finiteVolume::processes::SpeciesTransport>(eosIn, fluxCalculatorIn, transport, parameters),
-                                 std::make_shared<ablate::finiteVolume::processes::EVTransport>(eosIn, fluxCalculatorIn, evTransport ? evTransport : transport),
-                             },
-                             additionalProcesses),
-                         std::move(boundaryConditions)) {}
+                                                                     const std::shared_ptr<eos::transport::TransportModel>& evTransport, int compact)
+    : FiniteVolumeSolver(
+          std::move(solverId), std::move(region), std::move(options),
+          compact ? (compact == 1 ? utilities::VectorUtilities::Merge({std::make_shared<ablate::finiteVolume::processes::CompactCompressibleNSSpeciesTransport>(
+                                                                           parameters, eosIn, fluxCalculatorIn, transport,
+                                                                           utilities::VectorUtilities::Find<ablate::finiteVolume::processes::PressureGradientScaling>(additionalProcesses)),
+                                                                       std::make_shared<ablate::finiteVolume::processes::EVTransport>(eosIn, fluxCalculatorIn, evTransport ? evTransport : transport)},
+                                                                      additionalProcesses)
+                                  : utilities::VectorUtilities::Merge({std::make_shared<ablate::finiteVolume::processes::CompactCompressibleNSSpeciesSingleProgressTransport>(
+                                                                          parameters, eosIn, fluxCalculatorIn, transport, evTransport ? evTransport : transport,
+                                                                          utilities::VectorUtilities::Find<ablate::finiteVolume::processes::PressureGradientScaling>(additionalProcesses))},
+                                                                      additionalProcesses))
+                  : utilities::VectorUtilities::Merge(
+                        {
+                            // create assumed processes for compressible flow
+                            std::make_shared<ablate::finiteVolume::processes::NavierStokesTransport>(
+                                parameters, eosIn, fluxCalculatorIn, transport, utilities::VectorUtilities::Find<ablate::finiteVolume::processes::PressureGradientScaling>(additionalProcesses)),
+                            std::make_shared<ablate::finiteVolume::processes::SpeciesTransport>(eosIn, fluxCalculatorIn, transport, parameters),
+                            std::make_shared<ablate::finiteVolume::processes::EVTransport>(eosIn, fluxCalculatorIn, evTransport ? evTransport : transport),
+                        },
+                        additionalProcesses),
+          std::move(boundaryConditions)) {}
 
 ablate::finiteVolume::CompressibleFlowSolver::CompressibleFlowSolver(std::string solverId, std::shared_ptr<domain::Region> region, std::shared_ptr<parameters::Parameters> options,
                                                                      const std::shared_ptr<eos::EOS>& eosIn, const std::shared_ptr<parameters::Parameters>& parameters,
                                                                      const std::shared_ptr<eos::transport::TransportModel>& transport,
                                                                      const std::shared_ptr<fluxCalculator::FluxCalculator>& fluxCalculatorIn,
                                                                      std::vector<std::shared_ptr<boundaryConditions::BoundaryCondition>> boundaryConditions,
-                                                                     const std::shared_ptr<eos::transport::TransportModel>& evTransport)
-    : CompressibleFlowSolver(std::move(solverId), std::move(region), std::move(options), eosIn, parameters, transport, fluxCalculatorIn, {}, std::move(boundaryConditions), evTransport) {}
+                                                                     const std::shared_ptr<eos::transport::TransportModel>& evTransport, int compact)
+    : CompressibleFlowSolver(std::move(solverId), std::move(region), std::move(options), eosIn, parameters, transport, fluxCalculatorIn, {}, std::move(boundaryConditions), evTransport, compact) {}
 
 #include "registrar.hpp"
 REGISTER(ablate::solver::Solver, ablate::finiteVolume::CompressibleFlowSolver, "compressible finite volume flow", ARG(std::string, "id", "the name of the flow field"),
@@ -41,4 +59,5 @@ REGISTER(ablate::solver::Solver, ablate::finiteVolume::CompressibleFlowSolver, "
          OPT(ablate::finiteVolume::fluxCalculator::FluxCalculator, "fluxCalculator", "the flux calculators (defaults to none)"),
          OPT(std::vector<ablate::finiteVolume::processes::Process>, "additionalProcesses", "any additional processes besides euler/yi/ev transport"),
          OPT(std::vector<ablate::finiteVolume::boundaryConditions::BoundaryCondition>, "boundaryConditions", "the boundary conditions for the flow field"),
-         OPT(ablate::eos::transport::TransportModel, "evTransport", "when provided, this model will be used for ev transport instead of default"));
+         OPT(ablate::eos::transport::TransportModel, "evTransport", "when provided, this model will be used for ev transport instead of default"),
+         OPT(int, "compact", "Integer value describing whether to treat all the transport seperately, partially combined, or fully combined (see commented code above constructor for values)"));

src/finiteVolume/compressibleFlowSolver.hpp

@@ -31,7 +31,8 @@ class CompressibleFlowSolver : public FiniteVolumeSolver {
     CompressibleFlowSolver(std::string solverId, std::shared_ptr<domain::Region> region, std::shared_ptr<parameters::Parameters> options, const std::shared_ptr<eos::EOS>& eos,
                            const std::shared_ptr<parameters::Parameters>& parameters, const std::shared_ptr<eos::transport::TransportModel>& transport,
                            const std::shared_ptr<fluxCalculator::FluxCalculator>& = {}, std::vector<std::shared_ptr<processes::Process>> additionalProcesses = {},
-                           std::vector<std::shared_ptr<boundaryConditions::BoundaryCondition>> boundaryConditions = {}, const std::shared_ptr<eos::transport::TransportModel>& evTransport = {});
+                           std::vector<std::shared_ptr<boundaryConditions::BoundaryCondition>> boundaryConditions = {}, const std::shared_ptr<eos::transport::TransportModel>& evTransport = {},
+                           int compact = 0);
 
     /**
      * Constructor without ev or additional processes
@@ -48,7 +49,7 @@ class CompressibleFlowSolver : public FiniteVolumeSolver {
     CompressibleFlowSolver(std::string solverId, std::shared_ptr<domain::Region> region, std::shared_ptr<parameters::Parameters> options, const std::shared_ptr<eos::EOS>& eos,
                            const std::shared_ptr<parameters::Parameters>& parameters, const std::shared_ptr<eos::transport::TransportModel>& transport,
                            const std::shared_ptr<fluxCalculator::FluxCalculator>& = {}, std::vector<std::shared_ptr<boundaryConditions::BoundaryCondition>> boundaryConditions = {},
-                           const std::shared_ptr<eos::transport::TransportModel>& evTransport = {});
+                           const std::shared_ptr<eos::transport::TransportModel>& evTransport = {}, int compact = 0);
     ~CompressibleFlowSolver() override = default;
 };
 }  // namespace ablate::finiteVolume