Line data Source code
1 :
2 : /*
3 : * SPDX-FileCopyrightText: Copyright (c) 2014-2025: Goethe University Frankfurt
4 : * SPDX-License-Identifier: LicenseRef-UG4-LGPL-3.0
5 : *
6 : * Author: Arne Naegel
7 : *
8 : * This file is part of UG4.
9 : *
10 : * UG4 is free software: you can redistribute it and/or modify it under the
11 : * terms of the GNU Lesser General Public License version 3 (as published by the
12 : * Free Software Foundation) with the following additional attribution
13 : * requirements (according to LGPL/GPL v3 §7):
14 : *
15 : * (1) The following notice must be displayed in the Appropriate Legal Notices
16 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (2) The following notice must be displayed at a prominent place in the
19 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
20 : *
21 : * (3) The following bibliography is recommended for citation and must be
22 : * preserved in all covered files:
23 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
24 : * parallel geometric multigrid solver on hierarchically distributed grids.
25 : * Computing and visualization in science 16, 4 (2013), 151-164"
26 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
27 : * flexible software system for simulating pde based models on high performance
28 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
29 : *
30 : * This program is distributed in the hope that it will be useful,
31 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
32 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
33 : * GNU Lesser General Public License for more details.
34 : */
35 :
36 :
37 : // ug4 headers
38 : #include "common/assert.h"
39 : #include "lib_algebra/cpu_algebra_types.h"
40 : #include "lib_algebra/algebra_common/sparsematrix_util.h"
41 : #include "lib_disc/function_spaces/grid_function_util.h"
42 :
43 : // plugin headers
44 : #include "linear_implicit_timestep.h"
45 : #include "../limex_tools.h"
46 :
47 : namespace ug{
48 :
49 :
50 : template <typename TAlgebra>
51 0 : void LinearImplicitEuler<TAlgebra>::
52 : prepare_step(SmartPtr<VectorTimeSeries<vector_type> > prevSol,
53 : number dt)
54 : {
55 : UG_DLOG(LIB_LIMEX, 5, "LinearlyImplicitEuler:prepare_step" << std::endl);
56 :
57 : PROFILE_BEGIN_GROUP(LinearImplicitEuler_prepare_step, "discretization LinearImplicitEuler");
58 : // perform checks
59 0 : if(prevSol->size() < m_prevSteps)
60 0 : UG_THROW("LinearImplicitEuler::prepare_step:"
61 : " Number of previous solutions must be at least "<<
62 : m_prevSteps <<", but only "<< prevSol->size() << " passed.\n");
63 :
64 : // remember old values
65 0 : m_pPrevSol = prevSol;
66 :
67 : // remember time step size
68 0 : m_dt = dt;
69 :
70 : // update scalings
71 0 : m_futureTime = update_scaling(m_vScaleMass, m_vScaleStiff,
72 : m_dt, m_pPrevSol->time(0),
73 : m_pPrevSol);
74 :
75 0 : std::cout << "PREP: "<< m_vScaleMass[0] <<", " << m_vScaleStiff[0] << ", " <<m_dt << ", " << m_pPrevSol->time(0) << std::endl;
76 :
77 : // prepare time step (elemDisc-wise)
78 : try
79 : {
80 0 : this->m_spDomDisc->prepare_timestep(m_pPrevSol, m_futureTime);
81 0 : this->m_spMatrixJDisc->prepare_timestep(m_pPrevSol, m_futureTime);
82 :
83 0 : if (m_spGammaDisc.valid())
84 0 : { m_spGammaDisc->prepare_timestep(m_pPrevSol, m_futureTime);}
85 : }
86 0 : UG_CATCH_THROW("LinearlyImplicitEuler: Cannot prepare time step.");
87 :
88 : // create matrix J
89 0 : if (m_spMatrixJOp.invalid())
90 0 : { m_spMatrixJOp = make_sp(new AssembledLinearOperator<TAlgebra>(this->m_spMatrixJDisc)); }
91 :
92 : // create matrix Gamma
93 0 : if (m_spGammaDisc != SPNULL && m_spGammaOp == SPNULL)
94 0 : { m_spGammaOp = make_sp(new AssembledLinearOperator<TAlgebra>(this->m_spGammaDisc)); }
95 :
96 : /*{
97 : m_spMatrixJOp->init(*m_pPrevSol->oldest());
98 : }
99 : */
100 0 : }
101 :
102 : template <typename TAlgebra>
103 0 : void LinearImplicitEuler<TAlgebra>::
104 : prepare_step_elem(SmartPtr<VectorTimeSeries<vector_type> > prevSol,
105 : number dt, const GridLevel& gl)
106 : {
107 : UG_DLOG(LIB_LIMEX, 5, "LinearlyImplicitEuler:prepare_step_elem" << std::endl);
108 :
109 : PROFILE_BEGIN_GROUP(LinearImplicitEuler_step_elem, "discretization LinearImplicitEuler");
110 : // perform checks
111 0 : if(prevSol->size() < m_prevSteps)
112 0 : UG_THROW("LinearImplicitEuler::prepare_step_elem:"
113 : " Number of previous solutions must be at least "<<
114 : m_prevSteps <<", but only "<< prevSol->size() << " passed.\n");
115 :
116 : // remember old values
117 0 : m_pPrevSol = prevSol;
118 :
119 : // remember time step size
120 0 : m_dt = dt;
121 :
122 : // update scalings
123 0 : m_futureTime = update_scaling(m_vScaleMass, m_vScaleStiff,
124 : m_dt, m_pPrevSol->time(0),
125 : m_pPrevSol);
126 : // prepare timestep
127 : try{
128 0 : this->m_spDomDisc->prepare_timestep(m_pPrevSol, m_futureTime, gl);
129 0 : this->m_spMatrixJDisc->prepare_timestep(m_pPrevSol, m_futureTime, gl);
130 0 : } UG_CATCH_THROW("LinearImplicitEuler: Cannot prepare timestep.");
131 :
132 : // Aux linear operator
133 0 : if (m_spMatrixJOp.invalid()) {
134 0 : m_spMatrixJOp = make_sp(new AssembledLinearOperator<TAlgebra>(this->m_spMatrixJDisc));
135 : }
136 :
137 :
138 0 : if (m_spGammaDisc.valid() && m_spGammaOp == SPNULL) {
139 0 : m_spGammaOp = make_sp(new AssembledLinearOperator<TAlgebra>(this->m_spGammaDisc));
140 : }
141 :
142 : // m_spMatrixJOp->init(*m_pPrevSol->oldest());
143 : //std::cout << "PREPELEM: "<< m_vScaleMass[0] <<", " << m_vScaleStiff[0] << ", " <<m_dt << ", " << m_pPrevSol->time(0) << std::endl;
144 :
145 0 : }
146 :
147 : template <typename TAlgebra>
148 0 : void LinearImplicitEuler<TAlgebra>::
149 : adjust_solution(vector_type& u, const GridLevel& gl)
150 : {
151 : UG_DLOG(LIB_LIMEX, 5, "LinearlyImplicitEuler:adjust_solution" << &u << std::endl);
152 : PROFILE_BEGIN_GROUP(LinearImplicitEuler_adjust_solution, "discretization LinearImplicitEuler");
153 :
154 : // adjust solution
155 : try{
156 0 : this->m_spDomDisc->adjust_solution(u, m_futureTime, gl);
157 0 : }UG_CATCH_THROW("LinearImplicitEuler: Cannot adjust solution.");
158 :
159 0 : }
160 :
161 :
162 :
163 : template <typename TAlgebra>
164 0 : void LinearImplicitEuler<TAlgebra>::
165 : finish_step_elem(SmartPtr<VectorTimeSeries<vector_type> > currSol,
166 : const GridLevel& gl)
167 : {
168 : UG_DLOG(LIB_LIMEX, 5, "LinearlyImplicitEuler:finish_step_elem" << std::endl);
169 :
170 : // perform checks whether 'currSol' is a solutionTimeSeries only with the new values
171 0 : if(currSol->time(0) != m_futureTime)
172 0 : UG_THROW("LinearImplicitEuler::finish_step_elem:"
173 : " The solution of the SolutionTimeSeries used in this function"
174 : " does not coincide with the current solution! ");
175 :
176 : // finish timestep using the current solution
177 : try{
178 0 : this->m_spDomDisc->finish_timestep(currSol, gl);
179 0 : }UG_CATCH_THROW("LinearImplicitEuler: Cannot finish timestep.");
180 0 : }
181 :
182 :
183 : /*
184 : *
185 : * Non-linear system
186 : *
187 : * \partial m(u) = f(u)
188 : *
189 : * */
190 :
191 :
192 : /** WARNING: This function is abused
193 : * Must return: $Mk + \tau J$
194 : * */
195 : template <typename TAlgebra>
196 0 : void LinearImplicitEuler<TAlgebra>::
197 : assemble_jacobian(matrix_type& J_limex, const vector_type& u, const GridLevel& gl)
198 : {
199 : UG_DLOG(LIB_LIMEX, 5, "LinearlyImplicitEuler:assemble_jacobian" << std::endl);
200 :
201 : PROFILE_BEGIN_GROUP(LinearImplicitEuler_assemble_jacobian, "discretization LinearImplicitEuler");
202 :
203 : // perform checks
204 0 : if(m_pPrevSol->size() < m_prevSteps)
205 0 : UG_THROW("LinearImplicitEuler::assemble_jacobian:"
206 : " Number of previous solutions must be at least "<<
207 : m_prevSteps <<", but only "<< m_pPrevSol->size() << " passed.");
208 :
209 : // push unknown solution to solution time series
210 : // ATTENTION: Here, we must cast away the constness of the solution, but note,
211 : // that we pass pPrevSol as a const object in assemble_... Thus,
212 : // the solution will not be changed there and we pop it from the
213 : // Solution list afterwards, such that nothing happens to u
214 : // \todo: avoid this hack, use smart ptr properly
215 : int DummyRefCount = 2;
216 : SmartPtr<vector_type> pU(const_cast<vector_type*>(&u), &DummyRefCount);
217 0 : m_pPrevSol->push(pU, m_futureTime);
218 :
219 : // assemble "Jacobian" (M_{k-1} + \tau J) using current iterate
220 : try{
221 :
222 : // check assertions
223 : UG_ASSERT(m_spMatrixJDisc.valid(), "Huhh: Invalid matrix discretization")
224 :
225 0 : if (!m_useCachedMatrices)
226 : {
227 : // un-cached, (i.e., approximate Newton)
228 0 : this->m_spMatrixJDisc->assemble_jacobian(J_limex, m_pPrevSol, m_dt, gl);
229 : UG_DLOG(LIB_LIMEX, 5, "Computed J_limex =" << J_limex << " for dt=" << m_dt << std::endl);
230 0 : write_debug(J_limex, "myMatrixAssembled.mat");
231 : }
232 : else
233 : {
234 : /* Assemble (Mk + \tau J)
235 : * where Mk has Dirichlet rows (e.g., for inout bnd cond) */
236 :
237 : // First, compute J_limex as the derivative of the storage (mass) terms,
238 0 : this->m_spMatrixJDisc->assemble_jacobian(J_limex, m_pPrevSol, 0.0, gl);
239 : UG_DLOG(LIB_LIMEX, 3, "> Computed Mk (" << &J_limex << " : "<<J_limex <<
240 : " at " << m_pPrevSol->oldest_time() << ", " << GetNNZs(J_limex) << " nonzeros)" << std::endl);
241 0 : write_debug(J_limex, "myMk.mat");
242 :
243 : // Second, compute the derivative of the time dependent part.
244 : const double mydt = 1.0; // use time step 1.0 for assembly
245 0 : matrix_type &J_stiff = m_spMatrixJOp->get_matrix();
246 :
247 0 : if (m_bMatrixJNeedsUpdate)
248 : {
249 : // First part of J: -df/du = Jm + \tau Ja
250 0 : this->m_spMatrixJDisc->assemble_jacobian(J_stiff, m_pPrevSol, mydt, gl);
251 : UG_DLOG(LIB_LIMEX, 3, "> Cached J_0 = -df/du (" << &J_stiff << ":"<< J_stiff<< ")"<< std::endl);
252 0 : write_debug(J_stiff, "myStiff0.mat");
253 :
254 : // Subtracting mass matrix yields $\tau*Ja$.
255 0 : MatAddNonDirichlet<matrix_type>(J_stiff, 1.0, J_stiff, -1.0, J_limex);
256 0 : write_debug(J_stiff, "myStiff1.mat");
257 :
258 : // Second part of J: Gamma update (if applicable)
259 : // (Both time dependent and time-independent.)
260 0 : if (m_spGammaDisc.valid())
261 : {
262 : UG_ASSERT(m_spGammaOp != SPNULL, "Huhh: No operator??? ");
263 0 : matrix_type &myGamma = m_spGammaOp->get_matrix();
264 0 : this->m_spGammaDisc->assemble_jacobian(myGamma, m_pPrevSol, mydt, gl);
265 : UG_DLOG(LIB_LIMEX, 3, "Assembled Gamma0 ("<< myGamma);
266 : UG_DLOG(LIB_LIMEX, 3, " at " << m_pPrevSol->oldest_time() <<", " << GetNNZs(myGamma) << " nnz)" << std::endl);
267 0 : write_debug(myGamma, "myGamma.mat");
268 :
269 0 : MatAddNonDirichlet<matrix_type>(J_stiff, 1.0, J_stiff, 1.0, myGamma);
270 0 : write_debug(J_stiff, "myStiff2.mat");
271 : UG_DLOG(LIB_LIMEX, 3, "Added Gamma0"<< std::endl)
272 : }
273 :
274 : // do not need to recompute
275 0 : m_bMatrixJNeedsUpdate = false;
276 : } // m_bMatrixJNeedsUpdate == true;
277 :
278 0 : write_debug(J_stiff, "myStiffX.mat");
279 :
280 : // Updating Jtotal = Mk+ \tau Ja (for non-dirichlet)
281 : UG_ASSERT (J_limex.num_rows() == J_stiff.num_rows(), "Huhh: Row dimension does not match: "
282 : << J_limex.num_rows() <<"!=" << J_stiff.num_rows());
283 : UG_ASSERT (J_limex.num_cols() == J_stiff.num_cols(), "Huhh: Col dimension does not match:"
284 : << J_limex.num_cols() <<"!=" << J_stiff.num_cols());
285 :
286 : // Note: J_limex has Dirichlet values.
287 0 : MatAddNonDirichlet<matrix_type>(J_limex, 1.0, J_limex, m_dt, J_stiff);
288 : UG_DLOG(LIB_LIMEX, 3, "Computed $J=Mk+ tau J_a$ (" << J_limex << " with tau=" << m_dt);
289 : UG_DLOG(LIB_LIMEX, 3, " at " << m_pPrevSol->oldest_time() <<", " << GetNNZs(J_limex) << " nonzeros)" << std::endl);
290 :
291 0 : write_debug(J_limex, "myMatrixCached.mat");
292 : } // if(! m_useCachedMatrices)
293 :
294 : }
295 0 : UG_CATCH_THROW("LinearlyImplicitEuler: Cannot assemble jacobian(s).");
296 :
297 : //UG_ASSERT(0, "Crashing")
298 :
299 : // pop unknown solution to solution time series
300 : m_pPrevSol->remove_latest();
301 0 : }
302 :
303 : /** WARNING: This function is abused
304 : * Must return : d_A(k-1):= tau * F(k-1) - A(k-1) u(k-1) */
305 : template <typename TAlgebra>
306 0 : void LinearImplicitEuler<TAlgebra>::
307 : assemble_defect(vector_type& d, const vector_type& u, const GridLevel& gl)
308 : {
309 : UG_DLOG(LIB_LIMEX, 5, "LinearlyImplicitEuler:assemble_defect" << std::endl);
310 :
311 : PROFILE_BEGIN_GROUP(LinearImplicitEuler_assemble_defect, "discretization LinearImplicitEuler");
312 :
313 : // perform checks
314 0 : if(m_pPrevSol->size() < m_prevSteps)
315 0 : UG_THROW("LinearImplicitEuler::assemble_defect:"
316 : " Number of previous solutions must be at least "<<
317 : m_prevSteps <<", but only "<< m_pPrevSol->size() << " passed.");
318 :
319 : // push unknown solution to solution time series
320 : // ATTENTION: Here, we must cast away the constness of the solution, but note,
321 : // that we pass pPrevSol as a const object in assemble_... Thus,
322 : // the solution will not be changed there and we pop it from the
323 : // Solution list afterwards, such that nothing happens to u
324 : // \todo: avoid this hack, use smart ptr properly
325 : int DummyRefCount = 2;
326 : SmartPtr<vector_type> pU(const_cast<vector_type*>(&u), &DummyRefCount);
327 0 : m_pPrevSol->push(pU, m_futureTime);
328 :
329 : // future solution part
330 : try{
331 :
332 : // d:=\tau {f(k-1) - A(k-1) u(k-1)}
333 0 : std::vector<number> vScaleMass(1, 0.0); // 0
334 0 : std::vector<number> vScaleStiff(1, m_dt);
335 0 : this->m_spDomDisc->assemble_defect(d, m_pPrevSol, vScaleMass, vScaleStiff, gl);
336 :
337 : // d := d + (M - \tau J) (u(k-1)-u)
338 : /*
339 : vector_type deltau = *m_pPrevSol->oldest()->clone();
340 : deltau -= u;
341 :
342 : this->m_spDomDisc->assemble_jacobian(m_spMatrixJOp->get_matrix(), m_pPrevSol, m_dt, gl);
343 : m_spMatrixJOp->apply_sub(d, deltau);
344 : */
345 0 : }UG_CATCH_THROW("LinearlyImplicitEuler: Cannot assemble defect.");
346 :
347 : // pop unknown solution to solution time series
348 : m_pPrevSol->remove_latest();
349 :
350 0 : }
351 :
352 :
353 : /* RHS (non-linear system) */
354 : template <typename TAlgebra>
355 0 : void LinearImplicitEuler<TAlgebra>::
356 : assemble_rhs(vector_type& b, const vector_type& u, const GridLevel& gl)
357 : {
358 : UG_ASSERT(0, "Really wanna use me???");
359 : PROFILE_BEGIN_GROUP(LinearImplicitEuler_assemble_rhs, "discretization LinearImplicitEuler");
360 : // perform checks
361 0 : if(m_pPrevSol->size() < m_prevSteps)
362 0 : UG_THROW("LinearImplicitEuler::assemble_linear:"
363 : " Number of previous solutions must be at least "<<
364 : m_prevSteps <<", but only "<< m_pPrevSol->size() << " passed.");
365 :
366 : // assemble jacobian using current iterate
367 : try{
368 0 : this->m_spDomDisc->assemble_rhs(b, m_pPrevSol, m_vScaleMass, m_vScaleStiff, gl);
369 0 : }UG_CATCH_THROW("LinearImplicitEuler: Cannot assemble jacobian.");
370 :
371 0 : }
372 :
373 : /*
374 : *
375 : * Linear system
376 : *
377 : *
378 : */
379 :
380 : template <typename TAlgebra>
381 0 : void LinearImplicitEuler<TAlgebra>::
382 : assemble_linear(matrix_type& J, vector_type& b, const GridLevel& gl)
383 : {
384 : UG_DLOG(LIB_LIMEX, 5, "LinearlyImplicitEuler:assemble_linear" << std::endl);
385 :
386 : PROFILE_BEGIN_GROUP(LinearImplicitEuler_assemble_linear, "discretization LinearImplicitEuler");
387 :
388 : // perform checks
389 0 : if(m_pPrevSol->size() < m_prevSteps)
390 0 : UG_THROW("LinearImplicitEuler::assemble_defect:"
391 : " Number of previous solutions must be at least "<<
392 : m_prevSteps <<", but only "<< m_pPrevSol->size() << " passed.");
393 :
394 : // future solution part
395 : try {
396 : UG_ASSERT(m_pPrevSol->size()<2, "Only need one solution,,,");
397 :
398 : // J(u_{k+1} - u_k) = -\tau f_k
399 : // NOTE: both routines have invoked the contraints, so constraint adjustment are assumed to be correct!
400 0 : this->assemble_jacobian(J, *m_pPrevSol->oldest(), gl);
401 0 : this->assemble_defect(b, *m_pPrevSol->oldest(), gl);
402 :
403 : // J u_{k+1} = -\tau f_k + J u_k
404 : //AxpyCommonSparseMatrix(J, b, 0.0, *m_pPrevSol->oldest(), +1.0, *m_pPrevSol->oldest());
405 0 : J.axpy(b, -1.0, b, 1.0, *m_pPrevSol->oldest());
406 : }
407 0 : UG_CATCH_THROW("LinearImplicitEuler: Cannot assemble defect.");
408 0 : }
409 :
410 :
411 :
412 : /* RHS (linear system) */
413 : template <typename TAlgebra>
414 0 : void LinearImplicitEuler<TAlgebra>::
415 : assemble_rhs(vector_type& b, const GridLevel& gl)
416 : {
417 : UG_ASSERT(0, "Really wanna use me???");
418 : PROFILE_BEGIN_GROUP(LinearImplicitEuler_assemble_rhs, "discretization LinearImplicitEuler");
419 : // perform checks
420 0 : if(m_pPrevSol->size() < m_prevSteps)
421 0 : UG_THROW("LinearImplicitEuler::assemble_rhs:"
422 : " Number of previous solutions must be at least "<<
423 : m_prevSteps <<", but only "<< m_pPrevSol->size() << " passed.");
424 :
425 : // assemble jacobian using current iterate
426 : try{
427 0 : this->m_spDomDisc->assemble_rhs(b, m_pPrevSol, m_vScaleMass, m_vScaleStiff, gl);
428 0 : }UG_CATCH_THROW("LinearImplicitEuler: Cannot assemble jacobian.");
429 :
430 0 : }
431 :
432 :
433 : ////////////////////////////////////////////////////////////////////////
434 : // template instantiations for all current algebra types.
435 :
436 : UG_ALGEBRA_CPP_TEMPLATE_DEFINE_ALL(LinearImplicitEuler)
437 :
438 :
439 : }; // namespace ug
|
