Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #pragma once
11 :
12 : #include "DenseMatrix.h"
13 : #include "PertinentGeochemicalSystem.h"
14 : #include <libmesh/dense_vector.h>
15 :
16 : /**
17 : * Class to swap basis species with equilibrium species
18 : */
19 : class GeochemistrySpeciesSwapper
20 : {
21 : public:
22 : /**
23 : * @param basis_size Number of basis species in the geochemical model
24 : * @param stoi_tol Swapping involves inverting matrices via a singular value decomposition. During
25 : * this process: (1) if abs(singular value) < stoi_tol * L1norm(singular values), then the matrix
26 : * is deemed singular (so the basis swap is deemed invalid); (2) if abs(any stoichiometric
27 : * coefficient) < stoi_tol then it is set to zero. A suggested reasonable value is stoi_tol=1E-6
28 : */
29 : GeochemistrySpeciesSwapper(unsigned basis_size, Real stoi_tol);
30 :
31 : bool operator==(const GeochemistrySpeciesSwapper & rhs) const
32 : {
33 2 : return (_stoi_tol == rhs._stoi_tol) && (_swap_matrix.m() == rhs._swap_matrix.m());
34 : };
35 :
36 : /**
37 : * Check that replacing the named basis species with the named equilibrium species is valid. In
38 : * performing this check, the "swap" matrix is inverted, so the check is somewhat expensive.
39 : * Therefore, if the swap will be performed, assuming it's valid, the performSwap method should be
40 : * used instead.
41 : * A mooseError or mooseException will result if the swap is invalid
42 : * @param mgd The ModelGeochemicalDatabase that holds all information regarding basis
43 : * components, equilibrium species, stoichiometries, etc
44 : * @param replace_this The basis component that will be removed from the basis and added to the
45 : * equilibrium species list
46 : * @param with_this The equilibrium species that will be removed from the equilibrium species list
47 : * and added to the basis component list
48 : */
49 : void checkSwap(const ModelGeochemicalDatabase & mgd,
50 : const std::string & replace_this,
51 : const std::string & with_this);
52 :
53 : /**
54 : * Check that replacing the indexed basis species with the indexed equilibrium species is valid.
55 : * In performing this check, the "swap" matrix is inverted, so the check is somewhat expensive.
56 : * Therefore, if the swap will be performed, assuming it's valid, the performSwap method should be
57 : * used instead.
58 : * A mooseError or mooseException will result if the swap is invalid
59 : * @param mgd The ModelGeochemicalDatabase that holds all information regarding basis
60 : * components, equilibrium species, stoichiometries, etc
61 : * @param basis_index_to_replace The index of the basis component that will be removed from the
62 : * basis and added to the equilibrium species list
63 : * @param eqm_index_to_insert The index of the equilibrium species that will be removed from the
64 : * equilibrium species list and added to the basis component list
65 : */
66 : void checkSwap(const ModelGeochemicalDatabase & mgd,
67 : unsigned basis_index_to_replace,
68 : unsigned eqm_index_to_insert);
69 :
70 : /**
71 : * Check that replacing the named basis species with the named equilibrium species is valid, and
72 : * then perform this swap by altering the mgd data structure.
73 : * A mooseError or mooseException will result if the swap is invalid
74 : * @param mgd The ModelGeochemicalDatabase that holds all information regarding basis
75 : * components, equilibrium species, stoichiometries, etc
76 : * @param replace_this The basis component that will be removed from the basis and added to the
77 : * equilibrium species list
78 : * @param with_this The equilibrium species that will be removed from the equilibrium species list
79 : * and added to the basis component list
80 : */
81 : void performSwap(ModelGeochemicalDatabase & mgd,
82 : const std::string & replace_this,
83 : const std::string & with_this);
84 :
85 : /**
86 : * Check that replacing the indexed basis species with the indexed equilibrium species is valid,
87 : * and then perform this swap by altering the mgd data structure.
88 : * A mooseError or mooseException will result if the swap is invalid
89 : * @param mgd The ModelGeochemicalDatabase that holds all information regarding basis
90 : * components, equilibrium species, stoichiometries, etc
91 : * @param basis_index_to_replace The index of the basis component that will be removed from the
92 : * basis and added to the equilibrium species list
93 : * @param eqm_index_to_insert The index of the equilibrium species that will be removed from the
94 : * equilibrium species list and added to the basis component list
95 : */
96 : void performSwap(ModelGeochemicalDatabase & mgd,
97 : unsigned basis_index_to_replace,
98 : unsigned eqm_index_to_insert);
99 :
100 : /**
101 : * Check that replacing the named basis species with the named equilibrium species is valid, and
102 : * then perform this swap by altering the mgd data structure, and calculate the bulk composition
103 : * of the new system.
104 : * A mooseError or mooseException will result if the swap is invalid
105 : * @param mgd The ModelGeochemicalDatabase that holds all information regarding basis
106 : * components, equilibrium species, stoichiometries, etc
107 : * @param[in/out] bulk_composition Upon entry, the bulk composition in the original basis. Upon
108 : * exit, the bulk composition in the basis after the swap
109 : * @param replace_this The basis component that will be removed from the basis and added to the
110 : * equilibrium species list
111 : * @param with_this The equilibrium species that will be removed from the equilibrium species list
112 : * and added to the basis component list
113 : */
114 : void performSwap(ModelGeochemicalDatabase & mgd,
115 : DenseVector<Real> & bulk_composition,
116 : const std::string & replace_this,
117 : const std::string & with_this);
118 :
119 : /**
120 : * Check that replacing the indexed basis species with the indexed equilibrium species is valid,
121 : * and then perform this swap by altering the mgd data structure, and calculate the bulk
122 : * composition of the new system.
123 : * A mooseError or mooseException will result if the swap is invalid
124 : * @param mgd The ModelGeochemicalDatabase that holds all information regarding basis
125 : * components, equilibrium species, stoichiometries, etc
126 : * @param[in/out] bulk_composition Upon entry, the bulk composition in the original basis. Upon
127 : * exit, the bulk composition in the basis after the swap
128 : * @param basis_index_to_replace The index of the basis component that will be removed from the
129 : * basis and added to the equilibrium species list
130 : * @param eqm_index_to_insert The index of the equilibrium species that will be removed from the
131 : * equilibrium species list and added to the basis component list
132 : */
133 : void performSwap(ModelGeochemicalDatabase & mgd,
134 : DenseVector<Real> & bulk_composition,
135 : unsigned basis_index_to_replace,
136 : unsigned eqm_index_to_insert);
137 :
138 : /**
139 : * For the the given basis index, find the equilibrium index that it should be swapped with
140 : * @param basis_ind index of basis species in mgd that we'd like to swap out of the basis
141 : * @param mgd the model's geochemical databgase
142 : * @param eqm_molality the molality of all the equilibrium species
143 : * @param minerals_allowed if false, best_eqm_species will not correspond to a mineral. Since
144 : * equilibrium minerals have molality=0 it is very unlikely that best_eqm_species will correspond
145 : * to a mineral, even if this flag is true.
146 : * @param gases_allowed if false, best_eqm_species will not correspond to a gas. Since
147 : * equilibrium gases have molality=0 it is very unlikely that best_eqm_species will correspond to
148 : * a gas, even if this flag is true.
149 : * @param sorption_allowed if false, best_eqm_species will not be involved in sorption
150 : * @param[out] best_eqm_species the index of the equilibrium species that should be swapped with
151 : * basis_ind
152 : * @return true if a valid swap was found. if false, then best_eqm_species will be garbage
153 : */
154 : bool findBestEqmSwap(unsigned basis_ind,
155 : const ModelGeochemicalDatabase & mgd,
156 : const std::vector<Real> & eqm_molality,
157 : bool minerals_allowed,
158 : bool gas_allowed,
159 : bool sorption_allowed,
160 : unsigned & best_eqm_species) const;
161 :
162 : private:
163 : /**
164 : * Construct the swap matrix, and its inverse, that describes the swap between the indexed basis
165 : * species and the indexed equilibrium index. The inverse of the swap matrix is used in alterMGD
166 : * to modify its datastructures to implement the swap.
167 : * A mooseError or mooseException will result if the swap is invalid
168 : * @param mgd The ModelGeochemicalDatabase that holds all information regarding basis
169 : * components, equilibrium species, stoichiometries, etc
170 : * @param basis_index_to_replace The index of the basis component that will be removed from the
171 : * basis and added to the equilibrium species list
172 : * @param eqm_index_to_insert The index of the equilibrium species that will be removed from the
173 : * equilibrium species list and added to the basis component list
174 : */
175 : void constructInverseMatrix(const ModelGeochemicalDatabase & mgd,
176 : unsigned basis_index_to_replace,
177 : unsigned eqm_index_to_insert);
178 :
179 : /**
180 : * Modify the ModelGeochemicalDatabase mgd to swap the indexed basis
181 : * species and the indexed equilibrium index.
182 : * @param mgd The ModelGeochemicalDatabase that holds all information regarding basis
183 : * components, equilibrium species, stoichiometries, etc
184 : * @param basis_index_to_replace The index of the basis component that will be removed from the
185 : * basis and added to the equilibrium species list
186 : * @param eqm_index_to_insert The index of the equilibrium species that will be removed from the
187 : * equilibrium species list and added to the basis component list
188 : */
189 : void alterMGD(ModelGeochemicalDatabase & mgd,
190 : unsigned basis_index_to_replace,
191 : unsigned eqm_index_to_insert);
192 :
193 : /**
194 : * Calculate the bulk composition in the new basis based on the bulk composition in the original
195 : * basis
196 : * @param basis_size Size of basis, which is only used for error checking
197 : * @param[in/out] bulk_composition Upon entry, the bulk composition in the original basis. Upon
198 : * exit, the bulk composition in the basis after the swap
199 : */
200 : void alterBulkComposition(unsigned basis_size, DenseVector<Real> & bulk_composition) const;
201 :
202 : /// tolerance for matrix inversion and stoichiometric coefficients
203 : Real _stoi_tol;
204 :
205 : /// swap matrix
206 : DenseMatrix<Real> _swap_matrix;
207 :
208 : /**
209 : * Before _swap_matrix.svd (a non-const method) is called, we copy _true_swap_matrix =
210 : * _swap_matrix. It is necessary to record _swap_matrix in this fashion so that
211 : * mgd.swap_to_original_basis can be modified using it
212 : */
213 : DenseMatrix<Real> _true_swap_matrix;
214 :
215 : /// inverse of swap matrix
216 : DenseMatrix<Real> _inv_swap_matrix;
217 :
218 : /// used in SVD decomposition of swap matrix
219 : DenseVector<Real> _swap_sigma;
220 :
221 : /// used in SVD decomposition of swap matrix
222 : DenseMatrix<Real> _swap_U;
223 :
224 : /// used in SVD decomposition of swap matrix
225 : DenseMatrix<Real> _swap_VT;
226 : };
|
