Finite Element Domain Decomposition Library
FEDDLib
Loading...
Searching...
No Matches
NonLinearProblem_def.hpp
1#ifndef NONLINEARPROBLEM_DEF_hpp
2#define NONLINEARPROBLEM_DEF_hpp
3#include "NonLinearProblem_decl.hpp"
4
13
14namespace FEDD
15{
16
17 template <class SC, class LO, class GO, class NO>
18 NonLinearProblem<SC, LO, GO, NO>::NonLinearProblem(CommConstPtr_Type comm) : Problem<SC, LO, GO, NO>(comm),
19 previousSolution_(),
20 residualVec_(),
21 nonLinearTolerance_(1.e-6),
22 coeff_(0)
23 {
24 }
25
26 template <class SC, class LO, class GO, class NO>
27 NonLinearProblem<SC, LO, GO, NO>::NonLinearProblem(ParameterListPtr_Type &parameterList, CommConstPtr_Type comm) : Problem<SC, LO, GO, NO>(parameterList, comm),
28 previousSolution_(),
29 residualVec_(),
30 nonLinearTolerance_(1.e-6),
31 coeff_(0)
32 {
33 }
34 template <class SC, class LO, class GO, class NO>
35 NonLinearProblem<SC, LO, GO, NO>::~NonLinearProblem()
36 {
37 }
38
39 template <class SC, class LO, class GO, class NO>
40 void NonLinearProblem<SC, LO, GO, NO>::initializeProblem(int nmbVectors)
41 {
42 this->system_.reset(new BlockMatrix_Type(1));
43
44 this->initializeVectors(nmbVectors);
45
46 this->initializeVectorsNonLinear(nmbVectors);
47
48 // Init ThyraVectorSpcaes for NOX.
49 this->initVectorSpaces();
50 }
51
52 template <class SC, class LO, class GO, class NO>
53 void NonLinearProblem<SC, LO, GO, NO>::initializeVectorsNonLinear(int nmbVectors)
54 {
55
56 UN size = this->domainPtr_vec_.size();
57 this->previousSolution_.reset(new BlockMultiVector_Type(size));
58 this->residualVec_.reset(new BlockMultiVector_Type(size));
59
60 for (UN i = 0; i < size; i++)
61 {
62 if (this->dofsPerNode_vec_[i] > 1)
63 {
64 MapConstPtr_Type map = this->domainPtr_vec_[i]->getMapVecFieldUnique();
65 MultiVectorPtr_Type prevSolutionPart = Teuchos::rcp(new MultiVector_Type(map));
66 this->previousSolution_->addBlock(prevSolutionPart, i);
67 MultiVectorPtr_Type residualPart = Teuchos::rcp(new MultiVector_Type(map));
68 this->residualVec_->addBlock(residualPart, i);
69 }
70 else
71 {
72 MapConstPtr_Type map = this->domainPtr_vec_[i]->getMapUnique();
73 MultiVectorPtr_Type prevSolutionPart = Teuchos::rcp(new MultiVector_Type(map));
74 this->previousSolution_->addBlock(prevSolutionPart, i);
75 MultiVectorPtr_Type residualPart = Teuchos::rcp(new MultiVector_Type(map));
76 this->residualVec_->addBlock(residualPart, i);
77 }
78 }
79
80 this->residualVec_->putScalar(0.);
81 }
82
83 template <class SC, class LO, class GO, class NO>
84 void NonLinearProblem<SC, LO, GO, NO>::infoNonlinProblem()
85 {
86
87 bool verbose(this->comm_->getRank() == 0);
88 if (verbose)
89 {
90 std::cout << "\t ### ### ###" << std::endl;
91 std::cout << "\t ### Nonlinear Problem Information ###" << std::endl;
92 std::cout << "\t ### Linearization: " << this->parameterList_->sublist("General").get("Linearization", "default") << std::endl;
93 std::cout << "\t ### Relative tol: " << this->parameterList_->sublist("Parameter").get("relNonLinTol", 1.e-6) << "\t absolute tol: " << this->parameterList_->sublist("Parameter").get("absNonLinTol", 1.e-4) << "(not used for Newton or Fixed-Point)" << std::endl;
94 }
95 }
96
97 template <class SC, class LO, class GO, class NO>
98 void NonLinearProblem<SC, LO, GO, NO>::calculateNonLinResidualVec(SmallMatrix<double> &coeff, std::string type, double time)
99 {
100
101 coeff_ = coeff;
102 this->calculateNonLinResidualVec(type, time);
103 }
104
105 template <class SC, class LO, class GO, class NO>
106 void NonLinearProblem<SC, LO, GO, NO>::reAssembleAndFill(BlockMatrixPtr_Type bMat, std::string type)
107 {
108
109 this->assemble(type);
110 TEUCHOS_TEST_FOR_EXCEPTION(bMat->size() != this->system_->size(), std::logic_error, "Sizes of BlockMatrices are differen. reAssembleAndFill(...)");
111
112 for (int i = 0; i < this->system_->size(); i++)
113 {
114 for (int j = 0; j < this->system_->size(); j++)
115 {
116 if (this->system_->blockExists(i, j))
117 {
118 // MatrixPtr_Type mat = Teuchos::rcp(new Matrix_Type ( this->system_->getBlock( i, j ) ) );
119 bMat->addBlock(this->system_->getBlock(i, j), i, j);
120 }
121 }
122 }
123 }
124
125 template <class SC, class LO, class GO, class NO>
126 double NonLinearProblem<SC, LO, GO, NO>::calculateResidualNorm() const
127 {
128
129 Teuchos::Array<SC> residual(1);
130 residualVec_->norm2(residual());
131 TEUCHOS_TEST_FOR_EXCEPTION(residual.size() != 1, std::logic_error, "We need to change the code for numVectors>1.");
132 return residual[0];
133 }
134
135 template <class SC, class LO, class GO, class NO>
136 int NonLinearProblem<SC, LO, GO, NO>::solveUpdate()
137 {
138
139 // solution COPY!
140 *previousSolution_ = *this->solution_;
141
142 int its = this->solve(residualVec_);
143
144 return its;
145 }
146
147 template <class SC, class LO, class GO, class NO>
148 int NonLinearProblem<SC, LO, GO, NO>::solveAndUpdate(const std::string &criterion, double &criterionValue)
149 {
150 // BlockMatrixPtr_Type system
151 int its = solveUpdate();
152
153 if (criterion == "Update")
154 {
155 Teuchos::Array<SC> updateNorm(1);
156 this->solution_->norm2(updateNorm());
157 criterionValue = updateNorm[0];
158 }
159 this->solution_->update(1., *previousSolution_, 1.);
160
161 return its;
162 }
163
164 template <class SC, class LO, class GO, class NO>
165 typename NonLinearProblem<SC, LO, GO, NO>::BlockMultiVectorPtr_Type NonLinearProblem<SC, LO, GO, NO>::getResidualVector() const
166 {
167
168 return residualVec_;
169 }
170
171 template <class SC, class LO, class GO, class NO>
172 Thyra::ModelEvaluatorBase::InArgs<SC>
173 NonLinearProblem<SC, LO, GO, NO>::getNominalValues() const
174 {
175 return nominalValues_;
176 }
177
178 template <class SC, class LO, class GO, class NO>
179 Teuchos::RCP<const ::Thyra::VectorSpaceBase<SC>> NonLinearProblem<SC, LO, GO, NO>::get_x_space() const
180 {
181 return xSpace_;
182 }
183
184 template <class SC, class LO, class GO, class NO>
185 Teuchos::RCP<const ::Thyra::VectorSpaceBase<SC>> NonLinearProblem<SC, LO, GO, NO>::get_f_space() const
186 {
187 return fSpace_;
188 }
189
190 template <class SC, class LO, class GO, class NO>
191 Thyra::ModelEvaluatorBase::InArgs<SC>
192 NonLinearProblem<SC, LO, GO, NO>::createInArgs() const
193 {
194 return prototypeInArgs_;
195 }
196
197 // Private functions overridden from ModelEvaulatorDefaultBase
198
199 template <class SC, class LO, class GO, class NO>
200 Thyra::ModelEvaluatorBase::OutArgs<SC>
201 NonLinearProblem<SC, LO, GO, NO>::createOutArgsImpl() const
202 {
203 return prototypeOutArgs_;
204 }
205
206 template <class SC, class LO, class GO, class NO>
207 void NonLinearProblem<SC, LO, GO, NO>::initNOXParameters()
208 {
209
210 using Teuchos::RCP;
211 using Teuchos::rcp;
212 using ::Thyra::VectorBase;
213 typedef ::Thyra::ModelEvaluatorBase MEB;
214
215 MEB::InArgsSetup<SC> inArgs;
216 inArgs.setModelEvalDescription(this->description());
217 inArgs.setSupports(MEB::IN_ARG_x);
218 this->prototypeInArgs_ = inArgs;
219
220 MEB::OutArgsSetup<SC> outArgs;
221 outArgs.setModelEvalDescription(this->description());
222 outArgs.setSupports(MEB::OUT_ARG_f);
223 outArgs.setSupports(MEB::OUT_ARG_W_op);
224 outArgs.setSupports(MEB::OUT_ARG_W_prec);
225 this->prototypeOutArgs_ = outArgs;
226
227 this->nominalValues_ = inArgs;
228 }
229
230 template <class SC, class LO, class GO, class NO>
231 void NonLinearProblem<SC, LO, GO, NO>::initVectorSpaces()
232 {
233 this->initNOXParameters();
234
235 std::string type = this->parameterList_->sublist("General").get("Preconditioner Method", "Monolithic");
236 if (!type.compare("Monolithic"))
237 initVectorSpacesMonolithic();
238 else if (!type.compare("Teko") || type == "FaCSI" || type == "FaCSI-Teko" || type == "Diagonal")
239 initVectorSpacesBlock();
240 else
241 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Unkown preconditioner/solver type.");
242 }
243
244 template <class SC, class LO, class GO, class NO>
245 void NonLinearProblem<SC, LO, GO, NO>::initVectorSpacesMonolithic()
246 {
247
248 BlockMapPtr_Type map = Teuchos::rcp_const_cast<BlockMap_Type>(this->solution_->getMap());
249
250 // Teuchos::RCP<const XTpetra_Type> xTpetraMap = Teuchos::rcp_dynamic_cast<const XTpetra_Type>(map->getMergedMap()->getXpetraMap()->getMap());
251
252 TpetraMapConstPtr_Type tpetraMap = map->getMergedMap()->getTpetraMap();
253
254 this->xSpace_ = Thyra::createVectorSpace<SC, LO, GO, NO>(tpetraMap);
255 this->fSpace_ = Thyra::createVectorSpace<SC, LO, GO, NO>(tpetraMap);
256
257 typedef Teuchos::ScalarTraits<SC> ST;
258 x0_ = ::Thyra::createMember(this->xSpace_);
259 V_S(x0_.ptr(), ST::zero());
260
261 this->nominalValues_.set_x(x0_);
262 }
263
264 template <class SC, class LO, class GO, class NO>
265 void NonLinearProblem<SC, LO, GO, NO>::initVectorSpacesBlock()
266 {
267
268 BlockMapPtr_Type map = Teuchos::rcp_const_cast<BlockMap_Type>(this->solution_->getMap());
269
270 TpetraMap_Type tpetra_map;
271
272 Teuchos::Array<ThyraVecSpaceConstPtr_Type> vecSpaceArray(map->size());
273 for (int i = 0; i < map->size(); i++)
274 {
275 //Teuchos::RCP<const XTpetra_Type> xTpetraMap =
276 // Teuchos::rcp_dynamic_cast<const XTpetra_Type>(map->getBlock(i)->getXpetraMap()->getMap());
277 TpetraMapConstPtr_Type tpetraMap =map->getBlock(i)->getTpetraMap();
278 ThyraVecSpaceConstPtr_Type vecSpace = Thyra::createVectorSpace<SC, LO, GO, NO>(tpetraMap);
279 vecSpaceArray[i] = vecSpace;
280 }
281 this->xSpace_ = Teuchos::rcp(new Thyra::DefaultProductVectorSpace<SC>(vecSpaceArray()));
282 this->fSpace_ = Teuchos::rcp(new Thyra::DefaultProductVectorSpace<SC>(vecSpaceArray()));
283
284 typedef Teuchos::ScalarTraits<SC> ST;
285 x0_ = ::Thyra::createMember(this->xSpace_);
286 V_S(x0_.ptr(), ST::zero());
287
288 this->nominalValues_.set_x(x0_);
289 }
290
291}
292
293#endif
FEDD::Problem
Definition Problem_decl.hpp:38
FEDD
Adaptive Mesh Refinement.
Definition AdaptiveMeshRefinement.cpp:5