ug4
Loading...
Searching...
No Matches
nl_jacobi_impl.h
Go to the documentation of this file.
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
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
59namespace ug{
60
61template <typename TAlgebra>
64 m_spConvCheck(spConvCheck),
65 m_damp(1.0),
66 m_dgbCall(0)
67{};
68
69template <typename TAlgebra>
72 m_spConvCheck(new StdConvCheck<vector_type>(10, 1e-8, 1e-10, true)),
73 m_damp(1.0),
74 m_dgbCall(0)
75{};
76
77template <typename TAlgebra>
80{
81 m_spConvCheck = spConvCheck;
82 m_spConvCheck->set_offset(3);
83 m_spConvCheck->set_symbol('#');
84 m_spConvCheck->set_name("Nonlinear Jacobi Solver");
85}
86
87template <typename TAlgebra>
88bool
91{
92 NL_JACOBI_PROFILE_BEGIN(NL_JACOBISolver_init);
93 m_spAssOp = op.template cast_dynamic<AssembledOperator<TAlgebra> >();
94 if(m_spAssOp.invalid())
95 UG_THROW("NLJacobiSolver: currently only works for AssembledDiscreteOperator.");
96
97 m_spAss = m_spAssOp->discretization();
98 if(m_spAss.invalid())
99 UG_THROW("AssembledLinearOperator: Assembling routine not set.");
100
101 return true;
102}
103
104template <typename TAlgebra>
106{
107 return true;
108}
109
110
111template <typename TAlgebra>
113{
114 NL_JACOBI_PROFILE_BEGIN(NL_JACOBISolver_apply);
115 // increase call count
116 m_dgbCall++;
117
118// Jacobian
119 if(m_spJ.invalid() || m_spJ->discretization() != m_spAss) {
120 m_spJ = make_sp(new AssembledLinearOperator<TAlgebra>(m_spAss));
121 m_spJ->set_level(m_spAssOp->level());
122 }
123
124// create tmp vectors
125 SmartPtr<vector_type> spD = u.clone_without_values();
126 SmartPtr<vector_type> spC = u.clone_without_values();
127
128// Set dirichlet values
129 try{
130 m_spAssOp->prepare(u);
131 }
132 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 m_spAssOp->apply(*spD, u);
139 }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 std::string name("NLJacobi_Defect");
146 name.append(ext);
147 write_debug(*spD, name.c_str());
148 write_debug(u, "NLJacobi_StartSolution");
149
150// start convergence check
151 m_spConvCheck->start(*spD);
152
153 matrix_type& J = m_spJ->get_matrix();
154
155// loop iteration
156 while(!m_spConvCheck->iteration_ended())
157 {
158 // set correction c = 0
159 NL_JACOBI_PROFILE_BEGIN(NL_JACOBISetCorretionZero);
160 spC->set(0.0);
162
163 // Compute Jacobian J(u)
164 try{
165 NL_JACOBI_PROFILE_BEGIN(NL_JACOBIComputeJacobian);
166 m_spJ->init(u);
168 }UG_CATCH_THROW("NLJacobiSolver::apply: "
169 "Initialization of Jacobian failed.");
170
171 // Write Jacobian for debug
172 std::string matname("NLJacobi_Jacobian");
173 matname.append(ext);
174 write_debug(m_spJ->get_matrix(), matname.c_str());
175
176 NL_JACOBI_PROFILE_BEGIN(NL_JACOBIInvertBlocks);
177
178 // loop all indizes
179 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 InverseMatMult((*spC)[i], m_damp, J(i,i), (*spD)[i]);
184
185 // update solution
186 u[i] -= (*spC)[i];
187 }
189
190 // compute new Defect
191 NL_JACOBI_PROFILE_BEGIN(NL_JACOBIComputeDefect);
192 m_spAssOp->prepare(u);
193 m_spAssOp->apply(*spD, u);
195
196 // update counter
197 loopCnt++;
198 snprintf(ext, sizeof(ext),"_iter%03d", loopCnt);
199
200 // check convergence
201 m_spConvCheck->update(*spD);
202
203 // write defect for debug
204 std::string name("NLJacobi_Defect"); name.append(ext);
205 write_debug(*spD, name.c_str());
206 }
207
208 return m_spConvCheck->post();
209}
210
211template <typename TAlgebra>
212void NLJacobiSolver<TAlgebra>::write_debug(const vector_type& vec, const char* filename)
213{
214// add iter count to name
215 std::string name(filename);
216 char ext[20]; snprintf(ext, sizeof(ext),"_call%03d", m_dgbCall);
217 name.append(ext).append(".vec");
218
219// write
220 base_writer_type::write_debug(vec, name.c_str());
221}
222
223template <typename TAlgebra>
224void NLJacobiSolver<TAlgebra>::write_debug(const matrix_type& mat, const char* filename)
225{
226// add iter count to name
227 std::string name(filename);
228 char ext[20]; snprintf(ext, sizeof(ext),"_call%03d", m_dgbCall);
229 name.append(ext).append(".mat");
230
231// write
232 base_writer_type::write_debug(mat, name.c_str());
233}
234
235}
236
237#endif /* __H__UG__LIB_DISC__OPERATOR__NON_LINEAR_OPERATOR__NL_JACOBI__NL_JACOBIL_IMPL_H_ */
location name
Definition checkpoint_util.lua:128
Definition smart_pointer.h:108
matrix operator based on the assembling of a problem
Definition assembled_linear_operator.h:60
Definition convergence_check.h:72
describes a mapping X->Y
Definition operator.h:86
virtual bool apply(vector_type &u)
apply Operator, i.e. op^{-1}(0) = u
Definition nl_jacobi_impl.h:112
NLJacobiSolver()
default constructor
Definition nl_jacobi_impl.h:71
TAlgebra::matrix_type matrix_type
Matrix type.
Definition nl_jacobi.h:89
TAlgebra::vector_type vector_type
Vector type.
Definition nl_jacobi.h:86
void write_debug(const vector_type &vec, const char *filename)
Definition nl_jacobi_impl.h:212
virtual bool init(SmartPtr< IOperator< vector_type > > op)
This operator inverts the Operator op: Y -> X.
Definition nl_jacobi_impl.h:90
void set_convergence_check(SmartPtr< IConvergenceCheck< vector_type > > spConvCheck)
Definition nl_jacobi_impl.h:79
virtual bool prepare(vector_type &u)
prepare Operator
Definition nl_jacobi_impl.h:105
Definition convergence_check.h:179
#define UG_CATCH_THROW(msg)
Definition error.h:64
#define UG_THROW(msg)
Definition error.h:57
the ug namespace
bool InverseMatMult(number &dest, const double &beta, const TMat &mat, const TVec &vec)
you can implement this function with GetInverse and MatMult
#define NL_JACOBI_PROFILE_BEGIN(name)
Definition nl_jacobi_impl.h:51
#define NL_JACOBI_PROFILE_END()
Definition nl_jacobi_impl.h:52
SmartPtr< T, FreePolicy > make_sp(T *inst)
returns a SmartPtr for the passed raw pointer
Definition smart_pointer.h:836
function ProblemDisc new(problemDesc, dom)