Provide an interface to the micm State

[K20shores]

In #225, we started storing a single state. This is because large systems take a long time (~5 seconds) to initialize the state. Repeated calls to solve from python were constantly making and deleting a state which very drastically degraded solving time.

Additionally, fortran and python have to call solve with many arguments

musica/fortran/test/fetch_content_integration/test_micm_api.F90

Lines 122 to 123 in 954265b

	call micm%solve(time_step, temperature, pressure, air_density, concentrations, &
	user_defined_reaction_rates, solver_state, solver_stats, error)

musica/python/test/test_analytical.py

Lines 55 to 62 in 954265b

	musica.micm_solve(
	solver,
	time_step,
	temperature,
	pressure,
	air_density,
	ordered_concentrations,
	ordered_rate_constants)

And, when that data gets into C++, we have to copy it onto the state

musica/src/micm/micm.cpp

Lines 437 to 450 in 326b511

	std::size_t i_cond = 0;
	for (auto &cond : state.conditions_)
	{
	cond.temperature_ = temperature[i_cond];
	cond.pressure_ = pressure[i_cond];
	cond.air_density_ = air_density[i_cond++];
	}
	std::size_t i_species_elem = 0;
	for (auto &var : state.variables_.AsVector())
	var = concentrations[i_species_elem++];
	std::size_t i_custom_rate_elem = 0;
	for (auto &var : state.custom_rate_parameters_.AsVector())
	var = custom_rate_parameters[i_custom_rate_elem++];

musica/python/wrapper.cpp

Lines 48 to 122 in 326b511

	std::vector<double> temperature_cpp;
	if (py::isinstance<py::float_>(temperature))
	{
	temperature_cpp.push_back(temperature.cast<double>());
	}
	else if (py::isinstance<py::list>(temperature))
	{
	py::list temperature_list = temperature.cast<py::list>();
	temperature_cpp.reserve(len(temperature_list));
	for (auto item : temperature_list)
	{
	temperature_cpp.push_back(item.cast<double>());
	}
	}
	else
	{
	throw std::runtime_error(
	"Temperature must be a list or a double. Got " +
	std::string(py::str(temperature.get_type()).cast<std::string>()));
	}
	std::vector<double> pressure_cpp;
	if (py::isinstance<py::float_>(pressure))
	{
	pressure_cpp.push_back(pressure.cast<double>());
	}
	else if (py::isinstance<py::list>(pressure))
	{
	py::list pressure_list = pressure.cast<py::list>();
	pressure_cpp.reserve(len(pressure_list));
	for (auto item : pressure_list)
	{
	pressure_cpp.push_back(item.cast<double>());
	}
	}
	else
	{
	throw std::runtime_error(
	"Pressure must be a list or a double. Got " + std::string(py::str(pressure.get_type()).cast<std::string>()));
	}
	std::vector<double> air_density_cpp;
	if (py::isinstance<py::float_>(air_density))
	{
	air_density_cpp.push_back(air_density.cast<double>());
	}
	else if (py::isinstance<py::list>(air_density))
	{
	py::list air_density_list = air_density.cast<py::list>();
	air_density_cpp.reserve(len(air_density_list));
	for (auto item : air_density_list)
	{
	air_density_cpp.push_back(item.cast<double>());
	}
	}
	else
	{
	throw std::runtime_error(
	"Air density must be a list or a double. Got " +
	std::string(py::str(air_density.get_type()).cast<std::string>()));
	}
	std::vector<double> concentrations_cpp;
	concentrations_cpp.reserve(len(concentrations));
	for (auto item : concentrations)
	{
	concentrations_cpp.push_back(item.cast<double>());
	}
	std::vector<double> custom_rate_parameters_cpp;
	if (!custom_rate_parameters.is_none())
	{
	py::list parameters = custom_rate_parameters.cast<py::list>();
	custom_rate_parameters_cpp.reserve(len(parameters));
	for (auto item : parameters)
	{
	custom_rate_parameters_cpp.push_back(item.cast<double>());
	}
	}

This means that we are doing many copies of the data before we even get to solve.

Write a wrapper around the micm state class which allows for data to be accessed and set directly in both python and fortran.

Acceptance criteria

• A new State wrapper class is created and accessible everywhere
• A state is not stored alongside a solver

  • musica/include/musica/micm.hpp

    Lines 254 to 267 in 326b511

    	std::pair<std::unique_ptr<RosenbrockStandard>, StandardState> 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::pair<std::unique_ptr<BackwardEuler>, VectorState> 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::pair<std::unique_ptr<BackwardEulerStandard>, StandardState> backward_euler_standard_;

• All solving interfaces are updated to pass a state
  • You may be able to make a solver state pair and just pass that as we do now, but be able to remove the other arguments. Note that time_step is not a state parameter

  • musica/src/micm/micm.cpp

    Lines 418 to 428 in 326b511

    	void MICM::Solve(
    	auto &solver_state_pair,
    	double time_step,
    	double *temperature,
    	double *pressure,
    	double *air_density,
    	double *concentrations,
    	double *custom_rate_parameters,
    	String *solver_state,
    	SolverResultStats *solver_stats,
    	Error *error)

Ideas

• We need to figure out how this issue in micm gets resolved for us to really figure this out
  • That did get resolved and it was a bug in micm.
• Add a new State class in https://github.com/NCAR/musica/blob/main/include/musica/micm.hpp or in a new file
  • support operations to
    • get and set temperature
    • get and set pressure
    • get and set air density
    • get and set concentrations
    • set user defined reaction rates
• The implementation of all of the above can go in https://github.com/NCAR/musica/blob/main/src/micm/micm.cpp, or you can make a new file
• Add fortran interfaces to these new functions in https://github.com/NCAR/musica/blob/main/fortran/micm.F90
• Add unit tests for the state in https://github.com/NCAR/musica/tree/main/src/test/unit/micm
  • and update the existing tests so that this new state is used
• Update fortran tests
  • https://github.com/NCAR/musica/blob/main/fortran/test/fetch_content_integration/test_micm_api.F90
  • https://github.com/NCAR/musica/blob/main/fortran/test/fetch_content_integration/test_micm_box_model.F90