DOpE
instatreducedproblem.h
Go to the documentation of this file.
1 
24 #ifndef _INSTAT_REDUCED_PROBLEM_H_
25 #define _INSTAT_REDUCED_PROBLEM_H_
26 
27 #include "reducedprobleminterface.h"
28 #include "integrator.h"
29 #include "parameterreader.h"
30 #include "statevector.h"
31 #include "solutionextractor.h"
32 #include "pdeinterface.h"
33 #include "functionalinterface.h"
34 #include "dirichletdatainterface.h"
35 #include "dopeexception.h"
36 #include "instat_step_newtonsolver.h"
37 #include "fractional_step_theta_step_newtonsolver.h"
38 #include "newtonsolvermixeddims.h"
39 //#include "integratormixeddims.h"
40 #include "cglinearsolver.h"
41 #include "gmreslinearsolver.h"
42 #include "directlinearsolver.h"
43 #include "voidlinearsolver.h"
44 #include "constraintinterface.h"
45 #include "helper.h"
46 #include "dwrdatacontainer.h"
47 
48 #include <base/data_out_base.h>
49 #include <numerics/data_out.h>
50 #include <numerics/matrix_tools.h>
51 #include <numerics/vector_tools.h>
52 #include <base/function.h>
53 #include <lac/sparse_matrix.h>
54 #include <lac/compressed_simple_sparsity_pattern.h>
55 #include <lac/sparse_direct.h>
56 #include <lac/block_sparsity_pattern.h>
57 #include <lac/block_sparse_matrix.h>
58 #include <lac/vector.h>
59 
60 #include <fstream>
61 
62 namespace DOpE
63 {
79 template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER, typename CONTROLINTEGRATOR,
80  typename INTEGRATOR, typename PROBLEM, typename VECTOR, int dopedim, int dealdim>
81 class InstatReducedProblem: public ReducedProblemInterface<PROBLEM, VECTOR>
82 {
83  public:
96  template<typename INTEGRATORDATACONT>
97  InstatReducedProblem(PROBLEM *OP, std::string state_behavior,
98  ParameterReader &param_reader,
99  INTEGRATORDATACONT& idc,
100  int base_priority = 0);
101 
102 
116  template<typename STATEINTEGRATORDATACONT, typename CONTROLINTEGRATORCONT>
117  InstatReducedProblem(PROBLEM *OP, std::string state_behavior,
118  ParameterReader &param_reader, CONTROLINTEGRATORCONT& c_idc,
119  STATEINTEGRATORDATACONT & s_idc, int base_priority = 0);
120 
121  virtual ~InstatReducedProblem();
122 
123  /******************************************************/
124 
130  static void declare_params(ParameterReader &param_reader);
131 
132  /******************************************************/
133 
139  void ReInit();
140 
141  /******************************************************/
142 
148  bool ComputeReducedConstraints(const ControlVector<VECTOR>& q, ConstraintVector<VECTOR>& g);
149 
150  /******************************************************/
151 
157  void GetControlBoxConstraints(ControlVector<VECTOR>& lb, ControlVector<VECTOR>& ub);
158 
159 
160  /******************************************************/
161 
167  void ComputeReducedGradient(const ControlVector<VECTOR>& q, ControlVector<VECTOR>& gradient,
168  ControlVector<VECTOR>& gradient_transposed);
169 
170  /******************************************************/
171 
177  double ComputeReducedCostFunctional(const ControlVector<VECTOR>& q);
178 
179  /******************************************************/
180 
186  void ComputeReducedFunctionals(const ControlVector<VECTOR>& q);
187 
188  /******************************************************/
189 
195  void ComputeReducedHessianVector(const ControlVector<VECTOR>& q, const ControlVector<VECTOR>& direction,
196  ControlVector<VECTOR>& hessian_direction,
197  ControlVector<VECTOR>& hessian_direction_transposed);
198 
199  /******************************************************/
200 
220  template<class DWRC>
221  void
222  ComputeRefinementIndicators(DWRC&)
223  {
224  throw DOpEException("ExcNotImplemented",
225  "InstatReducedProblem::ComputeRefinementIndicators");
226  }
227 
228  /******************************************************/
229 
235  void StateSizeInfo(std::stringstream& out)
236  {
237  GetU().PrintInfos(out);
238  }
239 
240  /******************************************************/
241 
252  void WriteToFile(const VECTOR &v, std::string name, std::string outfile,
253  std::string dof_type, std::string filetype);
254 
255  /******************************************************/
256 
267  void WriteToFile(const ControlVector<VECTOR> &v, std::string name, std::string outfile,
268  std::string dof_type, std::string filetype);
269 
270  /******************************************************/
271 
279  void WriteToFile(const std::vector<double> &v, std::string outfile);
280 
281  protected:
282  const StateVector<VECTOR> & GetU() const
283  {
284  return _u;
285  }
286  StateVector<VECTOR> & GetU()
287  {
288  return _u;
289  }
290  StateVector<VECTOR> & GetZ()
291  {
292  return _z;
293  }
294  StateVector<VECTOR> & GetDU()
295  {
296  return _du;
297  }
298  StateVector<VECTOR> & GetDZ()
299  {
300  return _dz;
301  }
302 
303  NONLINEARSOLVER& GetNonlinearSolver(std::string type);
304  CONTROLNONLINEARSOLVER& GetControlNonlinearSolver();
305  INTEGRATOR& GetIntegrator()
306  {
307  return _integrator;
308  }
309  CONTROLINTEGRATOR& GetControlIntegrator()
310  {
311  return _control_integrator;
312  }
313 
314  /******************************************************/
315 
324  void ComputeTimeFunctionals(unsigned int step, unsigned int num_steps);
333  void ComputeReducedState(const ControlVector<VECTOR>& q);
334 
335  /******************************************************/
336 
346  void ComputeReducedAdjoint(const ControlVector<VECTOR>& q);
347 
348  /******************************************************/
349 
361  template<typename PDE>
362  void ForwardTimeLoop(PDE& problem, StateVector<VECTOR>& sol, std::string outname, bool eval_funcs);
363 
364  /******************************************************/
365 
374  template<typename PDE>
375  void BackwardTimeLoop(PDE& problem, StateVector<VECTOR>& sol, std::string outname);
376 
377  private:
378 
379  StateVector<VECTOR> _u;
380  StateVector<VECTOR> _z;
381  StateVector<VECTOR> _du;
382  StateVector<VECTOR> _dz;
383 
384  INTEGRATOR _integrator;
385  CONTROLINTEGRATOR _control_integrator;
386  NONLINEARSOLVER _nonlinear_state_solver;
387  NONLINEARSOLVER _nonlinear_adjoint_solver;
388  CONTROLNONLINEARSOLVER _nonlinear_gradient_solver;
389 
390  bool _build_state_matrix, _build_adjoint_matrix, _build_control_matrix;
391  bool _state_reinit, _adjoint_reinit, _gradient_reinit;
392 
393  bool _project_initial_data;
394 
395  friend class SolutionExtractor<InstatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
396  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR,dopedim, dealdim>, VECTOR > ;
397 };
398 
399 /*************************************************************************/
400 /*****************************IMPLEMENTATION******************************/
401 /*************************************************************************/
402 using namespace dealii;
403 
404 /******************************************************/
405 template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER, typename CONTROLINTEGRATOR,
406  typename INTEGRATOR, typename PROBLEM, typename VECTOR, int dopedim,
407  int dealdim>
408 void InstatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER, CONTROLINTEGRATOR, INTEGRATOR,
409  PROBLEM, VECTOR, dopedim, dealdim>::declare_params(
410  ParameterReader &param_reader)
411 {
412  NONLINEARSOLVER::declare_params(param_reader);
413 }
414 /******************************************************/
415 
416  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
417  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
418  typename VECTOR, int dopedim, int dealdim>
419  template<typename INTEGRATORDATACONT>
420  InstatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
421  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::InstatReducedProblem(
422  PROBLEM *OP,
423  std::string state_behavior,
424  ParameterReader &param_reader,
425  INTEGRATORDATACONT& idc,
426  int base_priority) :
427  ReducedProblemInterface<PROBLEM, VECTOR> (OP,
428  base_priority),
429  _u(OP->GetSpaceTimeHandler(), state_behavior, param_reader),
430  _z(OP->GetSpaceTimeHandler(), state_behavior, param_reader),
431  _du(OP->GetSpaceTimeHandler(), state_behavior, param_reader),
432  _dz(OP->GetSpaceTimeHandler(), state_behavior, param_reader),
433  _integrator(idc),
434  _control_integrator(idc),
435  _nonlinear_state_solver(_integrator, param_reader),
436  _nonlinear_adjoint_solver(_integrator, param_reader),
437  _nonlinear_gradient_solver(_control_integrator, param_reader)
438  {
439  //Solvers should be ReInited
440  {
441  _state_reinit = true;
442  _adjoint_reinit = true;
443  _gradient_reinit = true;
444  }
445  if(this->GetProblem()->GetSpaceTimeHandler()->GetControlType() == DOpEtypes::ControlType::stationary)
446  {
447  throw DOpEException("ControlType: stationary has no meeaning here." , "InstatReducedProblem::InstatReducedProblem");
448  }
449  }
450 
451 /******************************************************/
452  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
453  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
454  typename VECTOR, int dopedim, int dealdim>
455  template<typename STATEINTEGRATORDATACONT, typename CONTROLINTEGRATORCONT>
456  InstatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
457  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::InstatReducedProblem(
458  PROBLEM *OP, std::string state_behavior,
459  ParameterReader &param_reader,
460  CONTROLINTEGRATORCONT& c_idc,
461  STATEINTEGRATORDATACONT & s_idc,
462  int base_priority) :
463  ReducedProblemInterface<PROBLEM, VECTOR> (OP,
464  base_priority),
465  _u(OP->GetSpaceTimeHandler(), state_behavior, param_reader),
466  _z(OP->GetSpaceTimeHandler(), state_behavior, param_reader),
467  _du(OP->GetSpaceTimeHandler(), state_behavior, param_reader),
468  _dz(OP->GetSpaceTimeHandler(), state_behavior, param_reader),
469  _integrator(s_idc),
470  _control_integrator(c_idc),
471  _nonlinear_state_solver(_integrator, param_reader),
472  _nonlinear_adjoint_solver(_integrator, param_reader),
473  _nonlinear_gradient_solver(_control_integrator, param_reader)
474  {
475  //Solvers should be ReInited
476  {
477  _state_reinit = true;
478  _adjoint_reinit = true;
479  _gradient_reinit = true;
480  }
481  if(this->GetProblem()->GetSpaceTimeHandler()->GetControlType() == DOpEtypes::ControlType::stationary)
482  {
483  throw DOpEException("ControlType: stationary has no meeaning here." , "InstatReducedProblem::InstatReducedProblem");
484  }
485  }
486 
487 /******************************************************/
488 
489 template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER, typename CONTROLINTEGRATOR,
490  typename INTEGRATOR, typename PROBLEM, typename VECTOR, int dopedim,
491  int dealdim>
492 InstatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER, CONTROLINTEGRATOR, INTEGRATOR,
493  PROBLEM, VECTOR, dopedim, dealdim>::~InstatReducedProblem()
494 {
495 }
496 
497 /******************************************************/
498 
499 template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER, typename CONTROLINTEGRATOR,
500  typename INTEGRATOR, typename PROBLEM, typename VECTOR, int dopedim,
501  int dealdim>
502 NONLINEARSOLVER& InstatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER, CONTROLINTEGRATOR,
503  INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::GetNonlinearSolver(std::string type)
504 {
505  if ((type == "state") || (type == "tangent"))
506  {
507  return _nonlinear_state_solver;
508  }
509  else if ((type == "adjoint") || (type == "adjoint_hessian"))
510  {
511  return _nonlinear_adjoint_solver;
512  }
513  else
514  {
515  throw DOpEException("No Solver for Problem type:`" + type + "' found",
516  "InstatReducedProblem::GetNonlinearSolver");
517 
518  }
519 }
520 /******************************************************/
521 
522 template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER, typename CONTROLINTEGRATOR,
523  typename INTEGRATOR, typename PROBLEM, typename VECTOR, int dopedim,
524  int dealdim>
525 CONTROLNONLINEARSOLVER& InstatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
526  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::GetControlNonlinearSolver()
527 {
528  if ((this->GetProblem()->GetType() == "gradient") || (this->GetProblem()->GetType() == "hessian"))
529  {
530  return _nonlinear_gradient_solver;
531  }
532  else
533  {
534  throw DOpEException("No Solver for Problem type:`" + this->GetProblem()->GetType() + "' found",
535  "InstatReducedProblem::GetControlNonlinearSolver");
536 
537  }
538 }
539 /******************************************************/
540 
541 template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER, typename CONTROLINTEGRATOR,
542  typename INTEGRATOR, typename PROBLEM, typename VECTOR, int dopedim,
543  int dealdim>
544 void InstatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER, CONTROLINTEGRATOR, INTEGRATOR,
545  PROBLEM, VECTOR, dopedim, dealdim>::ReInit()
546 {
547  ReducedProblemInterface<PROBLEM, VECTOR>::ReInit();
548 
549  //Some Solvers must be reinited when called
550  // Better have subproblems, so that solver can be reinited here
551  {
552  _state_reinit = true;
553  _adjoint_reinit = true;
554  _gradient_reinit = true;
555  }
556 
557  _build_state_matrix = true;
558  _build_adjoint_matrix = true;
559 
560  GetU().ReInit();
561  GetZ().ReInit();
562  GetDU().ReInit();
563  GetDZ().ReInit();
564 
565  _build_control_matrix = true;
566 }
567 
568 /******************************************************/
569 
570 template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER, typename CONTROLINTEGRATOR,
571  typename INTEGRATOR, typename PROBLEM, typename VECTOR, int dopedim,
572  int dealdim>
573 void InstatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER, CONTROLINTEGRATOR, INTEGRATOR,
574  PROBLEM, VECTOR, dopedim, dealdim>::ComputeReducedState(const ControlVector<VECTOR>& q)
575 {
576  this->InitializeFunctionalValues(this->GetProblem()->GetNFunctionals() + 1);
577 
578  this->GetOutputHandler()->Write("Computing State Solution:", 4 + this->GetBasePriority());
579 
580  this->SetProblemType("state");
581  auto& problem = this->GetProblem()->GetStateProblem();
582 
583  if (_state_reinit == true)
584  {
585  GetNonlinearSolver("state").ReInit(problem);
586  _state_reinit = false;
587  }
588 
589  this->GetProblem()->AddAuxiliaryControl(&q,"control");
590  this->ForwardTimeLoop(problem,this->GetU(),"State",true);
591  this->GetProblem()->DeleteAuxiliaryControl("control");
592 }
593 /******************************************************/
594 
595 template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER, typename CONTROLINTEGRATOR,
596  typename INTEGRATOR, typename PROBLEM, typename VECTOR, int dopedim,
597  int dealdim>
598 bool InstatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER, CONTROLINTEGRATOR, INTEGRATOR,
599  PROBLEM, VECTOR, dopedim, dealdim>::ComputeReducedConstraints(
600  const ControlVector<VECTOR>& /*q*/,
601  ConstraintVector<VECTOR>& /*g*/)
602 {
603  abort();
604 }
605 
606 /******************************************************/
607 
608 template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER, typename CONTROLINTEGRATOR,
609  typename INTEGRATOR, typename PROBLEM, typename VECTOR, int dopedim,
610  int dealdim>
611 void InstatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER, CONTROLINTEGRATOR, INTEGRATOR,
612  PROBLEM, VECTOR, dopedim, dealdim>::GetControlBoxConstraints(ControlVector<VECTOR>& /*lb*/, ControlVector<VECTOR>& /*ub*/)
613 {
614  abort();
615 }
616 
617 /******************************************************/
618 
619 template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER, typename CONTROLINTEGRATOR,
620  typename INTEGRATOR, typename PROBLEM, typename VECTOR, int dopedim,
621  int dealdim>
622 void InstatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER, CONTROLINTEGRATOR, INTEGRATOR,
623  PROBLEM, VECTOR, dopedim, dealdim>::ComputeReducedAdjoint(
624  const ControlVector<VECTOR>& q)
625 {
626  this->GetOutputHandler()->Write("Computing Adjoint Solution:", 4 + this->GetBasePriority());
627 
628  this->SetProblemType("adjoint");
629  auto& problem = this->GetProblem()->GetAdjointProblem();
630  if (_adjoint_reinit == true)
631  {
632  GetNonlinearSolver("adjoint").ReInit(problem);
633  _adjoint_reinit = false;
634  }
635 
636  this->GetProblem()->AddAuxiliaryState(&(this->GetU()),"state");
637  this->GetProblem()->AddAuxiliaryControl(&q,"control");
638  this->BackwardTimeLoop(problem,this->GetZ(),"Adjoint");
639  this->GetProblem()->DeleteAuxiliaryControl("control");
640  this->GetProblem()->DeleteAuxiliaryState("state");
641 }
642 
643 /******************************************************/
644 
645 template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER, typename CONTROLINTEGRATOR,
646  typename INTEGRATOR, typename PROBLEM, typename VECTOR, int dopedim,
647  int dealdim>
648 void InstatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER, CONTROLINTEGRATOR, INTEGRATOR,
649  PROBLEM, VECTOR, dopedim, dealdim>::ComputeReducedGradient(
650  const ControlVector<VECTOR>& q,
651  ControlVector<VECTOR>& gradient,
652  ControlVector<VECTOR>& gradient_transposed)
653 {
654  this->ComputeReducedAdjoint(q);
655 
656  this->GetOutputHandler()->Write("Computing Reduced Gradient:",
657  4 + this->GetBasePriority());
658  if (this->GetProblem()->HasControlInDirichletData())
659  {
660  throw DOpEException("Control in Dirichlet Data for instationary problems not yet implemented!"
661  ,"InstatReducedProblem::ComputeReducedGradient");
662  }
663 
664  this->SetProblemType("gradient");
665  if (_gradient_reinit == true)
666  {
667  GetControlNonlinearSolver().ReInit(*(this->GetProblem()));
668  _state_reinit = false;
669  }
670 
671  if(this->GetProblem()->GetSpaceTimeHandler()->GetControlType() == DOpEtypes::ControlType::initial)
672  {
673  //Set time to initial
674  const std::vector<double> times =
675  this->GetProblem()->GetSpaceTimeHandler()->GetTimes();
676  const unsigned int
677  n_dofs_per_interval =
678  this->GetProblem()->GetSpaceTimeHandler()->GetTimeDoFHandler().GetLocalNbrOfDoFs();
679  std::vector<unsigned int> local_to_global(n_dofs_per_interval);
680  {
681  TimeIterator it =
682  this->GetProblem()->GetSpaceTimeHandler()->GetTimeDoFHandler().first_interval();
683  it.get_time_dof_indices(local_to_global);
684  this->GetProblem()->SetTime(times[local_to_global[0]], it);
685  GetU().SetTimeDoFNumber(local_to_global[0], it);
686  GetZ().SetTimeDoFNumber(local_to_global[0], it);
687  }
688 
689  this->GetProblem()->AddAuxiliaryToIntegrator(this->GetControlIntegrator());
690 
691  if (dopedim == dealdim)
692  {
693  this->GetControlIntegrator().AddDomainData("control",
694  &(q.GetSpacialVector()));
695  }
696  else if (dopedim == 0)
697  {
698  this->GetControlIntegrator().AddParamData("control",
699  &(q.GetSpacialVectorCopy()));
700  }
701  else
702  {
703  throw DOpEException("dopedim not implemented",
704  "InstatReducedProblem::ComputeReducedGradient");
705  }
706 
707  this->GetControlIntegrator().AddDomainData("state",
708  &(GetU().GetSpacialVector()));
709  this->GetControlIntegrator().AddDomainData("adjoint",
710  &(GetZ().GetSpacialVector()));
711  gradient_transposed = 0.;
712  if (dopedim == dealdim)
713  {
714  this->GetControlIntegrator().AddDomainData("last_newton_solution",
715  &(gradient_transposed.GetSpacialVector()));
716  this->GetControlIntegrator().ComputeNonlinearResidual(
717  *(this->GetProblem()), gradient.GetSpacialVector(), true);
718  this->GetControlIntegrator().DeleteDomainData("last_newton_solution");
719  }
720  else if (dopedim == 0)
721  {
722  this->GetControlIntegrator().AddParamData("last_newton_solution",
723  &(gradient_transposed.GetSpacialVectorCopy()));
724  this->GetControlIntegrator().ComputeNonlinearResidual(
725  *(this->GetProblem()), gradient.GetSpacialVector(), true);
726 
727  this->GetControlIntegrator().DeleteParamData("last_newton_solution");
728  gradient_transposed.UnLockCopy();
729  }
730 
731  gradient *= -1.;
732  gradient_transposed = gradient;
733 
734  //Compute l^2 representation of the Gradient
735 
736  _build_control_matrix = this->GetControlNonlinearSolver().NonlinearSolve(
737  *(this->GetProblem()), gradient_transposed.GetSpacialVector(), true,
738  _build_control_matrix);
739  if (dopedim == dealdim)
740  {
741  this->GetControlIntegrator().DeleteDomainData("control");
742  }
743  else if (dopedim == 0)
744  {
745  this->GetControlIntegrator().DeleteParamData("control");
746  q.UnLockCopy();
747  }
748  else
749  {
750  throw DOpEException("dopedim not implemented",
751  "InstatReducedProblem::ComputeReducedGradient");
752  }
753  this->GetControlIntegrator().DeleteDomainData("state");
754  this->GetControlIntegrator().DeleteDomainData("adjoint");
755 
756  this->GetProblem()->DeleteAuxiliaryFromIntegrator(
757  this->GetControlIntegrator());
758 
759  this->GetOutputHandler()->Write(gradient,
760  "Gradient" + this->GetPostIndex(), this->GetProblem()->GetDoFType());
761  this->GetOutputHandler()->Write(gradient_transposed,
762  "Gradient_Transposed" + this->GetPostIndex(),
763  this->GetProblem()->GetDoFType());
764  }
765  else
766  {
767  std::stringstream out;
768  out << "Unknown ControlType: "<<this->GetProblem()->GetSpaceTimeHandler()->GetControlType();
769  throw DOpEException(out.str(), "InstatReducedProblem::ComputeReducedGradient");
770  }
771 }
772 
773 /******************************************************/
774 
775 template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER, typename CONTROLINTEGRATOR,
776  typename INTEGRATOR, typename PROBLEM, typename VECTOR, int dopedim,
777  int dealdim>
778 double InstatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER, CONTROLINTEGRATOR, INTEGRATOR,
779  PROBLEM, VECTOR, dopedim, dealdim>::ComputeReducedCostFunctional(
780  const ControlVector<VECTOR>& q)
781 {
782  this->ComputeReducedState(q);
783 
784  if (this->GetFunctionalValues()[0].size() != 1)
785  {
786  if (this->GetFunctionalValues()[0].size() == 0)
787  throw DOpEException(
788  "Apparently the CostFunctional was never evaluated! \n\tCheck if the return value of `NeedTimes' is set correctly.",
789  "InstatReducedProblem::ComputeReducedCostFunctional");
790  else
791  throw DOpEException(
792  "The CostFunctional has been evaluated too many times! \n\tCheck if the return value of `NeedTimes' is set correctly.",
793  "InstatReducedProblem::ComputeReducedCostFunctional");
794  }
795  return this->GetFunctionalValues()[0][0];
796 }
797 
798 /******************************************************/
799 
800 template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER, typename CONTROLINTEGRATOR,
801  typename INTEGRATOR, typename PROBLEM, typename VECTOR, int dopedim,
802  int dealdim>
803 void InstatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER, CONTROLINTEGRATOR, INTEGRATOR,
804  PROBLEM, VECTOR, dopedim, dealdim>::ComputeReducedFunctionals(
805  const ControlVector<VECTOR>& /*q*/)
806 {
807  //We dont need q as the values are precomputed during Solve State...
808  this->GetOutputHandler()->Write("Computing Functionals:" + this->GetBasePriority(), 4);
809 
810  for (unsigned int i = 0; i < this->GetProblem()->GetNFunctionals(); i++)
811  {
812  this->SetProblemType("aux_functional", i);
813  if (this->GetProblem()->GetFunctionalType().find("timelocal"))
814  {
815  if (this->GetFunctionalValues()[i + 1].size() == 1)
816  {
817  std::stringstream out;
818  this->GetOutputHandler()->InitOut(out);
819  out << this->GetProblem()->GetFunctionalName() << ": " << this->GetFunctionalValues()[i + 1][0];
820  this->GetOutputHandler()->Write(out, 2 + this->GetBasePriority());
821  }
822  else if (this->GetFunctionalValues()[i + 1].size() > 1)
823  {
824  if (this->GetFunctionalValues()[i + 1].size()
825  == this->GetProblem()->GetSpaceTimeHandler()->GetMaxTimePoint() + 1)
826  {
827  std::stringstream out;
828  this->GetOutputHandler()->InitOut(out);
829  out << this->GetProblem()->GetFunctionalName() << " too large. Writing to file instead: ";
830  this->GetOutputHandler()->Write(out, 2 + this->GetBasePriority());
831  this->GetOutputHandler()->Write(this->GetFunctionalValues()[i + 1],
832  this->GetProblem()->GetFunctionalName()
833  + this->GetPostIndex(), "time");
834  }
835  else
836  {
837  std::stringstream out;
838  this->GetOutputHandler()->InitOut(out);
839  out << this->GetProblem()->GetFunctionalName() << ": ";
840  for (unsigned int k = 0; k < this->GetFunctionalValues()[i + 1].size(); k++)
841  out << this->GetFunctionalValues()[i + 1][k] << " ";
842  this->GetOutputHandler()->Write(out, 2 + this->GetBasePriority());
843  }
844  }
845  else
846  {
847  throw DOpEException("Functional: " + this->GetProblem()->GetFunctionalType()
848  + " was not evaluated ever!", "InstatReducedProblem::ComputeFunctionals");
849  }
850  }
851  else if (this->GetProblem()->GetFunctionalType().find("timedistributed"))
852  {
853  std::stringstream out;
854  this->GetOutputHandler()->InitOut(out);
855  out << this->GetProblem()->GetFunctionalName() << ": " << this->GetFunctionalValues()[i + 1][0];
856  this->GetOutputHandler()->Write(out, 2 + this->GetBasePriority());
857  }
858  else
859  {
860  throw DOpEException("Unknown Functional Type: " + this->GetProblem()->GetFunctionalType(),
861  "InstatReducedProblem::ComputeFunctionals");
862  }
863  }
864 }
865 
866 /******************************************************/
867 
868 template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER, typename CONTROLINTEGRATOR,
869  typename INTEGRATOR, typename PROBLEM, typename VECTOR, int dopedim,
870  int dealdim>
871 void InstatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER, CONTROLINTEGRATOR, INTEGRATOR,
872  PROBLEM, VECTOR, dopedim, dealdim>::ComputeReducedHessianVector(
873  const ControlVector<VECTOR>& q,
874  const ControlVector<VECTOR>& direction,
875  ControlVector<VECTOR>& hessian_direction,
876  ControlVector<VECTOR>& hessian_direction_transposed)
877 {
878  this->GetOutputHandler()->Write("Computing ReducedHessianVector:",
879  4 + this->GetBasePriority());
880 
881  //Solving the Tangent Problem
882  {
883  this->GetOutputHandler()->Write("\tSolving Tangent:",
884  5 + this->GetBasePriority());
885  this->SetProblemType("tangent");
886  auto& problem = this->GetProblem()->GetTangentProblem();
887 
888  this->GetProblem()->AddAuxiliaryControl(&q,"control");
889  this->GetProblem()->AddAuxiliaryControl(&direction,"dq");
890  this->GetProblem()->AddAuxiliaryState(&(this->GetU()),"state");
891  this->GetProblem()->AddAuxiliaryState(&(this->GetZ()),"adjoint");
892 
893  this->ForwardTimeLoop(problem,this->GetDU(),"Tangent",false);
894 
895  }
896  //Solving the Adjoint-Hessian Problem
897  {
898  this->GetOutputHandler()->Write("\tSolving Adjoint Hessian:",
899  5 + this->GetBasePriority());
900  this->SetProblemType("adjoint_hessian");
901  auto& problem = this->GetProblem()->GetAdjointHessianProblem();
902 
903  this->GetProblem()->AddAuxiliaryState(&(this->GetDU()),"tangent");
904 
905  this->BackwardTimeLoop(problem,this->GetDZ(),"Hessian");
906  }
907  if (this->GetProblem()->HasControlInDirichletData())
908  {
909  throw DOpEException("Control in Dirichlet Data for instationary problems not yet implemented!"
910  ,"InstatReducedProblem::ComputeReducedHessianVector");
911  }
912 
913  //Computing Hessian times Vector Representation
914  {
915  this->GetOutputHandler()->Write(
916  "\tComputing Representation of the Hessian:",
917  5 + this->GetBasePriority());
918  this->SetProblemType("hessian");
919 
920  this->GetProblem()->AddAuxiliaryState(&(this->GetDZ()),"adjoint_hessian");
921 
922  if(this->GetProblem()->GetSpaceTimeHandler()->GetControlType() == DOpEtypes::ControlType::initial)
923  {
924  //Set time to initial
925  const std::vector<double> times =
926  this->GetProblem()->GetSpaceTimeHandler()->GetTimes();
927  const unsigned int
928  n_dofs_per_interval =
929  this->GetProblem()->GetSpaceTimeHandler()->GetTimeDoFHandler().GetLocalNbrOfDoFs();
930  std::vector<unsigned int> local_to_global(n_dofs_per_interval);
931  {
932  TimeIterator it =
933  this->GetProblem()->GetSpaceTimeHandler()->GetTimeDoFHandler().first_interval();
934  it.get_time_dof_indices(local_to_global);
935  this->GetProblem()->SetTime(times[local_to_global[0]], it);
936  GetU().SetTimeDoFNumber(local_to_global[0], it);
937  GetZ().SetTimeDoFNumber(local_to_global[0], it);
938  GetDU().SetTimeDoFNumber(local_to_global[0], it);
939  GetDZ().SetTimeDoFNumber(local_to_global[0], it);
940  }
941 
942  this->GetProblem()->AddAuxiliaryToIntegrator(
943  this->GetControlIntegrator());
944 
945  hessian_direction_transposed = 0.;
946  if (dopedim == dealdim)
947  {
948  this->GetControlIntegrator().AddDomainData("last_newton_solution",
949  &(hessian_direction_transposed.GetSpacialVector()));
950  this->GetControlIntegrator().ComputeNonlinearResidual(
951  *(this->GetProblem()), hessian_direction.GetSpacialVector(),
952  true);
953  this->GetControlIntegrator().DeleteDomainData("last_newton_solution");
954  }
955  else if (dopedim == 0)
956  {
957  this->GetControlIntegrator().AddParamData("last_newton_solution",
958  &(hessian_direction_transposed.GetSpacialVectorCopy()));
959  this->GetControlIntegrator().ComputeNonlinearResidual(
960  *(this->GetProblem()), hessian_direction.GetSpacialVector(),
961  true);
962  this->GetControlIntegrator().DeleteParamData("last_newton_solution");
963  hessian_direction_transposed.UnLockCopy();
964  }
965 
966  hessian_direction *= -1.;
967  hessian_direction_transposed = hessian_direction;
968  //Compute l^2 representation of the HessianVector
969  //hessian Matrix is the same as control matrix
970  _build_control_matrix =
971  this->GetControlNonlinearSolver().NonlinearSolve(
972  *(this->GetProblem()),
973  hessian_direction_transposed.GetSpacialVector(), true,
974  _build_control_matrix);
975 
976  this->GetOutputHandler()->Write(hessian_direction,
977  "HessianDirection" + this->GetPostIndex(),
978  this->GetProblem()->GetDoFType());
979  this->GetOutputHandler()->Write(hessian_direction_transposed,
980  "HessianDirection_Transposed" + this->GetPostIndex(),
981  this->GetProblem()->GetDoFType());
982 
983  this->GetProblem()->DeleteAuxiliaryFromIntegrator(
984  this->GetControlIntegrator());
985  }//Endof the case of control in the initial values
986  else
987  {
988  std::stringstream out;
989  out << "Unknown ControlType: "<<this->GetProblem()->GetSpaceTimeHandler()->GetControlType();
990  throw DOpEException(out.str(), "InstatReducedProblem::ComputeReducedHessianVector");
991  }
992  }//End of HessianVector Repr.
993 
994  //Cleaning
995  this->GetProblem()->DeleteAuxiliaryControl("control");
996  this->GetProblem()->DeleteAuxiliaryControl("dq");
997  this->GetProblem()->DeleteAuxiliaryState("state");
998  this->GetProblem()->DeleteAuxiliaryState("adjoint");
999  this->GetProblem()->DeleteAuxiliaryState("tangent");
1000  this->GetProblem()->DeleteAuxiliaryState("adjoint_hessian");
1001 }
1002 
1003 /******************************************************/
1004 
1005 template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER, typename CONTROLINTEGRATOR,
1006  typename INTEGRATOR, typename PROBLEM, typename VECTOR, int dopedim,
1007  int dealdim>
1008 void InstatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER, CONTROLINTEGRATOR, INTEGRATOR,
1009  PROBLEM, VECTOR, dopedim, dealdim>::ComputeTimeFunctionals(unsigned int step, unsigned int num_steps)
1010 {
1011 
1012  this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
1013 
1014  this->GetIntegrator().AddDomainData("state", &(GetU().GetSpacialVector()));
1015  double ret = 0;
1016  bool found = false;
1017  {//CostFunctional
1018  this->SetProblemType("cost_functional");
1019  if (this->GetProblem()->NeedTimeFunctional())
1020  {
1021  if (this->GetProblem()->GetFunctionalType().find("domain") != std::string::npos)
1022  {
1023  found = true;
1024  ret += this->GetIntegrator().ComputeDomainScalar(*(this->GetProblem()));
1025  }
1026  if (this->GetProblem()->GetFunctionalType().find("point") != std::string::npos)
1027  {
1028  found = true;
1029  ret += this->GetIntegrator().ComputePointScalar(*(this->GetProblem()));
1030  }
1031  if (this->GetProblem()->GetFunctionalType().find("boundary") != std::string::npos)
1032  {
1033  found = true;
1034  ret += this->GetIntegrator().ComputeBoundaryScalar(*(this->GetProblem()));
1035  }
1036  if (this->GetProblem()->GetFunctionalType().find("face") != std::string::npos)
1037  {
1038  found = true;
1039  ret += this->GetIntegrator().ComputeFaceScalar(*(this->GetProblem()));
1040  }
1041 
1042  if (!found)
1043  {
1044  throw DOpEException("Unknown Functional Type: " + this->GetProblem()->GetFunctionalType(),
1045  "InstatReducedProblem::ComputeTimeFunctionals");
1046  }
1047  //Wert speichern
1048  if (this->GetProblem()->GetFunctionalType().find("timelocal"))
1049  {
1050  if (this->GetFunctionalValues()[0].size() != 1)
1051  {
1052  this->GetFunctionalValues()[0].resize(1);
1053  this->GetFunctionalValues()[0][0] = 0.;
1054  }
1055  this->GetFunctionalValues()[0][0] += ret;
1056  }
1057  else if (this->GetProblem()->GetFunctionalType().find("timedistributed"))
1058  {//TODO was passiert hier? Vermutlich sollte hier spaeter Zeitintegration durchgefuehrt werden?
1059  if (this->GetFunctionalValues()[0].size() != 1)
1060  {
1061  this->GetFunctionalValues()[0].resize(1);
1062  }
1063  double w = 0.;
1064  if ((step == 0))
1065  {
1066  w = 0.5 * (this->GetProblem()->GetSpaceTimeHandler()->GetTime(step + 1)
1067  - this->GetProblem()->GetSpaceTimeHandler()->GetTime(step));
1068  }
1069  else if (step + 1 == num_steps)
1070  {
1071  w = 0.5 * (this->GetProblem()->GetSpaceTimeHandler()->GetTime(step)
1072  - this->GetProblem()->GetSpaceTimeHandler()->GetTime(step - 1));
1073  }
1074  else
1075  {
1076  w = 0.5 * (this->GetProblem()->GetSpaceTimeHandler()->GetTime(step + 1)
1077  - this->GetProblem()->GetSpaceTimeHandler()->GetTime(step));
1078  w += 0.5 * (this->GetProblem()->GetSpaceTimeHandler()->GetTime(step)
1079  - this->GetProblem()->GetSpaceTimeHandler()->GetTime(step - 1));
1080  }
1081  this->GetFunctionalValues()[0][0] += w * ret;
1082  }
1083  else
1084  {
1085  throw DOpEException("Unknown Functional Type: " + this->GetProblem()->GetFunctionalType(),
1086  "InstatReducedProblem::ComputeTimeFunctionals");
1087  }
1088  }
1089  }
1090  {//Aux Functionals
1091  for (unsigned int i = 0; i < this->GetProblem()->GetNFunctionals(); i++)
1092  {
1093  ret = 0;
1094  found = false;
1095  this->SetProblemType("aux_functional", i);
1096  if (this->GetProblem()->NeedTimeFunctional())
1097  {
1098  if (this->GetProblem()->GetFunctionalType().find("domain") != std::string::npos)
1099  {
1100  found = true;
1101  ret += this->GetIntegrator().ComputeDomainScalar(*(this->GetProblem()));
1102  }
1103  if (this->GetProblem()->GetFunctionalType().find("point") != std::string::npos)
1104  {
1105  found = true;
1106  ret += this->GetIntegrator().ComputePointScalar(*(this->GetProblem()));
1107  }
1108  if (this->GetProblem()->GetFunctionalType().find("boundary") != std::string::npos)
1109  {
1110  found = true;
1111  ret += this->GetIntegrator().ComputeBoundaryScalar(*(this->GetProblem()));
1112  }
1113  if (this->GetProblem()->GetFunctionalType().find("face") != std::string::npos)
1114  {
1115  found = true;
1116  ret += this->GetIntegrator().ComputeFaceScalar(*(this->GetProblem()));
1117  }
1118 
1119  if (!found)
1120  {
1121  throw DOpEException(
1122  "Unknown Functional Type: " + this->GetProblem()->GetFunctionalType(),
1123  "InstatReducedProblem::ComputeTimeFunctionals");
1124  }
1125  //Wert speichern
1126  if (this->GetProblem()->GetFunctionalType().find("timelocal"))
1127  {
1128  std::stringstream out;
1129  this->GetOutputHandler()->InitOut(out);
1130  out << "\t" << this->GetProblem()->GetFunctionalName() << ": " << ret;
1131  this->GetOutputHandler()->Write(out, 5 + this->GetBasePriority());
1132  this->GetFunctionalValues()[i + 1].push_back(ret);
1133  }
1134  else if (this->GetProblem()->GetFunctionalType().find("timedistributed"))
1135  {
1136  if (this->GetFunctionalValues()[i + 1].size() != 1)
1137  {
1138  this->GetFunctionalValues()[i + 1].resize(1);
1139  }
1140  double w = 0.;
1141  if ((step == 0))
1142  {
1143  w = 0.5 * (this->GetProblem()->GetSpaceTimeHandler()->GetTime(step + 1)
1144  - this->GetProblem()->GetSpaceTimeHandler()->GetTime(step));
1145  }
1146  else if (step + 1 == num_steps)
1147  {
1148  w = 0.5 * (this->GetProblem()->GetSpaceTimeHandler()->GetTime(step)
1149  - this->GetProblem()->GetSpaceTimeHandler()->GetTime(step - 1));
1150  }
1151  else
1152  {
1153  w = 0.5 * (this->GetProblem()->GetSpaceTimeHandler()->GetTime(step + 1)
1154  - this->GetProblem()->GetSpaceTimeHandler()->GetTime(step));
1155  w += 0.5 * (this->GetProblem()->GetSpaceTimeHandler()->GetTime(step)
1156  - this->GetProblem()->GetSpaceTimeHandler()->GetTime(step - 1));
1157  }
1158  this->GetFunctionalValues()[i + 1][0] += w * ret;
1159  }
1160  else
1161  {
1162  throw DOpEException(
1163  "Unknown Functional Type: " + this->GetProblem()->GetFunctionalType(),
1164  "InstatReducedProblem::ComputeTimeFunctionals");
1165  }
1166  }
1167  }
1168  }
1169  this->GetIntegrator().DeleteDomainData("state");
1170  this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
1171 
1172 }
1173 
1174 /******************************************************/
1175 template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER, typename CONTROLINTEGRATOR,
1176 typename INTEGRATOR, typename PROBLEM, typename VECTOR, int dopedim,
1177 int dealdim>
1178 void InstatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER, CONTROLINTEGRATOR, INTEGRATOR,
1179 PROBLEM, VECTOR, dopedim, dealdim>::WriteToFile(const VECTOR &v, std::string name, std::string outfile, std::string dof_type, std::string filetype)
1180 {
1181  if (dof_type == "state")
1182  {
1183  auto& data_out = this->GetProblem()->GetSpaceTimeHandler()->GetDataOut();
1184  data_out.attach_dof_handler(this->GetProblem()->GetSpaceTimeHandler()->GetStateDoFHandler());
1185 
1186  data_out.add_data_vector(v, name);
1187  data_out.build_patches();
1188 
1189  std::ofstream output(outfile.c_str());
1190 
1191  if (filetype == ".vtk")
1192  {
1193  data_out.write_vtk(output);
1194  }
1195  else if (filetype == ".gpl")
1196  {
1197  data_out.write_gnuplot(output);
1198  }
1199  else
1200  {
1201  throw DOpEException("Don't know how to write filetype `" + filetype + "'!",
1202  "InstatReducedProblem::WriteToFile");
1203  }
1204  data_out.clear();
1205  }
1206  else if (dof_type == "control")
1207  {
1208 #if dope_dimension >0
1209  auto& data_out = this->GetProblem()->GetSpaceTimeHandler()->GetDataOut();
1210  data_out.attach_dof_handler (this->GetProblem()->GetSpaceTimeHandler()->GetControlDoFHandler());
1211 
1212  data_out.add_data_vector (v,name);
1213  data_out.build_patches ();
1214 
1215  std::ofstream output(outfile.c_str());
1216 
1217  if(filetype == ".vtk")
1218  {
1219  data_out.write_vtk (output);
1220  }
1221  else if(filetype == ".gpl")
1222  {
1223  data_out.write_gnuplot (output);
1224  }
1225  else
1226  {
1227  throw DOpEException("Don't know how to write filetype `" + filetype + "'!","InstatReducedProblem::WriteToFile");
1228  }
1229  data_out.clear();
1230 #else
1231  std::ofstream output(outfile.c_str());
1232  Vector<double> off;
1233  off = v;
1234  for (unsigned int i = 0; i < off.size(); i++)
1235  {
1236  output << off(i) << std::endl;
1237  }
1238 #endif
1239  }
1240  else
1241  {
1242  throw DOpEException("No such DoFHandler `" + dof_type + "'!",
1243  "InstatReducedProblem::WriteToFile");
1244  }
1245 }
1246 
1247 /******************************************************/
1248 
1249 template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER, typename CONTROLINTEGRATOR,
1250  typename INTEGRATOR, typename PROBLEM, typename VECTOR, int dopedim,
1251  int dealdim>
1252 void InstatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER, CONTROLINTEGRATOR, INTEGRATOR,
1253  PROBLEM, VECTOR, dopedim, dealdim>::WriteToFile(const ControlVector<VECTOR> &v,
1254  std::string name,
1255  std::string outfile,
1256  std::string dof_type,
1257  std::string filetype)
1258 {
1259  WriteToFile(v.GetSpacialVector(), name, outfile, dof_type, filetype);
1260 }
1261 
1262 /******************************************************/
1263 
1264  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
1265  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
1266  typename VECTOR, int dopedim, int dealdim>
1267  void
1268  InstatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
1269  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::WriteToFile(
1270  const std::vector<double> &v, std::string outfile)
1271  {
1272  //TODO This should get timedofhandler later on.
1273  const std::vector<double>& t =
1274  this->GetProblem()->GetSpaceTimeHandler()->GetTimes();
1275  std::ofstream out(outfile.c_str());
1276  assert( t.size() == v.size());
1277  assert(out.is_open());
1278 
1279  out << "#Time\tvalue" << std::endl;
1280  for (unsigned int i = 0; i < v.size(); i++)
1281  {
1282  out << t[i] << "\t" << v[i] << std::endl;
1283  }
1284  out.close();
1285  }
1286 
1287 /******************************************************/
1288 
1289 template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
1290  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
1291  typename VECTOR, int dopedim, int dealdim>
1292  template<typename PDE>
1293  void InstatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
1294  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::
1295  ForwardTimeLoop(PDE& problem, StateVector<VECTOR>& sol, std::string outname, bool eval_funcs)
1296  {
1297  VECTOR u_alt;
1298 
1299  unsigned int max_timestep =
1300  problem.GetSpaceTimeHandler()->GetMaxTimePoint();
1301  const std::vector<double> times =
1302  problem.GetSpaceTimeHandler()->GetTimes();
1303  const unsigned int
1304  n_dofs_per_interval =
1305  problem.GetSpaceTimeHandler()->GetTimeDoFHandler().GetLocalNbrOfDoFs();
1306  std::vector<unsigned int> local_to_global(n_dofs_per_interval);
1307  {
1308  TimeIterator it =
1309  problem.GetSpaceTimeHandler()->GetTimeDoFHandler().first_interval();
1310  it.get_time_dof_indices(local_to_global);
1311  problem.SetTime(times[local_to_global[0]], it);
1312  sol.SetTimeDoFNumber(local_to_global[0], it);
1313  }
1314  //u_alt auf initial_values setzen
1315  {
1316  //dazu erstmal gesamt-dof berechnen
1317  const std::vector<unsigned int>& dofs_per_block =
1318  this->GetProblem()->GetSpaceTimeHandler()->GetStateDoFsPerBlock();
1319  unsigned int n_dofs = 0;
1320  unsigned int n_blocks = dofs_per_block.size();
1321  for (unsigned int i = 0; i < n_blocks; i++)
1322  {
1323  n_dofs += dofs_per_block[i];
1324  }
1325  //und dann auf den Helper zuerueckgreifen (wegen Templateisierung)
1326  DOpEHelper::ReSizeVector(n_dofs, dofs_per_block, u_alt);
1327  }
1328 
1329  //Projection der Anfangsdaten
1330  this->GetOutputHandler()->SetIterationNumber(0, "Time");
1331  {
1332  this->GetOutputHandler()->Write("Computing Initial Values:",
1333  4 + this->GetBasePriority());
1334 
1335  auto& initial_problem = problem.GetInitialProblem();
1336  this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
1337 
1338  //TODO: Possibly another solver for the initial value than for the pde...
1339  _build_state_matrix = this->GetNonlinearSolver("state").NonlinearSolve_Initial(
1340  initial_problem, u_alt, true, true);
1341  _build_state_matrix = true;
1342 
1343  this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
1344 
1345  }
1346  sol.GetSpacialVector() = u_alt;
1347  this->GetOutputHandler()->Write(u_alt, outname + this->GetPostIndex(),
1348  problem.GetDoFType());
1349 
1350 
1351  if(eval_funcs)
1352  {//Funktional Auswertung in t_0
1353  ComputeTimeFunctionals(0,
1354  max_timestep);
1355  this->SetProblemType("state");
1356  }
1357 
1358 
1359  for (TimeIterator it =
1360  problem.GetSpaceTimeHandler()->GetTimeDoFHandler().first_interval(); it
1361  != problem.GetSpaceTimeHandler()->GetTimeDoFHandler().after_last_interval(); ++it)
1362  {
1363  it.get_time_dof_indices(local_to_global);
1364  problem.SetTime(times[local_to_global[0]], it);
1365  sol.SetTimeDoFNumber(local_to_global[0], it);
1366  //TODO Eventuell waere ein Test mit nicht-gleichmaessigen Zeitschritten sinnvoll!
1367 
1368  //we start here at i=1 because we assume that the most
1369  //left DoF in the actual interval is already computed!
1370  for (unsigned int i = 1; i < n_dofs_per_interval; i++)
1371  {
1372  this->GetOutputHandler()->SetIterationNumber(local_to_global[i],
1373  "Time");
1374  double time = times[local_to_global[i]];
1375 
1376  std::stringstream out;
1377  this->GetOutputHandler()->InitOut(out);
1378  out << "\t\t Timestep: " << local_to_global[i] << " ("
1379  << times[local_to_global[i - 1]] << " -> " << time
1380  << ") using " << problem.GetName();
1381  problem.GetOutputHandler()->Write(out,
1382  4 + this->GetBasePriority());
1383 
1384  sol.SetTimeDoFNumber(local_to_global[i], it);
1385  sol.GetSpacialVector() = 0;
1386 
1387  this->GetProblem()->AddAuxiliaryToIntegrator(
1388  this->GetIntegrator());
1389 
1390  this->GetNonlinearSolver("state").NonlinearLastTimeEvals(problem,
1391  u_alt, sol.GetSpacialVector());
1392 
1393  this->GetProblem()->DeleteAuxiliaryFromIntegrator(
1394  this->GetIntegrator());
1395 
1396  problem.SetTime(time, it);
1397  this->GetProblem()->AddAuxiliaryToIntegrator(
1398  this->GetIntegrator());
1399 
1400  _build_state_matrix
1401  = this->GetNonlinearSolver("state").NonlinearSolve(problem,
1402  u_alt, sol.GetSpacialVector(), true,
1403  _build_state_matrix);
1404 
1405  this->GetProblem()->DeleteAuxiliaryFromIntegrator(
1406  this->GetIntegrator());
1407 
1408  //TODO do a transfer to the next grid for changing spatial meshes!
1409  u_alt = sol.GetSpacialVector();
1410  this->GetOutputHandler()->Write(sol.GetSpacialVector(),
1411  outname + this->GetPostIndex(), problem.GetDoFType());
1412  if(eval_funcs)
1413  {//Funktional Auswertung in t_n//if abfrage, welcher typ
1414  ComputeTimeFunctionals(local_to_global[i], max_timestep);
1415  this->SetProblemType("state");
1416  }
1417  }
1418  }
1419  }
1420 
1421 /******************************************************/
1422 
1423 template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
1424  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
1425  typename VECTOR, int dopedim, int dealdim>
1426  template<typename PDE>
1427  void InstatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
1428  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::
1429  BackwardTimeLoop(PDE& problem, StateVector<VECTOR>& sol, std::string outname)
1430  {
1431  VECTOR u_alt;
1432 
1433  unsigned int max_timestep =
1434  problem.GetSpaceTimeHandler()->GetMaxTimePoint();
1435  const std::vector<double> times =
1436  problem.GetSpaceTimeHandler()->GetTimes();
1437  const unsigned int
1438  n_dofs_per_interval =
1439  problem.GetSpaceTimeHandler()->GetTimeDoFHandler().GetLocalNbrOfDoFs();
1440  std::vector<unsigned int> local_to_global(n_dofs_per_interval);
1441  {
1442  TimeIterator it =
1443  problem.GetSpaceTimeHandler()->GetTimeDoFHandler().last_interval();
1444  it.get_time_dof_indices(local_to_global);
1445  problem.SetTime(times[local_to_global[local_to_global.size()-1]], it);
1446  sol.SetTimeDoFNumber(local_to_global[local_to_global.size()-1], it);
1447  }
1448  //u_alt auf initial_values setzen
1449  {
1450  //dazu erstmal gesamt-dof berechnen
1451  const std::vector<unsigned int>& dofs_per_block =
1452  this->GetProblem()->GetSpaceTimeHandler()->GetStateDoFsPerBlock();
1453  unsigned int n_dofs = 0;
1454  unsigned int n_blocks = dofs_per_block.size();
1455  for (unsigned int i = 0; i < n_blocks; i++)
1456  {
1457  n_dofs += dofs_per_block[i];
1458  }
1459  //und dann auf den Helper zuerueckgreifen (wegen Templateisierung)
1460  DOpEHelper::ReSizeVector(n_dofs, dofs_per_block, u_alt);
1461  }
1462  //Projection der Anfangsdaten
1463  this->GetOutputHandler()->SetIterationNumber(max_timestep, "Time");
1464  {
1465  this->GetOutputHandler()->Write("Computing Initial Values:",
1466  4 + this->GetBasePriority());
1467 
1468  auto& initial_problem = problem.GetInitialProblem();
1469  this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
1470 
1471  //TODO: Possibly another solver for the initial value than for the pde...
1472  _build_state_matrix = this->GetNonlinearSolver("adjoint").NonlinearSolve_Initial(
1473  initial_problem, u_alt, true, true);
1474  _build_state_matrix = true;
1475 
1476  this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
1477 
1478  }
1479  sol.GetSpacialVector() = u_alt;
1480  this->GetOutputHandler()->Write(u_alt, outname + this->GetPostIndex(),
1481  problem.GetDoFType());
1482 
1483  //TODO: Maybe we should calculate the local gradient computations here
1484 
1485  for (TimeIterator it =
1486  problem.GetSpaceTimeHandler()->GetTimeDoFHandler().last_interval(); it
1487  != problem.GetSpaceTimeHandler()->GetTimeDoFHandler().before_first_interval(); --it)
1488  {
1489  it.get_time_dof_indices(local_to_global);
1490  problem.SetTime(times[local_to_global[local_to_global.size()-1]], it);
1491  sol.SetTimeDoFNumber(local_to_global[local_to_global.size()-1], it);
1492  //TODO Eventuell waere ein Test mit nicht-gleichmaessigen Zeitschritten sinnvoll!
1493 
1494  //we start here at i= 1 and transform i -> n_dofs_per_interval-1-i because we assume that the most
1495  //right DoF in the actual interval is already computed!
1496  for (unsigned int i = 1; i < n_dofs_per_interval; i++)
1497  {
1498  unsigned int j = n_dofs_per_interval-1-i;
1499  this->GetOutputHandler()->SetIterationNumber(local_to_global[j],
1500  "Time");
1501  double time = times[local_to_global[j]];
1502 
1503  std::stringstream out;
1504  this->GetOutputHandler()->InitOut(out);
1505  out << "\t\t Timestep: " << local_to_global[j+1] << " ("
1506  << times[local_to_global[j + 1]] << " -> " << time
1507  << ") using " << problem.GetName();
1508  problem.GetOutputHandler()->Write(out,
1509  4 + this->GetBasePriority());
1510 
1511  sol.SetTimeDoFNumber(local_to_global[j], it);
1512  sol.GetSpacialVector() = 0;
1513 
1514  this->GetProblem()->AddAuxiliaryToIntegrator(
1515  this->GetIntegrator());
1516 
1517  this->GetNonlinearSolver("adjoint").NonlinearLastTimeEvals(problem,
1518  u_alt, sol.GetSpacialVector());
1519 
1520  this->GetProblem()->DeleteAuxiliaryFromIntegrator(
1521  this->GetIntegrator());
1522 
1523  problem.SetTime(time, it);
1524  this->GetProblem()->AddAuxiliaryToIntegrator(
1525  this->GetIntegrator());
1526 
1527  _build_adjoint_matrix
1528  = this->GetNonlinearSolver("adjoint").NonlinearSolve(problem,
1529  u_alt, sol.GetSpacialVector(), true,
1530  _build_adjoint_matrix);
1531 
1532  this->GetProblem()->DeleteAuxiliaryFromIntegrator(
1533  this->GetIntegrator());
1534 
1535  //TODO do a transfer to the next grid for changing spatial meshes!
1536  u_alt = sol.GetSpacialVector();
1537  this->GetOutputHandler()->Write(sol.GetSpacialVector(),
1538  outname + this->GetPostIndex(), problem.GetDoFType());
1539 
1540  //Maybe build local gradient here
1541  }
1542  }
1543  }
1545 }
1546 #endif
instat_step_newtonsolver.h
DOpE::InstatReducedProblem::GetControlNonlinearSolver
CONTROLNONLINEARSOLVER & GetControlNonlinearSolver()
Definition: instatreducedproblem.h:526
DOpE::ControlVector::UnLockCopy
void UnLockCopy() const
Definition: controlvector.h:196
DOpE::ConstraintVector
Definition: constraintvector.h:47
DOpE::InstatReducedProblem::GetDU
StateVector< VECTOR > & GetDU()
Definition: instatreducedproblem.h:294
DOpE::ControlVector::GetSpacialVectorCopy
const dealii::Vector< double > & GetSpacialVectorCopy() const
Definition: controlvector.cc:132
DOpE::InstatReducedProblem::GetControlBoxConstraints
void GetControlBoxConstraints(ControlVector< VECTOR > &lb, ControlVector< VECTOR > &ub)
Definition: instatreducedproblem.h:612
DOpE::DOpEtypes::initial
Definition: dopetypes.h:106
DOpE::InstatReducedProblem::BackwardTimeLoop
void BackwardTimeLoop(PDE &problem, StateVector< VECTOR > &sol, std::string outname)
Definition: instatreducedproblem.h:1429
DOpE::StateVector::GetSpacialVector
VECTOR & GetSpacialVector()
Definition: statevector.cc:329
DOpEHelper::ReSizeVector
void ReSizeVector(unsigned int ndofs, const std::vector< unsigned int > &dofs_per_block, dealii::BlockVector< double > &vector)
Definition: helper.cc:30
DOpE::ParameterReader
Definition: parameterreader.h:36
DOpE::InstatReducedProblem::GetU
StateVector< VECTOR > & GetU()
Definition: instatreducedproblem.h:286
fractional_step_theta_step_newtonsolver.h
DOpE::TimeIterator
Definition: timeiterator.h:63
DOpE::DOpEtypes::stationary
Definition: dopetypes.h:105
DOpE::ControlVector
Definition: controlvector.h:48
DOpE::InstatReducedProblem::ComputeReducedHessianVector
void ComputeReducedHessianVector(const ControlVector< VECTOR > &q, const ControlVector< VECTOR > &direction, ControlVector< VECTOR > &hessian_direction, ControlVector< VECTOR > &hessian_direction_transposed)
Definition: instatreducedproblem.h:872
DOpE::InstatReducedProblem::GetZ
StateVector< VECTOR > & GetZ()
Definition: instatreducedproblem.h:290
DOpE::InstatReducedProblem::GetIntegrator
INTEGRATOR & GetIntegrator()
Definition: instatreducedproblem.h:305
DOpE::ReducedProblemInterface::GetProblem
PROBLEM * GetProblem()
Definition: reducedprobleminterface.h:493
DOpE::SolutionExtractor
Definition: solutionextractor.h:40
DOpE::InstatReducedProblem
Definition: instatreducedproblem.h:81
DOpE::InstatReducedProblem::WriteToFile
void WriteToFile(const VECTOR &v, std::string name, std::string outfile, std::string dof_type, std::string filetype)
Definition: instatreducedproblem.h:1179
DOpE::InstatReducedProblem::~InstatReducedProblem
virtual ~InstatReducedProblem()
Definition: instatreducedproblem.h:493
DOpE::InstatReducedProblem::ComputeReducedFunctionals
void ComputeReducedFunctionals(const ControlVector< VECTOR > &q)
Definition: instatreducedproblem.h:804
DOpE::InstatReducedProblem::StateSizeInfo
void StateSizeInfo(std::stringstream &out)
Definition: instatreducedproblem.h:235
DOpE::InstatReducedProblem::ReInit
void ReInit()
Definition: instatreducedproblem.h:545
DOpE::ControlVector::GetSpaceTimeHandler
const SpaceTimeHandlerBase< VECTOR > * GetSpaceTimeHandler() const
Definition: controlvector.h:211
DOpE::InstatReducedProblem::ForwardTimeLoop
void ForwardTimeLoop(PDE &problem, StateVector< VECTOR > &sol, std::string outname, bool eval_funcs)
Definition: instatreducedproblem.h:1295
DOpE::TimeIterator::get_time_dof_indices
void get_time_dof_indices(std::vector< unsigned int > &local_dof_indices) const
Definition: timeiterator.h:208
DOpE::InstatReducedProblem::GetU
const StateVector< VECTOR > & GetU() const
Definition: instatreducedproblem.h:282
DOpE::StateVector
Definition: statevector.h:49
DOpE::InstatReducedProblem::ComputeReducedAdjoint
void ComputeReducedAdjoint(const ControlVector< VECTOR > &q)
Definition: instatreducedproblem.h:623
DOpE::InstatReducedProblem::ComputeReducedGradient
void ComputeReducedGradient(const ControlVector< VECTOR > &q, ControlVector< VECTOR > &gradient, ControlVector< VECTOR > &gradient_transposed)
Definition: instatreducedproblem.h:649
DOpE::ControlVector::GetSpacialVector
VECTOR & GetSpacialVector()
Definition: controlvector.cc:118
DOpE::InstatReducedProblem::ComputeTimeFunctionals
void ComputeTimeFunctionals(unsigned int step, unsigned int num_steps)
Definition: instatreducedproblem.h:1009
DOpE::InstatReducedProblem::declare_params
static void declare_params(ParameterReader &param_reader)
Definition: instatreducedproblem.h:409
DOpE::InstatReducedProblem::GetNonlinearSolver
NONLINEARSOLVER & GetNonlinearSolver(std::string type)
Definition: instatreducedproblem.h:503
DOpE::InstatReducedProblem::InstatReducedProblem
InstatReducedProblem(PROBLEM *OP, std::string state_behavior, ParameterReader &param_reader, INTEGRATORDATACONT &idc, int base_priority=0)
Definition: instatreducedproblem.h:421
DOpE::ReducedProblemInterface
Definition: reducedprobleminterface.h:338
DOpE::InstatReducedProblem::ComputeReducedCostFunctional
double ComputeReducedCostFunctional(const ControlVector< VECTOR > &q)
Definition: instatreducedproblem.h:779
DOpE::InstatReducedProblem::ComputeReducedConstraints
bool ComputeReducedConstraints(const ControlVector< VECTOR > &q, ConstraintVector< VECTOR > &g)
Definition: instatreducedproblem.h:599
DOpE::StateVector::SetTimeDoFNumber
void SetTimeDoFNumber(unsigned int dof_number, const TimeIterator &interval) const
Definition: statevector.cc:231
DOpE::InstatReducedProblem::GetDZ
StateVector< VECTOR > & GetDZ()
Definition: instatreducedproblem.h:298
DOpE::InstatReducedProblem::ComputeRefinementIndicators
void ComputeRefinementIndicators(DWRC &)
Definition: instatreducedproblem.h:222
DOpE::DOpEException
Definition: dopeexception.h:35
DOpE::InstatReducedProblem::GetControlIntegrator
CONTROLINTEGRATOR & GetControlIntegrator()
Definition: instatreducedproblem.h:309
DOpE::ReducedProblemInterface::ReInit
virtual void ReInit()
Definition: reducedprobleminterface.h:360
DOpE::InstatReducedProblem::ComputeReducedState
void ComputeReducedState(const ControlVector< VECTOR > &q)
Definition: instatreducedproblem.h:574
Generated on Fri Oct 24 2014 17:56:40 for DOpE by   doxygen 1.8.6