Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 17906) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 17907) @@ -209,5 +209,4 @@ ./shared/Numerics/extrema.cpp\ ./shared/Numerics/XZvectorsToCoordinateSystem.cpp\ - ./shared/Numerics/OptArgs.h\ ./shared/Numerics/OptPars.h\ ./shared/Exceptions/exceptions.h\ @@ -420,9 +419,7 @@ ./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\ Index: /issm/trunk-jpl/src/c/analyses/Analysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/Analysis.h (revision 17906) +++ /issm/trunk-jpl/src/c/analyses/Analysis.h (revision 17907) @@ -17,4 +17,5 @@ class ElementMatrix; class Gauss; +class FemModel; class Analysis{ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 17906) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 17907) @@ -1174,5 +1174,5 @@ } /*}}}*/ -void FemModel::CostFunctionx(IssmDouble* pJ){/*{{{*/ +void FemModel::CostFunctionx(IssmDouble* pJ,IssmDouble** pJlist,int* pn){/*{{{*/ /*Intermediary*/ @@ -1188,13 +1188,12 @@ 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(num_responses); for(int i=0;iGetObjectByOffset(i); - J += reCast(result->GetValue()); - //_printf0_(J<<" "); - } - //_printf0_(" \n"); + Jlist[i] = reCast(result->GetValue()); + J += Jlist[i]; + } _assert_(cost_functions->Size()==num_responses); @@ -1202,5 +1201,8 @@ delete cost_functions; xDelete(responses); - *pJ=J; + if(pJ) *pJ = J; + if(pJlist) *pJlist = Jlist; + else xDelete(Jlist); + if(pn) *pn = num_responses; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 17906) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 17907) @@ -86,5 +86,5 @@ 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_ Index: /issm/trunk-jpl/src/c/cores/control_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/control_core.cpp (revision 17906) +++ /issm/trunk-jpl/src/c/cores/control_core.cpp (revision 17907) @@ -12,4 +12,5 @@ /*Local prototypes*/ bool controlconvergence(IssmDouble J, IssmDouble tol_cm); +IssmDouble objectivefunction(IssmDouble search_scalar,void* optargs); void control_core(FemModel* femmodel){ @@ -31,5 +32,4 @@ /*intermediary: */ IssmDouble search_scalar = 1.; - OptArgs optargs; OptPars optpars; @@ -65,5 +65,4 @@ /*Initialize some of the BrentSearch arguments: */ - optargs.femmodel=femmodel; optpars.xmin=0; optpars.xmax=1; @@ -84,5 +83,5 @@ if(VerboseControl()) _printf0_(" optimizing along gradient direction\n"); optpars.maxiter=reCast(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,2 +127,59 @@ 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; +} Index: /issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp (revision 17906) +++ /issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp (revision 17907) @@ -160,13 +160,12 @@ /*Compute objective function*/ - femmodel->CostFunctionx(pf); + IssmDouble* Jlist = NULL; + femmodel->CostFunctionx(pf,&Jlist,NULL); _printf0_("f(x) = "<Responsex(&f,EnumToStringx(responses[i])); - _printf0_(" "<(Jlist); if(indic==0){ Index: /issm/trunk-jpl/src/c/cores/controltao_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/controltao_core.cpp (revision 17906) +++ /issm/trunk-jpl/src/c/cores/controltao_core.cpp (revision 17907) @@ -135,5 +135,5 @@ /*Compute objective function*/ - femmodel->CostFunctionx(fcn); + femmodel->CostFunctionx(fcn,NULL,NULL); /*Compute gradient*/ Index: /issm/trunk-jpl/src/c/cores/cores.h =================================================================== --- /issm/trunk-jpl/src/c/cores/cores.h (revision 17906) +++ /issm/trunk-jpl/src/c/cores/cores.h (revision 17907) @@ -7,5 +7,4 @@ /*forward declarations: */ -struct OptArgs; class FemModel; class Parameters; @@ -45,5 +44,5 @@ 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 Index: sm/trunk-jpl/src/c/cores/objectivefunction.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/objectivefunction.cpp (revision 17906) +++ (revision ) @@ -1,77 +1,0 @@ -/*!\file: objectivefunction - * \brief objective function that returns a misfit, for a certain parameter. - */ - -/*include files: {{{*/ -#ifdef HAVE_CONFIG_H - #include -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif -#include "./cores.h" -#include "../toolkits/toolkits.h" -#include "../classes/classes.h" -#include "../shared/shared.h" -#include "../solutionsequences/solutionsequences.h" -#include "../modules/modules.h" -/*}}}*/ - -IssmDouble objectivefunction(IssmDouble search_scalar,OptArgs* optargs){ - - /*output: */ - IssmDouble J; - - /*parameters: */ - int solution_type,analysis_type; - bool isFS = false; - bool conserve_loads = true; - FemModel *femmodel = NULL; - - /*Recover finite element model: */ - femmodel=optargs->femmodel; - - /*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); - - /*Free ressources:*/ - return J; -} Index: /issm/trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp (revision 17906) +++ /issm/trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp (revision 17907) @@ -17,9 +17,8 @@ #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 Index: sm/trunk-jpl/src/c/shared/Numerics/OptArgs.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Numerics/OptArgs.h (revision 17906) +++ (revision ) @@ -1,14 +1,0 @@ -/*!\file: OptArgs.h - * \brief place holder for optimization function arguments - */ - -#ifndef _OPTARGS_H_ -#define _OPTARGS_H_ - -class FemModel; - -struct OptArgs{ - FemModel* femmodel; -}; - -#endif Index: sm/trunk-jpl/src/c/shared/Numerics/OptimalSearch.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Numerics/OptimalSearch.cpp (revision 17906) +++ (revision ) @@ -1,97 +1,0 @@ -/*!\file: OptimalSearch.cpp - * \brief optimization algorithm - */ - -#ifdef HAVE_CONFIG_H - #include -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif - -#include - -#include "../Exceptions/exceptions.h" -#include "../io/io.h" -#include "../MemOps/MemOps.h" -#include "./Verbosity.h" -#include "./OptPars.h" -#include "./OptArgs.h" -#include "./types.h" -#include "./isnan.h" - -void OptimalSearch(IssmDouble* psearch_scalar,IssmDouble* pJ,OptPars* optpars,IssmDouble (*f)(IssmDouble,OptArgs*), OptArgs* optargs){ - - /* This routine is optimizing a given function*/ - - /*function values: */ - IssmDouble fx1,fx2,fxbest; - IssmDouble x1,x2,xmin,xbest; - - /*tolerances: */ - IssmDouble seps; - IssmDouble tolerance=1.e-4; - int maxiter; - - /*counters: */ - int iter=0; - bool loop; - - /*Recover parameters:*/ - xmin =optpars->xmin; - x1 =optpars->xmin; - x2 =optpars->xmax; - maxiter =optpars->maxiter; - - //get the value of the function at the first boundary - fx1= (*f)(x1,optargs); - if (xIsNan(fx1)) _error_("Function evaluation returned NaN"); - cout<(fx2)) _error_("Function evaluation returned NaN"); - if(VerboseControl()) - _printf0_(" "<fabs(fx2-fx1)?fabs(fx2-fx1):fabs(x2-x1)) << "\n"); - - //Stop the optimization? - if ((fabs(x2-x1)+seps)=maxiter){ - if(VerboseControl()) _printf0_(" " << "exiting: Maximum number of iterations has been exceeded - increase 'maxiter'\n"); - loop=false; - } - else{ - //continue - loop=true; - } - - /*Compute x2*/ - if(fx2