Index: /issm/trunk-jpl/src/c/classes/Misfit.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Misfit.cpp (revision 22436) +++ /issm/trunk-jpl/src/c/classes/Misfit.cpp (revision 22437) @@ -12,4 +12,5 @@ #include "./classes.h" +#include "./FemModel.h" #include "./ExternalResults/ExternalResult.h" #include "./ExternalResults/Results.h" @@ -17,5 +18,4 @@ #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 22436) +++ /issm/trunk-jpl/src/c/classes/Misfit.h (revision 22437) @@ -8,12 +8,5 @@ /*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 22436) +++ /issm/trunk-jpl/src/c/classes/Nodalvalue.h (revision 22437) @@ -9,12 +9,5 @@ /*{{{*/ #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 22436) +++ /issm/trunk-jpl/src/c/classes/Numberedcostfunction.h (revision 22437) @@ -8,19 +8,5 @@ /*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" @@ -29,28 +15,28 @@ class Numberedcostfunction: public Object, public Definition{ - public: +public: - 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(); +int definitionenum; +char* name; +int number_cost_functions; +int* cost_functions_list; - /*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); +/*Numberedcostfunction constructors, destructors :*/ +Numberedcostfunction(); +Numberedcostfunction(char* in_name, int in_definitionenum,int number_cost_functions_in,int* cost_functions_list_in); +~Numberedcostfunction(); - /*Definition virtual function resolutoin: */ - int DefinitionEnum(); - char* Name(); - IssmDouble Response(FemModel* femmodel); +/*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); }; Index: /issm/trunk-jpl/src/c/classes/Regionaloutput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Regionaloutput.cpp (revision 22436) +++ /issm/trunk-jpl/src/c/classes/Regionaloutput.cpp (revision 22437) @@ -15,5 +15,4 @@ #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 22436) +++ /issm/trunk-jpl/src/c/classes/Regionaloutput.h (revision 22437) @@ -9,9 +9,5 @@ /*{{{*/ #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 22436) +++ /issm/trunk-jpl/src/c/cores/ad_core.cpp (revision 22437) @@ -26,5 +26,5 @@ int num_dependents=0; int num_independents=0; - bool isautodiff = false; + bool isautodiff,iscontrol; char *driver = NULL; size_t tape_stats[15]; @@ -37,6 +37,7 @@ /*AD mode on?: */ femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum); - - if(isautodiff){ + femmodel->parameters->FindParam(&iscontrol,InversionIscontrolEnum); + + if(isautodiff && !iscontrol){ #ifdef _HAVE_ADOLC_ Index: /issm/trunk-jpl/src/c/cores/controlADm1qn3_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/controlADm1qn3_core.cpp (revision 22437) +++ /issm/trunk-jpl/src/c/cores/controlADm1qn3_core.cpp (revision 22437) @@ -0,0 +1,470 @@ +/*!\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 22436) +++ /issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp (revision 22437) @@ -11,38 +11,319 @@ #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; - IssmDouble* Jlist; - int M; - int N; - int* i; + FemModel * femmodel; + IssmPDouble* Jlist; + int M; + int N; + int* i; } m1qn3_struct; -void controlm1qn3_core(FemModel* femmodel){ +/*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,dxmin,gttol; + double f; + double dxmin,gttol; int maxsteps,maxiter; int intn,numberofvertices,num_controls,num_cost_functions,solution_type; - IssmDouble *scaling_factors = NULL; - IssmDouble *X = NULL; - IssmDouble *G = NULL; + IssmDouble *scaling_factors = NULL; + double *X = NULL; + double *G = NULL; /*Recover some parameters*/ @@ -52,6 +333,6 @@ femmodel->parameters->FindParam(&maxsteps,InversionMaxstepsEnum); femmodel->parameters->FindParam(&maxiter,InversionMaxiterEnum); - femmodel->parameters->FindParam(&dxmin,InversionDxminEnum); - femmodel->parameters->FindParam(>tol,InversionGttolEnum); + femmodel->parameters->FindParamAndMakePassive(&dxmin,InversionDxminEnum); + femmodel->parameters->FindParamAndMakePassive(>tol,InversionGttolEnum); femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum); femmodel->parameters->SetParam(false,SaveResultsEnum); @@ -60,5 +341,5 @@ /*Initialize M1QN3 parameters*/ if(VerboseControl())_printf0_(" Initialize M1QN3 parameters\n"); - SimulFunc costfuncion = &simul; /*Cost function address*/ + 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*/ @@ -77,6 +358,6 @@ /*Get initial guess*/ - Vector *Xpetsc = NULL; - GetVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value"); + Vector *Xpetsc = NULL; + GetPassiveVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value"); X = Xpetsc->ToMPISerial(); Xpetsc->GetSize(&intn); @@ -92,5 +373,5 @@ for(int c=0;c(scaling_factors[c]); } } @@ -111,20 +392,17 @@ 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(&indic,&n,X,&f,G,izs,rzs,(void*)&mystruct); - + 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_(costfuncion,prosca,&ctonbe_,&ctcabe_, + 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; @@ -138,22 +416,34 @@ default: _printf0_(" Unknown end condition\n"); } - /*Constrain solution vector*/ - 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"); + 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]elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G); + + /*Set X as our new control*/ + IssmDouble* aX=xNew(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)); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,JEnum,mystruct.Jlist,(*mystruct.i),mystruct.N,0,0)); /*Finalize*/ @@ -164,4 +454,5 @@ solutioncore(femmodel); + /*Clean-up and return*/ xDelete(G); @@ -170,117 +461,10 @@ xDelete(XU); xDelete(XL); - xDelete(scaling_factors); - xDelete(mystruct.Jlist); + 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); -} +}/*}}}*/ #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 22436) +++ /issm/trunk-jpl/src/c/cores/stressbalance_core.cpp (revision 22437) @@ -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 22436) +++ /issm/trunk-jpl/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp (revision 22437) @@ -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,5 +98,7 @@ /*Add to cost function*/ - Jelem+=misfit*weight*Jdet*gauss->weight; + //Jelem+=misfit*weight*Jdet*gauss->weight; + Jelem+=misfit*Jdet*gauss->weight; + } Index: /issm/trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp (revision 22436) +++ /issm/trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp (revision 22437) @@ -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 22436) +++ /issm/trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp (revision 22437) @@ -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;