OGS
TimeDiscretizedODESystem.cpp
Go to the documentation of this file.
1
12
15#include "NumLib/Exceptions.h"
17
18namespace detail
19{
21template <typename Solutions, typename Vector>
22void applyKnownSolutions(std::vector<Solutions> const* const known_solutions,
23 Vector& x)
24{
25 if (!known_solutions)
26 {
27 return;
28 }
29
30 for (auto const& bc : *known_solutions)
31 {
32 for (std::size_t i = 0; i < bc.ids.size(); ++i)
33 {
34 // TODO that might have bad performance for some Vector types, e.g.,
35 // PETSc.
36 MathLib::setVector(x, bc.ids[i], bc.values[i]);
37 }
38 }
40}
41} // namespace detail
42
43namespace NumLib
44{
45TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
47 TimeDiscretizedODESystem(const int process_id, ODE& ode,
48 TimeDisc& time_discretization)
49 : _ode(ode),
50 _time_disc(time_discretization),
51 _mat_trans(createMatrixTranslator<ODETag>(time_discretization))
52{
54 _ode.getMatrixSpecifications(process_id), _Jac_id);
56 _ode.getMatrixSpecifications(process_id), _M_id);
58 _ode.getMatrixSpecifications(process_id), _K_id);
60 _ode.getMatrixSpecifications(process_id), _b_id);
61}
62
65 NonlinearSolverTag::Newton>::~TimeDiscretizedODESystem()
66{
71}
72
75 assemble(std::vector<GlobalVector*> const& x_new_timestep,
76 std::vector<GlobalVector*> const& x_prev,
77 int const process_id)
78{
79 namespace LinAlg = MathLib::LinAlg;
80
81 auto const t = _time_disc.getCurrentTime();
82 auto const dt = _time_disc.getCurrentTimeIncrement();
83 auto const& x_curr = *x_new_timestep[process_id];
84
85 std::vector<GlobalVector*> xdot(x_new_timestep.size());
86 _xdot_ids.resize(x_new_timestep.size());
87 for (std::size_t i = 0; i < xdot.size(); i++)
88 {
89 xdot[i] =
91 _time_disc.getXdot(*x_new_timestep[i], *x_prev[i], *xdot[i]);
92 }
93
94 _M->setZero();
95 _K->setZero();
96 _b->setZero();
97 _Jac->setZero();
98
99 _ode.preAssemble(t, dt, x_curr);
100 try
101 {
102 _ode.assembleWithJacobian(t, dt, x_new_timestep, xdot, process_id, *_M,
103 *_K, *_b, *_Jac);
104 }
105 catch (AssemblyException const&)
106 {
107 for (auto& v : xdot)
108 {
110 }
111 throw;
112 }
113
118
119 for (auto& v : xdot)
120 {
122 }
123}
124
127 NonlinearSolverTag::Newton>::getResidual(GlobalVector const& x_new_timestep,
128 GlobalVector const& x_prev,
129 GlobalVector& res) const
130{
131 // TODO Maybe the duplicate calculation of xdot here and in assembleJacobian
132 // can be optimized. However, that would make the interface a bit more
133 // fragile.
135 _time_disc.getXdot(x_new_timestep, x_prev, xdot);
136
137 _mat_trans->computeResidual(*_M, *_K, *_b, x_new_timestep, xdot, res);
138
140}
141
144 NonlinearSolverTag::Newton>::getJacobian(GlobalMatrix& Jac) const
145{
146 _mat_trans->computeJacobian(*_Jac, Jac);
147}
148
151 NonlinearSolverTag::Newton>::computeKnownSolutions(GlobalVector const& x,
152 int const process_id)
153{
154 _known_solutions =
155 _ode.getKnownSolutions(_time_disc.getCurrentTime(), x, process_id);
156}
157
161{
162 ::detail::applyKnownSolutions(_known_solutions, x);
163}
164
167 applyKnownSolutionsNewton(GlobalMatrix& Jac, GlobalVector& res,
168 GlobalVector& minus_delta_x) const
169{
170 if (!_known_solutions)
171 {
172 return;
173 }
174
176 std::vector<IndexType> ids;
177 for (auto const& bc : *_known_solutions)
178 {
179 ids.insert(end(ids), begin(bc.ids), end(bc.ids));
180 }
181
182 // For the Newton method the values must be zero
183 std::vector<double> values(ids.size(), 0);
184 MathLib::applyKnownSolution(Jac, res, minus_delta_x, ids, values);
185}
186
189 TimeDiscretizedODESystem(const int process_id, ODE& ode,
190 TimeDisc& time_discretization)
191 : _ode(ode),
192 _time_disc(time_discretization),
193 _mat_trans(createMatrixTranslator<ODETag>(time_discretization))
194{
196 ode.getMatrixSpecifications(process_id), _M_id);
198 ode.getMatrixSpecifications(process_id), _K_id);
200 ode.getMatrixSpecifications(process_id), _b_id);
201}
202
205 NonlinearSolverTag::Picard>::~TimeDiscretizedODESystem()
206{
210}
211
214 assemble(std::vector<GlobalVector*> const& x_new_timestep,
215 std::vector<GlobalVector*> const& x_prev,
216 int const process_id)
217{
218 namespace LinAlg = MathLib::LinAlg;
219
220 auto const t = _time_disc.getCurrentTime();
221 auto const dt = _time_disc.getCurrentTimeIncrement();
222 auto const& x_curr = *x_new_timestep[process_id];
223 std::vector<GlobalVector*> xdot(x_new_timestep.size());
224 _xdot_ids.resize(x_new_timestep.size());
225
226 for (std::size_t i = 0; i < xdot.size(); i++)
227 {
228 xdot[i] =
230 _time_disc.getXdot(*x_new_timestep[i], *x_prev[i], *xdot[i]);
231 }
232
233 _M->setZero();
234 _K->setZero();
235 _b->setZero();
236
237 _ode.preAssemble(t, dt, x_curr);
238 _ode.assemble(t, dt, x_new_timestep, xdot, process_id, *_M, *_K, *_b);
239
243
244 for (auto& v : xdot)
245 {
247 }
248}
249
252 NonlinearSolverTag::Picard>::computeKnownSolutions(GlobalVector const& x,
253 int const process_id)
254{
255 _known_solutions =
256 _ode.getKnownSolutions(_time_disc.getCurrentTime(), x, process_id);
257}
258
262{
263 ::detail::applyKnownSolutions(_known_solutions, x);
264}
265
268 applyKnownSolutionsPicard(GlobalMatrix& A,
269 GlobalVector& rhs,
270 GlobalVector& x) const
271{
272 if (!_known_solutions)
273 {
274 return;
275 }
276
278 std::vector<IndexType> ids;
279 std::vector<double> values;
280 for (auto const& bc : *_known_solutions)
281 {
282 ids.insert(end(ids), begin(bc.ids), end(bc.ids));
283 values.insert(end(values), begin(bc.values), end(bc.values));
284 }
285 MathLib::applyKnownSolution(A, rhs, x, ids, values);
286}
287
288} // namespace NumLib
Global vector based on Eigen vector.
Definition: EigenVector.h:28
virtual void releaseMatrix(GlobalMatrix const &A)=0
virtual GlobalMatrix & getMatrix(std::size_t &id)=0
Get an uninitialized matrix with the given id.
virtual GlobalVector & getVector(std::size_t &id)=0
Get an uninitialized vector with the given id.
virtual void releaseVector(GlobalVector const &x)=0
std::unique_ptr< MatrixTranslator< ODETag > > createMatrixTranslator(TimeDiscretization const &timeDisc)
std::array< double, 3 > Vector
Definition: VariableType.h:28
void finalizeAssembly(PETScMatrix &A)
Definition: LinAlg.cpp:163
void applyKnownSolution(EigenMatrix &A, EigenVector &b, EigenVector &, const std::vector< EigenMatrix::IndexType > &vec_knownX_id, const std::vector< double > &vec_knownX_x)
Definition: EigenTools.cpp:17
void setVector(PETScVector &v, std::initializer_list< double > values)
static const double v
static const double t
void applyKnownSolutions(std::vector< Solutions > const *const known_solutions, Vector &x)
Applies known solutions to the solution vector x.
static NUMLIB_EXPORT MatrixProvider & provider
static NUMLIB_EXPORT VectorProvider & provider