Changeset 17907

Timestamp:

05/01/14 13:40:47 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: simplifying inversions

Location:

issm/trunk-jpl/src/c

Files:

• Makefile.am (modified) (2 diffs)
• analyses/Analysis.h (modified) (1 diff)
• classes/FemModel.cpp (modified) (3 diffs)
• classes/FemModel.h (modified) (1 diff)
• cores/control_core.cpp (modified) (5 diffs)
• cores/controlm1qn3_core.cpp (modified) (1 diff)
• cores/controltao_core.cpp (modified) (1 diff)
• cores/cores.h (modified) (2 diffs)
• cores/objectivefunction.cpp (deleted)
• shared/Numerics/BrentSearch.cpp (modified) (1 diff)
• shared/Numerics/OptArgs.h (deleted)
• shared/Numerics/OptimalSearch.cpp (deleted)
• shared/Numerics/numerics.h (modified) (2 diffs)

issm/trunk-jpl/src/c/Makefile.am

@@ -209 +209 @@
                                         ./shared/Numerics/extrema.cpp\
                                         ./shared/Numerics/XZvectorsToCoordinateSystem.cpp\
-                                        ./shared/Numerics/OptArgs.h\
                                         ./shared/Numerics/OptPars.h\
                                         ./shared/Exceptions/exceptions.h\
@@ -420 +419 @@
                                         ./classes/Inputs/ControlInput.cpp\
                                         ./shared/Numerics/BrentSearch.cpp\
-                                        ./shared/Numerics/OptimalSearch.cpp \
                                         ./cores/control_core.cpp\
                                         ./cores/controltao_core.cpp\
                                         ./cores/controlm1qn3_core.cpp\
-                                        ./cores/objectivefunction.cpp\
                                         ./cores/gradient_core.cpp\
                                         ./cores/adjointstressbalance_core.cpp\

issm/trunk-jpl/src/c/analyses/Analysis.h

@@ -17 +17 @@
 class ElementMatrix;
 class Gauss;
+class FemModel;
 
 class Analysis{

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -1174 +1174 @@
 }
 /*}}}*/
-void FemModel::CostFunctionx(IssmDouble* pJ){/*{{{*/
+void FemModel::CostFunctionx(IssmDouble* pJ,IssmDouble** pJlist,int* pn){/*{{{*/
 
         /*Intermediary*/
@@ -1188 +1188 @@
         this->RequestedOutputsx(&cost_functions,responses,num_responses);
 
-        /*Add all contributions one by one*/
-        IssmDouble J=0;
-        //_printf0_("list of misfits: ");
+        /*Get and add all contributions one by one*/
+        IssmDouble  J=0;
+        IssmDouble* Jlist = xNew<IssmDouble>(num_responses);
         for(int i=0;i<num_responses;i++){
                 ExternalResult* result=(ExternalResult*)cost_functions->GetObjectByOffset(i);
-                J += reCast<IssmDouble>(result->GetValue());
-                //_printf0_(J<<" ");
-        }
-        //_printf0_(" \n");
+                Jlist[i] = reCast<IssmDouble>(result->GetValue());
+                J       += Jlist[i];
+        }
         _assert_(cost_functions->Size()==num_responses);
 
@@ -1202 +1201 @@
         delete cost_functions;
         xDelete<int>(responses);
-        *pJ=J;
+        if(pJ)     *pJ     = J;
+        if(pJlist) *pJlist = Jlist;
+        else        xDelete<IssmDouble>(Jlist);
+        if(pn)     *pn     = num_responses;
 }
 /*}}}*/

issm/trunk-jpl/src/c/classes/FemModel.h

@@ -86 +86 @@
                 void Responsex(IssmDouble* presponse,const char* response_descriptor);
                 void OutputControlsx(Results **presults);
-                void CostFunctionx( IssmDouble* pJ);
+                void CostFunctionx(IssmDouble* pJ,IssmDouble** pJlist,int* pn);
                 void ThicknessAbsGradientx( IssmDouble* pJ);
                 #ifdef _HAVE_GIA_

issm/trunk-jpl/src/c/cores/control_core.cpp

@@ -12 +12 @@
 /*Local prototypes*/
 bool controlconvergence(IssmDouble J, IssmDouble tol_cm);
+IssmDouble objectivefunction(IssmDouble search_scalar,void* optargs);
 
 void control_core(FemModel* femmodel){
@@ -31 +32 @@
         /*intermediary: */
         IssmDouble search_scalar = 1.;
-        OptArgs    optargs;
         OptPars    optpars;
 
@@ -65 +65 @@
 
         /*Initialize some of the BrentSearch arguments: */
-        optargs.femmodel=femmodel;
         optpars.xmin=0; optpars.xmax=1;
 
@@ -84 +83 @@
                 if(VerboseControl()) _printf0_("   optimizing along gradient direction\n");
                 optpars.maxiter=reCast<int,IssmDouble>(maxiter[n]); optpars.cm_jump=cm_jump[n];
-                BrentSearch(&search_scalar,J+n,&optpars,&objectivefunction,&optargs);
+                BrentSearch(&search_scalar,J+n,&optpars,&objectivefunction,(void*)femmodel);
 
                 if(VerboseControl()) _printf0_("   updating parameter using optimized search scalar\n"); //true means update save controls
@@ -128 +127 @@
         return converged;
 }
+
+IssmDouble objectivefunction(IssmDouble search_scalar,void* optargs){
+
+        /*output: */
+        IssmDouble J;
+
+        /*parameters: */
+        int        solution_type,analysis_type;
+        bool       isFS       = false;
+        bool       conserve_loads = true;
+        FemModel  *femmodel       = (FemModel*)optargs;
+
+        /*Recover parameters: */
+        femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
+        femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+
+        /*set analysis type to compute velocity: */
+        if (solution_type==SteadystateSolutionEnum || solution_type==StressbalanceSolutionEnum){
+                femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
+        }
+        else if (solution_type==BalancethicknessSolutionEnum){
+                femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
+        }
+        else if (solution_type==BalancethicknessSoftSolutionEnum){
+                femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
+        }
+        else{
+                _error_("Solution " << EnumToStringx(solution_type) << " not implemented yet");
+        }
+
+        /*update parameter according to scalar: */ //false means: do not save control
+        InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,false);
+
+        /*Run stressbalance with updated inputs: */
+        if (solution_type==SteadystateSolutionEnum){
+                stressbalance_core(femmodel);   //We need a 3D velocity!! (vz is required for the next thermal run)
+        }
+        else if (solution_type==StressbalanceSolutionEnum){
+                solutionsequence_nonlinear(femmodel,conserve_loads); 
+        }
+        else if (solution_type==BalancethicknessSolutionEnum){
+                solutionsequence_linear(femmodel); 
+        }
+        else if (solution_type==BalancethicknessSoftSolutionEnum){
+                /*Don't do anything*/
+        }
+        else{
+                _error_("Solution " << EnumToStringx(solution_type) << " not implemented yet");
+        }
+
+        /*Compute misfit for this velocity field.*/
+        femmodel->CostFunctionx(&J,NULL,NULL);
+
+        /*Free ressources:*/
+        return J;
+}

issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp

@@ -160 +160 @@
 
         /*Compute objective function*/
-        femmodel->CostFunctionx(pf);
+        IssmDouble* Jlist = NULL;
+        femmodel->CostFunctionx(pf,&Jlist,NULL);
         _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
 
         /*Retrieve objective functions independently*/
-        for(int i=0;i<num_responses;i++){
-                femmodel->Responsex(&f,EnumToStringx(responses[i]));
-                _printf0_(" "<<setw(12)<<setprecision(7)<<f);
-        }
+        for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[i]);
         _printf0_("\n");
+        xDelete<IssmDouble>(Jlist);
 
         if(indic==0){

issm/trunk-jpl/src/c/cores/controltao_core.cpp

@@ -135 +135 @@
 
         /*Compute objective function*/
-        femmodel->CostFunctionx(fcn);
+        femmodel->CostFunctionx(fcn,NULL,NULL);
 
         /*Compute gradient*/

issm/trunk-jpl/src/c/cores/cores.h

@@ -7 +7 @@
 
 /*forward declarations: */
-struct OptArgs;
 class FemModel;
 class Parameters;
@@ -45 +44 @@
 void gia_core(FemModel* femmodel);
 void damage_core(FemModel* femmodel);
-IssmDouble objectivefunction(IssmDouble search_scalar,OptArgs* optargs);
+IssmDouble objectivefunction(IssmDouble search_scalar,FemModel* femmodel);
 
 //optimization

issm/trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp

@@ -17 +17 @@
 #include "./Verbosity.h"
 #include "./OptPars.h"
-#include "./OptArgs.h"
 #include "./types.h"
 #include "./isnan.h"
 
-void BrentSearch(IssmDouble* psearch_scalar,IssmDouble* pJ,OptPars* optpars,IssmDouble (*f)(IssmDouble,OptArgs*), OptArgs* optargs){
+void BrentSearch(IssmDouble* psearch_scalar,IssmDouble* pJ,OptPars* optpars,IssmDouble (*f)(IssmDouble,void*),void* optargs){
 
         /* This routine is optimizing a given function using Brent's method

issm/trunk-jpl/src/c/shared/Numerics/numerics.h

@@ -18 +18 @@
 #include "./types.h"
 #include "./constants.h"
-#include "./OptArgs.h"
 #include "./OptPars.h"
 
@@ -31 +30 @@
 int         min(int a,int b);
 int         max(int a,int b);
-IssmDouble  OptFunc(IssmDouble scalar, OptArgs *optargs);
-void        BrentSearch(IssmDouble *psearch_scalar,IssmDouble*pJ,OptPars*optpars,IssmDouble (*f)(IssmDouble,OptArgs*), OptArgs*optargs);
-void        OptimalSearch(IssmDouble *psearch_scalar,IssmDouble*pJ,OptPars*optpars,IssmDouble (*f)(IssmDouble,OptArgs*), OptArgs*optargs);
+void        BrentSearch(IssmDouble *psearch_scalar,IssmDouble*pJ,OptPars*optpars,IssmDouble (*f)(IssmDouble,void*),void* optargs);
 void        cross(IssmDouble *result,IssmDouble*vector1,IssmDouble*vector2);
 void        XZvectorsToCoordinateSystem(IssmDouble *T,IssmDouble*xzvectors);