controlvalidation_core.cpp File Reference

: core of the control solution More...

#include "./cores.h"
#include "../toolkits/toolkits.h"
#include "../classes/classes.h"
#include "../shared/shared.h"
#include "../modules/modules.h"

Go to the source code of this file.

Functions
void	controlvalidation_core (FemModel *femmodel)

Detailed Description

: core of the control solution

Definition in file controlvalidation_core.cpp.

Function Documentation

controlvalidation_core()

void controlvalidation_core

(

FemModel *

femmodel

)

Definition at line 155 of file controlvalidation_core.cpp.

                                               {
 
   int         solution_type,n;
   int         num_responses;
   IssmDouble  j0,j;
   IssmDouble  Ialpha,exponent,alpha;
   IssmDouble* scaling_factors = NULL;
   IssmDouble* jlist = NULL;
   int my_rank=IssmComm::GetRank();
 
   /*Recover parameters used throughout the solution*/
   femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
   femmodel->parameters->SetParam(false,SaveResultsEnum);
   femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
   femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
 
   /*Get initial guess*/
   IssmPDouble* X0 = NULL;
   GetPassiveVectorFromControlInputsx(&X0,&n,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
 
   /*Allocate vectors*/
   IssmDouble*  X = xNew<IssmDouble>(n);
   IssmPDouble* G = xNew<IssmPDouble>(n);
 
   /*out of solution_type, figure out solution core and adjoint function pointer*/
   void (*solutioncore)(FemModel*)=NULL;
   CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
 
   #if defined(_HAVE_ADOLC_)
   /*{{{*/
   IssmDouble* aX=xNew<IssmDouble>(n);
   if(my_rank==0){
      for(int i=0;i<n;i++){
         aX[i]<<=X0[i];
      }
   }
   _error_("not implemented yet...");
   /*}}}*/
   #elif defined(_HAVE_CODIPACK_)
   simul_starttrace2(femmodel);
   IssmDouble* aX=xNew<IssmDouble>(n);
   auto& tape_codi = IssmDouble::getGlobalTape();
   codi_global.input_indices.clear();
   if(my_rank==0){
      for (int i=0;i<n;i++) {
         aX[i]=X0[i];
         tape_codi.registerInput(aX[i]);
         codi_global.input_indices.push_back(aX[i].getGradientData());
      }
   }
   SetControlInputsFromVectorx(femmodel,aX);
   xDelete(aX);
 
   if(VerboseControl()) _printf0_("   Compute Initial cost function\n");
   solutioncore(femmodel);
 
   /*Get Dependents*/
   IssmDouble  output_value;
   int         num_dependents;
   IssmPDouble *dependents;
   DataSet*    dependent_objects=NULL;
   IssmDouble  J=0.;
   femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
   femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
 
   /*Go through our dependent variables, and compute the response:*/
   dependents=xNew<IssmPDouble>(num_dependents);
   codi_global.output_indices.clear();
   for(int i=0;i<dependent_objects->Size();i++){
      DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
      if(solution_type==TransientSolutionEnum){
         output_value = dep->GetValue();
      }
      else{
         dep->Responsex(&output_value,femmodel);
      }
      _printf0_("=== output ="<<output_value<<" \n");
      if(my_rank==0) {
         tape_codi.registerOutput(output_value);
         dependents[i] = output_value.getValue();
         codi_global.output_indices.push_back(output_value.getGradientData());
         J+=output_value;
      }
   }
   j0 = J;
   _printf0_("Initial cost function J(x) = "<<setw(12)<<setprecision(7)<<j0<<"\n");
   _assert_(j0>0.); 
   simul_stoptrace2();
   /*initialize direction index in the weights vector: */
   if(my_rank==0){
      tape_codi.setGradient(codi_global.output_indices[0],1.0);
   }
   tape_codi.evaluate();
 
   /*Get gradient for this dependent */
   auto in_size = codi_global.input_indices.size();
   for(size_t i = 0; i < in_size; ++i) {
      G[i] = tape_codi.getGradient(codi_global.input_indices[i]);
   }
   #else
   /*{{{*/
   void (*adjointcore)(FemModel*)  = NULL;
   AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type);
 
   if(VerboseControl()) _printf0_("   Compute Initial solution\n");
   solutioncore(femmodel);
   if(VerboseControl()) _printf0_("   Compute Adjoint\n");
   adjointcore(femmodel);
 
   if(VerboseControl()) _printf0_("   Compute Initial cost function\n");
   femmodel->CostFunctionx(&j0,&jlist,NULL);
   _printf0_("Initial cost function J(x) = "<<setw(12)<<setprecision(7)<<j0<<"\n");
   xDelete<IssmDouble>(jlist);
 
   if(VerboseControl()) _printf0_("   Compute Gradient\n");
   Gradjx(&G,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
   for(int i=0;i<n;i++) G[i] = -G[i];
   /*}}}*/
   #endif
 
   /*Allocate output*/
   int num = 26;
   IssmDouble* output = xNew<IssmDouble>(2*num);
 
   /*Start loop*/
   _printf0_("       alpha      Ialpha \n");
   _printf0_("_________________________\n");
   for(int m=0;m<num;m++){
 
      /*Create new vector*/
      alpha    = pow(2.,-m);
      for(int i=0;i<n;i++) X[i] = X0[i] + alpha*scaling_factors[0];
 
      /*Calculate j(k+alpha delta k) */
      SetControlInputsFromVectorx(femmodel,X);
      solutioncore(femmodel);
 
      #if defined(_HAVE_CODIPACK_)
      j=0.;
      for(int i=0;i<dependent_objects->Size();i++){
         DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
         if(solution_type==TransientSolutionEnum){
            output_value = dep->GetValue();
         }
         else{
            dep->Responsex(&output_value,femmodel);
         }
         j+=output_value;
      }
      #else
      femmodel->CostFunctionx(&j,NULL,NULL);
      #endif
 
      IssmDouble Den = 0.;
      for(int i=0;i<n;i++) Den += alpha* G[i] * scaling_factors[0];
      Ialpha = fabs((j - j0)/Den - 1.);
      _assert_(fabs(Den)>0.); 
 
      _printf0_(" " << setw(11) << setprecision (5)<<alpha<<" " << setw(11) << setprecision (5)<<Ialpha<<"\n");
      output[m*2+0] = alpha;
      output[m*2+1] = Ialpha;
   }
 
   /*output*/
   #ifdef _HAVE_AD_
   IssmPDouble* J_passive=xNew<IssmPDouble>(2*num);
   for(int i=0;i<2*num;i++) J_passive[i]=reCast<IssmPDouble>(output[i]);
   femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,J_passive,num,2,0,0));
   xDelete<IssmPDouble>(J_passive);
   IssmDouble* aG=xNew<IssmDouble>(n);
   for(int i=0;i<n;i++) aG[i] = G[i];
   ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG);
   xDelete<IssmDouble>(aG);
   #else
   femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,output,num,2,0,0));
   ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
   #endif
   femmodel->OutputControlsx(&femmodel->results);
 
   /*Clean up and return*/
   xDelete<IssmDouble>(output);
   xDelete<IssmPDouble>(G);
   xDelete<IssmDouble>(X);
   xDelete<double>(X0);
   xDelete<IssmDouble>(scaling_factors);
}