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