DOpE
statreducedproblem.h
Go to the documentation of this file.
1 
24 #ifndef STAT_REDUCED_PROBLEM_H_
25 #define STAT_REDUCED_PROBLEM_H_
26 
27 #include "reducedprobleminterface.h"
28 #include "integrator.h"
29 #include "parameterreader.h"
30 #include "statevector.h"
31 #include "stateproblem.h"
32 
33 #include <deal.II/lac/vector.h>
34 #include <deal.II/lac/block_sparsity_pattern.h>
35 #include <deal.II/lac/block_sparse_matrix.h>
36 
37 #include "optproblemcontainer.h"
38 #include "pdeinterface.h"
39 #include "functionalinterface.h"
40 #include "dirichletdatainterface.h"
41 #include "dopeexception.h"
42 #include "newtonsolver.h"
43 #include "newtonsolvermixeddims.h"
44 #include "cglinearsolver.h"
45 #include "gmreslinearsolver.h"
46 #include "directlinearsolver.h"
47 #include "voidlinearsolver.h"
48 #include "constraintinterface.h"
49 #include "solutionextractor.h"
50 
51 #include <deal.II/base/data_out_base.h>
52 #include <deal.II/numerics/data_out.h>
53 #include <deal.II/numerics/matrix_tools.h>
54 #include <deal.II/numerics/vector_tools.h>
55 #include <deal.II/base/function.h>
56 #include <deal.II/lac/sparse_matrix.h>
57 #include <deal.II/lac/compressed_simple_sparsity_pattern.h>
58 #include <deal.II/lac/block_sparsity_pattern.h>
59 #include <deal.II/lac/sparse_direct.h>
60 
61 #include <fstream>
62 namespace DOpE
63 {
79  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
80  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
81  typename VECTOR, int dopedim, int dealdim>
82  class StatReducedProblem : public ReducedProblemInterface<PROBLEM, VECTOR>
83  {
84  public:
97  template<typename INTEGRATORDATACONT>
98  StatReducedProblem(PROBLEM *OP, DOpEtypes::VectorStorageType state_behavior,
99  ParameterReader &param_reader, INTEGRATORDATACONT& idc,
100  int base_priority = 0);
101 
115  template<typename STATEINTEGRATORDATACONT,
116  typename CONTROLINTEGRATORCONT>
117  StatReducedProblem(PROBLEM *OP, DOpEtypes::VectorStorageType state_behavior,
118  ParameterReader &param_reader, CONTROLINTEGRATORCONT& c_idc,
119  STATEINTEGRATORDATACONT & s_idc, int base_priority = 0);
120 
121  ~StatReducedProblem();
122 
123  /******************************************************/
124 
130  static void
131  declare_params(ParameterReader &param_reader);
132 
133  /******************************************************/
134 
140  void
141  ReInit();
142 
143  /******************************************************/
144 
150  bool
151  ComputeReducedConstraints(const ControlVector<VECTOR>& q,
152  ConstraintVector<VECTOR>& g);
153 
154  /******************************************************/
155 
161  void
162  GetControlBoxConstraints(ControlVector<VECTOR>& lb,
163  ControlVector<VECTOR>& ub);
164 
165  /******************************************************/
166 
172  void
173  ComputeReducedGradient(const ControlVector<VECTOR>& q,
174  ControlVector<VECTOR>& gradient,
175  ControlVector<VECTOR>& gradient_transposed);
176 
177  /******************************************************/
178 
184  double
185  ComputeReducedCostFunctional(const ControlVector<VECTOR>& q);
186 
187  /******************************************************/
188 
194  void
195  ComputeReducedFunctionals(const ControlVector<VECTOR>& q);
196 
197  /******************************************************/
198 
218  template<class DWRC,class PDE>
219  void
220  ComputeRefinementIndicators(const ControlVector<VECTOR>& q,
221  DWRC& dwrc, PDE& pde);
222 
223  /******************************************************/
224 
230  void
231  ComputeReducedHessianVector(const ControlVector<VECTOR>& q,
232  const ControlVector<VECTOR>& direction,
233  ControlVector<VECTOR>& hessian_direction,
234  ControlVector<VECTOR>& hessian_direction_transposed);
235 
236  /******************************************************/
237 
243  void
244  ComputeReducedGradientOfGlobalConstraints(unsigned int num,
245  const ControlVector<VECTOR>& q, const ConstraintVector<VECTOR>& g,
246  ControlVector<VECTOR>& gradient,
247  ControlVector<VECTOR>& gradient_transposed);
248 
249  /******************************************************/
250 
256  void
257  StateSizeInfo(std::stringstream& out)
258  {
259  GetU().PrintInfos(out);
260  }
261 
262  /******************************************************/
263 
274  void
275  WriteToFile(const VECTOR &v, std::string name, std::string outfile,
276  std::string dof_type, std::string filetype);
277 
278  /******************************************************/
279 
288  void
289  WriteToFile(const ControlVector<VECTOR> &v, std::string name, std::string dof_type);
290 
298  void
299  WriteToFile(const std::vector<double> &/*v*/,
300  std::string /*outfile*/)
301  {
302  abort();
303  }
304 
305  protected:
313  void
314  ComputeReducedState(const ControlVector<VECTOR>& q);
315 
316  /******************************************************/
325  void
326  ComputeReducedAdjoint(const ControlVector<VECTOR>& q);
327 
328  /******************************************************/
329 
340  void
341  ComputeDualForErrorEstimation(const ControlVector<VECTOR>& q,
342  DOpEtypes::WeightComputation weight_comp);
343 
344  const StateVector<VECTOR> &
345  GetU() const
346  {
347  return u_;
348  }
349  StateVector<VECTOR> &
350  GetU()
351  {
352  return u_;
353  }
354  StateVector<VECTOR> &
355  GetZ()
356  {
357  return z_;
358  }
359  StateVector<VECTOR> &
360  GetDU()
361  {
362  return du_;
363  }
364  StateVector<VECTOR> &
365  GetDZ()
366  {
367  return dz_;
368  }
372  const StateVector<VECTOR> &
373  GetZForEE() const
374  {
375  return z_for_ee_;
376  }
377  StateVector<VECTOR> &
378  GetZForEE()
379  {
380  return z_for_ee_;
381  }
382 
383  NONLINEARSOLVER&
384  GetNonlinearSolver(std::string type);
385  CONTROLNONLINEARSOLVER&
386  GetControlNonlinearSolver();
387  INTEGRATOR&
388  GetIntegrator()
389  {
390  return integrator_;
391  }
392  CONTROLINTEGRATOR&
393  GetControlIntegrator()
394  {
395  return control_integrator_;
396  }
397 
398  private:
399  StateVector<VECTOR> u_;
400  StateVector<VECTOR> z_;
401  StateVector<VECTOR> du_;
402  StateVector<VECTOR> dz_;
403  StateVector<VECTOR> z_for_ee_;
404 
405  INTEGRATOR integrator_;
406  CONTROLINTEGRATOR control_integrator_;
407  NONLINEARSOLVER nonlinear_state_solver_;
408  NONLINEARSOLVER nonlinear_adjoint_solver_;
409  CONTROLNONLINEARSOLVER nonlinear_gradient_solver_;
410 
411  bool build_state_matrix_, build_adjoint_matrix_, build_control_matrix_;
412  bool state_reinit_, adjoint_reinit_, gradient_reinit_;
413 
414  friend class SolutionExtractor<
415  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
416  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>,
417  VECTOR> ;
418  };
419 
420  /*************************************************************************/
421  /*****************************IMPLEMENTATION******************************/
422  /*************************************************************************/
423  using namespace dealii;
424 
425  /******************************************************/
426  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
427  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
428  typename VECTOR, int dopedim, int dealdim>
429  void
430  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
431  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::declare_params(
432  ParameterReader &param_reader)
433  {
434  NONLINEARSOLVER::declare_params(param_reader);
435  }
436  /******************************************************/
437 
438  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
439  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
440  typename VECTOR, int dopedim, int dealdim>
441  template<typename INTEGRATORDATACONT>
442  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
443  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::StatReducedProblem(
444  PROBLEM *OP, DOpEtypes::VectorStorageType state_behavior,
445  ParameterReader &param_reader, INTEGRATORDATACONT& idc,
446  int base_priority)
447  : ReducedProblemInterface<PROBLEM, VECTOR>(OP,
448  base_priority), u_(OP->GetSpaceTimeHandler(), state_behavior,
449  param_reader), z_(OP->GetSpaceTimeHandler(), state_behavior,
450  param_reader), du_(OP->GetSpaceTimeHandler(), state_behavior,
451  param_reader), dz_(OP->GetSpaceTimeHandler(), state_behavior,
452  param_reader), z_for_ee_(OP->GetSpaceTimeHandler(),
453  state_behavior, param_reader), integrator_(idc), control_integrator_(
454  idc), nonlinear_state_solver_(integrator_, param_reader), nonlinear_adjoint_solver_(
455  integrator_, param_reader), nonlinear_gradient_solver_(
456  control_integrator_, param_reader)
457 
458  {
459  //ReducedProblems should be ReInited
460  {
461  state_reinit_ = true;
462  adjoint_reinit_ = true;
463  gradient_reinit_ = true;
464  }
465  }
466 
467  /******************************************************/
468 
469  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
470  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
471  typename VECTOR, int dopedim, int dealdim>
472  template<typename STATEINTEGRATORDATACONT, typename CONTROLINTEGRATORCONT>
473  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
474  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::StatReducedProblem(
475  PROBLEM *OP, DOpEtypes::VectorStorageType state_behavior,
476  ParameterReader &param_reader, CONTROLINTEGRATORCONT& c_idc,
477  STATEINTEGRATORDATACONT & s_idc, int base_priority)
478  : ReducedProblemInterface<PROBLEM, VECTOR>(OP,
479  base_priority), u_(OP->GetSpaceTimeHandler(), state_behavior,
480  param_reader), z_(OP->GetSpaceTimeHandler(), state_behavior,
481  param_reader), du_(OP->GetSpaceTimeHandler(), state_behavior,
482  param_reader), dz_(OP->GetSpaceTimeHandler(), state_behavior,
483  param_reader), z_for_ee_(OP->GetSpaceTimeHandler(),
484  state_behavior, param_reader), integrator_(s_idc), control_integrator_(
485  c_idc), nonlinear_state_solver_(integrator_, param_reader), nonlinear_adjoint_solver_(
486  integrator_, param_reader), nonlinear_gradient_solver_(
487  control_integrator_, param_reader)
488 
489  {
490  //ReducedProblems should be ReInited
491  {
492  state_reinit_ = true;
493  adjoint_reinit_ = true;
494  gradient_reinit_ = true;
495  }
496  }
497 
498  /******************************************************/
499 
500  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
501  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
502  typename VECTOR, int dopedim, int dealdim>
503  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
504  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::~StatReducedProblem()
505  {
506  }
507 
508  /******************************************************/
509 
510  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
511  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
512  typename VECTOR, int dopedim, int dealdim>
513  NONLINEARSOLVER&
514  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
515  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::GetNonlinearSolver(
516  std::string type)
517  {
518  if ((type == "state") || (type == "tangent"))
519  {
520  return nonlinear_state_solver_;
521  }
522  else if ((type == "adjoint") || (type == "adjoint_hessian")
523  || (type == "adjoint_for_ee"))
524  {
525  return nonlinear_adjoint_solver_;
526  }
527  else
528  {
529  throw DOpEException("No Solver for Problem type:`" + type + "' found",
530  "StatReducedProblem::GetNonlinearSolver");
531 
532  }
533  }
534  /******************************************************/
535 
536  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
537  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
538  typename VECTOR, int dopedim, int dealdim>
539  CONTROLNONLINEARSOLVER&
540  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
541  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::GetControlNonlinearSolver()
542  {
543  if ((this->GetProblem()->GetType() == "gradient")
544  || (this->GetProblem()->GetType() == "hessian"))
545  {
546  return nonlinear_gradient_solver_;
547  }
548  else
549  {
550  throw DOpEException(
551  "No Solver for Problem type:`" + this->GetProblem()->GetType()
552  + "' found", "StatReducedProblem::GetControlNonlinearSolver");
553 
554  }
555  }
556  /******************************************************/
557 
558  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
559  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
560  typename VECTOR, int dopedim, int dealdim>
561  void
562  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
563  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ReInit()
564  {
565  ReducedProblemInterface<PROBLEM, VECTOR>::ReInit();
566 
567  //Some Solvers must be reinited when called
568  // Better have subproblems, so that solver can be reinited here
569  {
570  state_reinit_ = true;
571  adjoint_reinit_ = true;
572  gradient_reinit_ = true;
573  }
574 
575  build_state_matrix_ = true;
576  build_adjoint_matrix_ = true;
577 
578  GetU().ReInit();
579  GetZ().ReInit();
580  GetDU().ReInit();
581  GetDZ().ReInit();
582  GetZForEE().ReInit();
583 
584  build_control_matrix_ = true;
585  }
586 
587  /******************************************************/
588 
589  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
590  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
591  typename VECTOR, int dopedim, int dealdim>
592  void
593  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
594  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ComputeReducedState(
595  const ControlVector<VECTOR>& q)
596  {
597  this->InitializeFunctionalValues(
598  this->GetProblem()->GetNFunctionals() + 1);
599 
600  this->GetOutputHandler()->Write("Computing State Solution:",
601  4 + this->GetBasePriority());
602 
603  this->SetProblemType("state");
604  auto& problem = this->GetProblem()->GetStateProblem();
605  if (state_reinit_ == true)
606  {
607  GetNonlinearSolver("state").ReInit(problem);
608  state_reinit_ = false;
609  }
610 
611  this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
612  if (dopedim == dealdim)
613  {
614  this->GetIntegrator().AddDomainData("control", &(q.GetSpacialVector()));
615  }
616  else if (dopedim == 0)
617  {
618  this->GetIntegrator().AddParamData("control",
619  &(q.GetSpacialVectorCopy()));
620  }
621  else
622  {
623  throw DOpEException("dopedim not implemented",
624  "StatReducedProblem::ComputeReducedState");
625  }
626 
627  build_state_matrix_ = this->GetNonlinearSolver("state").NonlinearSolve(
628  problem, (GetU().GetSpacialVector()), true, build_state_matrix_);
629 
630  if (dopedim == dealdim)
631  {
632  this->GetIntegrator().DeleteDomainData("control");
633  }
634  else if (dopedim == 0)
635  {
636  this->GetIntegrator().DeleteParamData("control");
637  q.UnLockCopy();
638  }
639  else
640  {
641  throw DOpEException("dopedim not implemented",
642  "StatReducedProblem::ComputeReducedState");
643  }
644  this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
645 
646  this->GetOutputHandler()->Write((GetU().GetSpacialVector()),
647  "State" + this->GetPostIndex(), problem.GetDoFType());
648 
649  }
650  /******************************************************/
651 
652  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
653  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
654  typename VECTOR, int dopedim, int dealdim>
655  bool
656  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
657  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ComputeReducedConstraints(
658  const ControlVector<VECTOR>& q, ConstraintVector<VECTOR>& g)
659  {
660  this->GetOutputHandler()->Write("Evaluating Constraints:",
661  4 + this->GetBasePriority());
662 
663  this->SetProblemType("constraints");
664 
665  g = 0;
666  //Local constraints
667  // this->GetProblem()->ComputeLocalConstraints(q.GetSpacialVector(), GetU().GetSpacialVector(),
668  // g.GetSpacialVector("local"));
669  if (dopedim == dealdim)
670  {
671  this->GetControlIntegrator().AddDomainData("control",
672  &(q.GetSpacialVector()));
673  }
674  else if (dopedim == 0)
675  {
676  this->GetControlIntegrator().AddParamData("control",
677  &(q.GetSpacialVectorCopy()));
678  }
679  else
680  {
681  throw DOpEException("dopedim not implemented",
682  "StatReducedProblem::ComputeReducedConstraints");
683  }
684  this->GetControlIntegrator().ComputeLocalControlConstraints(
685  *(this->GetProblem()), g.GetSpacialVector("local"));
686  if (dopedim == dealdim)
687  {
688  this->GetControlIntegrator().DeleteDomainData("control");
689  }
690  else if (dopedim == 0)
691  {
692  this->GetControlIntegrator().DeleteParamData("control");
693  q.UnLockCopy();
694  }
695  else
696  {
697  throw DOpEException("dopedim not implemented",
698  "StatReducedProblem::ComputeReducedConstraints");
699  }
700  //Global in Space-Time Constraints
701  dealii::Vector<double>& gc = g.GetGlobalConstraints();
702  //dealii::Vector<double> global_values(gc.size());
703 
704  unsigned int nglobal = gc.size(); //global_values.size();
705 
706  if (nglobal > 0)
707  {
708  this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
709 
710  for (unsigned int i = 0; i < nglobal; i++)
711  {
712  //this->SetProblemType("local_global_constraints", i);
713  this->SetProblemType("global_constraints", i);
714  this->GetIntegrator().AddDomainData("state",
715  &(GetU().GetSpacialVector()));
716  if (dopedim == dealdim)
717  {
718  this->GetIntegrator().AddDomainData("control",
719  &(q.GetSpacialVector()));
720  }
721  else if (dopedim == 0)
722  {
723  this->GetIntegrator().AddParamData("control",
724  &(q.GetSpacialVectorCopy()));
725  }
726  else
727  {
728  throw DOpEException("dopedim not implemented",
729  "StatReducedProblem::ComputeReducedConstraints");
730  }
731 
732  double ret = 0;
733  bool found = false;
734 
735  if (this->GetProblem()->GetConstraintType().find("domain")
736  != std::string::npos)
737  {
738  found = true;
739  ret += this->GetIntegrator().ComputeDomainScalar(
740  *(this->GetProblem()));
741  }
742  if (this->GetProblem()->GetConstraintType().find("point")
743  != std::string::npos)
744  {
745  found = true;
746  ret += this->GetIntegrator().ComputePointScalar(
747  *(this->GetProblem()));
748  }
749  if (this->GetProblem()->GetConstraintType().find("boundary")
750  != std::string::npos)
751  {
752  found = true;
753  ret += this->GetIntegrator().ComputeBoundaryScalar(
754  *(this->GetProblem()));
755  }
756  if (this->GetProblem()->GetConstraintType().find("face")
757  != std::string::npos)
758  {
759  found = true;
760  ret += this->GetIntegrator().ComputeFaceScalar(
761  *(this->GetProblem()));
762  }
763 
764  if (!found)
765  {
766  throw DOpEException(
767  "Unknown Constraint Type: "
768  + this->GetProblem()->GetConstraintType(),
769  "StatReducedProblem::ComputeReducedConstraints");
770  }
771  // global_values(i) = ret;
772  gc(i) = ret;
773 
774  if (dopedim == dealdim)
775  {
776  this->GetIntegrator().DeleteDomainData("control");
777  }
778  else if (dopedim == 0)
779  {
780  this->GetIntegrator().DeleteParamData("control");
781  q.UnLockCopy();
782  }
783  else
784  {
785  throw DOpEException("dopedim not implemented",
786  "StatReducedProblem::ComputeReducedConstraints");
787  }
788  this->GetIntegrator().DeleteDomainData("state");
789  }
790 
791  this->GetProblem()->DeleteAuxiliaryFromIntegrator(
792  this->GetIntegrator());
793  //gc = global_values;
794  }
795 
796  //Check that no global in space, local in time constraints are given!
797  if (g.HasType("local_global_control") || g.HasType("local_global_state"))
798  {
799  throw DOpEException(
800  "There are global in space, local in time constraints given. In the stationary case they should be moved to global in space and time!",
801  "StatReducedProblem::ComputeReducedConstraints");
802  }
803 
804  //this->GetProblem()->PostProcessConstraints(g, true);
805  this->GetProblem()->PostProcessConstraints(g);
806 
807 
808  return g.IsFeasible();
809  }
810 
811  /******************************************************/
812 
813  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
814  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
815  typename VECTOR, int dopedim, int dealdim>
816  void
817  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
818  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::GetControlBoxConstraints(
819  ControlVector<VECTOR>& lb, ControlVector<VECTOR>& ub)
820  {
821  this->GetProblem()->GetControlBoxConstraints(lb.GetSpacialVector(),
822  ub.GetSpacialVector());
823  }
824 
825  /******************************************************/
826 
827  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
828  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
829  typename VECTOR, int dopedim, int dealdim>
830  void
831  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
832  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ComputeReducedAdjoint(
833  const ControlVector<VECTOR>& q)
834  {
835  this->GetOutputHandler()->Write("Computing Reduced Adjoint:",
836  4 + this->GetBasePriority());
837 
838  this->SetProblemType("adjoint");
839  if (adjoint_reinit_ == true)
840  {
841  GetNonlinearSolver("adjoint").ReInit(*(this->GetProblem()));
842  adjoint_reinit_ = false;
843  }
844 
845  this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
846 
847  this->GetIntegrator().AddDomainData("state",
848  &(GetU().GetSpacialVector())); //&(GetU().GetSpacialVector())
849 
850  if (dopedim == dealdim)
851  {
852  this->GetIntegrator().AddDomainData("control", &(q.GetSpacialVector()));
853  }
854  else if (dopedim == 0)
855  {
856  this->GetIntegrator().AddParamData("control",
857  &(q.GetSpacialVectorCopy()));
858  }
859  else
860  {
861  throw DOpEException("dopedim not implemented",
862  "StatReducedProblem::ComputeReducedAdjoint");
863  }
864 
865  build_adjoint_matrix_ =
866  this->GetNonlinearSolver("adjoint").NonlinearSolve(
867  *(this->GetProblem()), (GetZ().GetSpacialVector()), true,
868  build_adjoint_matrix_);
869 
870  if (dopedim == dealdim)
871  {
872  this->GetIntegrator().DeleteDomainData("control");
873  }
874  else if (dopedim == 0)
875  {
876  this->GetIntegrator().DeleteParamData("control");
877  q.UnLockCopy();
878  }
879  else
880  {
881  throw DOpEException("dopedim not implemented",
882  "StatReducedProblem::ComputeReducedAdjoint");
883  }
884  this->GetIntegrator().DeleteDomainData("state");
885 
886  this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
887 
888  this->GetOutputHandler()->Write((GetZ().GetSpacialVector()),
889  "Adjoint" + this->GetPostIndex(), this->GetProblem()->GetDoFType());
890  }
891 
892  /******************************************************/
893 
894  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
895  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
896  typename VECTOR, int dopedim, int dealdim>
897  void
898  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
899  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ComputeDualForErrorEstimation(
900  const ControlVector<VECTOR>& q,
901  DOpEtypes::WeightComputation weight_comp)
902  {
903  this->GetOutputHandler()->Write("Computing Dual for Error Estimation:",
904  4 + this->GetBasePriority());
905 
906  if (weight_comp == DOpEtypes::higher_order_interpolation)
907  {
908  this->SetProblemType("adjoint_for_ee");
909  }
910  else
911  {
912  throw DOpEException("Unknown WeightComputation",
913  "StatPDEProblem::ComputeDualForErrorEstimation");
914  }
915 
916  // auto& problem = this->GetProblem()->GetStateProblem();//Hier ist adjoint problem einzufuegen
917  auto& problem = *(this->GetProblem());
918  if (adjoint_reinit_ == true)
919  {
920  GetNonlinearSolver("adjoint_for_ee").ReInit(problem);
921  adjoint_reinit_ = false;
922  }
923 
924  this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
925 
926  this->GetIntegrator().AddDomainData("state",
927  &(GetU().GetSpacialVector())); //&(GetU().GetSpacialVector())
928 
929  if (dopedim == dealdim)
930  {
931  this->GetIntegrator().AddDomainData("control", &(q.GetSpacialVector()));
932  }
933  else if (dopedim == 0)
934  {
935  this->GetIntegrator().AddParamData("control",
936  &(q.GetSpacialVectorCopy()));
937  }
938  else
939  {
940  throw DOpEException("dopedim not implemented",
941  "StatReducedProblem::ComputeReducedAdjoint");
942  }
943 
944  build_adjoint_matrix_ =
945  this->GetNonlinearSolver("adjoint_for_ee").NonlinearSolve(problem,
946  (GetZForEE().GetSpacialVector()), true, build_adjoint_matrix_);
947 
948  if (dopedim == dealdim)
949  {
950  this->GetIntegrator().DeleteDomainData("control");
951  }
952  else if (dopedim == 0)
953  {
954  this->GetIntegrator().DeleteParamData("control");
955  q.UnLockCopy();
956  }
957  else
958  {
959  throw DOpEException("dopedim not implemented",
960  "StatReducedProblem::ComputeReducedAdjoint");
961  }
962 
963  this->GetIntegrator().DeleteDomainData("state");
964 
965  this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
966 
967  this->GetOutputHandler()->Write((GetZForEE().GetSpacialVector()),
968  "Adjoint_for_ee" + this->GetPostIndex(), problem.GetDoFType());
969 
970  }
971 
972  /******************************************************/
973 
974  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
975  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
976  typename VECTOR, int dopedim, int dealdim>
977  void
978  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
979  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ComputeReducedGradient(
980  const ControlVector<VECTOR>& q, ControlVector<VECTOR>& gradient,
981  ControlVector<VECTOR>& gradient_transposed)
982  {
983  this->ComputeReducedAdjoint(q);
984 
985  this->GetOutputHandler()->Write("Computing Reduced Gradient:",
986  4 + this->GetBasePriority());
987 
988  //Preparations for ControlInTheDirichletData
989  VECTOR tmp;
990  if (this->GetProblem()->HasControlInDirichletData())
991  {
992  tmp.reinit(GetU().GetSpacialVector());
993  this->SetProblemType("adjoint");
994  this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
995 
996  if (dopedim == dealdim)
997  {
998  this->GetIntegrator().AddDomainData("control",
999  &(q.GetSpacialVector()));
1000  }
1001  else if (dopedim == 0)
1002  {
1003  this->GetIntegrator().AddParamData("control",
1004  &(q.GetSpacialVectorCopy()));
1005  }
1006  else
1007  {
1008  throw DOpEException("dopedim not implemented",
1009  "StatReducedProblem::ComputeReducedGradient");
1010  }
1011 
1012  this->GetIntegrator().AddDomainData("state",
1013  &(GetU().GetSpacialVector()));
1014  this->GetIntegrator().AddDomainData("last_newton_solution",
1015  &(GetZ().GetSpacialVector()));
1016 
1017  this->GetIntegrator().ComputeNonlinearResidual(*(this->GetProblem()),
1018  tmp, false);
1019  tmp *= -1.;
1020 
1021  if (dopedim == dealdim)
1022  {
1023  this->GetIntegrator().DeleteDomainData("control");
1024  }
1025  else if (dopedim == 0)
1026  {
1027  this->GetIntegrator().DeleteParamData("control");
1028  q.UnLockCopy();
1029  }
1030  else
1031  {
1032  throw DOpEException("dopedim not implemented",
1033  "StatReducedProblem::ComputeReducedGradient");
1034  }
1035 
1036  this->GetIntegrator().DeleteDomainData("state");
1037  this->GetIntegrator().DeleteDomainData("last_newton_solution");
1038  this->GetProblem()->DeleteAuxiliaryFromIntegrator(
1039  this->GetIntegrator());
1040  }
1041  //Endof Dirichletdata Preparations
1042  this->SetProblemType("gradient");
1043  if (gradient_reinit_ == true)
1044  {
1045  GetControlNonlinearSolver().ReInit(*(this->GetProblem()));
1046  gradient_reinit_ = false;
1047  }
1048 
1049  this->GetProblem()->AddAuxiliaryToIntegrator(
1050  this->GetControlIntegrator());
1051 
1052  if (dopedim == dealdim)
1053  {
1054  this->GetControlIntegrator().AddDomainData("control",
1055  &(q.GetSpacialVector()));
1056  }
1057  else if (dopedim == 0)
1058  {
1059  this->GetControlIntegrator().AddParamData("control",
1060  &(q.GetSpacialVectorCopy()));
1061  }
1062  else
1063  {
1064  throw DOpEException("dopedim not implemented",
1065  "StatReducedProblem::ComputeReducedGradient");
1066  }
1067 
1068  this->GetControlIntegrator().AddDomainData("state",
1069  &(GetU().GetSpacialVector()));
1070  this->GetControlIntegrator().AddDomainData("adjoint",
1071  &(GetZ().GetSpacialVector()));
1072  if (this->GetProblem()->HasControlInDirichletData())
1073  this->GetControlIntegrator().AddDomainData("adjoint_residual", &tmp);
1074 
1075  gradient_transposed = 0.;
1076  if (dopedim == dealdim)
1077  {
1078  this->GetControlIntegrator().AddDomainData("last_newton_solution",
1079  &(gradient_transposed.GetSpacialVector()));
1080  this->GetControlIntegrator().ComputeNonlinearResidual(
1081  *(this->GetProblem()), gradient.GetSpacialVector(), true);
1082  this->GetControlIntegrator().DeleteDomainData("last_newton_solution");
1083  }
1084  else if (dopedim == 0)
1085  {
1086  this->GetControlIntegrator().AddParamData("last_newton_solution",
1087  &(gradient_transposed.GetSpacialVectorCopy()));
1088  this->GetControlIntegrator().ComputeNonlinearResidual(
1089  *(this->GetProblem()), gradient.GetSpacialVector(), true);
1090 
1091  this->GetControlIntegrator().DeleteParamData("last_newton_solution");
1092  gradient_transposed.UnLockCopy();
1093 
1094  }
1095 
1096  gradient *= -1.;
1097  gradient_transposed = gradient;
1098 
1099  //Compute l^2 representation of the Gradient
1100 
1101  build_control_matrix_ = this->GetControlNonlinearSolver().NonlinearSolve(
1102  *(this->GetProblem()), gradient_transposed.GetSpacialVector(), true,
1103  build_control_matrix_);
1104  if (dopedim == dealdim)
1105  {
1106  this->GetControlIntegrator().DeleteDomainData("control");
1107  }
1108  else if (dopedim == 0)
1109  {
1110  this->GetControlIntegrator().DeleteParamData("control");
1111  q.UnLockCopy();
1112  }
1113  else
1114  {
1115  throw DOpEException("dopedim not implemented",
1116  "StatReducedProblem::ComputeReducedGradient");
1117  }
1118  this->GetControlIntegrator().DeleteDomainData("state");
1119  this->GetControlIntegrator().DeleteDomainData("adjoint");
1120  if (this->GetProblem()->HasControlInDirichletData())
1121  this->GetControlIntegrator().DeleteDomainData("adjoint_residual");
1122 
1123  this->GetProblem()->DeleteAuxiliaryFromIntegrator(
1124  this->GetControlIntegrator());
1125 
1126  this->GetOutputHandler()->Write(gradient,
1127  "Gradient" + this->GetPostIndex(), this->GetProblem()->GetDoFType());
1128  this->GetOutputHandler()->Write(gradient_transposed,
1129  "Gradient_Transposed" + this->GetPostIndex(),
1130  this->GetProblem()->GetDoFType());
1131 
1132  }
1133 
1134  /******************************************************/
1135 
1136  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
1137  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
1138  typename VECTOR, int dopedim, int dealdim>
1139  double
1140  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
1141  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ComputeReducedCostFunctional(
1142  const ControlVector<VECTOR>& q)
1143  {
1144  this->ComputeReducedState(q);
1145 
1146  double ret = 0;
1147  bool found = false;
1148 
1149  this->GetOutputHandler()->Write("Computing Cost Functional:",
1150  4 + this->GetBasePriority());
1151 
1152  this->SetProblemType("cost_functional");
1153 
1154  this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
1155 
1156  if (dopedim == dealdim)
1157  {
1158  this->GetIntegrator().AddDomainData("control", &(q.GetSpacialVector()));
1159  }
1160  else if (dopedim == 0)
1161  {
1162  this->GetIntegrator().AddParamData("control",
1163  &(q.GetSpacialVectorCopy()));
1164  }
1165  else
1166  {
1167  throw DOpEException("dopedim not implemented",
1168  "StatReducedProblem::ComputeReducedCostFunctional");
1169  }
1170  this->GetIntegrator().AddDomainData("state",
1171  &(GetU().GetSpacialVector()));
1172 
1173  if (this->GetProblem()->GetFunctionalType().find("domain")
1174  != std::string::npos)
1175  {
1176  found = true;
1177  ret += this->GetIntegrator().ComputeDomainScalar(*(this->GetProblem()));
1178  }
1179  if (this->GetProblem()->GetFunctionalType().find("point")
1180  != std::string::npos)
1181  {
1182  found = true;
1183  ret += this->GetIntegrator().ComputePointScalar(*(this->GetProblem()));
1184  }
1185  if (this->GetProblem()->GetFunctionalType().find("boundary")
1186  != std::string::npos)
1187  {
1188  found = true;
1189  ret += this->GetIntegrator().ComputeBoundaryScalar(
1190  *(this->GetProblem()));
1191  }
1192  if (this->GetProblem()->GetFunctionalType().find("face")
1193  != std::string::npos)
1194  {
1195  found = true;
1196  ret += this->GetIntegrator().ComputeFaceScalar(*(this->GetProblem()));
1197  }
1198 
1199  if (!found)
1200  {
1201  throw DOpEException(
1202  "Unknown Functional Type: "
1203  + this->GetProblem()->GetFunctionalType(),
1204  "StatReducedProblem::ComputeReducedCostFunctional");
1205  }
1206 
1207  if (dopedim == dealdim)
1208  {
1209  this->GetIntegrator().DeleteDomainData("control");
1210  }
1211  else if (dopedim == 0)
1212  {
1213  this->GetIntegrator().DeleteParamData("control");
1214  q.UnLockCopy();
1215  }
1216  else
1217  {
1218  throw DOpEException("dopedim not implemented",
1219  "StatReducedProblem::ComputeReducedCostFunctional");
1220  }
1221  this->GetIntegrator().DeleteDomainData("state");
1222  this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
1223 
1224  this->GetFunctionalValues()[0].push_back(ret);
1225  return ret;
1226  }
1227 
1228  /******************************************************/
1229 
1230  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
1231  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
1232  typename VECTOR, int dopedim, int dealdim>
1233  void
1234  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
1235  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ComputeReducedFunctionals(
1236  const ControlVector<VECTOR>& q)
1237  {
1238  this->GetOutputHandler()->Write("Computing Functionals:",
1239  4 + this->GetBasePriority());
1240 
1241  this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
1242 
1243  if (dopedim == dealdim)
1244  {
1245  this->GetIntegrator().AddDomainData("control", &(q.GetSpacialVector()));
1246  }
1247  else if (dopedim == 0)
1248  {
1249  this->GetIntegrator().AddParamData("control",
1250  &(q.GetSpacialVectorCopy()));
1251  }
1252  else
1253  {
1254  throw DOpEException("dopedim not implemented",
1255  "StatReducedProblem::ComputeReducedFunctionals");
1256  }
1257  this->GetIntegrator().AddDomainData("state",
1258  &(GetU().GetSpacialVector()));
1259 
1260  for (unsigned int i = 0; i < this->GetProblem()->GetNFunctionals(); i++)
1261  {
1262  double ret = 0;
1263  bool found = false;
1264 
1265  this->SetProblemType("aux_functional", i);
1266  if (this->GetProblem()->GetFunctionalType().find("domain")
1267  != std::string::npos)
1268  {
1269  found = true;
1270  ret += this->GetIntegrator().ComputeDomainScalar(
1271  *(this->GetProblem()));
1272  }
1273  if (this->GetProblem()->GetFunctionalType().find("point")
1274  != std::string::npos)
1275  {
1276  found = true;
1277  ret += this->GetIntegrator().ComputePointScalar(
1278  *(this->GetProblem()));
1279  }
1280  if (this->GetProblem()->GetFunctionalType().find("boundary")
1281  != std::string::npos)
1282  {
1283  found = true;
1284  ret += this->GetIntegrator().ComputeBoundaryScalar(
1285  *(this->GetProblem()));
1286  }
1287  if (this->GetProblem()->GetFunctionalType().find("face")
1288  != std::string::npos)
1289  {
1290  found = true;
1291  ret += this->GetIntegrator().ComputeFaceScalar(*(this->GetProblem()));
1292  }
1293 
1294  if (!found)
1295  {
1296  throw DOpEException(
1297  "Unknown Functional Type: "
1298  + this->GetProblem()->GetFunctionalType(),
1299  "StatReducedProblem::ComputeReducedFunctionals");
1300  }
1301  this->GetFunctionalValues()[i + 1].push_back(ret);
1302  std::stringstream out;
1303  this->GetOutputHandler()->InitOut(out);
1304  out << this->GetProblem()->GetFunctionalName() << ": " << ret;
1305  this->GetOutputHandler()->Write(out, 2 + this->GetBasePriority());
1306  }
1307 
1308  if (dopedim == dealdim)
1309  {
1310  this->GetIntegrator().DeleteDomainData("control");
1311  }
1312  else if (dopedim == 0)
1313  {
1314  this->GetIntegrator().DeleteParamData("control");
1315  q.UnLockCopy();
1316  }
1317  else
1318  {
1319  throw DOpEException("dopedim not implemented",
1320  "StatReducedProblem::ComputeReducedFunctionals");
1321  }
1322  this->GetIntegrator().DeleteDomainData("state");
1323  this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
1324 
1325  }
1326 
1327  /******************************************************/
1328 
1329  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
1330  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
1331  typename VECTOR, int dopedim, int dealdim>
1332  template<class DWRC,class PDE>
1333  void
1334  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
1335  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ComputeRefinementIndicators(
1336  const ControlVector<VECTOR>& q, DWRC& dwrc, PDE& pde)
1337 // StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
1338 // CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ComputeRefinementIndicators(
1339 // const ControlVector<VECTOR>& q, DWRDataContainerBase<VECTOR>& dwrc)
1340  {
1341  //Attach the ResidualModifier to the PDE.
1342  pde.ResidualModifier = boost::bind<void>(boost::mem_fn(&DWRC::ResidualModifier),boost::ref(dwrc),_1);
1343  pde.VectorResidualModifier = boost::bind<void>(boost::mem_fn(&DWRC::VectorResidualModifier),boost::ref(dwrc),_1);
1344 
1345  //first we reinit the dwrdatacontainer (this
1346  //sets the weight-vectors to their correct length)
1347  const unsigned int n_elements =
1348  this->GetProblem()->GetSpaceTimeHandler()->GetStateDoFHandler().get_tria().n_active_cells();
1349  dwrc.ReInit(n_elements);
1350 
1351  //Estimation for Costfunctional or if no dual is needed
1352  if(this->GetProblem()->EEFunctionalIsCost() || !dwrc.NeedDual())
1353  {
1354  this->GetOutputHandler()->Write("Computing Error Indicators:",
1355  4 + this->GetBasePriority());
1356  this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
1357 
1358  //add the primal and dual solution to the integrator
1359  this->GetIntegrator().AddDomainData("state",
1360  &(GetU().GetSpacialVector()));
1361  this->GetIntegrator().AddDomainData("adjoint_for_ee",
1362  &(GetZ().GetSpacialVector()));
1363 
1364  if (dopedim == dealdim)
1365  {
1366  this->GetIntegrator().AddDomainData("control",
1367  &(q.GetSpacialVector()));
1368  }
1369  else if (dopedim == 0)
1370  {
1371  this->GetIntegrator().AddParamData("control",
1372  &(q.GetSpacialVectorCopy()));
1373  }
1374  else
1375  {
1376  throw DOpEException("dopedim not implemented",
1377  "StatReducedProblem::ComputeRefinementIndicators");
1378  }
1379 
1380  this->SetProblemType("error_evaluation");
1381 
1382  //prepare the weights...
1383  dwrc.PrepareWeights(GetU(), GetZ());
1384  dwrc.PrepareWeights(q);
1385 
1386  //now we finally compute the refinement indicators
1387  this->GetIntegrator().ComputeRefinementIndicators(*this->GetProblem(),
1388  dwrc);
1389  // release the lock on the refinement indicators (see dwrcontainer.h)
1390  dwrc.ReleaseLock();
1391  dwrc.ClearWeightData();
1392 
1393  // clear the data
1394  if (dopedim == dealdim)
1395  {
1396  this->GetIntegrator().DeleteDomainData("control");
1397  }
1398  else if (dopedim == 0)
1399  {
1400  this->GetIntegrator().DeleteParamData("control");
1401  q.UnLockCopy();
1402  }
1403  else
1404  {
1405  throw DOpEException("dopedim not implemented",
1406  "StatReducedProblem::ComputeRefinementIndicators");
1407  }
1408  this->GetIntegrator().DeleteDomainData("state");
1409  this->GetIntegrator().DeleteDomainData("adjoint_for_ee");
1410  this->GetProblem()->DeleteAuxiliaryFromIntegrator(
1411  this->GetIntegrator());
1412  }
1413  else //Estimation for other (not the cost) functional
1414  {
1415  throw DOpEException("Estimating the error in other functionals than cost is not implemented",
1416  "StatReducedProblem::ComputeRefinementIndicators");
1417 
1418  }
1419 
1420  std::stringstream out;
1421  this->GetOutputHandler()->InitOut(out);
1422  out << "Error estimate using "<<dwrc.GetName();
1423  if(dwrc.NeedDual())
1424  out<<" For the computation of "<<this->GetProblem()->GetFunctionalName();
1425  out<< ": "<< dwrc.GetError();
1426  this->GetOutputHandler()->Write(out, 2 + this->GetBasePriority());
1427  }
1428 
1429  /******************************************************/
1430 
1431  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
1432  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
1433  typename VECTOR, int dopedim, int dealdim>
1434  void
1435  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
1436  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ComputeReducedHessianVector(
1437  const ControlVector<VECTOR>& q, const ControlVector<VECTOR>& direction,
1438  ControlVector<VECTOR>& hessian_direction,
1439  ControlVector<VECTOR>& hessian_direction_transposed)
1440  {
1441  this->GetOutputHandler()->Write("Computing ReducedHessianVector:",
1442  4 + this->GetBasePriority());
1443  this->GetOutputHandler()->Write("\tSolving Tangent:",
1444  5 + this->GetBasePriority());
1445 
1446  this->SetProblemType("tangent");
1447 
1448  this->GetProblem()->AddAuxiliaryToIntegrator(this->GetIntegrator());
1449 
1450  this->GetIntegrator().AddDomainData("state",
1451  &(GetU().GetSpacialVector()));
1452  this->GetControlIntegrator().AddDomainData("state",
1453  &(GetU().GetSpacialVector()));
1454 
1455  if (dopedim == dealdim)
1456  {
1457  this->GetIntegrator().AddDomainData("dq",
1458  &(direction.GetSpacialVector()));
1459  this->GetIntegrator().AddDomainData("control", &(q.GetSpacialVector()));
1460  }
1461  else if (dopedim == 0)
1462  {
1463  this->GetIntegrator().AddParamData("dq",
1464  &(direction.GetSpacialVectorCopy()));
1465  this->GetIntegrator().AddParamData("control",
1466  &(q.GetSpacialVectorCopy()));
1467  }
1468  else
1469  {
1470  throw DOpEException("dopedim not implemented",
1471  "StatReducedProblem::ComputeReducedHessianVector");
1472  }
1473 
1474  //tangent Matrix is the same as state matrix
1475  build_state_matrix_ = this->GetNonlinearSolver("tangent").NonlinearSolve(
1476  *(this->GetProblem()), (GetDU().GetSpacialVector()), true,
1477  build_state_matrix_);
1478 
1479  this->GetOutputHandler()->Write((GetDU().GetSpacialVector()),
1480  "Tangent" + this->GetPostIndex(), this->GetProblem()->GetDoFType());
1481 
1482  this->GetOutputHandler()->Write("\tSolving Adjoint Hessian:",
1483  5 + this->GetBasePriority());
1484  this->SetProblemType("adjoint_hessian");
1485  this->GetIntegrator().AddDomainData("adjoint",
1486  &(GetZ().GetSpacialVector()));
1487  this->GetIntegrator().AddDomainData("tangent",
1488  &(GetDU().GetSpacialVector()));
1489  this->GetControlIntegrator().AddDomainData("adjoint",
1490  &(GetZ().GetSpacialVector()));
1491  this->GetControlIntegrator().AddDomainData("tangent",
1492  &(GetDU().GetSpacialVector()));
1493 
1494  //adjoint_hessian Matrix is the same as adjoint matrix
1495  build_adjoint_matrix_ =
1496  this->GetNonlinearSolver("adjoint_hessian").NonlinearSolve(
1497  *(this->GetProblem()), (GetDZ().GetSpacialVector()), true,
1498  build_adjoint_matrix_);
1499 
1500  this->GetOutputHandler()->Write((GetDZ().GetSpacialVector()),
1501  "Hessian" + this->GetPostIndex(), this->GetProblem()->GetDoFType());
1502 
1503  this->GetIntegrator().AddDomainData("adjoint_hessian",
1504  &(GetDZ().GetSpacialVector()));
1505  this->GetControlIntegrator().AddDomainData("adjoint_hessian",
1506  &(GetDZ().GetSpacialVector()));
1507 
1508  this->GetOutputHandler()->Write(
1509  "\tComputing Representation of the Hessian:",
1510  5 + this->GetBasePriority());
1511  //Preparations for Control In The Dirichlet Data
1512  VECTOR tmp;
1513  VECTOR tmp_second;
1514  if (this->GetProblem()->HasControlInDirichletData())
1515  {
1516  tmp.reinit(GetU().GetSpacialVector());
1517  tmp_second.reinit(GetU().GetSpacialVector());
1518  this->SetProblemType("adjoint");
1519  this->GetIntegrator().AddDomainData("last_newton_solution",
1520  &(GetZ().GetSpacialVector()));
1521 
1522  this->GetIntegrator().ComputeNonlinearResidual(*(this->GetProblem()),
1523  tmp_second, false);
1524  tmp_second *= -1.;
1525 
1526  this->GetIntegrator().DeleteDomainData("last_newton_solution");
1527  this->SetProblemType("adjoint_hessian");
1528  this->GetIntegrator().AddDomainData("last_newton_solution",
1529  &(GetDZ().GetSpacialVector()));
1530 
1531  this->GetIntegrator().ComputeNonlinearResidual(*(this->GetProblem()),
1532  tmp, false);
1533  tmp *= -1.;
1534 
1535  this->GetIntegrator().DeleteDomainData("last_newton_solution");
1536  }
1537  //Endof Dirichletdata Preparations
1538  this->GetProblem()->DeleteAuxiliaryFromIntegrator(this->GetIntegrator());
1539 
1540  this->SetProblemType("hessian");
1541  this->GetProblem()->AddAuxiliaryToIntegrator(
1542  this->GetControlIntegrator());
1543  if (dopedim == dealdim)
1544  {
1545  this->GetIntegrator().DeleteDomainData("dq");
1546  this->GetIntegrator().DeleteDomainData("control");
1547  this->GetControlIntegrator().AddDomainData("dq",
1548  &(direction.GetSpacialVector()));
1549  this->GetControlIntegrator().AddDomainData("control",
1550  &(q.GetSpacialVector()));
1551  }
1552  else if (dopedim == 0)
1553  {
1554  this->GetIntegrator().DeleteParamData("dq");
1555  this->GetIntegrator().DeleteParamData("control");
1556  direction.UnLockCopy();
1557  q.UnLockCopy();
1558  this->GetControlIntegrator().AddParamData("dq",
1559  &(direction.GetSpacialVectorCopy()));
1560  this->GetControlIntegrator().AddParamData("control",
1561  &(q.GetSpacialVectorCopy()));
1562  }
1563  else
1564  {
1565  throw DOpEException("dopedim not implemented",
1566  "StatReducedProblem::ComputeReducedHessianVector");
1567  }
1568  if (this->GetProblem()->HasControlInDirichletData())
1569  {
1570  this->GetControlIntegrator().AddDomainData("adjoint_residual", &tmp);
1571  this->GetControlIntegrator().AddDomainData("hessian_residual",
1572  &tmp_second);
1573  }
1574 
1575  {
1576  hessian_direction_transposed = 0.;
1577  if (dopedim == dealdim)
1578  {
1579  this->GetControlIntegrator().AddDomainData("last_newton_solution",
1580  &(hessian_direction_transposed.GetSpacialVector()));
1581  this->GetControlIntegrator().ComputeNonlinearResidual(
1582  *(this->GetProblem()), hessian_direction.GetSpacialVector(),
1583  true);
1584  this->GetControlIntegrator().DeleteDomainData("last_newton_solution");
1585  }
1586  else if (dopedim == 0)
1587  {
1588  this->GetControlIntegrator().AddParamData("last_newton_solution",
1589  &(hessian_direction_transposed.GetSpacialVectorCopy()));
1590  this->GetControlIntegrator().ComputeNonlinearResidual(
1591  *(this->GetProblem()), hessian_direction.GetSpacialVector(),
1592  true);
1593  this->GetControlIntegrator().DeleteParamData("last_newton_solution");
1594  hessian_direction_transposed.UnLockCopy();
1595  }
1596  hessian_direction *= -1.;
1597  hessian_direction_transposed = hessian_direction;
1598  //Compute l^2 representation of the HessianVector
1599  //hessian Matrix is the same as control matrix
1600  build_control_matrix_ =
1601  this->GetControlNonlinearSolver().NonlinearSolve(
1602  *(this->GetProblem()),
1603  hessian_direction_transposed.GetSpacialVector(), true,
1604  build_control_matrix_);
1605 
1606  this->GetOutputHandler()->Write(hessian_direction,
1607  "HessianDirection" + this->GetPostIndex(),
1608  this->GetProblem()->GetDoFType());
1609  this->GetOutputHandler()->Write(hessian_direction_transposed,
1610  "HessianDirection_Transposed" + this->GetPostIndex(),
1611  this->GetProblem()->GetDoFType());
1612  }
1613 
1614  if (dopedim == dealdim)
1615  {
1616  this->GetControlIntegrator().DeleteDomainData("dq");
1617  this->GetControlIntegrator().DeleteDomainData("control");
1618  }
1619  else if (dopedim == 0)
1620  {
1621  this->GetControlIntegrator().DeleteParamData("dq");
1622  this->GetControlIntegrator().DeleteParamData("control");
1623  direction.UnLockCopy();
1624  q.UnLockCopy();
1625  }
1626  else
1627  {
1628  throw DOpEException("dopedim not implemented",
1629  "StatReducedProblem::ComputeReducedHessianVector");
1630  }
1631  this->GetIntegrator().DeleteDomainData("state");
1632  this->GetIntegrator().DeleteDomainData("adjoint");
1633  this->GetIntegrator().DeleteDomainData("tangent");
1634  this->GetIntegrator().DeleteDomainData("adjoint_hessian");
1635  this->GetControlIntegrator().DeleteDomainData("state");
1636  this->GetControlIntegrator().DeleteDomainData("adjoint");
1637  this->GetControlIntegrator().DeleteDomainData("tangent");
1638  this->GetControlIntegrator().DeleteDomainData("adjoint_hessian");
1639  if (this->GetProblem()->HasControlInDirichletData())
1640  {
1641  this->GetControlIntegrator().DeleteDomainData("adjoint_residual");
1642  this->GetControlIntegrator().DeleteDomainData("hessian_residual");
1643  }
1644  this->GetProblem()->DeleteAuxiliaryFromIntegrator(
1645  this->GetControlIntegrator());
1646 
1647  }
1648 
1649  /******************************************************/
1650 
1651  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
1652  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
1653  typename VECTOR, int dopedim, int dealdim>
1654  void
1655  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
1656  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::ComputeReducedGradientOfGlobalConstraints(
1657  unsigned int num, const ControlVector<VECTOR>& q,
1658  const ConstraintVector<VECTOR>& g, ControlVector<VECTOR>& gradient,
1659  ControlVector<VECTOR>& gradient_transposed)
1660  {
1661  //FIXME: If the global constraints depend on u we need to calculate a corresponding
1662  // dual solution before we can calculate the gradient.
1663  std::stringstream out;
1664  out << "Computing Reduced Gradient of global constraint " << num << " :";
1665  this->GetOutputHandler()->Write(out, 4 + this->GetBasePriority());
1666  //Compute derivatives of global constraints
1667  this->SetProblemType("global_constraint_gradient", num);
1668 
1669  if (dopedim == dealdim)
1670  {
1671  this->GetControlIntegrator().AddDomainData("control",
1672  &(q.GetSpacialVector()));
1673  }
1674  else if (dopedim == 0)
1675  {
1676  this->GetControlIntegrator().AddParamData("control",
1677  &(q.GetSpacialVectorCopy()));
1678  }
1679  else
1680  {
1681  throw DOpEException("dopedim not implemented",
1682  "StatReducedProblem::ComputeReducedGradient");
1683  }
1684  this->GetProblem()->AddAuxiliaryToIntegrator(
1685  this->GetControlIntegrator());
1686  this->GetControlIntegrator().AddDomainData("constraints_local",
1687  &g.GetSpacialVector("local"));
1688  this->GetControlIntegrator().AddParamData("constraints_global",
1689  &g.GetGlobalConstraints());
1690 
1691  //Compute
1692  this->GetControlIntegrator().ComputeNonlinearRhs(*(this->GetProblem()),
1693  gradient.GetSpacialVector(), true);
1694  gradient_transposed = gradient;
1695 
1696  this->GetControlIntegrator().DeleteDomainData("constraints_local");
1697  this->GetControlIntegrator().DeleteParamData("constraints_global");
1698  this->GetProblem()->DeleteAuxiliaryFromIntegrator(
1699  this->GetControlIntegrator());
1700  if (dopedim == dealdim)
1701  {
1702  this->GetControlIntegrator().DeleteDomainData("control");
1703  }
1704  else if (dopedim == 0)
1705  {
1706  this->GetControlIntegrator().DeleteParamData("control");
1707  q.UnLockCopy();
1708  }
1709  else
1710  {
1711  throw DOpEException("dopedim not implemented",
1712  "StatReducedProblem::ComputeReducedGradient");
1713  }
1714  }
1715 
1716  /******************************************************/
1717 
1718  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
1719  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
1720  typename VECTOR, int dopedim, int dealdim>
1721  void
1722  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
1723  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::WriteToFile(
1724  const VECTOR &v, std::string name, std::string outfile,
1725  std::string dof_type, std::string filetype)
1726  {
1727  if (dof_type == "state")
1728  {
1729  auto& data_out =
1730  this->GetProblem()->GetSpaceTimeHandler()->GetDataOut();
1731  data_out.attach_dof_handler(
1732  this->GetProblem()->GetSpaceTimeHandler()->GetStateDoFHandler());
1733 
1734  data_out.add_data_vector(v, name);
1735  data_out.build_patches();
1736 
1737  std::ofstream output(outfile.c_str());
1738 
1739  if (filetype == ".vtk")
1740  {
1741  data_out.write_vtk(output);
1742  }
1743  else if (filetype == ".gpl")
1744  {
1745  data_out.write_gnuplot(output);
1746  }
1747  else
1748  {
1749  throw DOpEException(
1750  "Don't know how to write filetype `" + filetype + "'!",
1751  "StatReducedProblem::WriteToFile");
1752  }
1753  data_out.clear();
1754  }
1755  else if (dof_type == "control")
1756  {
1757 #if dope_dimension >0
1758  auto& data_out = this->GetProblem()->GetSpaceTimeHandler()->GetDataOut();
1759  data_out.attach_dof_handler (this->GetProblem()->GetSpaceTimeHandler()->GetControlDoFHandler());
1760 
1761  data_out.add_data_vector (v,name);
1762  data_out.build_patches ();
1763 
1764  std::ofstream output(outfile.c_str());
1765 
1766  if(filetype == ".vtk")
1767  {
1768  data_out.write_vtk (output);
1769  }
1770  else if (filetype == ".gpl")
1771  {
1772  data_out.write_gnuplot(output);
1773  }
1774  else
1775  {
1776  throw DOpEException("Don't know how to write filetype `" + filetype + "'!",
1777  "StatReducedProblem::WriteToFile");
1778  }
1779  data_out.clear();
1780 #else
1781  if (filetype == ".txt")
1782  {
1783  std::ofstream output(outfile.c_str());
1784  Vector<double> off;
1785  off = v;
1786  for (unsigned int i = 0; i < off.size(); i++)
1787  {
1788  output << off(i) << std::endl;
1789  }
1790  }
1791  else
1792  {
1793  throw DOpEException(
1794  "Don't know how to write filetype `" + filetype + "'!",
1795  "StatReducedProblem::WriteToFile");
1796  }
1797 #endif
1798  }
1799  else
1800  {
1801  throw DOpEException("No such DoFHandler `" + dof_type + "'!",
1802  "StatReducedProblem::WriteToFile");
1803  }
1804  }
1805 
1806  /******************************************************/
1807 
1808  template<typename CONTROLNONLINEARSOLVER, typename NONLINEARSOLVER,
1809  typename CONTROLINTEGRATOR, typename INTEGRATOR, typename PROBLEM,
1810  typename VECTOR, int dopedim, int dealdim>
1811  void
1812  StatReducedProblem<CONTROLNONLINEARSOLVER, NONLINEARSOLVER,
1813  CONTROLINTEGRATOR, INTEGRATOR, PROBLEM, VECTOR, dopedim, dealdim>::WriteToFile(
1814  const ControlVector<VECTOR> &v, std::string name,
1815  std::string dof_type)
1816  {
1817  this->GetOutputHandler()->Write(v.GetSpacialVector(), name, dof_type);
1818  }
1819 
1820 /******************************************************/
1821 }
1822 #endif
DOpE::StatReducedProblem::ComputeRefinementIndicators
void ComputeRefinementIndicators(const ControlVector< VECTOR > &q, DWRC &dwrc, PDE &pde)
Definition: statreducedproblem.h:1335
DOpE::StatReducedProblem::GetControlIntegrator
CONTROLINTEGRATOR & GetControlIntegrator()
Definition: statreducedproblem.h:393
DOpE::ControlVector::UnLockCopy
void UnLockCopy() const
Definition: controlvector.h:199
DOpE::ConstraintVector
Definition: constraintvector.h:48
DOpE::ControlVector::GetSpacialVectorCopy
const dealii::Vector< double > & GetSpacialVectorCopy() const
Definition: controlvector.cc:276
DOpE::StatReducedProblem::GetNonlinearSolver
NONLINEARSOLVER & GetNonlinearSolver(std::string type)
Definition: statreducedproblem.h:515
DOpE::StatReducedProblem::GetControlNonlinearSolver
CONTROLNONLINEARSOLVER & GetControlNonlinearSolver()
Definition: statreducedproblem.h:541
DOpE::StatReducedProblem::ComputeReducedFunctionals
void ComputeReducedFunctionals(const ControlVector< VECTOR > &q)
Definition: statreducedproblem.h:1235
DOpE::StatReducedProblem::ReInit
void ReInit()
Definition: statreducedproblem.h:563
DOpE::StatReducedProblem::GetIntegrator
INTEGRATOR & GetIntegrator()
Definition: statreducedproblem.h:388
DOpE::StatReducedProblem::ComputeDualForErrorEstimation
void ComputeDualForErrorEstimation(const ControlVector< VECTOR > &q, DOpEtypes::WeightComputation weight_comp)
Definition: statreducedproblem.h:899
DOpE::StatReducedProblem::GetU
StateVector< VECTOR > & GetU()
Definition: statreducedproblem.h:350
DOpE::ParameterReader
Definition: parameterreader.h:36
DOpE::StatReducedProblem::WriteToFile
void WriteToFile(const VECTOR &v, std::string name, std::string outfile, std::string dof_type, std::string filetype)
Definition: statreducedproblem.h:1723
DOpE::StatReducedProblem::ComputeReducedAdjoint
void ComputeReducedAdjoint(const ControlVector< VECTOR > &q)
Definition: statreducedproblem.h:832
DOpE::ConstraintVector::GetGlobalConstraints
const dealii::Vector< double > & GetGlobalConstraints() const
Definition: constraintvector.cc:174
DOpE::StatReducedProblem
Definition: statreducedproblem.h:82
DOpE::StatReducedProblem::GetZ
StateVector< VECTOR > & GetZ()
Definition: statreducedproblem.h:355
DOpE::StatReducedProblem::GetZForEE
StateVector< VECTOR > & GetZForEE()
Definition: statreducedproblem.h:378
DOpE::ControlVector
Definition: controlvector.h:49
DOpE::StatReducedProblem::declare_params
static void declare_params(ParameterReader &param_reader)
Definition: statreducedproblem.h:431
DOpE::StatReducedProblem::GetControlBoxConstraints
void GetControlBoxConstraints(ControlVector< VECTOR > &lb, ControlVector< VECTOR > &ub)
Definition: statreducedproblem.h:818
DOpE::StatReducedProblem::ComputeReducedCostFunctional
double ComputeReducedCostFunctional(const ControlVector< VECTOR > &q)
Definition: statreducedproblem.h:1141
DOpE::SolutionExtractor
Definition: solutionextractor.h:40
DOpE::StatReducedProblem::GetZForEE
const StateVector< VECTOR > & GetZForEE() const
Definition: statreducedproblem.h:373
DOpE::StatReducedProblem::GetDZ
StateVector< VECTOR > & GetDZ()
Definition: statreducedproblem.h:365
DOpE::StatReducedProblem::ComputeReducedGradient
void ComputeReducedGradient(const ControlVector< VECTOR > &q, ControlVector< VECTOR > &gradient, ControlVector< VECTOR > &gradient_transposed)
Definition: statreducedproblem.h:979
DOpE::StateVector
Definition: statevector.h:50
DOpE::ControlVector::GetSpacialVector
VECTOR & GetSpacialVector()
******************************************************/
Definition: controlvector.cc:204
DOpE::ConstraintVector::HasType
bool HasType(std::string name) const
******************************************************/
Definition: constraintvector.cc:119
DOpE::DOpEtypes::VectorStorageType
VectorStorageType
Definition: dopetypes.h:120
DOpE::DOpEtypes::WeightComputation
WeightComputation
Definition: dopetypes.h:80
DOpE::StatReducedProblem::ComputeReducedState
void ComputeReducedState(const ControlVector< VECTOR > &q)
Definition: statreducedproblem.h:594
DOpE::ReducedProblemInterface
Definition: reducedprobleminterface.h:335
DOpE::StatReducedProblem::StatReducedProblem
StatReducedProblem(PROBLEM *OP, DOpEtypes::VectorStorageType state_behavior, ParameterReader &param_reader, INTEGRATORDATACONT &idc, int base_priority=0)
Definition: statreducedproblem.h:443
DOpE::StatReducedProblem::ComputeReducedHessianVector
void ComputeReducedHessianVector(const ControlVector< VECTOR > &q, const ControlVector< VECTOR > &direction, ControlVector< VECTOR > &hessian_direction, ControlVector< VECTOR > &hessian_direction_transposed)
Definition: statreducedproblem.h:1436
DOpE::StatReducedProblem::ComputeReducedConstraints
bool ComputeReducedConstraints(const ControlVector< VECTOR > &q, ConstraintVector< VECTOR > &g)
Definition: statreducedproblem.h:657
DOpE::StatReducedProblem::StateSizeInfo
void StateSizeInfo(std::stringstream &out)
Definition: statreducedproblem.h:257
DOpE::ConstraintVector::GetSpacialVector
VECTOR & GetSpacialVector(std::string name)
Definition: constraintvector.cc:130
DOpE::DOpEException
Definition: dopeexception.h:35
DOpE::StatReducedProblem::WriteToFile
void WriteToFile(const std::vector< double > &, std::string)
Definition: statreducedproblem.h:299
DOpE::StatReducedProblem::GetU
const StateVector< VECTOR > & GetU() const
Definition: statreducedproblem.h:345
DOpE::ConstraintVector::IsFeasible
virtual bool IsFeasible() const
Definition: constraintvector.cc:574
DOpE::StatReducedProblem::GetDU
StateVector< VECTOR > & GetDU()
Definition: statreducedproblem.h:360
DOpE::ReducedProblemInterface::ReInit
virtual void ReInit()
Definition: reducedprobleminterface.h:357
DOpE::StatReducedProblem::~StatReducedProblem
~StatReducedProblem()
Definition: statreducedproblem.h:504
DOpE::StatReducedProblem::ComputeReducedGradientOfGlobalConstraints
void ComputeReducedGradientOfGlobalConstraints(unsigned int num, const ControlVector< VECTOR > &q, const ConstraintVector< VECTOR > &g, ControlVector< VECTOR > &gradient, ControlVector< VECTOR > &gradient_transposed)
Definition: statreducedproblem.h:1656
Generated on Mon Jan 25 2016 11:32:19 for DOpE by   doxygen 1.8.6