Hermes2D  3.0
nox_solver.cpp
Go to the documentation of this file.
1 // This file is part of HermesCommon
2 //
3 // Hermes2D is free software: you can redistribute it and/or modify
4 // it under the terms of the GNU General Public License as published by
5 // the Free Software Foundation, either version 2 of the License, or
6 // (at your option) any later version.
7 //
8 // Hermes2D is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 // GNU General Public License for more details.
12 //
13 // You should have received a copy of the GNU General Public License
14 // along with Hermes2D. If not, see <http://www.gnu.org/licenses/>.
18 #include "solver/nox_solver.h"
19 #include "projections/ogprojection.h"
20 
21 #if(defined HAVE_NOX && defined HAVE_EPETRA && defined HAVE_TEUCHOS)
22 
23 namespace Hermes
24 {
25  namespace Hermes2D
26  {
27  static Epetra_SerialComm seq_comm;
28 
29  template<typename Scalar>
30  DiscreteProblemNOX<Scalar>::DiscreteProblemNOX() : DiscreteProblem<Scalar>(), jacobian(nullptr)
31  {
32  this->precond = Teuchos::null;
33  }
34 
35  template<typename Scalar>
36  DiscreteProblemNOX<Scalar>::DiscreteProblemNOX(WeakFormSharedPtr<Scalar> wf, std::vector<SpaceSharedPtr<Scalar> > spaces) : DiscreteProblem<Scalar>(wf, spaces), jacobian(new EpetraMatrix<Scalar>())
37  {
38  this->precond = Teuchos::null;
39  this->jacobian->prealloc(Space<Scalar>::get_num_dofs(spaces));
40  }
41 
42  template<typename Scalar>
43  DiscreteProblemNOX<Scalar>::DiscreteProblemNOX(WeakFormSharedPtr<Scalar> wf, SpaceSharedPtr<Scalar> space) : DiscreteProblem<Scalar>(wf, space), jacobian(new EpetraMatrix<Scalar>())
44  {
45  this->precond = Teuchos::null;
46  this->jacobian->prealloc(Space<Scalar>::get_num_dofs(space));
47  }
48 
49  template<typename Scalar>
50  DiscreteProblemNOX<Scalar>::~DiscreteProblemNOX()
51  {
52  delete this->jacobian;
53  }
54 
55  template<typename Scalar>
56  bool DiscreteProblemNOX<Scalar>::computeF(const Epetra_Vector &x, Epetra_Vector &f, FillType flag)
57  {
58  // wrap our structures around core Epetra objects
59  EpetraVector<Scalar> xx(x);
60  EpetraVector<Scalar> rhs(f);
61 
62  rhs.zero();
63 
64  Scalar* coeff_vec = malloc_with_check<Scalar>(xx.get_size());
65  xx.extract(coeff_vec);
66  // nullptr is for the global matrix.
67  this->assemble(coeff_vec, nullptr, &rhs);
68  free_with_check(coeff_vec);
69 
70  return true;
71  }
72 
73  template<typename Scalar>
74  bool DiscreteProblemNOX<Scalar>::computeJacobian(const Epetra_Vector &x, Epetra_Operator &op)
75  {
76  Epetra_RowMatrix *jac = dynamic_cast<Epetra_RowMatrix *>(&op);
77  assert(jac != nullptr);
78 
79  // wrap our structures around core Epetra objects
80  EpetraVector<Scalar> xx(x);
81  EpetraMatrix<Scalar> jacob(*jac);
82 
83  jacob.zero();
84 
85  Scalar* coeff_vec = malloc_with_check<Scalar>(xx.get_size());
86  xx.extract(coeff_vec);
87  // nullptr is for the right-hand side.
88  this->assemble(coeff_vec, &jacob, nullptr);
89  free_with_check(coeff_vec);
90  //jacob.finish();
91 
92  return true;
93  }
94 
95  template<typename Scalar>
96  bool DiscreteProblemNOX<Scalar>::computePreconditioner(const Epetra_Vector &x, Epetra_Operator &m,
97  Teuchos::ParameterList *precParams)
98  {
99  // wrap our structures around core Epetra objects
100  EpetraVector<Scalar> xx(x);
101 
102  Scalar* coeff_vec = malloc_with_check<Scalar>(xx.get_size());
103  xx.extract(coeff_vec);
104  this->assemble(coeff_vec, jacobian);
105  free_with_check(coeff_vec);
106 
107  precond->create(jacobian);
108  precond->compute();
109  m = *precond->get_obj();
110 
111  return true;
112  }
113 
114  template<typename Scalar>
115  Teuchos::RCP<Hermes::Preconditioners::EpetraPrecond<Scalar> > DiscreteProblemNOX<Scalar>::get_precond()
116  {
117  return precond;
118  }
119 
120  template<typename Scalar>
121  EpetraMatrix<Scalar> *DiscreteProblemNOX<Scalar>::get_jacobian()
122  {
123  return jacobian;
124  }
125 
126  template<typename Scalar>
127  Scalar* NewtonSolverNOX<Scalar>::get_sln_vector()
128  {
129  return this->sln_vector;
130  }
131 
132  template<typename Scalar>
133  NewtonSolverNOX<Scalar>::NewtonSolverNOX(DiscreteProblemNOX<Scalar>* problem) : dp(problem), sln_vector(nullptr), own_dp(false)
134  {
135  init();
136  }
137 
138  template<typename Scalar>
139  NewtonSolverNOX<Scalar>::NewtonSolverNOX(WeakFormSharedPtr<Scalar> wf, std::vector<SpaceSharedPtr<Scalar> > spaces) : dp(new DiscreteProblemNOX<Scalar>(wf, spaces)), sln_vector(nullptr), own_dp(true)
140  {
141  init();
142  }
143 
144  template<typename Scalar>
145  NewtonSolverNOX<Scalar>::NewtonSolverNOX(WeakFormSharedPtr<Scalar> wf, SpaceSharedPtr<Scalar> space) : dp(new DiscreteProblemNOX<Scalar>(wf, space)), sln_vector(nullptr), own_dp(true)
146  {
147  init();
148  }
149 
150  template<typename Scalar>
151  void NewtonSolverNOX<Scalar>::init()
152  {
153  // default values
154  // convergence test
155  conv.max_iters = 10;
156  conv.abs_resid = 1.0e-6;
157  conv.rel_resid = 1.0e-2;
158  conv.norm_type = NOX::Abstract::Vector::TwoNorm;
159  conv.stype = NOX::StatusTest::NormF::Scaled;
160  conv.update = 1.0e-5;
161  conv.wrms_rtol = 1.0e-2;
162  conv.wrms_atol = 1.0e-8;
163 
164  conv_flag.absresid = 1;
165  conv_flag.relresid = 0;
166  conv_flag.update = 0;
167  conv_flag.wrms = 0;
168 
169  nl_pars = Teuchos::rcp(new Teuchos::ParameterList);
170  // Set the nonlinear solver method
171  nl_pars->set("Nonlinear Solver", "Line Search Based");
172  nl_pars->sublist("Printing").set("Output Information", NOX::Utils::Error);
173 
174  // Sublist for line search
175  Teuchos::ParameterList &search_pars = nl_pars->sublist("Line Search");
176  search_pars.set("Method", "Full Step");
177 
178  // Sublist for direction
179  Teuchos::ParameterList &dir_pars = nl_pars->sublist("Direction");
180  dir_pars.set("Method", "Newton");
181  Teuchos::ParameterList &newton_pars = dir_pars.sublist("Newton");
182 
183  newton_pars.set("Forcing Term Method", "Constant");
184 
185  // Sublist for linear solver for the Newton method
186  Teuchos::ParameterList &ls_pars = newton_pars.sublist("Linear Solver");
187  ls_pars.set("Aztec Solver", "GMRES");
188  ls_pars.set("Max Iterations", 800);
189  ls_pars.set("Tolerance", 1e-8);
190  ls_pars.set("Size of Krylov Subspace", 50);
191  ls_pars.set("Preconditioner Reuse Policy", "Recompute");
192  ls_pars.set("Output Frequency", AZ_all);
193  ls_pars.set("Max Age Of Prec", 999);
194  }
195 
196  template<typename Scalar>
197  void NewtonSolverNOX<Scalar>::set_time(double time)
198  {
199  }
200 
201  template<typename Scalar>
202  void NewtonSolverNOX<Scalar>::set_time_step(double time_step)
203  {
204  }
205 
206  template<typename Scalar>
207  NewtonSolverNOX<Scalar>::~NewtonSolverNOX()
208  {
209  // FIXME: this does not destroy the "interface_", and Trilinos
210  // complains at closing main.cpp.
211  if (this->own_dp)
212  delete this->dp;
213  }
214 
215  template<typename Scalar>
216  void DiscreteProblemNOX<Scalar>::set_precond(Teuchos::RCP<Hermes::Preconditioners::EpetraPrecond<Scalar> > &pc)
217  {
218  precond = pc;
219  }
220 
221  template<typename Scalar>
222  void NewtonSolverNOX<Scalar>::set_precond(Hermes::Preconditioners::EpetraPrecond<Scalar> &pc)
223  {
224  Teuchos::RCP<Hermes::Preconditioners::EpetraPrecond<Scalar> > tpc = Teuchos::rcpFromRef(pc);
225  dp->set_precond(tpc);
226  nl_pars->sublist("Direction").sublist("Newton").sublist("Linear Solver").set("Preconditioner", "User Defined");
227  }
228 
229  template<typename Scalar>
230  void NewtonSolverNOX<Scalar>::set_precond(const char *pc)
231  {
232  nl_pars->sublist("Direction").sublist("Newton").sublist("Linear Solver").set("Preconditioner", pc);
233  }
234 
235  template<typename Scalar>
236  void NewtonSolverNOX<Scalar>::set_output_flags(int flags)
237  {
238  nl_pars->sublist("Printing").set("Output Information", flags);
239  }
240 
241  template<typename Scalar>
242  void NewtonSolverNOX<Scalar>::set_ls_type(const char *type)
243  {
244  nl_pars->sublist("Direction").sublist("Newton").sublist("Linear Solver").set("Aztec Solver", type);
245  }
246 
247  template<typename Scalar>
248  void NewtonSolverNOX<Scalar>::set_ls_max_iters(int iters)
249  {
250  nl_pars->sublist("Direction").sublist("Newton").sublist("Linear Solver").set("Max Iterations", iters);
251  }
252 
253  template<typename Scalar>
254  void NewtonSolverNOX<Scalar>::set_ls_tolerance(double tolerance)
255  {
256  nl_pars->sublist("Direction").sublist("Newton").sublist("Linear Solver").set("Tolerance", tolerance);
257  }
258 
259  template<typename Scalar>
260  void NewtonSolverNOX<Scalar>::set_ls_sizeof_krylov_subspace(int size)
261  {
262  nl_pars->sublist("Direction").sublist("Newton").sublist("Linear Solver").set("Size of Krylov Subspace", size);
263  }
264 
265  template<typename Scalar>
266  void NewtonSolverNOX<Scalar>::set_precond_reuse(const char * pc_reuse)
267  {
268  nl_pars->sublist("Direction").sublist("Newton").sublist("Linear Solver").set("Preconditioner Reuse Policy", pc_reuse);
269  }
270  template<typename Scalar>
271  void NewtonSolverNOX<Scalar>::set_precond_max_age(int max_age)
272  {
273  nl_pars->sublist("Direction").sublist("Newton").sublist("Linear Solver").set("Max Age Of Prec", max_age);
274  }
275 
276  template<typename Scalar>
277  void NewtonSolverNOX<Scalar>::solve(std::vector<MeshFunctionSharedPtr<Scalar> > initial_guess)
278  {
279  Scalar* coeff_vec = malloc_with_check<Scalar>(Space<Scalar>::get_num_dofs(this->dp->spaces));
280  OGProjection<Scalar>::project_global(this->dp->spaces, initial_guess, coeff_vec);
281  this->solve(coeff_vec);
282  free_with_check(coeff_vec);
283  }
284 
285  template<typename Scalar>
286  void NewtonSolverNOX<Scalar>::solve(MeshFunctionSharedPtr<Scalar> initial_guess)
287  {
288  Scalar* coeff_vec = malloc_with_check<Scalar>(Space<Scalar>::get_num_dofs(this->dp->spaces));
289  OGProjection<Scalar>::project_global(this->dp->spaces[0], initial_guess, coeff_vec);
290  this->solve(coeff_vec);
291  free_with_check(coeff_vec);
292  }
293 
294  template<typename Scalar>
295  void NewtonSolverNOX<Scalar>::solve(Scalar* coeff_vec)
296  {
297  // Put the initial coeff_vec into the inner structure for the initial guess.
298  Hermes::Algebra::EpetraVector<Scalar> temp_init_sln;
299  int ndofs = Space<Scalar>::get_num_dofs(this->dp->get_spaces());
300  temp_init_sln.alloc(ndofs);
301  if (!coeff_vec)
302  temp_init_sln.zero();
303  else
304  {
305  for (int i = 0; i < ndofs; i++)
306  temp_init_sln.set(i, coeff_vec[i]);
307  }
308 
309  NOX::Epetra::Vector init_sln(*temp_init_sln.vec);
310 
311  // prepare variables
312  // print settings
313  Teuchos::ParameterList &print_pars = nl_pars->sublist("Printing");
314  // linear system settings
315  Teuchos::ParameterList &ls_pars = nl_pars->sublist("Direction").sublist("Newton").sublist("Linear Solver");
316  // preconditioner
317  Teuchos::RCP<Hermes::Preconditioners::EpetraPrecond<Scalar> > precond = dp->get_precond();
318  //linear system
319  Teuchos::RCP<NOX::Epetra::LinearSystemAztecOO> lin_sys;
320  // problem
321  Teuchos::RCP<NOX::Epetra::Interface::Required> i_req = Teuchos::rcpFromRef(*dp);
322  // jacobian matrix
323  Teuchos::RCP<Epetra_RowMatrix> jac_mat;
324 
325  // Create linear system
326  if (this->dp->get_weak_formulation()->is_matrix_free())
327  {
328  if (precond == Teuchos::null)
329  { //Matrix free without preconditioner
330  lin_sys = Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(print_pars, ls_pars, i_req, init_sln));
331  }
332  else
333  { //Matrix free with preconditioner
334  Teuchos::RCP<NOX::Epetra::Interface::Preconditioner> i_prec = Teuchos::rcpFromRef(*dp);
335  lin_sys = Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(print_pars, ls_pars, i_req,
336  i_prec, precond, init_sln));
337  }
338  }
339  else
340  { // not matrix Free (with jacobian)
341  Teuchos::RCP<NOX::Epetra::Interface::Jacobian> i_jac = Teuchos::rcpFromRef(*dp);
342  jac_mat = Teuchos::rcp(dp->get_jacobian()->mat);
343  if (precond == Teuchos::null)
344  { //Matrix without preconditioner
345  lin_sys = Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(print_pars, ls_pars, i_req,
346  i_jac, jac_mat, init_sln));
347  }
348  else
349  { //Matrix with preconditioner
350  Teuchos::RCP<NOX::Epetra::Interface::Preconditioner> i_prec = Teuchos::rcpFromRef(*dp);
351  lin_sys = Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(print_pars, ls_pars,
352  i_jac, jac_mat, i_prec, precond, init_sln));
353  }
354  }
355 
356  // Create the Group
357  Teuchos::RCP<NOX::Epetra::Group> grp =
358  Teuchos::rcp(new NOX::Epetra::Group(print_pars, i_req, init_sln, lin_sys));
359 
360  // Create convergence tests
361  Teuchos::RCP<NOX::StatusTest::Combo> converged =
362  Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::AND));
363 
364  if (conv_flag.absresid)
365  {
366  Teuchos::RCP<NOX::StatusTest::NormF> absresid =
367  Teuchos::rcp(new NOX::StatusTest::NormF(conv.abs_resid, conv.norm_type, conv.stype));
368  converged->addStatusTest(absresid);
369  }
370 
371  if (conv_flag.relresid)
372  {
373  Teuchos::RCP<NOX::StatusTest::NormF> relresid =
374  Teuchos::rcp(new NOX::StatusTest::NormF(*grp.get(), conv.rel_resid));
375  converged->addStatusTest(relresid);
376  }
377 
378  if (conv_flag.update)
379  {
380  Teuchos::RCP<NOX::StatusTest::NormUpdate> update =
381  Teuchos::rcp(new NOX::StatusTest::NormUpdate(conv.update));
382  converged->addStatusTest(update);
383  }
384 
385  if (conv_flag.wrms)
386  {
387  Teuchos::RCP<NOX::StatusTest::NormWRMS> wrms =
388  Teuchos::rcp(new NOX::StatusTest::NormWRMS(conv.wrms_rtol, conv.wrms_atol));
389  converged->addStatusTest(wrms);
390  }
391 
392  Teuchos::RCP<NOX::StatusTest::MaxIters> maxiters =
393  Teuchos::rcp(new NOX::StatusTest::MaxIters(conv.max_iters));
394 
395  Teuchos::RCP<NOX::StatusTest::FiniteValue> fv =
396  Teuchos::rcp(new NOX::StatusTest::FiniteValue);
397 
398  Teuchos::RCP<NOX::StatusTest::Combo> cmb =
399  Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR));
400 
401  cmb->addStatusTest(fv);
402  cmb->addStatusTest(converged);
403  cmb->addStatusTest(maxiters);
404 
405  // Output parameters
406  Teuchos::RCP<Teuchos::ParameterList> final_pars = nl_pars;
407 
408  // Create the solver
409  Teuchos::RCP<NOX::Solver::Generic> solver = NOX::Solver::buildSolver(grp, cmb, final_pars);
410 
411  // Solve.
412  NOX::StatusTest::StatusType status = solver->solve();
413 
414  if (!this->dp->get_weak_formulation()->is_matrix_free())
415  // release the ownership (we take care of jac_mat by ourselves)
416  jac_mat.release();
417 
418  if (status == NOX::StatusTest::Converged)
419  {
420  // get result informations
421  num_iters = solver->getNumIterations();
422  residual = solver->getSolutionGroup().getNormF();
423  num_lin_iters = final_pars->sublist("Direction").sublist("Newton").sublist("Linear Solver").sublist("Output").get("Total Number of Linear Iterations", -1);
424  achieved_tol = final_pars->sublist("Direction").sublist("Newton").sublist("Linear Solver").sublist("Output").get("Achieved Tolerance", 0.0);
425 
426  // Get the Epetra_Vector with the final solution from the solver
427  const NOX::Epetra::Group &f_grp =
428  dynamic_cast<const NOX::Epetra::Group &>(solver->getSolutionGroup());
429  const Epetra_Vector &f_sln =
430  (dynamic_cast<const NOX::Epetra::Vector &>(f_grp.getX())).getEpetraVector();
431  // extract solution
432  if (this->sln_vector)
433  free_with_check(this->sln_vector);
434  this->sln_vector = (Scalar*)calloc(ndofs, sizeof(Scalar));
435 
436  f_sln.ExtractCopy(this->sln_vector);
437  }
438  else // not converged
439  {
440  num_iters = -1;
441  throw Exceptions::Exception("Nox solver did not converge");
442  }
443  }
444 
445  template<typename Scalar>
446  int NewtonSolverNOX<Scalar>::get_num_iters()
447  {
448  return num_iters;
449  }
450 
451  template<typename Scalar>
452  double NewtonSolverNOX<Scalar>::get_residual()
453  {
454  return residual;
455  }
456 
457  template<typename Scalar>
458  int NewtonSolverNOX<Scalar>::get_num_lin_iters()
459  {
460  return num_lin_iters;
461  }
462 
463  template<typename Scalar>
464  double NewtonSolverNOX<Scalar>::get_achieved_tol()
465  {
466  return achieved_tol;
467  }
468 
469  template<typename Scalar>
470  void NewtonSolverNOX<Scalar>::set_norm_type(NOX::Abstract::Vector::NormType type)
471  {
472  conv.norm_type = type;
473  }
474 
475  template<typename Scalar>
476  void NewtonSolverNOX<Scalar>::set_scale_type(NOX::StatusTest::NormF::ScaleType type)
477  {
478  conv.stype = type;
479  }
480 
481  template<typename Scalar>
482  void NewtonSolverNOX<Scalar>::set_conv_iters(int iters)
483  {
484  conv.max_iters = iters;
485  }
486 
487  template<typename Scalar>
488  void NewtonSolverNOX<Scalar>::set_conv_abs_resid(double resid)
489  {
490  conv_flag.absresid = 1;
491  conv.abs_resid = resid;
492  }
493 
494  template<typename Scalar>
495  void NewtonSolverNOX<Scalar>::set_conv_rel_resid(double resid)
496  {
497  conv_flag.relresid = 1;
498  conv.rel_resid = resid;
499  }
500 
501  template<typename Scalar>
502  void NewtonSolverNOX<Scalar>::disable_abs_resid()
503  {
504  conv_flag.absresid = 0;
505  }
506 
507  template<typename Scalar>
508  void NewtonSolverNOX<Scalar>::disable_rel_resid()
509  {
510  conv_flag.relresid = 0;
511  }
512 
513  template<typename Scalar>
514  void NewtonSolverNOX<Scalar>::set_conv_update(double update)
515  {
516  conv_flag.update = 1;
517  conv.update = update;
518  }
519 
520  template<typename Scalar>
521  void NewtonSolverNOX<Scalar>::set_conv_wrms(double rtol, double atol)
522  {
523  conv_flag.wrms = 1;
524  conv.wrms_rtol = rtol;
525  conv.wrms_atol = atol;
526  }
527 
528  template class HERMES_API DiscreteProblemNOX < double > ;
529  //complex version of nox solver is not implemented
530  // template class HERMES_API DiscreteProblemNOX<std::complex<double> >;
531  template class HERMES_API NewtonSolverNOX < double > ;
532  //complex version of nox solver is not implemented
533  // template class HERMES_API NewtonSolverNOX<std::complex<double> >;
534  }
535 }
536 #endif
Definition: adapt.h:24
::xsd::cxx::tree::flags flags
Parsing and serialization flags.
::xsd::cxx::tree::time< char, simple_type > time
C++ type corresponding to the time XML Schema built-in type.
::xsd::cxx::tree::type type
C++ type corresponding to the anyType XML Schema built-in type.
int get_num_dofs() const
Returns the number of basis functions contained in the space.
Definition: space.cpp:286
static void project_global(SpaceSharedPtr< Scalar > space, MatrixFormVol< Scalar > *custom_projection_jacobian, VectorFormVol< Scalar > *custom_projection_residual, Scalar *target_vec)
This method allows to specify your own OG-projection form.