GeochemicalSolver.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#pragma once

#include "GeochemicalSystem.h"
#include "GeochemistryIonicStrength.h"

class GeochemicalSolver
{
public:
  GeochemicalSolver(unsigned num_basis,
                    unsigned num_kin,
                    GeochemistryIonicStrength & is,
                    Real abs_tol,
                    Real rel_tol,
                    unsigned max_iter,
                    Real max_initial_residual,
                    Real swap_threshold,
                    unsigned max_swaps_allowed,
                    const std::vector<std::string> & prevent_precipitation,
                    Real max_ionic_strength,
                    unsigned ramp_max_ionic_strength,
                    bool evaluate_kin_always);

  void solveSystem(GeochemicalSystem & egs,
                   std::stringstream & ss,
                   unsigned & tot_iter,
                   Real & abs_residual,
                   Real dt,
                   DenseVector<Real> & mole_additions,
                   DenseMatrix<Real> & dmole_additions);

  void setMaxInitialResidual(Real max_initial_residual);

  Real getMaxInitialResidual() const;

  void setRampMaxIonicStrength(unsigned ramp_max_ionic_strength);

  unsigned getRampMaxIonicStrength() const;

private:
  GeochemistryIonicStrength & _is;
  const unsigned _num_basis;
  const unsigned _num_kin;
  unsigned _num_basis_in_algebraic_system;
  unsigned _num_in_algebraic_system;
  DenseVector<Real> _residual;
  Real _abs_residual;
  DenseMatrix<Real> _jacobian;
  DenseVector<Real> _new_mol;
  const Real _abs_tol;
  const Real _rel_tol;
  Real _res0_times_rel;
  const unsigned _max_iter;
  Real _max_initial_residual;
  const Real _swap_threshold;
  const unsigned _max_swaps_allowed;
  const std::vector<std::string> _prevent_precipitation;
  const Real _max_ionic_strength;
  unsigned _ramp_max_ionic_strength;
  bool _evaluate_kin_always;
  DenseVector<Real> _input_mole_additions;
  DenseMatrix<Real> _input_dmole_additions;

  Real computeResidual(const GeochemicalSystem & egs,
                       DenseVector<Real> & residual,
                       const DenseVector<Real> & mole_additions) const;

  void solveAndUnderrelax(const GeochemicalSystem & egs,
                          DenseMatrix<Real> & jacobian,
                          DenseVector<Real> & new_mol) const;

  bool swapNeeded(const GeochemicalSystem & egs,
                  unsigned & swap_out_of_basis,
                  unsigned & swap_into_basis,
                  std::stringstream & ss) const;

  bool reduceInitialResidual(GeochemicalSystem & egs,
                             Real dt,
                             DenseVector<Real> & mole_additions,
                             DenseMatrix<Real> & dmole_additions);
};