TransportSolver.hpp

// Reaktoro is a unified framework for modeling chemically reactive systems.
//
// Copyright (C) 2014-2018 Allan Leal
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with this library. If not, see <http://www.gnu.org/licenses/>.
 
#pragma once
 
// C++ includes
#include <memory>
#include <vector>
 
// Reaktoro includes
#include <Reaktoro/Common/Index.hpp>
#include <Reaktoro/Common/StringList.hpp>
#include <Reaktoro/Core/ChemicalOutput.hpp>
#include <Reaktoro/Core/ChemicalProperties.hpp>
#include <Reaktoro/Core/ChemicalState.hpp>
#include <Reaktoro/Core/ChemicalSystem.hpp>
#include <Reaktoro/Equilibrium/EquilibriumSolver.hpp>
#include <Reaktoro/Math/Matrix.hpp>
 
namespace Reaktoro {
 
// Forward declarations
//class ChemicalProperties;
//class ChemicalState;
//class ChemicalSystem;
 
//class BoundaryState
//{
//public:
//    BoundaryState(const ChemicalState& state);
//
//private:
//};
//
//class BoundaryFlux
//{
//public:
//    BoundaryFlux(const ChemicalState& state);
//
//private:
//};
 
 
class ChemicalField
{
public:
    using Iterator = std::vector<ChemicalState>::iterator;
 
    using ConstIterator = std::vector<ChemicalState>::const_iterator;
 
    ChemicalField(Index size, const ChemicalSystem& system);
 
    ChemicalField(Index size, const ChemicalState& state);
 
    auto size() const -> Index { return m_size; }
 
    auto begin() const -> ConstIterator { return m_states.cbegin(); }
 
    auto begin() -> Iterator { return m_states.begin(); }
 
    auto end() const -> ConstIterator { return m_states.cend(); }
 
    auto end() -> Iterator { return m_states.end(); }
 
    auto operator[](Index index) const -> const ChemicalState& { return m_states[index]; }
 
    auto operator[](Index index) -> ChemicalState& { return m_states[index]; }
 
    auto set(const ChemicalState& state) -> void;
 
    auto temperature(VectorRef values) -> void;
 
    auto pressure(VectorRef values) -> void;
 
    auto elementAmounts(VectorRef values) -> void;
 
    auto output(std::string filename, StringList quantities) -> void;
 
private:
    Index m_size;
 
//    Vector temperatures;
//
//    Vector pressures;
//
//    /// The matrix of amounts for every element (
//    Matrix element_amounts;
 
    ChemicalSystem m_system;
 
    std::vector<ChemicalState> m_states;
 
    std::vector<ChemicalProperties> m_properties;
};
 
class TridiagonalMatrix
{
public:
    TridiagonalMatrix() : TridiagonalMatrix(0) {}
 
    TridiagonalMatrix(Index size) : m_size(size), m_data(size * 3) {}
 
    auto size() const -> Index { return m_size; }
 
    auto data() -> VectorRef { return m_data; }
 
    auto data() const -> VectorConstRef { return m_data; }
 
    auto row(Index index) -> VectorRef { return m_data.segment(3 * index, 3); }
 
    auto row(Index index) const -> VectorConstRef { return m_data.segment(3 * index, 3); }
 
    auto a() -> VectorStridedRef { return Vector::Map(m_data.data() + 3, size() - 1, Eigen::InnerStride<3>()); }
 
    auto a() const -> VectorConstRef { return Vector::Map(m_data.data() + 3, size() - 1, Eigen::InnerStride<3>()); }
 
    auto b() -> VectorStridedRef { return Vector::Map(m_data.data() + 1, size(), Eigen::InnerStride<3>()); }
 
    auto b() const -> VectorConstRef { return Vector::Map(m_data.data() + 1, size(), Eigen::InnerStride<3>()); }
 
    auto c() -> VectorStridedRef { return Vector::Map(m_data.data() + 2, size() - 1, Eigen::InnerStride<3>()); }
 
    auto c() const -> VectorConstRef { return Vector::Map(m_data.data() + 2, size() - 1, Eigen::InnerStride<3>()); }
 
    auto resize(Index size) -> void;
 
    auto factorize() -> void;
 
    auto solve(VectorRef x, VectorConstRef b) const -> void;
 
    auto solve(VectorRef x) const -> void;
 
    operator Matrix() const;
 
private:
    Index m_size;
 
    Vector m_data;
};
 
class Mesh
{
public:
    Mesh();
 
    Mesh(Index num_cells, double xl = 0.0, double xr = 1.0);
 
    auto setDiscretization(Index num_cells, double xl = 0.0, double xr = 1.0) -> void;
 
    auto numCells() const -> Index { return m_num_cells; }
 
    auto xl() const -> double { return m_xl; }
 
    auto xr() const -> double { return m_xr; }
 
    auto dx() const -> double { return m_dx; }
 
    auto xcells() const -> VectorConstRef { return m_xcells; }
 
private:
    Index m_num_cells = 10;
 
    double m_xl = 0.0;
 
    double m_xr = 1.0;
 
    double m_dx = 0.1;
 
    Vector m_xcells;
};
 
class TransportSolver
{
public:
    TransportSolver();
 
    auto setMesh(const Mesh& mesh) -> void { mesh_ = mesh; }
 
    auto setVelocity(double val) -> void { velocity = val; }
 
    auto setDiffusionCoeff(double val) -> void { diffusion = val; }
 
    auto setBoundaryValue(double val) -> void { ul = val; };
 
    auto setTimeStep(double val) -> void { dt = val; }
 
    auto mesh() const -> const Mesh& { return mesh_; }
 
    auto initialize() -> void;
 
    auto step(VectorRef u, VectorConstRef q) -> void;
 
    auto step(VectorRef u) -> void;
 
private:
    Mesh mesh_;
 
    double dt = 0.0;
 
    double velocity = 0.0;
 
    double diffusion = 0.0;
 
    double ul;
 
    TridiagonalMatrix A;
 
    Vector phi;
 
    Vector u0;
};
 
class ReactiveTransportSolver
{
public:
    ReactiveTransportSolver(const ChemicalSystem& system);
 
    auto setMesh(const Mesh& mesh) -> void;
 
    auto setVelocity(double val) -> void;
 
    auto setDiffusionCoeff(double val) -> void;
 
    auto setBoundaryState(const ChemicalState& state) -> void;
 
    auto setTimeStep(double val) -> void;
 
    auto system() const -> const ChemicalSystem& { return system_; }
 
    auto output() -> ChemicalOutput;
 
    auto initialize(const ChemicalField& field) -> void;
 
    auto step(ChemicalField& field) -> void;
 
private:
    ChemicalSystem system_;
 
    TransportSolver transportsolver;
 
    EquilibriumSolver equilibriumsolver;
 
    std::vector<ChemicalOutput> outputs;
 
    Vector bbc;
 
    Matrix bf;
 
    Matrix bs;
 
    Matrix b;
 
    Index steps = 0;
};
 
} // namespace Reaktoro