Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

162 euler option redo #204

Merged
merged 5 commits into from
Aug 23, 2024
Merged

162 euler option redo #204

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
9fb15b6
Added Backward Euler definitions to musics/mich.hpp.
dwfncar Aug 17, 2024
33e55a4
Added Backward Euler definitions micm.cpp.
dwfncar Aug 17, 2024
dea6f4d
Added unit test BackwardEulerStandardOrder.
dwfncar Aug 17, 2024
5f08aec
Added unit test BackwardEulerVectorOrder.
dwfncar Aug 17, 2024
c6fcf7a
Save.
dwfncar Aug 18, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion fortran/micm.F90
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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);
}
}
Loading