Finite Element Domain Decomposition Library
FEDDLib
Loading...
Searching...
No Matches
NonLinearSolver_def.hpp
1#ifndef NONLINEARSOLVER_DEF_hpp
2#define NONLINEARSOLVER_DEF_hpp
3
12
13namespace FEDD {
14
15
16template<class SC,class LO,class GO,class NO>
17NonLinearSolver<SC,LO,GO,NO>::NonLinearSolver():
18type_("")
19{}
20
21
22template<class SC,class LO,class GO,class NO>
23NonLinearSolver<SC,LO,GO,NO>::NonLinearSolver(std::string type):
24type_(type)
25{}
26
27template<class SC,class LO,class GO,class NO>
28NonLinearSolver<SC,LO,GO,NO>::~NonLinearSolver(){
29
30}
31
32template<class SC,class LO,class GO,class NO>
33void NonLinearSolver<SC,LO,GO,NO>::solve(NonLinearProblem_Type &problem,vec_dbl_ptr_Type valuesForExport){
34
35 if (!type_.compare("FixedPoint")) {
36 solveFixedPoint(problem,valuesForExport);
37 }
38 else if(!type_.compare("Newton")){
39 solveNewton(problem,valuesForExport);
40 }
41 else if(!type_.compare("FixedPointNewton")){
42 solveFixedPointNewton(problem, valuesForExport);
43 }
44 else if(!type_.compare("NOX")){
45#ifdef FEDD_HAVE_NOX
46 solveNOX(problem,valuesForExport);
47#endif
48 }
49
50}
51
52template<class SC,class LO,class GO,class NO>
53void NonLinearSolver<SC,LO,GO,NO>::solve(TimeProblem_Type &problem, double time, vec_dbl_ptr_Type valuesForExport){
54
55 if (!type_.compare("FixedPoint")) {
56 solveFixedPoint(problem,time);
57 }
58 else if(!type_.compare("Newton")){
59 solveNewton(problem,time, valuesForExport);
60 }
61 else if(!type_.compare("NOX")){
62#ifdef FEDD_HAVE_NOX
63 solveNOX(problem, valuesForExport);
64#endif
65 }
66 else if(!type_.compare("Extrapolation")){
67 solveExtrapolation(problem, time);
68 }
69}
70
71#ifdef FEDD_HAVE_NOX
72template<class SC,class LO,class GO,class NO>
73void NonLinearSolver<SC,LO,GO,NO>::solveNOX(NonLinearProblem_Type &problem,vec_dbl_ptr_Type valuesForExport){
74
75 bool verbose = problem.getVerbose();
76 Teuchos::RCP<NonLinearProblem<SC,LO,GO,NO> > problemPtr = Teuchos::rcpFromRef(problem);
77 Teuchos::RCP<Teuchos::ParameterList> p = sublist(problemPtr->getParameterList(),"ThyraSolver");
78// p->set("Preconditioner Type", "None"); // CH 16.04.19: preconditioner will be built seperately
79// sublist( sublist(p, "Linear Solver Types") , "Belos")->set("Left Preconditioner If Unspecified",true);
80 problemPtr->getLinearSolverBuilder()->setParameterList(p);
81
82 Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<SC> > lowsFactory = problemPtr->getLinearSolverBuilder()->createLinearSolveStrategy("");
83
84 //problemPtr->set_W_factory(lowsFactory);
85
86 // Create the initial guess
87 Teuchos::RCP<Thyra::VectorBase<SC> > initial_guess = problemPtr->getNominalValues().get_x()->clone_v();
88 Thyra::V_S(initial_guess.ptr(),Teuchos::ScalarTraits<SC>::zero());
89
90 Teuchos::RCP<Thyra::LinearOpBase<SC> > W_op = problemPtr->create_W_op();
91 Teuchos::RCP<Thyra::PreconditionerBase<SC> > W_prec = problemPtr->create_W_prec();
92
93 // cout << " NonLinearSolver<SC,LO,GO,NO>::solveNOX " << endl;
94
95 // problemPtr->create_W_op();
96
97 Teuchos::RCP<NOX::Thyra::Group> nox_group(new NOX::Thyra::Group(initial_guess,
98 problemPtr.getConst(),
99 W_op,
100 lowsFactory.getConst(),
101 W_prec,
102 Teuchos::null,
103 Teuchos::null,
104 Teuchos::null));
105
106 nox_group->computeF();
107
108 // Create the NOX status tests and the solver
109 // Create the convergence tests
110
111 Teuchos::RCP<NOX::StatusTest::NormUpdate> updateTol =
112 Teuchos::rcp(new NOX::StatusTest::NormUpdate( problemPtr->getParameterList()->sublist("Parameter").get("updateTol",1.0e-6) ) );
113
114 Teuchos::RCP<NOX::StatusTest::RelativeNormF> relresid =
115 Teuchos::rcp(new NOX::StatusTest::RelativeNormF( problemPtr->getParameterList()->sublist("Parameter").get("relNonLinTol",1.0e-4) ) );
116
117 Teuchos::RCP<NOX::StatusTest::NormWRMS> wrms =
118 Teuchos::rcp(new NOX::StatusTest::NormWRMS(problemPtr->getParameterList()->sublist("Parameter").get("relNonLinTol",1.0e-4), problemPtr->getParameterList()->sublist("Parameter").get("absNonLinTol",1.0e-6)));
119
120 Teuchos::RCP<NOX::StatusTest::NormF> absRes =
121 Teuchos::rcp(new NOX::StatusTest::NormF( problemPtr->getParameterList()->sublist("Parameter").get("absNonLinTol",1.0e-6) ) );
122
123 Teuchos::RCP<NOX::StatusTest::Combo> converged;
124
125 if ( !problemPtr->getParameterList()->sublist("Parameter").get("Combo","AND").compare("AND") )
126 converged = Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::AND));
127 else if (!problemPtr->getParameterList()->sublist("Parameter").get("Combo","AND").compare("OR") )
128 converged = Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR));
129
130 if ( problemPtr->getParameterList()->sublist("Parameter").get("Use rel tol",true) )
131 converged->addStatusTest(relresid);
132 if ( problemPtr->getParameterList()->sublist("Parameter").get("Use update tol",false) )
133 converged->addStatusTest(updateTol);
134 if (problemPtr->getParameterList()->sublist("Parameter").get("Use WRMS",false))
135 converged->addStatusTest(wrms);
136 if (problemPtr->getParameterList()->sublist("Parameter").get("Use abs tol",false))
137 converged->addStatusTest(absRes);
138
139 Teuchos::RCP<NOX::StatusTest::MaxIters> maxiters =
140 Teuchos::rcp(new NOX::StatusTest::MaxIters(problemPtr->getParameterList()->sublist("Parameter").get("MaxNonLinIts",10)));
141 Teuchos::RCP<NOX::StatusTest::FiniteValue> fv =
142 Teuchos::rcp(new NOX::StatusTest::FiniteValue);
143 Teuchos::RCP<NOX::StatusTest::Combo> combo =
144 Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR));
145 combo->addStatusTest(fv);
146 combo->addStatusTest(converged);
147 combo->addStatusTest(maxiters);
148
149 // Create nox parameter list
150 Teuchos::RCP<Teuchos::ParameterList> nl_params = sublist(problemPtr->getParameterList(),"NOXSolver");
151
152 // Create the solver
153 Teuchos::RCP<NOX::Solver::Generic> solver = NOX::Solver::buildSolver(nox_group, combo, nl_params);
154 NOX::StatusTest::StatusType solveStatus = solver->solve();
155 double nonLinearIts = solver->getSolverStatistics()->linearSolve.allNonlinearSolves_NumLinearSolves;
156 double linearIts = solver->getSolverStatistics()->linearSolve.allNonlinearSolves_NumLinearIterations;
157
158 linearIts/=nonLinearIts;
159 if (verbose){
160 std::cout << "############################################################" << std::endl;
161 std::cout << "### Total nonlinear iterations : " << nonLinearIts << " with an average of " << linearIts << " linear iterations ###" << std::endl;
162 std::cout << "############################################################" << std::endl;
163 }
164
165 if ( problemPtr->getParameterList()->sublist("Parameter").get("Cancel MaxNonLinIts",false) ) {
166 TEUCHOS_TEST_FOR_EXCEPTION((int)nonLinearIts == problemPtr->getParameterList()->sublist("Parameter").get("MaxNonLinIts",10) ,std::runtime_error,"Maximum nonlinear Iterations reached. Problem might have converged in the last step. Still we cancel here.");
167 }
168 nonLinearIts_ = nonLinearIts;
169
170}
171
172template<class SC,class LO,class GO,class NO>
173void NonLinearSolver<SC,LO,GO,NO>::solveNOX(TimeProblem_Type &problem, vec_dbl_ptr_Type valuesForExport){
174
175 bool verbose = problem.getVerbose();
176 Teuchos::RCP<TimeProblem_Type> problemPtr = Teuchos::rcpFromRef(problem);
177 Teuchos::RCP<Teuchos::ParameterList> p = sublist(problemPtr->getParameterList(),"ThyraSolver");
178
179 problemPtr->getLinearSolverBuilder()->setParameterList(p);
180
181 Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<SC> >
182 lowsFactory = problemPtr->getLinearSolverBuilder()->createLinearSolveStrategy("");
183
184 TEUCHOS_TEST_FOR_EXCEPTION(problemPtr->getSolution()->getNumVectors()>1, std::runtime_error, "With the current implementation NOX can only be used with 1 MultiVector column.");
185 // Create the initial guess and fill with last solution
186 Teuchos::RCP<Thyra::VectorBase<SC> > initialGuess = problemPtr->getNominalValues().get_x()->clone_v();
187 // Try to convert to a ProductVB. If resulting pointer is not null we need to use the ProductMV below, otherwise it is a monolithic vector.
188 Teuchos::RCP<Thyra::ProductVectorBase<SC> > initialGuessProd = Teuchos::rcp_dynamic_cast<Thyra::ProductVectorBase<SC> >(initialGuess);
189 Teuchos::RCP<Thyra::MultiVectorBase<SC> > solMV;
190 if (!initialGuessProd.is_null())
191 solMV = problemPtr->getSolution()->getProdThyraMultiVector();
192 else
193 solMV = problemPtr->getSolution()->getThyraMultiVector();
194
195 Thyra::assign(initialGuess.ptr(), *solMV->col(0));
196 //Thyra::V_S(initialGuess.ptr(),Teuchos::ScalarTraits<SC>::zero());
197 Teuchos::RCP<Thyra::LinearOpBase<SC> > W_op = problemPtr->create_W_op();
198 Teuchos::RCP<Thyra::PreconditionerBase<SC> > W_prec = problemPtr->create_W_prec();
199 // problemPtr->create_W_op();
200
201 Teuchos::RCP<NOX::Thyra::Group> nox_group(new NOX::Thyra::Group(initialGuess,
202 problemPtr.getConst(),
203 W_op,
204 lowsFactory.getConst(),
205 W_prec,
206 Teuchos::null,
207 Teuchos::null,
208 Teuchos::null));
209
210 nox_group->computeF();
211
212 // Create the NOX status tests and the solver
213 // Create the convergence tests
214 Teuchos::RCP<NOX::StatusTest::NormUpdate> updateTol =
215 Teuchos::rcp(new NOX::StatusTest::NormUpdate( problemPtr->getParameterList()->sublist("Parameter").get("updateTol",1.0e-6) ) );
216
217 Teuchos::RCP<NOX::StatusTest::RelativeNormF> relresid =
218 Teuchos::rcp(new NOX::StatusTest::RelativeNormF( problemPtr->getParameterList()->sublist("Parameter").get("relNonLinTol",1.0e-4) ) );
219
220 Teuchos::RCP<NOX::StatusTest::NormWRMS> wrms =
221 Teuchos::rcp(new NOX::StatusTest::NormWRMS(problemPtr->getParameterList()->sublist("Parameter").get("relNonLinTol",1.0e-4), problemPtr->getParameterList()->sublist("Parameter").get("absNonLinTol",1.0e-6)));
222
223 Teuchos::RCP<NOX::StatusTest::NormF> absRes =
224 Teuchos::rcp(new NOX::StatusTest::NormF( problemPtr->getParameterList()->sublist("Parameter").get("absNonLinTol",1.0e-6) ) );
225
226 Teuchos::RCP<NOX::StatusTest::Combo> converged;
227
228 if ( !problemPtr->getParameterList()->sublist("Parameter").get("Combo","AND").compare("AND") )
229 converged = Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::AND));
230 else if (!problemPtr->getParameterList()->sublist("Parameter").get("Combo","AND").compare("OR") )
231 converged = Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR));
232
233 if ( problemPtr->getParameterList()->sublist("Parameter").get("Use rel tol",true) )
234 converged->addStatusTest(relresid);
235 if ( problemPtr->getParameterList()->sublist("Parameter").get("Use update tol",false) )
236 converged->addStatusTest(updateTol);
237 if (problemPtr->getParameterList()->sublist("Parameter").get("Use WRMS",false))
238 converged->addStatusTest(wrms);
239 if (problemPtr->getParameterList()->sublist("Parameter").get("Use abs tol",false))
240 converged->addStatusTest(absRes);
241
242 Teuchos::RCP<NOX::StatusTest::MaxIters> maxiters =
243 Teuchos::rcp(new NOX::StatusTest::MaxIters(problemPtr->getParameterList()->sublist("Parameter").get("MaxNonLinIts",10)));
244 Teuchos::RCP<NOX::StatusTest::FiniteValue> fv =
245 Teuchos::rcp(new NOX::StatusTest::FiniteValue);
246 Teuchos::RCP<NOX::StatusTest::Combo> combo =
247 Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR));
248 combo->addStatusTest(fv);
249 combo->addStatusTest(converged);
250 combo->addStatusTest(maxiters);
251
252 // Create nox parameter list
253 Teuchos::RCP<Teuchos::ParameterList> nl_params = sublist(problemPtr->getParameterList(),"NOXSolver");
254 // Create the solver
255 Teuchos::RCP<NOX::Solver::Generic> solver =
256 NOX::Solver::buildSolver(nox_group, combo, nl_params);
257 NOX::StatusTest::StatusType solveStatus = solver->solve();
258
259 double nonLinearIts = solver->getSolverStatistics()->linearSolve.allNonlinearSolves_NumLinearSolves;
260 double linearIts = solver->getSolverStatistics()->linearSolve.allNonlinearSolves_NumLinearIterations;
261
262 linearIts/=nonLinearIts;
263 if (verbose){
264 std::cout << "############################################################" << std::endl;
265 std::cout << "### Total nonlinear iterations : " << nonLinearIts << " with an average of " << linearIts << " linear iterations ###" << std::endl;
266 std::cout << "############################################################" << std::endl;
267 }
268
269 if ( problemPtr->getParameterList()->sublist("Parameter").get("Cancel MaxNonLinIts",false) ) {
270 TEUCHOS_TEST_FOR_EXCEPTION((int)nonLinearIts == problemPtr->getParameterList()->sublist("Parameter").get("MaxNonLinIts",10) ,std::runtime_error,"Maximum nonlinear Iterations reached. Problem might have converged in the last step. Still we cancel here.");
271 }
272
273 if (!valuesForExport.is_null()) {
274 if (valuesForExport->size() == 2){
275 (*valuesForExport)[0] = linearIts;
276 (*valuesForExport)[1] = nonLinearIts;
277 }
278 }
279}
280#endif
281
282template<class SC,class LO,class GO,class NO>
283void NonLinearSolver<SC,LO,GO,NO>::solveFixedPoint(NonLinearProblem_Type &problem,vec_dbl_ptr_Type valuesForExport){
284
285 bool verbose = problem.getVerbose();
286 TEUCHOS_TEST_FOR_EXCEPTION(problem.getRhs()->getNumVectors()!=1,std::logic_error,"We need to change the code for numVectors>1.");
287 // -------
288 // fix point iteration
289 // -------
290 double gmresIts = 0.;
291 double residual0 = 1.;
292 double residual = 1.;
293
294 double tol = problem.getParameterList()->sublist("Parameter").get("relNonLinTol",1.0e-6);
295 int maxNonLinIts = problem.getParameterList()->sublist("Parameter").get("MaxNonLinIts",10);
296 int nlIts=0;
297
298 double criterionValue = 1.;
299 std::string criterion = problem.getParameterList()->sublist("Parameter").get("Criterion","Residual");
300
301
302 while ( nlIts < maxNonLinIts ) {
303
304 // This will compute the residual vector
305 problem.calculateNonLinResidualVec("reverse");
306
307 // Analogous to Newton but we assemble here only the parts which are needed for Fixed Point Iteration
308 problem.assemble("FixedPoint");
309
310 problem.setBoundariesSystem();
311
312 if (criterion=="Residual")
313 residual = problem.calculateResidualNorm();
314
315 if (nlIts==0)
316 residual0 = residual;
317
318 if (criterion=="Residual"){
319 criterionValue = residual/residual0;
320 if (verbose)
321 std::cout << "### Fixed Point iteration : " << nlIts << " relative nonlinear residual : " << criterionValue << std::endl;
322 if ( criterionValue < tol )
323 break;
324 }
325
326 // Set the current Nonlinear Step in the Problem class
327 problem.setNonlinearIterationStep(nlIts);
328 gmresIts += problem.solveAndUpdate( criterion, criterionValue );
329 nlIts++;
330 if(criterion=="Update"){
331 if (verbose)
332 std::cout << "### Fixed Point iteration : " << nlIts << " residual of update : " << criterionValue << std::endl;
333 if ( criterionValue < tol )
334 break;
335 }
336 // ####### end FPI #######
337 }
338
339 gmresIts/=nlIts;
340 if (verbose)
341 std::cout << "### Total FPI : " << nlIts << " with average gmres its : " << gmresIts << std::endl;
342 if ( problem.getParameterList()->sublist("Parameter").get("Cancel MaxNonLinIts",false) ) {
343 TEUCHOS_TEST_FOR_EXCEPTION( nlIts == maxNonLinIts ,std::runtime_error,"Maximum nonlinear Iterations reached. Problem might have converged in the last step. Still we cancel here.");
344 }
345
346}
347
348template<class SC,class LO,class GO,class NO>
349void NonLinearSolver<SC,LO,GO,NO>::solveNewton( NonLinearProblem_Type &problem, vec_dbl_ptr_Type valuesForExport ){
350
351 bool verbose = problem.getVerbose();
352
353 TEUCHOS_TEST_FOR_EXCEPTION(problem.getRhs()->getNumVectors()!=1,std::logic_error,"We need to change the code for numVectors>1.")
354 // -------
355 // Newton
356 // -------
357 double gmresIts = 0.;
358 double residual0 = 1.;
359 double residual = 1.;
360 double tol = problem.getParameterList()->sublist("Parameter").get("relNonLinTol",1.0e-6);
361 int nlIts=0;
362 int maxNonLinIts = problem.getParameterList()->sublist("Parameter").get("MaxNonLinIts",10);
363 double criterionValue = 1.;
364 std::string criterion = problem.getParameterList()->sublist("Parameter").get("Criterion","Residual");
365
366 while ( nlIts < maxNonLinIts ) {
367 //this makes only sense for Navier-Stokes/Stokes, for other problems, e.g., non linear elasticity, it should do nothing.
368
369 problem.calculateNonLinResidualVec("reverse");
370
371 if (criterion=="Residual")
372 residual = problem.calculateResidualNorm();
373
374 problem.assemble("Newton");
375
376 problem.setBoundariesSystem();
377
378 if (nlIts==0)
379 residual0 = residual;
380
381 if (criterion=="Residual"){
382 criterionValue = residual/residual0;
383 if (verbose)
384 std::cout << "### Newton iteration : " << nlIts << " relative nonlinear residual : " << criterionValue << std::endl;
385 if ( criterionValue < tol )
386 break;
387 }
388 // PRINT INFOS
389 // Set the current Nonlinear Step in the Problem class
390 problem.setNonlinearIterationStep(nlIts);
391 gmresIts += problem.solveAndUpdate( criterion, criterionValue );
392 nlIts++;
393 if(criterion=="Update"){
394 if (verbose)
395 std::cout << "### Newton iteration : " << nlIts << " residual of update : " << criterionValue << std::endl;
396 if ( criterionValue < tol )
397 break;
398 }
399
400 // ####### end FPI #######
401 }
402
403 gmresIts/=nlIts;
404 if (verbose)
405 std::cout << "### Total Newton iterations : " << nlIts << " with average gmres its : " << gmresIts << std::endl;
406 if ( problem.getParameterList()->sublist("Parameter").get("Cancel MaxNonLinIts",false) ) {
407 TEUCHOS_TEST_FOR_EXCEPTION(nlIts == maxNonLinIts ,std::runtime_error,"Maximum nonlinear Iterations reached. Problem might have converged in the last step. Still we cancel here.");
408 }
409
410 if (!valuesForExport.is_null()) {
411 if (valuesForExport->size() == 2){
412 (*valuesForExport)[0] = gmresIts;
413 (*valuesForExport)[1] = nlIts;
414 }
415 }
416}
417
418
419/*
420@TODO: In future user would like to switch between linearization based on the fact that the criteriumValue is not
421 decreasing "enough" or that the it is increasing -> For this one would need to store the previous criterionValues in a vector
422*/
423template<class SC,class LO,class GO,class NO>
424void NonLinearSolver<SC,LO,GO,NO>::solveFixedPointNewton( NonLinearProblem_Type &problem, vec_dbl_ptr_Type valuesForExport ){
425
426 bool verbose = problem.getVerbose();
427
428 TEUCHOS_TEST_FOR_EXCEPTION(problem.getRhs()->getNumVectors()!=1,std::logic_error,"We need to change the code for numVectors>1.")
429 // -------
430 // FixedPoint-Newton
431 // -------
432 double gmresIts = 0.;
433 double residual0 = 1.;
434 double residual = 1.;
435 double tol = problem.getParameterList()->sublist("Parameter").get("relNonLinTol",1.0e-6);
436 int nlIts=0;
437 int maxNonLinIts = problem.getParameterList()->sublist("Parameter").get("MaxNonLinIts",10);
438 double criterionValue = 1.;
439 std::string criterion = problem.getParameterList()->sublist("Parameter").get("Criterion","Residual");
440
441 std::string linearization = problem.getParameterList()->sublist("General").get("Initial_Linearization", "FixedPoint"); // Default is FixedPoint for first iteration
442 Teuchos::RCP<NonLinearProblem<SC,LO,GO,NO> > problemPtr = Teuchos::rcpFromRef(problem);
443 Teuchos::RCP<Teuchos::ParameterList> parameterListGeneral = sublist(problemPtr->getParameterList(),"General");
444 bool linearizationChanged = false;
445
446
447 if (verbose)
448 std::cout << "####### Run FixedPoint-Newton Combination with Initial Linearization: " << linearization << " ####### " << std::endl;
449
450 while ( nlIts < maxNonLinIts ) { // User still can decinde on a maximum number of nonlinear iterations
451
452 // Compute residual vector
453 problem.calculateNonLinResidualVec("reverse"); // Based on current solution (which is in first iteration the initial guess) compute the nonlinear residual vector
454
455 if (criterion=="Residual")
456 residual = problem.calculateResidualNorm(); // Compute norm of nonlinear residual vector
457
458
459 if (nlIts==0)
460 residual0 = residual; // Store initial residual for relative convergence check
461
462 if (criterion=="Residual"){ // Check convergence based on relative nonlinear residual
463 criterionValue = residual/residual0;
464 if (verbose)
465 std::cout << "### " << linearization << " iteration : " << nlIts << " relative nonlinear residual : " << criterionValue << std::endl;
466 // One could add here a check for if nlIts=0 then initial linearization is shown which is not used for assemble?
467 if ( criterionValue < tol )
468 break;
469 }
470
471 // Use the user-defined switching strategy to decide whether to switch linearization and which one to use
472 linearizationChanged = this->switchingStrategy_(linearization, nlIts, criterionValue, parameterListGeneral); // Get the next linearization type
473 if (linearizationChanged)
474 {
475 if (verbose)
476 std::cout << "### Switching linearization to " << linearization << " in iteration " << nlIts << std::endl;
477 problem.changeAssFELinearization(linearization); // Change linearization in all elements based on chosen strategy (Picard/ Newton) - But only if different from previous one
478 }
479
480 problem.assemble(linearization);
481
482
483 problem.setBoundariesSystem();
484
485
486 // PRINT INFOS
487 // Set the current Nonlinear Step in the Problem class
488 problem.setNonlinearIterationStep(nlIts);
489 gmresIts += problem.solveAndUpdate( criterion, criterionValue );
490 nlIts++;
491 if(criterion=="Update"){
492 if (verbose)
493 std::cout << "### " << linearization << " iteration : " << nlIts << " residual of update : " << criterionValue << std::endl;
494 if ( criterionValue < tol )
495 break;
496 }
497
498 // ####### end FPI #######
499 }
500
501 gmresIts/=nlIts;
502 if (verbose)
503 std::cout << "### Total " << "nonlinear" << " iterations : " << nlIts << " with average gmres its : " << gmresIts << std::endl;
504 if ( problem.getParameterList()->sublist("Parameter").get("Cancel MaxNonLinIts",false) ) {
505 TEUCHOS_TEST_FOR_EXCEPTION(nlIts == maxNonLinIts ,std::runtime_error,"Maximum nonlinear Iterations reached. Problem might have converged in the last step. Still we cancel here.");
506 }
507
508 if (!valuesForExport.is_null()) {
509 if (valuesForExport->size() == 2){
510 (*valuesForExport)[0] = gmresIts;
511 (*valuesForExport)[1] = nlIts;
512 }
513 }
514}
515
516
517
518
519template<class SC,class LO,class GO,class NO>
520void NonLinearSolver<SC,LO,GO,NO>::solveFixedPoint(TimeProblem_Type &problem, double time){
521
522 bool verbose = problem.getVerbose();
523 problem.setBoundariesRHS(time);
524 TEUCHOS_TEST_FOR_EXCEPTION(problem.getRhs()->getNumVectors()!=1,std::logic_error,"We need to change the code for numVectors>1.")
525
526 // -------
527 // fix point iteration
528 // -------
529 double gmresIts = 0.;
530 double residual0 = 1.;
531 double residual = 1.;
532 double tol = problem.getParameterList()->sublist("Parameter").get("relNonLinTol",1.0e-6);
533 int nlIts=0;
534 int maxNonLinIts = problem.getParameterList()->sublist("Parameter").get("MaxNonLinIts",10);
535 double criterionValue = 1.;
536 std::string criterion = problem.getParameterList()->sublist("Parameter").get("Criterion","Residual");
537
538 while ( nlIts < maxNonLinIts ) {
539
540 problem.calculateNonLinResidualVec("reverse", time);
541
542 if (criterion=="Residual")
543 residual = problem.calculateResidualNorm();
544
545 if (nlIts==0)
546 residual0 = residual;
547
548 // Linearization of system matrix is done in calculateNonLinResidualVec
549 // Now we need to combine it with the mass matrix
550 problem.combineSystems();
551
552 problem.setBoundariesSystem();
553
554 if (criterion=="Residual"){
555 criterionValue = residual/residual0;
556 if (verbose)
557 std::cout << "### Fixed Point iteration : " << nlIts << " relative nonlinear residual : " << criterionValue << std::endl;
558 if ( criterionValue < tol )
559 break;
560 }
561
562 gmresIts += problem.solveAndUpdate( criterion, criterionValue );
563
564 nlIts++;
565 if(criterion=="Update"){
566 if (verbose)
567 std::cout << "### Fixed Point iteration : " << nlIts << " residual of update : " << criterionValue << std::endl;
568 if ( criterionValue < tol )
569 break;
570 }
571 // ####### end FPI #######
572 }
573
574 gmresIts/=nlIts;
575 if (verbose)
576 std::cout << "### Total FPI : " << nlIts << " with average gmres its : " << gmresIts << std::endl;
577 if ( problem.getParameterList()->sublist("Parameter").get("Cancel MaxNonLinIts",false) ) {
578 TEUCHOS_TEST_FOR_EXCEPTION( nlIts == maxNonLinIts ,std::runtime_error,"Maximum nonlinear Iterations reached. Problem might have converged in the last step. Still we cancel here.");
579 }
580}
581
582
583
584template<class SC,class LO,class GO,class NO>
585void NonLinearSolver<SC,LO,GO,NO>::solveNewton(TimeProblem_Type &problem, double time, vec_dbl_ptr_Type valuesForExport ){
586
587 bool verbose = problem.getVerbose();
588 problem.setBoundariesRHS(time);
589
590
591 TEUCHOS_TEST_FOR_EXCEPTION(problem.getRhs()->getNumVectors()!=1,std::logic_error,"We need to change the code for numVectors>1.")
592
593 // -------
594 // Newton iteration
595 // -------
596 double gmresIts = 0.;
597 double residual0 = 1.;
598 double residual = 1.;
599 double tol = problem.getParameterList()->sublist("Parameter").get("relNonLinTol",1.0e-6);
600 int nlIts=0;
601 int maxNonLinIts = problem.getParameterList()->sublist("Parameter").get("MaxNonLinIts",10);
602 double criterionValue = 1.;
603 std::string criterion = problem.getParameterList()->sublist("Parameter").get("Criterion","Residual");
604 std::string timestepping = problem.getParameterList()->sublist("Timestepping Parameter").get("Class","Singlestep");
605
606 while ( nlIts < maxNonLinIts ) {
607 if (timestepping == "External")
608 problem.calculateNonLinResidualVec("external", time);
609 else
610 problem.calculateNonLinResidualVec("reverse", time);
611 if (criterion=="Residual")
612 residual = problem.calculateResidualNorm();
613
614 if (nlIts==0)
615 residual0 = residual;
616
617 if (criterion=="Residual"){
618 criterionValue = residual/residual0;
619// exporterTxt->exportData( criterionValue );
620 if (verbose)
621 std::cout << "### Newton iteration : " << nlIts << " relative nonlinear residual : " << criterionValue << std::endl;
622 if ( criterionValue < tol )
623 break;
624 }
625
626 // Systems are combined in timeProblem.assemble("Newton") and then combined
627 problem.assemble("Newton");
628
629 problem.setBoundariesSystem();
630
631
632 if (timestepping == "External"){//AceGen
633 gmresIts += problem.solveAndUpdate( "ResidualAceGen", criterionValue );
634 // exporterTxt->exportData( criterionValue );
635
636 //problem.assembleExternal( "OnlyUpdate" );// update AceGEN internal variables
637 }
638 else
639 gmresIts += problem.solveAndUpdate( criterion, criterionValue );
640
641 nlIts++;
642
643 //problem.getSolution()->getBlock(0)->print();
644 if(criterion=="Update"){
645 if (verbose)
646 std::cout << "### Newton iteration : " << nlIts << " residual of update : " << criterionValue << std::endl;
647 if ( criterionValue < tol )
648 break;
649 }
650
651 // ####### end FPI #######
652 }
653
654 gmresIts/=nlIts;
655 if (verbose)
656 std::cout << "### Total Newton iteration : " << nlIts << " with average gmres its : " << gmresIts << std::endl;
657 if ( problem.getParameterList()->sublist("Parameter").get("Cancel MaxNonLinIts",false) ) {
658 TEUCHOS_TEST_FOR_EXCEPTION(nlIts == maxNonLinIts ,std::runtime_error,"Maximum nonlinear Iterations reached. Problem might have converged in the last step. Still we cancel here.");
659 }
660 if (!valuesForExport.is_null()) {
661 if (valuesForExport->size() == 2){
662 (*valuesForExport)[0] = gmresIts;
663 (*valuesForExport)[1] = nlIts;
664 }
665
666
667 }
668
669}
670
671
672template<class SC,class LO,class GO,class NO>
673void NonLinearSolver<SC,LO,GO,NO>::solveExtrapolation(TimeProblem<SC,LO,GO,NO> &problem, double time){
674
675 bool verbose = problem.getVerbose();
676
677 problem.assemble("Extrapolation");
678
679 problem.setBoundaries(time); // Setting boundaries to system rhs. The rest of the rhs (e.g. M*u_t) must/should be implemented in DAESolver
680
681 int gmresIts = problem.solve( );
682
683 if (verbose) {
684 std::cout << "### GMRES Its : " << gmresIts << std::endl;
685 }
686}
687}
688#endif
FEDD::NonLinearSolver
Definition NonLinearSolver_decl.hpp:36
FEDD::NonLinearSolver::solve
void solve(NonLinearProblem_Type &problem, vec_dbl_ptr_Type valuesForExport=Teuchos::null)
Call for solving a nonlinear problem, depending on linearization solveNOX/solveFixedPoint/solveNewton...
Definition NonLinearSolver_def.hpp:33
FEDD
Adaptive Mesh Refinement.
Definition AdaptiveMeshRefinement_decl.hpp:36