Index: /issm/trunk-jpl/src/c/classes/Misfit.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Misfit.cpp (revision 22437) +++ /issm/trunk-jpl/src/c/classes/Misfit.cpp (revision 22438) @@ -12,5 +12,4 @@ #include "./classes.h" -#include "./FemModel.h" #include "./ExternalResults/ExternalResult.h" #include "./ExternalResults/Results.h" @@ -18,4 +17,5 @@ #include "./Elements/Element.h" #include "./Elements/Elements.h" +#include "./FemModel.h" #include "../modules/SurfaceAreax/SurfaceAreax.h" #include "../classes/Params/Parameters.h" Index: /issm/trunk-jpl/src/c/classes/Misfit.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Misfit.h (revision 22437) +++ /issm/trunk-jpl/src/c/classes/Misfit.h (revision 22438) @@ -8,5 +8,12 @@ /*Headers:*/ #include "./Definition.h" +#include "../datastructures/datastructures.h" +#include "./Elements/Element.h" +#include "./Elements/Elements.h" #include "./FemModel.h" +#include "../modules/SurfaceAreax/SurfaceAreax.h" +#include "../classes/Params/Parameters.h" +#include "../classes/Inputs/Input.h" +#include "../classes/gauss/Gauss.h" IssmDouble OutputDefinitionsResponsex(FemModel* femmodel,int output_enum); Index: /issm/trunk-jpl/src/c/classes/Nodalvalue.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Nodalvalue.h (revision 22437) +++ /issm/trunk-jpl/src/c/classes/Nodalvalue.h (revision 22438) @@ -9,5 +9,12 @@ /*{{{*/ #include "./Definition.h" +#include "../datastructures/datastructures.h" +#include "./Elements/Element.h" +#include "./Elements/Elements.h" #include "./FemModel.h" +#include "../modules/SurfaceAreax/SurfaceAreax.h" +#include "../classes/Params/Parameters.h" +#include "../classes/Inputs/Input.h" +#include "../classes/gauss/Gauss.h" /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Numberedcostfunction.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Numberedcostfunction.h (revision 22437) +++ /issm/trunk-jpl/src/c/classes/Numberedcostfunction.h (revision 22438) @@ -8,5 +8,19 @@ /*Headers:*/ #include "./Definition.h" +#include "../datastructures/datastructures.h" +#include "./Elements/Element.h" +#include "./Elements/Elements.h" #include "./FemModel.h" +#include "./ExternalResults/ExternalResult.h" +#include "./ExternalResults/Results.h" +#include "../modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.h" +#include "../modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.h" +#include "../modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.h" +#include "../modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.h" +#include "../modules/ThicknessAbsMisfitx/ThicknessAbsMisfitx.h" +#include "../modules/ThicknessAlongGradientx/ThicknessAlongGradientx.h" +#include "../modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.h" +#include "../modules/RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.h" +#include "../modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h" @@ -15,28 +29,28 @@ class Numberedcostfunction: public Object, public Definition{ -public: + public: -int definitionenum; -char* name; -int number_cost_functions; -int* cost_functions_list; + int definitionenum; + char* name; + int number_cost_functions; + int* cost_functions_list; + + /*Numberedcostfunction constructors, destructors :*/ + Numberedcostfunction(); + Numberedcostfunction(char* in_name, int in_definitionenum,int number_cost_functions_in,int* cost_functions_list_in); + ~Numberedcostfunction(); -/*Numberedcostfunction constructors, destructors :*/ -Numberedcostfunction(); -Numberedcostfunction(char* in_name, int in_definitionenum,int number_cost_functions_in,int* cost_functions_list_in); -~Numberedcostfunction(); + /*Object virtual function resolutoin: */ + Object* copy(); + void DeepEcho(void); + void Echo(void); + int Id(void); + void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction); + int ObjectEnum(void); -/*Object virtual function resolutoin: */ -Object* copy(); -void DeepEcho(void); -void Echo(void); -int Id(void); -void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction); -int ObjectEnum(void); - -/*Definition virtual function resolutoin: */ -int DefinitionEnum(); -char* Name(); -IssmDouble Response(FemModel* femmodel); + /*Definition virtual function resolutoin: */ + int DefinitionEnum(); + char* Name(); + IssmDouble Response(FemModel* femmodel); }; Index: /issm/trunk-jpl/src/c/classes/Regionaloutput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Regionaloutput.cpp (revision 22437) +++ /issm/trunk-jpl/src/c/classes/Regionaloutput.cpp (revision 22438) @@ -15,4 +15,5 @@ #include "./Elements/Element.h" #include "./Elements/Elements.h" +#include "./FemModel.h" #include "../classes/Params/Parameters.h" Index: /issm/trunk-jpl/src/c/classes/Regionaloutput.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Regionaloutput.h (revision 22437) +++ /issm/trunk-jpl/src/c/classes/Regionaloutput.h (revision 22438) @@ -9,5 +9,9 @@ /*{{{*/ #include "./Definition.h" +#include "../datastructures/datastructures.h" +#include "./Elements/Element.h" +#include "./Elements/Elements.h" #include "./FemModel.h" +#include "../classes/Params/Parameters.h" /*}}}*/ Index: /issm/trunk-jpl/src/c/cores/ad_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/ad_core.cpp (revision 22437) +++ /issm/trunk-jpl/src/c/cores/ad_core.cpp (revision 22438) @@ -26,5 +26,5 @@ int num_dependents=0; int num_independents=0; - bool isautodiff,iscontrol; + bool isautodiff = false; char *driver = NULL; size_t tape_stats[15]; @@ -37,7 +37,6 @@ /*AD mode on?: */ femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum); - femmodel->parameters->FindParam(&iscontrol,InversionIscontrolEnum); - - if(isautodiff && !iscontrol){ + + if(isautodiff){ #ifdef _HAVE_ADOLC_ Index: sm/trunk-jpl/src/c/cores/controlADm1qn3_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/controlADm1qn3_core.cpp (revision 22437) +++ (revision ) @@ -1,470 +1,0 @@ -/*!\file: controlm1qn3_core.cpp - * \brief: core of the control solution - */ - -#include -#include "./cores.h" -#include "../toolkits/toolkits.h" -#include "../classes/classes.h" -#include "../shared/shared.h" -#include "../modules/modules.h" -#include "../solutionsequences/solutionsequences.h" - -#if defined (_HAVE_M1QN3_) -//& !defined(_HAVE_ADOLC_) -/*m1qn3 prototypes*/ -extern "C" void *ctonbe_; // DIS mode : Conversion -extern "C" void *ctcabe_; // DIS mode : Conversion -extern "C" void *euclid_; // Scalar product -typedef void (*SimulFunc) (long* indic,long* n, double* x,double* pf,double* g,long [],float [],void* dzs); -extern "C" void m1qn3_ (void f(long* indic,long* n, double* x,double* pf,double* g,long [],float [],void* dzs), - void **, void **, void **, - long *, double [],double *, double[], double*, double *, - double *, char [], long *, long *, long *, long *, long *, long *, long [],double [], long *, - long *, long *, long [], float [],void* ); - -/*Use struct to provide arguments*/ -typedef struct{ - FemModel * femmodel; - IssmPDouble* Jlist; - int M; - int N; - int* i; -} m1qn3_struct; - -/*m1qm3 functions*/ -void simul_starttrace(FemModel* femmodel){/*{{{*/ - - /*Retrive ADOLC parameters*/ - IssmDouble gcTriggerRatio; - IssmDouble gcTriggerMaxSize; - IssmDouble obufsize; - IssmDouble lbufsize; - IssmDouble cbufsize; - IssmDouble tbufsize; - femmodel->parameters->FindParam(&gcTriggerRatio,AutodiffGcTriggerRatioEnum); - femmodel->parameters->FindParam(&gcTriggerMaxSize,AutodiffGcTriggerMaxSizeEnum); - femmodel->parameters->FindParam(&obufsize,AutodiffObufsizeEnum); - femmodel->parameters->FindParam(&lbufsize,AutodiffLbufsizeEnum); - femmodel->parameters->FindParam(&cbufsize,AutodiffCbufsizeEnum); - femmodel->parameters->FindParam(&tbufsize,AutodiffTbufsizeEnum); - - /*Set garbage collection parameters: */ - setStoreManagerControl(reCast(gcTriggerRatio),reCast(gcTriggerMaxSize)); - - /*Start trace: */ - int skipFileDeletion=1; - int keepTaylors=1; - int my_rank=IssmComm::GetRank(); - trace_on(my_rank,keepTaylors,reCast(obufsize),reCast(lbufsize),reCast(cbufsize),reCast(tbufsize),skipFileDeletion); -}/*}}}*/ -void simul_ad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){/*{{{*/ - - /*Get rank*/ - int my_rank=IssmComm::GetRank(); - - /*Recover Arguments*/ - m1qn3_struct *input_struct = (m1qn3_struct*)dzs; - FemModel *femmodel = input_struct->femmodel; - IssmPDouble *Jlist = input_struct->Jlist; - int JlistM = input_struct->M; - int JlistN = input_struct->N; - int *Jlisti = input_struct->i; - int intn = (int)*n; - - /*Recover some parameters*/ - int num_responses,num_controls,numberofvertices,solution_type; - IssmDouble* scaling_factors = NULL; - femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); - femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum); - femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum); - femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); - numberofvertices=femmodel->vertices->NumberOfVertices(); - - /*Constrain input vector and update controls*/ - double *XL = NULL; - double *XU = NULL; - GetPassiveVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound"); - GetPassiveVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound"); - for(int i=0;i(scaling_factors[c]); - if(X[index]>XU[index]) X[index]=XU[index]; - if(X[index](num_controls*numberofvertices,"t"); -#else - IssmDouble* aX=xNew(num_controls*numberofvertices); -#endif - if(my_rank==0){ - for(int i=0;i(aX); - - /*Compute solution (forward)*/ - void (*solutioncore)(FemModel*)=NULL; - CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type); - solutioncore(femmodel); - - /*Get Dependents*/ - IssmDouble output_value; - int num_dependents; - IssmPDouble *dependents; - DataSet* dependent_objects=NULL; - - femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum); - femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum); - - /*Go through our dependent variables, and compute the response:*/ - dependents=xNew(num_dependents); - for(int i=0;iSize();i++){ - DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i); - dep->Responsex(&output_value,femmodel); - if (my_rank==0) { - output_value>>=dependents[i]; - } - } - - /*Turning off trace tape*/ - trace_off(); - - /*Print tape statistics so that user can kill this run if something is off already:*/ - if(VerboseAutodiff()){ /*{{{*/ - size_t tape_stats[15]; - tapestats(my_rank,tape_stats); //reading of tape statistics - int commSize=IssmComm::GetSize(); - int *sstats=new int[7]; - sstats[0]=tape_stats[NUM_OPERATIONS]; - sstats[1]=tape_stats[OP_FILE_ACCESS]; - sstats[2]=tape_stats[NUM_LOCATIONS]; - sstats[3]=tape_stats[LOC_FILE_ACCESS]; - sstats[4]=tape_stats[NUM_VALUES]; - sstats[5]=tape_stats[VAL_FILE_ACCESS]; - sstats[6]=tape_stats[TAY_STACK_SIZE]; - int *rstats=NULL; - if (my_rank==0) rstats=new int[commSize*7]; - ISSM_MPI_Gather(sstats,7,ISSM_MPI_INT,rstats,7,ISSM_MPI_INT,0,IssmComm::GetComm()); - if (my_rank==0) { - int offset=50; - int rOffset=(commSize/10)+1; - _printf_(" ADOLC statistics: \n"); - _printf_(" "<file":"") << "\n"); - _printf_(" locations: entry size " << sizeof(locint) << " Bytes\n"); - _printf_(" "<file":"") << "\n"); - _printf_(" constant values: entry size " << sizeof(double) << " Bytes\n"); - _printf_(" "<file":"") << "\n"); - _printf_(" Taylor stack: entry size " << sizeof(revreal) << " Bytes\n"); - _printf_(" "<tape_stats[TAY_BUFFER_SIZE]?" ->file":"") << "\n"); - delete []rstats; - } - delete [] sstats; - } /*}}}*/ - - /*diverse: */ - int dummy; - int num_independents=0; - - /*intermediary: */ - IssmPDouble *aWeightVector=NULL; - IssmPDouble *weightVectorTimesJac=NULL; - - /*output: */ - IssmPDouble *totalgradient=NULL; - - /*retrieve parameters: */ - num_independents = intn; - - /*if no dependents, no point in running a driver: */ - if(!(num_dependents*num_independents)) _error_("this is not allowed"); - - /*for adolc to run in parallel, we 0 out on rank~=0. But we still keep track of num_dependents:*/ - int num_dependents_old = num_dependents; - int num_independents_old = num_independents; - if(my_rank!=0){ - num_dependents = 0; - num_independents = 0; - } - - /*get the EDF pointer:*/ - ext_diff_fct *anEDF_for_solverx_p=xDynamicCast * >(femmodel->parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p; - - /* these are always needed regardless of the interpreter */ - anEDF_for_solverx_p->dp_x=xNew(anEDF_for_solverx_p->max_n); - anEDF_for_solverx_p->dp_y=xNew(anEDF_for_solverx_p->max_m); - - /* Ok, now we are going to call the fos_reverse in a loop on the index, from 0 to num_dependents, so - * as to generate num_dependents gradients: */ - totalgradient=xNewZeroInit(num_independents); - - for(int aDepIndex=0;aDepIndex(num_dependents); - if (my_rank==0) aWeightVector[aDepIndex]=1.; - - /*initialize output gradient: */ - weightVectorTimesJac=xNew(num_independents); - - /*set the forward method function pointer: */ - #ifdef _HAVE_GSL_ - anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx; - #endif - #ifdef _HAVE_MUMPS_ - anEDF_for_solverx_p->fos_reverse_iArr=fos_reverse_mumpsSolveEDF; - #endif - - anEDF_for_solverx_p->dp_U=xNew(anEDF_for_solverx_p->max_m); - anEDF_for_solverx_p->dp_Z=xNew(anEDF_for_solverx_p->max_n); - - /*call driver: */ - fos_reverse(my_rank,num_dependents,num_independents, aWeightVector, weightVectorTimesJac ); - - /*Add to totalgradient: */ - if(my_rank==0) for(int i=0;iCostFunctionx(&J,&Jtemp,NULL); - *pf = reCast(J); - _printf0_("f(x) = "<(Jtemp[i]); - Jlist[(*Jlisti)*JlistN+num_responses] = reCast(J); - xDelete(Jtemp); - - if(*indic==0){ - /*dry run, no gradient required*/ - - /*Retrieve objective functions independently*/ - _printf0_(" N/A |\n"); - for(int i=0;i(XU); - xDelete(XL); - return; - } - - /*Compute gradient*/ - for(long i=0;i=XU[index]) G[index]=0.; - if(X[index]<=XL[index]) G[index]=0.; - G[index] = G[index]*reCast(scaling_factors[c]); - X[index] = X[index]/reCast(scaling_factors[c]); - Gnorm += G[index]*G[index]; - } - } - Gnorm = sqrt(Gnorm); - - /*Print info*/ - _printf0_(" "<(XU); - xDelete(XL); - xDelete(scaling_factors); - xDelete(totalgradient); - -}/*}}}*/ -void controlm1qn3_core(FemModel* femmodel){/*{{{*/ - - /*Intermediaries*/ - long omode; - double f; - double dxmin,gttol; - int maxsteps,maxiter; - int intn,numberofvertices,num_controls,num_cost_functions,solution_type; - IssmDouble *scaling_factors = NULL; - double *X = NULL; - double *G = NULL; - - /*Recover some parameters*/ - femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); - femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum); - femmodel->parameters->FindParam(&num_cost_functions,InversionNumCostFunctionsEnum); - femmodel->parameters->FindParam(&maxsteps,InversionMaxstepsEnum); - femmodel->parameters->FindParam(&maxiter,InversionMaxiterEnum); - femmodel->parameters->FindParamAndMakePassive(&dxmin,InversionDxminEnum); - femmodel->parameters->FindParamAndMakePassive(>tol,InversionGttolEnum); - femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum); - femmodel->parameters->SetParam(false,SaveResultsEnum); - numberofvertices=femmodel->vertices->NumberOfVertices(); - - /*Initialize M1QN3 parameters*/ - if(VerboseControl())_printf0_(" Initialize M1QN3 parameters\n"); - SimulFunc simul_ptr = &simul_ad; /*Cost function address*/ - void** prosca = &euclid_; /*Dot product function (euclid is the default)*/ - char normtype[] = "dfn"; /*Norm type: dfn = scalar product defined by prosca*/ - long izs[5]; /*Arrays used by m1qn3 subroutines*/ - long iz[5]; /*Integer m1qn3 working array of size 5*/ - float rzs[1]; /*Arrays used by m1qn3 subroutines*/ - long impres = 0; /*verbosity level*/ - long imode[3] = {0}; /*scaling and starting mode, 0 by default*/ - long indic = 4; /*compute f and g*/ - long reverse = 0; /*reverse or direct mode*/ - long io = 6; /*Channel number for the output*/ - - /*Optimization criterions*/ - long niter = long(maxsteps); /*Maximum number of iterations*/ - long nsim = long(maxiter);/*Maximum number of function calls*/ - - /*Get initial guess*/ - Vector *Xpetsc = NULL; - GetPassiveVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value"); - X = Xpetsc->ToMPISerial(); - Xpetsc->GetSize(&intn); - delete Xpetsc; - _assert_(intn==numberofvertices*num_controls); - - /*Get problem dimension and initialize gradient and initial guess*/ - long n = long(intn); - G = xNew(n); - - /*Scale control for M1QN3*/ - for(int i=0;i(scaling_factors[c]); - } - } - - /*Allocate m1qn3 working arrays (see documentation)*/ - long m = 100; - long ndz = 4*n+m*(2*n+1); - double* dz = xNew(ndz); - - if(VerboseControl())_printf0_(" Computing initial solution\n"); - _printf0_("\n"); - _printf0_("Cost function f(x) | Gradient norm |g(x)| | List of contributions\n"); - _printf0_("____________________________________________________________________\n"); - - /*Prepare structure for m1qn3*/ - m1qn3_struct mystruct; - mystruct.femmodel = femmodel; - mystruct.M = maxiter; - mystruct.N = num_cost_functions+1; - mystruct.Jlist = xNewZeroInit(mystruct.M*mystruct.N); - mystruct.i = xNewZeroInit(1); - /*Initialize Gradient and cost function of M1QN3*/ - indic = 4; /*gradient required*/ - simul_ad(&indic,&n,X,&f,G,izs,rzs,(void*)&mystruct); - /*Estimation of the expected decrease in f during the first iteration*/ - double df1=f; - - /*Call M1QN3 solver*/ - m1qn3_(simul_ptr,prosca,&ctonbe_,&ctcabe_, - &n,X,&f,G,&dxmin,&df1, - >tol,normtype,&impres,&io,imode,&omode,&niter,&nsim,iz,dz,&ndz, - &reverse,&indic,izs,rzs,(void*)&mystruct); - switch(int(omode)){ - case 0: _printf0_(" Stop requested (indic = 0)\n"); break; - case 1: _printf0_(" Convergence reached (gradient satisfies stopping criterion)\n"); break; - case 2: _printf0_(" Bad initialization\n"); break; - case 3: _printf0_(" Line search failure\n"); break; - case 4: _printf0_(" Maximum number of iterations exceeded\n");break; - case 5: _printf0_(" Maximum number of function calls exceeded\n"); break; - case 6: _printf0_(" stopped on dxmin during line search\n"); break; - case 7: _printf0_(" > 0 or <0\n"); break; - default: _printf0_(" Unknown end condition\n"); - } - /*Constrain solution vector*/ - double *XL = NULL; - double *XU = NULL; - GetPassiveVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound"); - GetPassiveVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound"); - - for(int i=0;i(scaling_factors[c]); - if(X[index]>XU[index]) X[index]=XU[index]; - if(X[index](intn); - for(int i=0;i(X[i]); - SetControlInputsFromVectorx(femmodel,aX); - xDelete(aX); - - /*Set final gradient in inputs*/ - IssmDouble* aG=xNew(intn); - for(int i=0;i(G[i]); - ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG); - xDelete(aG); - - /*Add last cost function to results*/ - femmodel->OutputControlsx(&femmodel->results); - femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,JEnum,mystruct.Jlist,(*mystruct.i),mystruct.N,0,0)); - - /*Finalize*/ - if(VerboseControl()) _printf0_(" preparing final solution\n"); - femmodel->parameters->SetParam(true,SaveResultsEnum); - void (*solutioncore)(FemModel*)=NULL; - CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type); - solutioncore(femmodel); - - - /*Clean-up and return*/ - xDelete(G); - xDelete(X); - xDelete(dz); - xDelete(XU); - xDelete(XL); - xDelete(scaling_factors); - xDelete(mystruct.Jlist); - xDelete(mystruct.i); -}/*}}}*/ - -#else -void controlm1qn3_core(FemModel* femmodel){_error_("M1QN3 not installed");} -#endif //_HAVE_M1QN3_ Index: /issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp (revision 22437) +++ /issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp (revision 22438) @@ -11,319 +11,38 @@ #include "../solutionsequences/solutionsequences.h" -#if defined (_HAVE_M1QN3_) -//& !defined(_HAVE_ADOLC_) +#if defined (_HAVE_M1QN3_) & !defined(_HAVE_ADOLC_) /*m1qn3 prototypes*/ extern "C" void *ctonbe_; // DIS mode : Conversion extern "C" void *ctcabe_; // DIS mode : Conversion extern "C" void *euclid_; // Scalar product -typedef void (*SimulFunc) (long* indic,long* n, double* x,double* pf,double* g,long [],float [],void* dzs); -extern "C" void m1qn3_ (void f(long* indic,long* n, double* x,double* pf,double* g,long [],float [],void* dzs), +typedef void (*SimulFunc) (long* indic,long* n, double* x, double* pf,double* g,long [],float [],void* dzs); +extern "C" void m1qn3_ (void f(long* indic,long* n, double* x, double* pf,double* g,long [],float [],void* dzs), void **, void **, void **, - long *, double [],double *, double[], double*, double *, - double *, char [], long *, long *, long *, long *, long *, long *, long [],double [], long *, + long *, double [], double *, double [], double*, double *, + double *, char [], long *, long *, long *, long *, long *, long *, long [], double [], long *, long *, long *, long [], float [],void* ); + +/*Cost function prototype*/ +void simul(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs); /*Use struct to provide arguments*/ typedef struct{ - FemModel * femmodel; - IssmPDouble* Jlist; - int M; - int N; - int* i; + FemModel * femmodel; + IssmDouble* Jlist; + int M; + int N; + int* i; } m1qn3_struct; -/*m1qm3 functions*/ -void simul_starttrace(FemModel* femmodel){/*{{{*/ - - /*Retrive ADOLC parameters*/ - IssmDouble gcTriggerRatio; - IssmDouble gcTriggerMaxSize; - IssmDouble obufsize; - IssmDouble lbufsize; - IssmDouble cbufsize; - IssmDouble tbufsize; - femmodel->parameters->FindParam(&gcTriggerRatio,AutodiffGcTriggerRatioEnum); - femmodel->parameters->FindParam(&gcTriggerMaxSize,AutodiffGcTriggerMaxSizeEnum); - femmodel->parameters->FindParam(&obufsize,AutodiffObufsizeEnum); - femmodel->parameters->FindParam(&lbufsize,AutodiffLbufsizeEnum); - femmodel->parameters->FindParam(&cbufsize,AutodiffCbufsizeEnum); - femmodel->parameters->FindParam(&tbufsize,AutodiffTbufsizeEnum); - - /*Set garbage collection parameters: */ - setStoreManagerControl(reCast(gcTriggerRatio),reCast(gcTriggerMaxSize)); - - /*Start trace: */ - int skipFileDeletion=1; - int keepTaylors=1; - int my_rank=IssmComm::GetRank(); - trace_on(my_rank,keepTaylors,reCast(obufsize),reCast(lbufsize),reCast(cbufsize),reCast(tbufsize),skipFileDeletion); -}/*}}}*/ -void simul_ad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){/*{{{*/ - - /*Get rank*/ - int my_rank=IssmComm::GetRank(); - - /*Recover Arguments*/ - m1qn3_struct *input_struct = (m1qn3_struct*)dzs; - FemModel *femmodel = input_struct->femmodel; - IssmPDouble *Jlist = input_struct->Jlist; - int JlistM = input_struct->M; - int JlistN = input_struct->N; - int *Jlisti = input_struct->i; - int intn = (int)*n; - - /*Recover some parameters*/ - int num_responses,num_controls,numberofvertices,solution_type; - IssmDouble* scaling_factors = NULL; - femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); - femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum); - femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum); - femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); - numberofvertices=femmodel->vertices->NumberOfVertices(); - - /*Constrain input vector and update controls*/ - double *XL = NULL; - double *XU = NULL; - GetPassiveVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound"); - GetPassiveVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound"); - for(int i=0;i(scaling_factors[c]); - if(X[index]>XU[index]) X[index]=XU[index]; - if(X[index](num_controls*numberofvertices,"t"); -#else - IssmDouble* aX=xNew(num_controls*numberofvertices); -#endif - if(my_rank==0){ - for(int i=0;i(aX); - - /*Compute solution (forward)*/ - void (*solutioncore)(FemModel*)=NULL; - CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type); - solutioncore(femmodel); - - /*Get Dependents*/ - IssmDouble output_value; - int num_dependents; - IssmPDouble *dependents; - DataSet* dependent_objects=NULL; - - femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum); - femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum); - - /*Go through our dependent variables, and compute the response:*/ - dependents=xNew(num_dependents); - for(int i=0;iSize();i++){ - DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i); - dep->Responsex(&output_value,femmodel); - if (my_rank==0) { - output_value>>=dependents[i]; - } - } - - /*Turning off trace tape*/ - trace_off(); - - /*Print tape statistics so that user can kill this run if something is off already:*/ - if(VerboseAutodiff()){ /*{{{*/ - size_t tape_stats[15]; - tapestats(my_rank,tape_stats); //reading of tape statistics - int commSize=IssmComm::GetSize(); - int *sstats=new int[7]; - sstats[0]=tape_stats[NUM_OPERATIONS]; - sstats[1]=tape_stats[OP_FILE_ACCESS]; - sstats[2]=tape_stats[NUM_LOCATIONS]; - sstats[3]=tape_stats[LOC_FILE_ACCESS]; - sstats[4]=tape_stats[NUM_VALUES]; - sstats[5]=tape_stats[VAL_FILE_ACCESS]; - sstats[6]=tape_stats[TAY_STACK_SIZE]; - int *rstats=NULL; - if (my_rank==0) rstats=new int[commSize*7]; - ISSM_MPI_Gather(sstats,7,ISSM_MPI_INT,rstats,7,ISSM_MPI_INT,0,IssmComm::GetComm()); - if (my_rank==0) { - int offset=50; - int rOffset=(commSize/10)+1; - _printf_(" ADOLC statistics: \n"); - _printf_(" "<file":"") << "\n"); - _printf_(" locations: entry size " << sizeof(locint) << " Bytes\n"); - _printf_(" "<file":"") << "\n"); - _printf_(" constant values: entry size " << sizeof(double) << " Bytes\n"); - _printf_(" "<file":"") << "\n"); - _printf_(" Taylor stack: entry size " << sizeof(revreal) << " Bytes\n"); - _printf_(" "<tape_stats[TAY_BUFFER_SIZE]?" ->file":"") << "\n"); - delete []rstats; - } - delete [] sstats; - } /*}}}*/ - - /*diverse: */ - int dummy; - int num_independents=0; - - /*intermediary: */ - IssmPDouble *aWeightVector=NULL; - IssmPDouble *weightVectorTimesJac=NULL; - - /*output: */ - IssmPDouble *totalgradient=NULL; - - /*retrieve parameters: */ - num_independents = intn; - - /*if no dependents, no point in running a driver: */ - if(!(num_dependents*num_independents)) _error_("this is not allowed"); - - /*for adolc to run in parallel, we 0 out on rank~=0. But we still keep track of num_dependents:*/ - int num_dependents_old = num_dependents; - int num_independents_old = num_independents; - if(my_rank!=0){ - num_dependents = 0; - num_independents = 0; - } - - /*get the EDF pointer:*/ - ext_diff_fct *anEDF_for_solverx_p=xDynamicCast * >(femmodel->parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p; - - /* these are always needed regardless of the interpreter */ - anEDF_for_solverx_p->dp_x=xNew(anEDF_for_solverx_p->max_n); - anEDF_for_solverx_p->dp_y=xNew(anEDF_for_solverx_p->max_m); - - /* Ok, now we are going to call the fos_reverse in a loop on the index, from 0 to num_dependents, so - * as to generate num_dependents gradients: */ - totalgradient=xNewZeroInit(num_independents); - - for(int aDepIndex=0;aDepIndex(num_dependents); - if (my_rank==0) aWeightVector[aDepIndex]=1.; - - /*initialize output gradient: */ - weightVectorTimesJac=xNew(num_independents); - - /*set the forward method function pointer: */ - #ifdef _HAVE_GSL_ - anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx; - #endif - #ifdef _HAVE_MUMPS_ - anEDF_for_solverx_p->fos_reverse_iArr=fos_reverse_mumpsSolveEDF; - #endif - - anEDF_for_solverx_p->dp_U=xNew(anEDF_for_solverx_p->max_m); - anEDF_for_solverx_p->dp_Z=xNew(anEDF_for_solverx_p->max_n); - - /*call driver: */ - fos_reverse(my_rank,num_dependents,num_independents, aWeightVector, weightVectorTimesJac ); - - /*Add to totalgradient: */ - if(my_rank==0) for(int i=0;iCostFunctionx(&J,&Jtemp,NULL); - *pf = reCast(J); - _printf0_("f(x) = "<(Jtemp[i]); - Jlist[(*Jlisti)*JlistN+num_responses] = reCast(J); - xDelete(Jtemp); - - if(*indic==0){ - /*dry run, no gradient required*/ - - /*Retrieve objective functions independently*/ - _printf0_(" N/A |\n"); - for(int i=0;i(XU); - xDelete(XL); - return; - } - - /*Compute gradient*/ - for(long i=0;i=XU[index]) G[index]=0.; - if(X[index]<=XL[index]) G[index]=0.; - G[index] = G[index]*reCast(scaling_factors[c]); - X[index] = X[index]/reCast(scaling_factors[c]); - Gnorm += G[index]*G[index]; - } - } - Gnorm = sqrt(Gnorm); - - /*Print info*/ - _printf0_(" "<(XU); - xDelete(XL); - xDelete(scaling_factors); - xDelete(totalgradient); - -}/*}}}*/ -void controlm1qn3_core(FemModel* femmodel){/*{{{*/ +void controlm1qn3_core(FemModel* femmodel){ /*Intermediaries*/ long omode; - double f; - double dxmin,gttol; + double f,dxmin,gttol; int maxsteps,maxiter; int intn,numberofvertices,num_controls,num_cost_functions,solution_type; - IssmDouble *scaling_factors = NULL; - double *X = NULL; - double *G = NULL; + IssmDouble *scaling_factors = NULL; + IssmDouble *X = NULL; + IssmDouble *G = NULL; /*Recover some parameters*/ @@ -333,6 +52,6 @@ femmodel->parameters->FindParam(&maxsteps,InversionMaxstepsEnum); femmodel->parameters->FindParam(&maxiter,InversionMaxiterEnum); - femmodel->parameters->FindParamAndMakePassive(&dxmin,InversionDxminEnum); - femmodel->parameters->FindParamAndMakePassive(>tol,InversionGttolEnum); + femmodel->parameters->FindParam(&dxmin,InversionDxminEnum); + femmodel->parameters->FindParam(>tol,InversionGttolEnum); femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum); femmodel->parameters->SetParam(false,SaveResultsEnum); @@ -341,5 +60,5 @@ /*Initialize M1QN3 parameters*/ if(VerboseControl())_printf0_(" Initialize M1QN3 parameters\n"); - SimulFunc simul_ptr = &simul_ad; /*Cost function address*/ + SimulFunc costfuncion = &simul; /*Cost function address*/ void** prosca = &euclid_; /*Dot product function (euclid is the default)*/ char normtype[] = "dfn"; /*Norm type: dfn = scalar product defined by prosca*/ @@ -358,6 +77,6 @@ /*Get initial guess*/ - Vector *Xpetsc = NULL; - GetPassiveVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value"); + Vector *Xpetsc = NULL; + GetVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value"); X = Xpetsc->ToMPISerial(); Xpetsc->GetSize(&intn); @@ -373,5 +92,5 @@ for(int c=0;c(scaling_factors[c]); + X[index] = X[index]/scaling_factors[c]; } } @@ -392,17 +111,20 @@ mystruct.M = maxiter; mystruct.N = num_cost_functions+1; - mystruct.Jlist = xNewZeroInit(mystruct.M*mystruct.N); + mystruct.Jlist = xNewZeroInit(mystruct.M*mystruct.N); mystruct.i = xNewZeroInit(1); + /*Initialize Gradient and cost function of M1QN3*/ indic = 4; /*gradient required*/ - simul_ad(&indic,&n,X,&f,G,izs,rzs,(void*)&mystruct); + simul(&indic,&n,X,&f,G,izs,rzs,(void*)&mystruct); + /*Estimation of the expected decrease in f during the first iteration*/ double df1=f; /*Call M1QN3 solver*/ - m1qn3_(simul_ptr,prosca,&ctonbe_,&ctcabe_, + m1qn3_(costfuncion,prosca,&ctonbe_,&ctcabe_, &n,X,&f,G,&dxmin,&df1, >tol,normtype,&impres,&io,imode,&omode,&niter,&nsim,iz,dz,&ndz, &reverse,&indic,izs,rzs,(void*)&mystruct); + switch(int(omode)){ case 0: _printf0_(" Stop requested (indic = 0)\n"); break; @@ -416,34 +138,22 @@ default: _printf0_(" Unknown end condition\n"); } + /*Constrain solution vector*/ - double *XL = NULL; - double *XU = NULL; - GetPassiveVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound"); - GetPassiveVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound"); - - for(int i=0;i(scaling_factors[c]); + IssmDouble *XL = NULL; + IssmDouble *XU = NULL; + GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound"); + GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound"); + for(int i=0;iXU[index]) X[index]=XU[index]; if(X[index](intn); - for(int i=0;i(X[i]); - SetControlInputsFromVectorx(femmodel,aX); - xDelete(aX); - - /*Set final gradient in inputs*/ - IssmDouble* aG=xNew(intn); - for(int i=0;i(G[i]); - ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG); - xDelete(aG); - - /*Add last cost function to results*/ + SetControlInputsFromVectorx(femmodel,X); + ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G); femmodel->OutputControlsx(&femmodel->results); - femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,JEnum,mystruct.Jlist,(*mystruct.i),mystruct.N,0,0)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,JEnum,mystruct.Jlist,(*mystruct.i),mystruct.N,0,0)); /*Finalize*/ @@ -454,5 +164,4 @@ solutioncore(femmodel); - /*Clean-up and return*/ xDelete(G); @@ -461,10 +170,117 @@ xDelete(XU); xDelete(XL); + xDelete(scaling_factors); + xDelete(mystruct.Jlist); + xDelete(mystruct.i); +} + +/*Cost function definition*/ +void simul(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){ + + /*Recover Arguments*/ + m1qn3_struct *input_struct = (m1qn3_struct*)dzs; + FemModel *femmodel = input_struct->femmodel; + IssmDouble *Jlist = input_struct->Jlist; + int JlistM = input_struct->M; + int JlistN = input_struct->N; + int *Jlisti = input_struct->i; + + /*Recover some parameters*/ + int num_responses,num_controls,numberofvertices,solution_type; + IssmDouble* scaling_factors = NULL; + femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); + femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum); + femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + numberofvertices=femmodel->vertices->NumberOfVertices(); + + /*Constrain input vector and update controls*/ + IssmDouble *XL = NULL; + IssmDouble *XU = NULL; + GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound"); + GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound"); + for(int i=0;iXU[index]) X[index]=XU[index]; + if(X[index]parameters,solution_type); + solutioncore(femmodel); + + /*Check size of Jlist to avoid crashes*/ + _assert_((*Jlisti)CostFunctionx(pf,&Jtemp,NULL); + _printf0_("f(x) = "<(Jtemp); + + if(*indic==0){ + /*dry run, no gradient required*/ + + /*Retrieve objective functions independently*/ + _printf0_(" N/A |\n"); + for(int i=0;i(XU); + xDelete(XL); + return; + } + + /*Compute Adjoint*/ + AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type); + adjointcore(femmodel); + + /*Compute gradient*/ + IssmDouble* G2 = NULL; + Gradjx(&G2,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + for(long i=0;i<*n;i++) G[i] = -G2[i]; + xDelete(G2); + + /*Constrain Gradient*/ + IssmDouble Gnorm = 0.; + for(int i=0;i=XU[index]) G[index]=0.; + if(X[index]<=XL[index]) G[index]=0.; + G[index] = G[index]*scaling_factors[c]; + X[index] = X[index]/scaling_factors[c]; + Gnorm += G[index]*G[index]; + } + } + Gnorm = sqrt(Gnorm); + + /*Print info*/ + _printf0_(" "<(XU); + xDelete(XL); xDelete(scaling_factors); - xDelete(mystruct.Jlist); - xDelete(mystruct.i); -}/*}}}*/ +} #else -void controlm1qn3_core(FemModel* femmodel){_error_("M1QN3 not installed");} +void controlm1qn3_core(FemModel* femmodel){ + _error_("M1QN3 not installed"); +} #endif //_HAVE_M1QN3_ Index: /issm/trunk-jpl/src/c/cores/stressbalance_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/stressbalance_core.cpp (revision 22437) +++ /issm/trunk-jpl/src/c/cores/stressbalance_core.cpp (revision 22438) @@ -79,6 +79,6 @@ femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); } - -// if(solution_type==StressbalanceSolutionEnum)femmodel->RequestedDependentsx(); + + if(solution_type==StressbalanceSolutionEnum)femmodel->RequestedDependentsx(); /*Free ressources:*/ Index: /issm/trunk-jpl/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp (revision 22437) +++ /issm/trunk-jpl/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp (revision 22438) @@ -59,5 +59,5 @@ /*Retrieve all inputs we will be needing: */ -// Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); Input* vx_input =topelement->GetInput(VxEnum); _assert_(vx_input); Input* vxobs_input =topelement->GetInput(InversionVxObsEnum); _assert_(vxobs_input); @@ -79,5 +79,5 @@ /*Get all parameters at gaussian point*/ - //weights_input->GetInputValue(&weight,gauss,SurfaceAbsVelMisfitEnum); + weights_input->GetInputValue(&weight,gauss,SurfaceAbsVelMisfitEnum); vx_input->GetInputValue(&vx,gauss); vxobs_input->GetInputValue(&vxobs,gauss); @@ -98,7 +98,5 @@ /*Add to cost function*/ - //Jelem+=misfit*weight*Jdet*gauss->weight; - Jelem+=misfit*Jdet*gauss->weight; - + Jelem+=misfit*weight*Jdet*gauss->weight; } Index: /issm/trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp (revision 22437) +++ /issm/trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp (revision 22438) @@ -62,5 +62,5 @@ /*Retrieve all inputs we will be needing: */ -// Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); Input* vx_input =topelement->GetInput(VxEnum); _assert_(vx_input); Input* vxobs_input =topelement->GetInput(InversionVxObsEnum); _assert_(vxobs_input); @@ -82,5 +82,5 @@ /*Get all parameters at gaussian point*/ - //weights_input->GetInputValue(&weight,gauss,SurfaceLogVelMisfitEnum); + weights_input->GetInputValue(&weight,gauss,SurfaceLogVelMisfitEnum); vx_input->GetInputValue(&vx,gauss); vxobs_input->GetInputValue(&vxobs,gauss); @@ -106,5 +106,5 @@ misfit=4*pow(meanvel,2)*pow(log(velocity_mag/obs_velocity_mag),2); - weight = 1.; + /*Add to cost function*/ Jelem+=misfit*weight*Jdet*gauss->weight; Index: /issm/trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp (revision 22437) +++ /issm/trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp (revision 22438) @@ -61,5 +61,5 @@ /*Retrieve all inputs we will be needing: */ -// Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); Input* vx_input =topelement->GetInput(VxEnum); _assert_(vx_input); Input* vxobs_input =topelement->GetInput(InversionVxObsEnum); _assert_(vxobs_input); @@ -81,5 +81,5 @@ /*Get all parameters at gaussian point*/ - //weights_input->GetInputValue(&weight,gauss,SurfaceRelVelMisfitEnum); + weights_input->GetInputValue(&weight,gauss,SurfaceRelVelMisfitEnum); vx_input->GetInputValue(&vx,gauss); vxobs_input->GetInputValue(&vxobs,gauss); @@ -106,5 +106,5 @@ misfit=0.5*(scalex*pow((vx-vxobs),2)); } - weight = 1.; + /*Add to cost function*/ Jelem+=misfit*weight*Jdet*gauss->weight;