Skip to content

Commit

Permalink
Merge pull request #251 from NCAR/242-add-tutorial-for-rate-constant-…
Browse files Browse the repository at this point in the history
…types-except-custom-ones

242 add tutorial for rate constant types except custom ones
  • Loading branch information
K20shores committed Sep 26, 2023
2 parents 89ebafe + 2202a2a commit 6e7e991
Show file tree
Hide file tree
Showing 22 changed files with 764 additions and 79 deletions.
12 changes: 4 additions & 8 deletions docs/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -41,21 +41,17 @@ set(SPHINX_BUILD ${CMAKE_CURRENT_BINARY_DIR}/sphinx)
set(SPHINX_INDEX_FILE ${SPHINX_BUILD}/index.html)

add_custom_command(
OUTPUT ${SPHINX_INDEX_FILE}
OUTPUT ${SPHINX_INDEX_FILE} __fake_file_to_ensure_this_always_run
COMMAND
${SPHINX_EXECUTABLE} -b html
# Tell Breathe where to find the Doxygen output
-Dbreathe_projects.micm=${DOXYGEN_OUTPUT_DIR}/xml
${SPHINX_SOURCE} ${SPHINX_BUILD}
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
DEPENDS
# Other docs files you want to track should go here (or in some variable)
${CMAKE_CURRENT_SOURCE_DIR}/source/index.rst
${CMAKE_CURRENT_SOURCE_DIR}/source/getting_started.rst
${CMAKE_CURRENT_SOURCE_DIR}/source/user_guide.rst
${DOXYGEN_INDEX_FILE}
MAIN_DEPENDENCY ${SPHINX_SOURCE}/conf.py
COMMENT "Generating documentation with Sphinx"
# run this command each time. This way we don't have to keep track of each individual rst file
BuildDocs
)

add_custom_target(docs ALL DEPENDS ${SPHINX_INDEX_FILE})
add_custom_target(docs ALL DEPENDS ${SPHINX_INDEX_FILE} Doxygen)
4 changes: 3 additions & 1 deletion docs/requirements.txt
Original file line number Diff line number Diff line change
@@ -1,4 +1,6 @@
breathe
sphinx
sphinx-book-theme
sphinx-copybutton
sphinx-design
breathe
sphinx-tabs
30 changes: 30 additions & 0 deletions docs/source/_static/custom.css
Original file line number Diff line number Diff line change
Expand Up @@ -4,4 +4,34 @@
margin-left: auto;
margin-right: auto;
margin-top: 10px;
}

.download-div {
display: flex;
flex-direction: row;
justify-content: space-around;
margin-bottom: 1rem;
}

.download-button {
display: inline-block;
padding: 10px 20px;
background-color: #3498db;
color: white;
border: none;
border-radius: 5px;
text-align: center;
text-decoration: none;
font-size: 16px;
cursor: pointer;
transition: background-color 0.3s;
}

.download-button:hover {
background-color: #2980b9;
}

blockquote {
padding: unset;
border-left: unset;
}
Binary file not shown.
2 changes: 2 additions & 0 deletions docs/source/conf.py
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,9 @@
# ones.
extensions = [
'breathe',
'sphinx_copybutton',
'sphinx_design',
'sphinx_tabs.tabs',
]

breathe_default_project = "micm"
Expand Down
56 changes: 2 additions & 54 deletions docs/source/getting_started.rst
Original file line number Diff line number Diff line change
Expand Up @@ -78,60 +78,8 @@ See the [MICM documentation](https://ncar.github.io/micm/)
for details on the types of reactions available in MICM and how to configure them.
To solve this system save the following code in a file named `foo_chem.cpp`

.. code-block:: C++

#include <iomanip>
#include <iostream>
#include <micm/process/arrhenius_rate_constant.hpp>
#include <micm/solver/rosenbrock.hpp>

using namespace micm;

int main(const int argc, const char *argv[])
{
auto foo = Species{ "Foo" };
auto bar = Species{ "Bar" };
auto baz = Species{ "Baz" };
Phase gas_phase{ std::vector<Species>{ foo, bar, baz } };

System chemical_system{ SystemParameters{ .gas_phase_ = gas_phase } };

Process r1 = Process::create()
.reactants({ foo })
.products({ Yield(bar, 0.8), Yield(baz, 0.2) })
.rate_constant(ArrheniusRateConstant({ .A_ = 1.0e-3 }))
.phase(gas_phase);

Process r2 = Process::create()
.reactants({ foo, bar })
.products({ Yield(baz, 1) })
.rate_constant(ArrheniusRateConstant({ .A_ = 1.0e-5, .C_ = 110.0 }))
.phase(gas_phase);

std::vector<Process> reactions{ r1, r2 };

RosenbrockSolver solver{ chemical_system, reactions, RosenbrockSolverParameters::three_stage_rosenbrock_parameters() };

State state = solver.GetState();

state.conditions_[0].temperature_ = 287.45; // K
state.conditions_[0].pressure_ = 101319.9; // Pa
state.SetConcentration(foo, 20.0); // mol m-3

std::cout << "foo, bar, baz" << std::endl;
for (int i = 0; i < 10; ++i)
{
auto result = solver.Solve(500.0, state);
state.variables_ = result.result_;
std::cout << std::fixed << std::setprecision(6)
<< state.variables_[0][state.variable_map_["Foo"]] << ", "
<< state.variables_[0][state.variable_map_["Bar"]] << ", "
<< state.variables_[0][state.variable_map_["Baz"]] << std::endl;
}

return 0;
}
.. literalinclude:: ../../test/tutorial/test_README_example.cpp
:language: cpp

To build and run the example using GNU::

Expand Down
5 changes: 3 additions & 2 deletions docs/source/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
.. ~ for subsubsubsections
.. " for paragraphs
================================
Welcome to MICM's documentation!
================================

Expand All @@ -22,7 +23,7 @@ Welcome to MICM's documentation!

.. grid-item-card:: User guide
:img-top: _static/index_user_guide.svg
:link: user_guide
:link: user_guide/index
:link-type: doc

Learn how to configure micm for your mechanisms here!
Expand All @@ -48,7 +49,7 @@ Welcome to MICM's documentation!
:caption: Contents:

getting_started
user_guide
user_guide/index
api/index
contributing/index
citation
Expand Down
5 changes: 0 additions & 5 deletions docs/source/user_guide.rst

This file was deleted.

22 changes: 22 additions & 0 deletions docs/source/user_guide/index.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,22 @@
##########
User Guide
##########

MICM is quite expansive in its capabilities. We've written examples showing many different use cases.
Hopefully our examples are exhaustive enough to provide you with the tools you need to solve some atmospheric chemistry.

If you happen to find our examples are lacking for your needs, please,
`fill out an issue <https://github.com/NCAR/micm/issues/new>`_ and request the kind of example you'd like.



1. :ref:`Rate constants`
2. :ref:`User defined rate constants`


.. toctree::
:maxdepth: 1
:caption: Contents:

rate_constant_tutorial
user_defined_rate_constant_tutorial
178 changes: 178 additions & 0 deletions docs/source/user_guide/rate_constant_tutorial.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,178 @@
.. _Rate constants:

Rate Constants (except user defined ones)
=========================================

MICM supports a subset of the rate constants defined as part of the
`OpenAtmos Mechanism Configuration <https://open-atmos.github.io/MechanismConfiguration/reactions/index.html>`_
We will be adding more in the future.
The links to the ``micm`` classes below detail the class and methods. Please check the OpenAtmos standard for
specific on the algorithm implemented for each rate constant type.
At present, supported rate constants are:

- :cpp:class:`micm::ArrheniusRateConstant`, `OpenAtmos Arrhenius description <https://open-atmos.github.io/MechanismConfiguration/reactions/arrhenius.html>`_
- :cpp:class:`micm::BranchedRateConstant`, `OpenAtmos Branched description <https://open-atmos.github.io/MechanismConfiguration/reactions/branched.html>`_
- :cpp:class:`micm::SurfaceRateConstant`, `OpenAtmos Surface description <https://open-atmos.github.io/MechanismConfiguration/reactions/surface.html>`_
- :cpp:class:`micm::TernaryChemicalActivationRateConstant`, `OpenAtmos Ternay chemical activiation description <https://open-atmos.github.io/MechanismConfiguration/reactions/ternary_chemical_activation.html>`_
- :cpp:class:`micm::TroeRateConstant`, `OpenAtmos Troe description <https://open-atmos.github.io/MechanismConfiguration/reactions/troe.html>`_
- :cpp:class:`micm::TunnelingRateConstant`, `OpenAtmos Tunneling description <https://open-atmos.github.io/MechanismConfiguration/reactions/tunneling.html>`_
- :cpp:class:`micm::UserDefinedRateConstant`, this is a custom type, but this is how photolysis rates are setup `OpenAtmos Photolysis description <https://open-atmos.github.io/MechanismConfiguration/reactions/photolysis.html>`_

This tutorial covers all but the last one. See the :ref:`User defined rate constants` tutorial for examples and use
cases on that.

We'll setup and solve a fake chemical system with 7 species and 6 reactions,

.. math::
A &\longrightarrow B, &k_{1, \mathrm{arrhenius}} \\
B &\longrightarrow C (\mathrm{alkoxy\ products}) + D (\mathrm{nitrate\ products}), &k_{2, \mathrm{branched}} \\
C &\longrightarrow E, &k_{3, \mathrm{surface}} \\
D &\longrightarrow 2F, &k_{4, \mathrm{ternary\ chemical\ activation}} \\
2E &\longrightarrow G, &k_{5, \mathrm{troe}} \\
F &\longrightarrow G, &k_{6, \mathrm{tunneling}} \\
MICM can be configured in two ways. We can either build the mechanism up by hand with the ``micm`` API,
or parse a valid mechanism Configuration
in the OpenAtmos format. In this tutorial, we will do both. If you're looking for a copy and paste, choose
the appropriate tab below and be on your way! Otherwise, stick around for a line by line explanation.


.. tabs::

.. tab:: Build the Mechanism with the API

.. literalinclude:: ../../../test/tutorial/test_rate_constants_no_user_defined_by_hand.cpp
:language: cpp

.. tab:: OpenAtmos Configuration reading

.. raw:: html

<div class="download-div">
<a href="../_static/tutorials/rate_constants_no_user_defined.zip" download>
<button class="download-button">Download zip configuration</button>
</a>
</div>

.. literalinclude:: ../../../test/tutorial/test_rate_constants_no_user_defined_with_config.cpp
:language: cpp

Line-by-line explanation
------------------------

To get started, we'll need to include each rate constant type and the
rosenbrock solver at the top of the file.

.. literalinclude:: ../../../test/tutorial/test_rate_constants_no_user_defined_by_hand.cpp
:language: cpp
:lines: 1-13

After that, we'll use the ``micm`` namespace and setup a template alias so that we can instantiate the
rosenbrock solver.

.. literalinclude:: ../../../test/tutorial/test_rate_constants_no_user_defined_by_hand.cpp
:language: cpp
:lines: 15-22

To create a :cpp:class:`micm::RosenbrockSolver`, we have to define a chemical system (:cpp:class:`micm::System`)
and our reactions, which will be a vector of :cpp:class:`micm::Process` We will use the species to define these.

.. tabs::

.. tab:: Build the Mechanism with the API

To do this by hand, we have to define all of the chemical species in the system. This allows us to set
any properties of the species that may be necessary for rate constanta calculations, like molecular weights
and diffusion coefficients for the surface reaction. We will also put these species into the gas phase.

.. literalinclude:: ../../../test/tutorial/test_rate_constants_no_user_defined_by_hand.cpp
:language: cpp
:lines: 56-67

Now that we have a gas phase and our species, we can start building the reactions. Two things to note are that
stoichiemtric coefficients for reactants are represented by repeating that product as many times as you need.
To specify the yield of a product, we've created a typedef :cpp:type:`micm::Yield`
and a function :cpp:func:`micm::yields` that produces these.

.. literalinclude:: ../../../test/tutorial/test_rate_constants_no_user_defined_by_hand.cpp
:language: cpp
:lines: 69-133

And finally we define our chemical system and reactions

.. literalinclude:: ../../../test/tutorial/test_rate_constants_no_user_defined_by_hand.cpp
:language: cpp
:lines: 135-136

.. tab:: OpenAtmos Configuration reading

After defining a valid OpenAtmos configuration with reactions that ``micm`` supports, configuring the chemical
system and the processes is as simple as using the :cpp:class:`micm::SolverConfig` class

.. literalinclude:: ../../../test/tutorial/test_rate_constants_no_user_defined_with_config.cpp
:language: cpp
:lines: 57-69

Now that we have a chemical system and a list of reactions, we can create the RosenbrockSolver.
There are several ways to configure the solver. Here we are using a three stage solver. More options
can be found in the :cpp:class:`micm::RosenbrockSolverParameters`

.. literalinclude:: ../../../test/tutorial/test_rate_constants_no_user_defined_by_hand.cpp
:language: cpp
:lines: 138-140

The rosenbrock solver will provide us a state, which we can use to set the concentrations,
custom rate parameters, and temperature and pressure

.. tabs::

.. tab:: Build the Mechanism with the API

.. literalinclude:: ../../../test/tutorial/test_rate_constants_no_user_defined_by_hand.cpp
:language: cpp
:lines: 141-155

.. tab:: OpenAtmos Configuration reading

.. literalinclude:: ../../../test/tutorial/test_rate_constants_no_user_defined_with_config.cpp
:language: cpp
:lines: 75-93


Finally, we are ready to pick a timestep ans solve the system.

.. literalinclude:: ../../../test/tutorial/test_rate_constants_no_user_defined_by_hand.cpp
:language: cpp
:lines: 157-183


This is the output:


+-------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
| time | A | B | C | D | E | F | G |
+=======+=============+=============+=============+=============+=============+=============+=============+
| 0 | 1.00e+00 | 0.00e+00 | 0.00e+00 | 0.00e+00 | 0.00e+00 | 0.00e+00 | 0.00e+00 |
+-------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
| 500 | 3.18e-09 | 3.66e-09 | 9.83e-01 | 3.88e-14 | 1.41e-03 | 2.02e-13 | 7.92e-03 |
+-------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
| 1000 | 1.14e-14 | 1.31e-14 | 9.66e-01 | 1.39e-19 | 1.40e-03 | 7.24e-19 | 1.64e-02 |
+-------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
| 1500 | 4.09e-20 | 4.71e-20 | 9.49e-01 | 4.98e-25 | 1.39e-03 | 2.59e-24 | 2.48e-02 |
+-------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
| 2000 | 1.47e-25 | 1.69e-25 | 9.33e-01 | 1.79e-30 | 1.38e-03 | 9.30e-30 | 3.30e-02 |
+-------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
| 2500 | 5.26e-31 | 6.05e-31 | 9.17e-01 | 6.40e-36 | 1.37e-03 | 3.33e-35 | 4.11e-02 |
+-------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
| 3000 | 1.89e-36 | 2.17e-36 | 9.01e-01 | 2.30e-41 | 1.36e-03 | 1.20e-40 | 4.90e-02 |
+-------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
| 3500 | 6.77e-42 | 7.78e-42 | 8.85e-01 | 8.23e-47 | 1.34e-03 | 4.29e-46 | 5.68e-02 |
+-------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
| 4000 | 2.43e-47 | 2.79e-47 | 8.70e-01 | 2.95e-52 | 1.33e-03 | 1.54e-51 | 6.44e-02 |
+-------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
| 4500 | 8.70e-53 | 1.00e-52 | 8.55e-01 | 1.06e-57 | 1.32e-03 | 5.51e-57 | 7.20e-02 |
+-------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
| 5000 | 3.12e-58 | 3.59e-58 | 8.40e-01 | 3.80e-63 | 1.31e-03 | 1.98e-62 | 7.94e-02 |
+-------+-------------+-------------+-------------+-------------+-------------+-------------+-------------+
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
.. _User defined rate constants:

User Defined Rate Constants
###########################
1 change: 1 addition & 0 deletions include/micm/process/process.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
namespace micm
{

/// @brief An alias that allows making products easily
using Yield = std::pair<micm::Species, double>;

inline Yield yields(micm::Species species, double yield)
Expand Down
Loading

0 comments on commit 6e7e991

Please sign in to comment.