Skip to content

Commit

Permalink
162 euler option redo (#204)
Browse files Browse the repository at this point in the history 
* Added Backward Euler definitions to musics/mich.hpp.

* Added Backward Euler definitions micm.cpp.

* Added unit test BackwardEulerStandardOrder.

* Added unit test BackwardEulerVectorOrder.

* Save.
  • Loading branch information
@dwfncar
dwfncar authored Aug 23, 2024
1 parent d7d70e3 commit 8353ff7
Show file tree
Hide file tree
Showing 4 changed files with 175 additions and 6 deletions.
2 changes: 1 addition & 1 deletion fortran/micm.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -452,4 +452,4 @@ subroutine finalize(this)
ASSERT(error%is_success())
end subroutine finalize

end module musica_micm
end module musica_micm
32 changes: 28 additions & 4 deletions include/musica/micm.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,8 +39,10 @@ namespace musica
/// @brief Types of MICM solver
enum MICMSolver
{
Rosenbrock = 1, // Vector-ordered Rosenbrock solver
RosenbrockStandardOrder, // Standard-ordered Rosenbrock solver
Rosenbrock = 1, // Vector-ordered Rosenbrock solver
RosenbrockStandardOrder, // Standard-ordered Rosenbrock solver
BackwardEuler, // Vector-ordered BackwardEuler solver
BackwardEulerStandardOrder, // Standard-ordered BackwardEuler solver
};

struct SolverResultStats
Expand Down Expand Up @@ -181,6 +183,16 @@ namespace musica
/// @param error Error struct to indicate success or failure
void CreateRosenbrockStandardOrder(const std::string &config_path, Error *error);

/// @brief Create a BackwardEuler solver of vector-ordered matrix type by reading and parsing configuration file
/// @param config_path Path to configuration file or directory containing configuration file
/// @param error Error struct to indicate success or failure
void CreateBackwardEuler(const std::string &config_path, Error *error);

/// @brief Create a BackwardEuler solver of standard-ordered matrix type by reading and parsing configuration file
/// @param config_path Path to configuration file or directory containing configuration file
/// @param error Error struct to indicate success or failure
void CreateBackwardEulerStandardOrder(const std::string &config_path, Error *error);

/// @brief Solve the system
/// @param solver Pointer to solver
/// @param time_step Time [s] to advance the state by
Expand Down Expand Up @@ -242,7 +254,7 @@ namespace musica
public:
MICMSolver solver_type_;

/// @brief Vector-ordered Rosenbrock solver type
/// @brief Vector-ordered Rosenbrock
using DenseMatrixVector = micm::VectorMatrix<double, MICM_VECTOR_MATRIX_SIZE>;
using SparseMatrixVector = micm::SparseMatrix<double, micm::SparseMatrixVectorOrdering<MICM_VECTOR_MATRIX_SIZE>>;
using RosenbrockVectorType = typename micm::RosenbrockSolverParameters::
Expand All @@ -260,6 +272,18 @@ namespace musica
using StandardState = micm::State<DenseMatrixStandard, SparseMatrixStandard>;
std::unique_ptr<RosenbrockStandard> rosenbrock_standard_;

/// @brief Vector-ordered Backward Euler
using BackwardEulerVectorType = typename micm::BackwardEulerSolverParameters::
template SolverType<micm::ProcessSet, micm::LinearSolver<SparseMatrixVector, micm::LuDecomposition>>;
using BackwardEuler = micm::Solver<BackwardEulerVectorType, micm::State<DenseMatrixVector, SparseMatrixVector>>;
std::unique_ptr<BackwardEuler> backward_euler_;

/// @brief Standard-ordered Backward Euler
using BackwardEulerStandardType = typename micm::BackwardEulerSolverParameters::
template SolverType<micm::ProcessSet, micm::LinearSolver<SparseMatrixStandard, micm::LuDecomposition>>;
using BackwardEulerStandard = micm::Solver<BackwardEulerStandardType, micm::State<DenseMatrixStandard, SparseMatrixStandard>>;
std::unique_ptr<BackwardEulerStandard> backward_euler_standard_;

/// @brief Returns the number of grid cells
/// @return Number of grid cells
int NumGridCells() const
Expand Down Expand Up @@ -333,4 +357,4 @@ namespace musica
*error = ToError(MUSICA_ERROR_CATEGORY, MUSICA_ERROR_CODE_SPECIES_NOT_FOUND, msg.c_str());
return T();
}
} // namespace musica
} // namespace musica
115 changes: 115 additions & 0 deletions src/micm/micm.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,16 @@ namespace musica
micm->SetSolverType(MICMSolver::RosenbrockStandardOrder);
micm->CreateRosenbrockStandardOrder(std::string(config_path), error);
}
else if (solver_type == MICMSolver::BackwardEuler)
{
micm->SetSolverType(MICMSolver::BackwardEuler);
micm->CreateBackwardEuler(std::string(config_path), error);
}
else if (solver_type == MICMSolver::BackwardEulerStandardOrder)
{
micm->SetSolverType(MICMSolver::BackwardEulerStandardOrder);
micm->CreateBackwardEulerStandardOrder(std::string(config_path), error);
}
else
{
std::string msg = "Solver type '" + std::to_string(solver_type) + "' not found";
Expand Down Expand Up @@ -113,6 +123,34 @@ namespace musica
solver_stats,
error);
}
else if (micm->solver_type_ == MICMSolver::BackwardEuler)
{
micm->Solve(
micm->backward_euler_,
time_step,
temperature,
pressure,
air_density,
concentrations,
custom_rate_parameters,
solver_state,
solver_stats,
error);
}
else if (micm->solver_type_ == MICMSolver::BackwardEulerStandardOrder)
{
micm->Solve(
micm->backward_euler_standard_,
time_step,
temperature,
pressure,
air_density,
concentrations,
custom_rate_parameters,
solver_state,
solver_stats,
error);
}
};

String MicmVersion()
Expand All @@ -134,6 +172,14 @@ namespace musica
{
map = micm->GetSpeciesOrdering(micm->rosenbrock_standard_, error);
}
else if (micm->solver_type_ == MICMSolver::BackwardEuler)
{
map = micm->GetSpeciesOrdering(micm->backward_euler_, error);
}
else if (micm->solver_type_ == MICMSolver::BackwardEulerStandardOrder)
{
map = micm->GetSpeciesOrdering(micm->backward_euler_standard_, error);
}
if (!IsSuccess(*error))
{
return nullptr;
Expand Down Expand Up @@ -168,6 +214,14 @@ namespace musica
{
map = micm->GetUserDefinedReactionRatesOrdering(micm->rosenbrock_standard_, error);
}
else if (micm->solver_type_ == MICMSolver::BackwardEuler)
{
map = micm->GetUserDefinedReactionRatesOrdering(micm->backward_euler_, error);
}
else if (micm->solver_type_ == MICMSolver::BackwardEulerStandardOrder)
{
map = micm->GetUserDefinedReactionRatesOrdering(micm->backward_euler_standard_, error);
}
if (!IsSuccess(*error))
{
return nullptr;
Expand Down Expand Up @@ -282,6 +336,67 @@ namespace musica
}
}

void MICM::CreateBackwardEuler(const std::string &config_path, Error *error)
{
try
{
micm::SolverConfig<> solver_config;
solver_config.ReadAndParse(std::filesystem::path(config_path));
solver_parameters_ = std::make_unique<micm::SolverParameters>(solver_config.GetSolverParams());

backward_euler_ = std::make_unique<BackwardEuler>(
micm::SolverBuilder<
micm::BackwardEulerSolverParameters,
micm::VectorMatrix<double, MICM_VECTOR_MATRIX_SIZE>,
micm::SparseMatrix<double, micm::SparseMatrixVectorOrdering<MICM_VECTOR_MATRIX_SIZE>>,
micm::ProcessSet,
micm::LinearSolver<
micm::SparseMatrix<double, micm::SparseMatrixVectorOrdering<MICM_VECTOR_MATRIX_SIZE>>,
micm::LuDecomposition>,
VectorState>(micm::BackwardEulerSolverParameters())
.SetSystem(solver_parameters_->system_)
.SetReactions(solver_parameters_->processes_)
.SetNumberOfGridCells(num_grid_cells_)
.SetIgnoreUnusedSpecies(true)
.Build());

DeleteError(error);
*error = NoError();
}
catch (const std::system_error &e)
{
DeleteError(error);
*error = ToError(e);
}
}

void MICM::CreateBackwardEulerStandardOrder(const std::string &config_path, Error *error)
{
try
{
micm::SolverConfig<> solver_config;
solver_config.ReadAndParse(std::filesystem::path(config_path));
solver_parameters_ = std::make_unique<micm::SolverParameters>(solver_config.GetSolverParams());

backward_euler_standard_ =
std::make_unique<BackwardEulerStandard>(micm::CpuSolverBuilder<micm::BackwardEulerSolverParameters>(
micm::BackwardEulerSolverParameters())
.SetSystem(solver_parameters_->system_)
.SetReactions(solver_parameters_->processes_)
.SetNumberOfGridCells(num_grid_cells_)
.SetIgnoreUnusedSpecies(true)
.Build());

DeleteError(error);
*error = NoError();
}
catch (const std::system_error &e)
{
DeleteError(error);
*error = ToError(e);
}
}

void MICM::Solve(
auto &solver,
double time_step,
Expand Down
32 changes: 31 additions & 1 deletion src/test/unit/micm/micm_c_api.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -295,6 +295,36 @@ TEST(RosenbrockStandardOrder, SolveUsingStandardOrderedRosenbrock)
DeleteError(&error);
}

// Test case for solving system using vector-ordered Backward Euler solver
TEST(BackwardEulerVectorOrder, SolveUsingVectordOrderedBackwardEuler)
{
const char* config_path = "configs/chapman";
int num_grid_cells = 1;
Error error;
MICM* micm = CreateMicm(config_path, MICMSolver::BackwardEuler, num_grid_cells, &error);

TestSingleGridCell(micm);

DeleteMicm(micm, &error);
ASSERT_TRUE(IsSuccess(error));
DeleteError(&error);
}

// Test case for solving system using standard-ordered Backward Euler solver
TEST(BackwardEulerStandardOrder, SolveUsingStandardOrderedBackwardEuler)
{
const char* config_path = "configs/chapman";
int num_grid_cells = 1;
Error error;
MICM* micm = CreateMicm(config_path, MICMSolver::BackwardEulerStandardOrder, num_grid_cells, &error);

TestSingleGridCell(micm);

DeleteMicm(micm, &error);
ASSERT_TRUE(IsSuccess(error));
DeleteError(&error);
}

struct ArrheniusReaction
{
double A_{ 1 };
Expand Down Expand Up @@ -479,4 +509,4 @@ TEST_F(MicmCApiTest, GetSpeciesProperty)
GetSpeciesPropertyDouble(micm, "O3", "bad property", &error);
ASSERT_TRUE(IsError(error, MICM_ERROR_CATEGORY_SPECIES, MICM_SPECIES_ERROR_CODE_PROPERTY_NOT_FOUND));
DeleteError(&error);
}
}

0 comments on commit 8353ff7

Please sign in to comment.