Line data Source code
1 : /*
2 : * SPDX-FileCopyrightText: Copyright (c) 2014-2025: G-CSC, Goethe University Frankfurt
3 : * SPDX-License-Identifier: LicenseRef-UG4-LGPL-3.0
4 : *
5 : * Author: Arne Naegel
6 : *
7 : * This file is part of UG4.
8 : *
9 : * UG4 is free software: you can redistribute it and/or modify it under the
10 : * terms of the GNU Lesser General Public License version 3 (as published by the
11 : * Free Software Foundation) with the following additional attribution
12 : * requirements (according to LGPL/GPL v3 §7):
13 : *
14 : * (1) The following notice must be displayed in the Appropriate Legal Notices
15 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
16 : *
17 : * (2) The following notice must be displayed at a prominent place in the
18 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
19 : *
20 : * (3) The following bibliography is recommended for citation and must be
21 : * preserved in all covered files:
22 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
23 : * parallel geometric multigrid solver on hierarchically distributed grids.
24 : * Computing and visualization in science 16, 4 (2013), 151-164"
25 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
26 : * flexible software system for simulating pde based models on high performance
27 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
28 : *
29 : * This program is distributed in the hope that it will be useful,
30 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
31 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
32 : * GNU Lesser General Public License for more details.
33 : */
34 :
35 : #ifndef LIMEX_H_
36 : #define LIMEX_H_
37 :
38 : // extern libraries
39 : #include <vector>
40 : #include <cmath>
41 :
42 : // ug libraries
43 : #include "common/common.h"
44 : #include "common/util/smart_pointer.h"
45 :
46 : #include "lib_disc/time_disc/time_disc_interface.h"
47 : #include "lib_disc/assemble_interface.h"
48 : #include "lib_disc/operator/non_linear_operator/assembled_non_linear_operator.h"
49 : #include "lib_disc/operator/linear_operator/assembled_linear_operator.h"
50 :
51 : #include "lib_algebra/algebra_template_define_helper.h"
52 : #include "lib_algebra/operator/debug_writer.h"
53 :
54 : namespace ug{
55 :
56 : /// \addtogroup limex
57 : /// \{
58 : /// linear implicit time stepping scheme
59 : /**
60 : * This time stepping scheme discretizes equations of the form
61 : * \f[
62 : * M \partial_t u(t) = f(t)
63 : * \f]
64 : * as
65 : * \f[
66 : * (M - \Delta t J) \left( u(t^{k+1}) - u(t^k) \right) = \Delta t \cdot f(t^{k})
67 : * \f]
68 : *
69 : * Thus, for \f$\theta = 1 \f$ this is the Backward-Euler time stepping.
70 : */
71 : template <class TAlgebra>
72 : class LinearImplicitEuler
73 : : public ITimeDiscretization<TAlgebra>,
74 : public DebugWritingObject<TAlgebra>
75 :
76 : {
77 : public:
78 : /// Type of base class
79 : typedef ITimeDiscretization<TAlgebra> base_type;
80 :
81 : /// Type of algebra
82 : typedef TAlgebra algebra_type;
83 :
84 : /// Type of algebra matrix
85 : typedef typename algebra_type::matrix_type matrix_type;
86 :
87 : /// Type of algebra vector
88 : typedef typename algebra_type::vector_type vector_type;
89 :
90 : /// Domain Discretization type
91 : typedef IDomainDiscretization<algebra_type> domain_discretization_type;
92 :
93 : public:
94 : /// CTOR
95 0 : LinearImplicitEuler(SmartPtr<IDomainDiscretization<algebra_type> > spDD)
96 : : ITimeDiscretization<TAlgebra>(spDD),
97 : m_pPrevSol(NULL),
98 0 : m_dt(0.0),
99 0 : m_futureTime(0.0),
100 0 : m_spMatrixJDisc(spDD), m_spMatrixJOp(SPNULL), m_bMatrixJNeedsUpdate(true), m_useLinearMode(false),
101 0 : m_spGammaDisc(SPNULL), m_spGammaOp(SPNULL), m_bGammaNeedsUpdate(true),
102 0 : m_useCachedMatrices(true)
103 0 : {}
104 :
105 : /// CTOR
106 0 : LinearImplicitEuler(SmartPtr<IDomainDiscretization<algebra_type> > spDefectDisc,
107 : SmartPtr<IDomainDiscretization<algebra_type> > spMatrixJDisc)
108 : : ITimeDiscretization<TAlgebra>(spDefectDisc),
109 : m_pPrevSol(NULL),
110 0 : m_dt(0.0),
111 0 : m_futureTime(0.0),
112 0 : m_spMatrixJDisc(spMatrixJDisc), m_spMatrixJOp(SPNULL), m_bMatrixJNeedsUpdate(true), m_useLinearMode(false),
113 0 : m_spGammaDisc(SPNULL), m_spGammaOp(SPNULL), m_bGammaNeedsUpdate(true),
114 0 : m_useCachedMatrices(true)
115 0 : {}
116 :
117 : /// CTOR
118 0 : LinearImplicitEuler(SmartPtr<IDomainDiscretization<algebra_type> > spDefectDisc,
119 : SmartPtr<IDomainDiscretization<algebra_type> > spMatrixJDisc,
120 : SmartPtr<IDomainDiscretization<algebra_type> > spGammaDisc)
121 : : ITimeDiscretization<TAlgebra>(spDefectDisc),
122 : m_pPrevSol(NULL),
123 0 : m_dt(0.0),
124 0 : m_futureTime(0.0),
125 0 : m_spMatrixJDisc(spMatrixJDisc), m_spMatrixJOp(SPNULL), m_bMatrixJNeedsUpdate(true), m_useLinearMode(false),
126 0 : m_spGammaDisc(spGammaDisc), m_spGammaOp(SPNULL), m_bGammaNeedsUpdate(true),
127 0 : m_useCachedMatrices(true)
128 0 : {}
129 :
130 : /// DTOR
131 0 : virtual ~LinearImplicitEuler(){};
132 :
133 : public:
134 : /// \copydoc ITimeDiscretization::num_prev_steps()
135 0 : virtual size_t num_prev_steps() const {return m_prevSteps;}
136 :
137 : /// \copydoc ITimeDiscretization::prepare_step()
138 : virtual void prepare_step(SmartPtr<VectorTimeSeries<vector_type> > prevSol,
139 : number dt);
140 :
141 : /// \copydoc ITimeDiscretization::prepare_step_elem()
142 : virtual void prepare_step_elem(SmartPtr<VectorTimeSeries<vector_type> > prevSol,
143 : number dt, const GridLevel& gl);
144 :
145 : /// \copydoc ITimeDiscretization::finish_step()
146 0 : virtual void finish_step(SmartPtr<VectorTimeSeries<vector_type> > currSol) {};
147 :
148 : /// \copydoc ITimeDiscretization::finish_step_elem()
149 : virtual void finish_step_elem(SmartPtr<VectorTimeSeries<vector_type> > currSol,
150 : const GridLevel& gl);
151 :
152 0 : virtual number future_time() const {return m_futureTime;}
153 :
154 : public:
155 :
156 : /// Meant to assemble J(u) c = d(u), but abused here... (u not used!)
157 : void assemble_jacobian(matrix_type& J, const vector_type& u, const GridLevel& gl);
158 :
159 : /// Meant to assemble d(u), but abused here... (u not used!)
160 : void assemble_defect(vector_type& d, const vector_type& u, const GridLevel& gl);
161 :
162 : /// Should return (M+tau A) delta = tau f
163 : void assemble_linear(matrix_type& A, vector_type& b, const GridLevel& gl);
164 :
165 : void assemble_rhs(vector_type& b, const vector_type& u, const GridLevel& gl);
166 :
167 : void assemble_rhs(vector_type& b, const GridLevel& gl);
168 :
169 : void adjust_solution(vector_type& u, const GridLevel& gl);
170 :
171 0 : virtual size_t num_stages() const {return 1;};
172 0 : virtual void set_stage(size_t stage) {};
173 :
174 :
175 : /// Some simplifications for linear systems. (In this case, the mass matrix is not re-assembled.)
176 0 : void enable_linear_mode() { m_useLinearMode = true; }
177 0 : void disable_linear_mode() { m_useLinearMode = false; }
178 0 : bool use_linear_mode() const { return m_useLinearMode; }
179 :
180 0 : void enable_matrix_cache() { m_useCachedMatrices = true; }
181 0 : void disable_matrix_cache() { m_useCachedMatrices = false; }
182 0 : void set_matrix_cache(bool useCache) { m_useCachedMatrices = useCache; }
183 :
184 : protected:
185 :
186 0 : virtual number update_scaling(std::vector<number>& vSM,
187 : std::vector<number>& vSA,
188 : number dt, number currentTime,
189 : ConstSmartPtr<VectorTimeSeries<vector_type> > prevSol)
190 :
191 : {
192 : // resize scaling factors
193 0 : vSM.resize(1);
194 0 : vSM[0] = 1.0;
195 :
196 0 : vSA.resize(1);
197 0 : vSA[0] = dt;
198 :
199 0 : return currentTime + dt;
200 : }
201 :
202 : static const size_t m_prevSteps=1; ///< number of previous steps needed.
203 : std::vector<number> m_vScaleMass; ///< Scaling for mass part
204 : std::vector<number> m_vScaleStiff; ///< Scaling for stiffness part
205 :
206 : SmartPtr<VectorTimeSeries<vector_type> > m_pPrevSol; ///< Previous solutions
207 : number m_dt; ///< Time Step size
208 : number m_futureTime; ///< Future Time
209 :
210 : // discretization for defect
211 : using base_type::m_spDomDisc;
212 :
213 :
214 : // constant matrix $$ \tau J$$
215 : SmartPtr<IDomainDiscretization<algebra_type> > m_spMatrixJDisc;
216 : SmartPtr<AssembledLinearOperator<algebra_type> > m_spMatrixJOp; ///< Operator
217 : bool m_bMatrixJNeedsUpdate;
218 : bool m_useLinearMode;
219 :
220 : // Matrix: $\Gamma[u0, u0']$
221 : SmartPtr<IDomainDiscretization<algebra_type> > m_spGammaDisc; ///< Gamma disc
222 : SmartPtr<AssembledLinearOperator<algebra_type> > m_spGammaOp; ///< Gamma operator
223 : bool m_bGammaNeedsUpdate;
224 :
225 : SmartPtr<matrix_type> m_spMatrixCacheMk;
226 :
227 : bool m_useCachedMatrices;
228 :
229 : public:
230 :
231 : /// Invalidate all cached operators
232 0 : void invalidate()
233 : {
234 0 : m_spMatrixCacheMk = SPNULL;
235 0 : m_bMatrixJNeedsUpdate = true;
236 : invalidate_gamma();
237 0 : }
238 :
239 : /// Invalidate Gamma operator
240 0 : void invalidate_gamma()
241 0 : { m_spGammaOp = SPNULL; m_bGammaNeedsUpdate = true; }
242 :
243 0 : void set_gamma_disc(SmartPtr<IDomainDiscretization<algebra_type> > spGammaDisc)
244 0 : {m_spGammaDisc = spGammaDisc;}
245 :
246 : using DebugWritingObject<TAlgebra>::set_debug;
247 :
248 : protected:
249 : using DebugWritingObject<TAlgebra>::write_debug;
250 : using DebugWritingObject<TAlgebra>::debug_writer;
251 :
252 : };
253 :
254 : // end group limex
255 : /// \}
256 :
257 : } // namespace ug
258 : #endif
|
