Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Raphael Prohl
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : /*
34 : * (main parts are based on the structure of
35 : * newton_impl.h by Andreas Vogel)
36 : */
37 :
38 : #ifndef __H__UG__LIB_DISC__OPERATOR__NON_LINEAR_OPERATOR__NL_JACOBI__NL_JACOBIL_IMPL_H_
39 : #define __H__UG__LIB_DISC__OPERATOR__NON_LINEAR_OPERATOR__NL_JACOBI__NL_JACOBIL_IMPL_H_
40 :
41 : // extern includes
42 : #include <iostream>
43 :
44 : #include "lib_disc/function_spaces/grid_function_util.h"
45 : #include "lib_disc/common/local_algebra.h"
46 : #include "nl_jacobi.h"
47 :
48 : #define PROFILE_NL_JACOBI
49 : #ifdef PROFILE_NL_JACOBI
50 : #define NL_JACOBI_PROFILE_FUNC() PROFILE_FUNC_GROUP("NL Jacobi")
51 : #define NL_JACOBI_PROFILE_BEGIN(name) PROFILE_BEGIN_GROUP(name, "NL Jacobi")
52 : #define NL_JACOBI_PROFILE_END() PROFILE_END()
53 : #else
54 : #define NL_JACOBI_PROFILE_FUNC()
55 : #define NL_JACOBI_PROFILE_BEGIN(name)
56 : #define NL_JACOBI_PROFILE_END()
57 : #endif
58 :
59 : namespace ug{
60 :
61 : template <typename TAlgebra>
62 0 : NLJacobiSolver<TAlgebra>::
63 : NLJacobiSolver(SmartPtr<IConvergenceCheck<vector_type> > spConvCheck) :
64 : m_spConvCheck(spConvCheck),
65 0 : m_damp(1.0),
66 0 : m_dgbCall(0)
67 : {};
68 :
69 : template <typename TAlgebra>
70 0 : NLJacobiSolver<TAlgebra>::
71 : NLJacobiSolver() :
72 0 : m_spConvCheck(new StdConvCheck<vector_type>(10, 1e-8, 1e-10, true)),
73 0 : m_damp(1.0),
74 0 : m_dgbCall(0)
75 0 : {};
76 :
77 : template <typename TAlgebra>
78 0 : void NLJacobiSolver<TAlgebra>::
79 : set_convergence_check(SmartPtr<IConvergenceCheck<vector_type> > spConvCheck)
80 : {
81 0 : m_spConvCheck = spConvCheck;
82 0 : m_spConvCheck->set_offset(3);
83 0 : m_spConvCheck->set_symbol('#');
84 0 : m_spConvCheck->set_name("Nonlinear Jacobi Solver");
85 0 : }
86 :
87 : template <typename TAlgebra>
88 : bool
89 0 : NLJacobiSolver<TAlgebra>::
90 : init(SmartPtr<IOperator<vector_type> > op)
91 : {
92 : NL_JACOBI_PROFILE_BEGIN(NL_JACOBISolver_init);
93 0 : m_spAssOp = op.template cast_dynamic<AssembledOperator<TAlgebra> >();
94 0 : if(m_spAssOp.invalid())
95 0 : UG_THROW("NLJacobiSolver: currently only works for AssembledDiscreteOperator.");
96 :
97 0 : m_spAss = m_spAssOp->discretization();
98 0 : if(m_spAss.invalid())
99 0 : UG_THROW("AssembledLinearOperator: Assembling routine not set.");
100 :
101 0 : return true;
102 : }
103 :
104 : template <typename TAlgebra>
105 0 : bool NLJacobiSolver<TAlgebra>::prepare(vector_type& u)
106 : {
107 0 : return true;
108 : }
109 :
110 :
111 : template <typename TAlgebra>
112 0 : bool NLJacobiSolver<TAlgebra>::apply(vector_type& u)
113 : {
114 : NL_JACOBI_PROFILE_BEGIN(NL_JACOBISolver_apply);
115 : // increase call count
116 0 : m_dgbCall++;
117 :
118 : // Jacobian
119 0 : if(m_spJ.invalid() || m_spJ->discretization() != m_spAss) {
120 0 : m_spJ = make_sp(new AssembledLinearOperator<TAlgebra>(m_spAss));
121 : m_spJ->set_level(m_spAssOp->level());
122 : }
123 :
124 : // create tmp vectors
125 0 : SmartPtr<vector_type> spD = u.clone_without_values();
126 0 : SmartPtr<vector_type> spC = u.clone_without_values();
127 :
128 : // Set dirichlet values
129 : try{
130 0 : m_spAssOp->prepare(u);
131 : }
132 0 : UG_CATCH_THROW("NLJacobiSolver::apply: Prepare of Operator failed.");
133 :
134 : // Compute first Defect d = L(u)
135 : try{
136 : NL_JACOBI_PROFILE_BEGIN(NL_JACOBIComputeDefect1);
137 0 : m_spAssOp->apply(*spD, u);
138 : NL_JACOBI_PROFILE_END();
139 0 : }UG_CATCH_THROW("NLJacobiSolver::apply: "
140 : "Computation of Start-Defect failed.");
141 :
142 : // write start defect for debug
143 : int loopCnt = 0;
144 : char ext[20]; snprintf(ext, sizeof(ext),"_iter%03d", loopCnt);
145 0 : std::string name("NLJacobi_Defect");
146 0 : name.append(ext);
147 0 : write_debug(*spD, name.c_str());
148 0 : write_debug(u, "NLJacobi_StartSolution");
149 :
150 : // start convergence check
151 0 : m_spConvCheck->start(*spD);
152 :
153 0 : matrix_type& J = m_spJ->get_matrix();
154 :
155 : // loop iteration
156 0 : while(!m_spConvCheck->iteration_ended())
157 : {
158 : // set correction c = 0
159 : NL_JACOBI_PROFILE_BEGIN(NL_JACOBISetCorretionZero);
160 : spC->set(0.0);
161 : NL_JACOBI_PROFILE_END();
162 :
163 : // Compute Jacobian J(u)
164 : try{
165 : NL_JACOBI_PROFILE_BEGIN(NL_JACOBIComputeJacobian);
166 0 : m_spJ->init(u);
167 : NL_JACOBI_PROFILE_END();
168 0 : }UG_CATCH_THROW("NLJacobiSolver::apply: "
169 : "Initialization of Jacobian failed.");
170 :
171 : // Write Jacobian for debug
172 0 : std::string matname("NLJacobi_Jacobian");
173 0 : matname.append(ext);
174 0 : write_debug(m_spJ->get_matrix(), matname.c_str());
175 :
176 : NL_JACOBI_PROFILE_BEGIN(NL_JACOBIInvertBlocks);
177 :
178 : // loop all indizes
179 0 : for (size_t i = 0; i < u.size(); i++)
180 : {
181 : // get i,i-th block of J: J(i,i)
182 : // m_c_i = m_damp * d_i /J_ii
183 0 : InverseMatMult((*spC)[i], m_damp, J(i,i), (*spD)[i]);
184 :
185 : // update solution
186 0 : u[i] -= (*spC)[i];
187 : }
188 : NL_JACOBI_PROFILE_END();
189 :
190 : // compute new Defect
191 : NL_JACOBI_PROFILE_BEGIN(NL_JACOBIComputeDefect);
192 0 : m_spAssOp->prepare(u);
193 0 : m_spAssOp->apply(*spD, u);
194 : NL_JACOBI_PROFILE_END();
195 :
196 : // update counter
197 0 : loopCnt++;
198 : snprintf(ext, sizeof(ext),"_iter%03d", loopCnt);
199 :
200 : // check convergence
201 0 : m_spConvCheck->update(*spD);
202 :
203 : // write defect for debug
204 0 : std::string name("NLJacobi_Defect"); name.append(ext);
205 0 : write_debug(*spD, name.c_str());
206 : }
207 :
208 0 : return m_spConvCheck->post();
209 : }
210 :
211 : template <typename TAlgebra>
212 0 : void NLJacobiSolver<TAlgebra>::write_debug(const vector_type& vec, const char* filename)
213 : {
214 : // add iter count to name
215 0 : std::string name(filename);
216 0 : char ext[20]; snprintf(ext, sizeof(ext),"_call%03d", m_dgbCall);
217 0 : name.append(ext).append(".vec");
218 :
219 : // write
220 0 : base_writer_type::write_debug(vec, name.c_str());
221 0 : }
222 :
223 : template <typename TAlgebra>
224 0 : void NLJacobiSolver<TAlgebra>::write_debug(const matrix_type& mat, const char* filename)
225 : {
226 : // add iter count to name
227 0 : std::string name(filename);
228 0 : char ext[20]; snprintf(ext, sizeof(ext),"_call%03d", m_dgbCall);
229 0 : name.append(ext).append(".mat");
230 :
231 : // write
232 0 : base_writer_type::write_debug(mat, name.c_str());
233 0 : }
234 :
235 : }
236 :
237 : #endif /* __H__UG__LIB_DISC__OPERATOR__NON_LINEAR_OPERATOR__NL_JACOBI__NL_JACOBIL_IMPL_H_ */
|
