Changeset 22438


Ignore:
Timestamp:
02/21/18 10:47:55 (7 years ago)
Author:
erobo
Message:

CHG: reverting previous commit

Location:
issm/trunk-jpl/src/c
Files:
1 deleted
12 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Misfit.cpp

    r22437 r22438  
    1212
    1313#include "./classes.h"
    14 #include "./FemModel.h"
    1514#include "./ExternalResults/ExternalResult.h"
    1615#include "./ExternalResults/Results.h"
    … …  
    1817#include "./Elements/Element.h"
    1918#include "./Elements/Elements.h"
     19#include "./FemModel.h"
    2020#include "../modules/SurfaceAreax/SurfaceAreax.h"
    2121#include "../classes/Params/Parameters.h"

  • issm/trunk-jpl/src/c/classes/Misfit.h

    r22437 r22438  
    88/*Headers:*/
    99#include "./Definition.h"
     10#include "../datastructures/datastructures.h"
     11#include "./Elements/Element.h"
     12#include "./Elements/Elements.h"
    1013#include "./FemModel.h"
     14#include "../modules/SurfaceAreax/SurfaceAreax.h"
     15#include "../classes/Params/Parameters.h"
     16#include "../classes/Inputs/Input.h"
     17#include "../classes/gauss/Gauss.h"
    1118
    1219IssmDouble OutputDefinitionsResponsex(FemModel* femmodel,int output_enum);

  • issm/trunk-jpl/src/c/classes/Nodalvalue.h

    r22437 r22438  
    99/*{{{*/
    1010#include "./Definition.h"
     11#include "../datastructures/datastructures.h"
     12#include "./Elements/Element.h"
     13#include "./Elements/Elements.h"
    1114#include "./FemModel.h"
     15#include "../modules/SurfaceAreax/SurfaceAreax.h"
     16#include "../classes/Params/Parameters.h"
     17#include "../classes/Inputs/Input.h"
     18#include "../classes/gauss/Gauss.h"
    1219/*}}}*/
    1320

  • issm/trunk-jpl/src/c/classes/Numberedcostfunction.h

    r22437 r22438  
    88/*Headers:*/
    99#include "./Definition.h"
     10#include "../datastructures/datastructures.h"
     11#include "./Elements/Element.h"
     12#include "./Elements/Elements.h"
    1013#include "./FemModel.h"
     14#include "./ExternalResults/ExternalResult.h"
     15#include "./ExternalResults/Results.h"
     16#include "../modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.h"
     17#include "../modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.h"
     18#include "../modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.h"
     19#include "../modules/SurfaceLogVxVyMisfitx/SurfaceLogVxVyMisfitx.h"
     20#include "../modules/ThicknessAbsMisfitx/ThicknessAbsMisfitx.h"
     21#include "../modules/ThicknessAlongGradientx/ThicknessAlongGradientx.h"
     22#include "../modules/ThicknessAcrossGradientx/ThicknessAcrossGradientx.h"
     23#include "../modules/RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.h"
     24#include "../modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h"
    1125
    1226
    … …  
    1529class Numberedcostfunction: public Object, public Definition{
    1630
    17 public:
     31        public:
    1832
    19 int   definitionenum;
    20 char* name;
    21 int   number_cost_functions;
    22 int*  cost_functions_list;
     33                int   definitionenum;
     34                char* name;
     35                int   number_cost_functions;
     36                int*  cost_functions_list;
     37               
     38                /*Numberedcostfunction constructors, destructors :*/
     39                Numberedcostfunction();
     40                Numberedcostfunction(char* in_name, int in_definitionenum,int number_cost_functions_in,int* cost_functions_list_in);
     41                ~Numberedcostfunction();
    2342
    24 /*Numberedcostfunction constructors, destructors :*/
    25 Numberedcostfunction();
    26 Numberedcostfunction(char* in_name, int in_definitionenum,int number_cost_functions_in,int* cost_functions_list_in);
    27 ~Numberedcostfunction();
     43                /*Object virtual function resolutoin: */
     44                Object* copy();
     45                void            DeepEcho(void);
     46                void            Echo(void);
     47                int             Id(void);
     48                void            Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction);
     49                int             ObjectEnum(void);
    2850
    29 /*Object virtual function resolutoin: */
    30 Object* copy();
    31 void            DeepEcho(void);
    32 void            Echo(void);
    33 int             Id(void);
    34 void            Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction);
    35 int             ObjectEnum(void);
    36 
    37 /*Definition virtual function resolutoin: */
    38 int             DefinitionEnum();
    39 char*           Name();
    40 IssmDouble Response(FemModel* femmodel);
     51                /*Definition virtual function resolutoin: */
     52                int             DefinitionEnum();
     53                char*           Name();
     54                IssmDouble Response(FemModel* femmodel);
    4155};
    4256

  • issm/trunk-jpl/src/c/classes/Regionaloutput.cpp

    r22437 r22438  
    1515#include "./Elements/Element.h"
    1616#include "./Elements/Elements.h"
     17#include "./FemModel.h"
    1718#include "../classes/Params/Parameters.h"
    1819

  • issm/trunk-jpl/src/c/classes/Regionaloutput.h

    r22437 r22438  
    99/*{{{*/
    1010#include "./Definition.h"
     11#include "../datastructures/datastructures.h"
     12#include "./Elements/Element.h"
     13#include "./Elements/Elements.h"
    1114#include "./FemModel.h"
     15#include "../classes/Params/Parameters.h"
    1216
    1317/*}}}*/

  • issm/trunk-jpl/src/c/cores/ad_core.cpp

    r22437 r22438  
    2626        int     num_dependents=0;
    2727        int     num_independents=0;
    28         bool    isautodiff,iscontrol;
     28        bool    isautodiff       = false;
    2929        char   *driver           = NULL;
    3030        size_t  tape_stats[15];
    … …  
    3737        /*AD mode on?: */
    3838        femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
    39         femmodel->parameters->FindParam(&iscontrol,InversionIscontrolEnum);
    40 
    41         if(isautodiff && !iscontrol){
     39
     40        if(isautodiff){
    4241
    4342                #ifdef _HAVE_ADOLC_

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

    r22437 r22438  
    1111#include "../solutionsequences/solutionsequences.h"
    1212
    13 #if defined (_HAVE_M1QN3_)
    14 //& !defined(_HAVE_ADOLC_)
     13#if defined (_HAVE_M1QN3_) & !defined(_HAVE_ADOLC_)
    1514/*m1qn3 prototypes*/
    1615extern "C" void *ctonbe_; // DIS mode : Conversion
    1716extern "C" void *ctcabe_; // DIS mode : Conversion
    1817extern "C" void *euclid_; // Scalar product
    19 typedef void (*SimulFunc) (long* indic,long* n, double* x,double* pf,double* g,long [],float [],void* dzs);
    20 extern "C" void m1qn3_ (void f(long* indic,long* n, double* x,double* pf,double* g,long [],float [],void* dzs),
     18typedef void (*SimulFunc) (long* indic,long* n, double* x, double* pf,double* g,long [],float [],void* dzs);
     19extern "C" void m1qn3_ (void f(long* indic,long* n, double* x, double* pf,double* g,long [],float [],void* dzs),
    2120                        void **, void **, void **,
    22                         long *, double [],double *, double[], double*, double *,
    23                         double *, char [], long *, long *, long *, long *, long *, long *, long [],double [], long *,
     21                        long *, double [], double *, double [], double*, double *,
     22                        double *, char [], long *, long *, long *, long *, long *, long *, long [], double [], long *,
    2423                        long *, long *, long [], float [],void* );
     24
     25/*Cost function prototype*/
     26void simul(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs);
    2527
    2628/*Use struct to provide arguments*/
    2729typedef struct{
    28         FemModel   * femmodel;
    29         IssmPDouble* Jlist;
    30         int          M;
    31         int          N;
    32         int*         i;
     30        FemModel  * femmodel;
     31        IssmDouble* Jlist;
     32        int         M;
     33        int         N;
     34        int*        i;
    3335} m1qn3_struct;
    3436
    35 /*m1qm3 functions*/
    36 void simul_starttrace(FemModel* femmodel){/*{{{*/
    37 
    38         /*Retrive ADOLC parameters*/
    39         IssmDouble gcTriggerRatio;
    40         IssmDouble gcTriggerMaxSize;
    41         IssmDouble obufsize;
    42         IssmDouble lbufsize;
    43         IssmDouble cbufsize;
    44         IssmDouble tbufsize;
    45         femmodel->parameters->FindParam(&gcTriggerRatio,AutodiffGcTriggerRatioEnum);
    46         femmodel->parameters->FindParam(&gcTriggerMaxSize,AutodiffGcTriggerMaxSizeEnum);
    47         femmodel->parameters->FindParam(&obufsize,AutodiffObufsizeEnum);
    48         femmodel->parameters->FindParam(&lbufsize,AutodiffLbufsizeEnum);
    49         femmodel->parameters->FindParam(&cbufsize,AutodiffCbufsizeEnum);
    50         femmodel->parameters->FindParam(&tbufsize,AutodiffTbufsizeEnum);
    51 
    52         /*Set garbage collection parameters: */
    53         setStoreManagerControl(reCast<IssmPDouble>(gcTriggerRatio),reCast<size_t>(gcTriggerMaxSize));
    54 
    55         /*Start trace: */
    56         int skipFileDeletion=1;
    57         int keepTaylors=1;
    58         int my_rank=IssmComm::GetRank();
    59         trace_on(my_rank,keepTaylors,reCast<size_t>(obufsize),reCast<size_t>(lbufsize),reCast<size_t>(cbufsize),reCast<size_t>(tbufsize),skipFileDeletion);
    60 }/*}}}*/
    61 void simul_ad(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){/*{{{*/
    62 
    63         /*Get rank*/
    64         int my_rank=IssmComm::GetRank();
    65 
    66         /*Recover Arguments*/
    67         m1qn3_struct *input_struct = (m1qn3_struct*)dzs;
    68         FemModel     *femmodel     = input_struct->femmodel;
    69         IssmPDouble  *Jlist        = input_struct->Jlist;
    70         int           JlistM       = input_struct->M;
    71         int           JlistN       = input_struct->N;
    72         int          *Jlisti       = input_struct->i;
    73         int           intn         = (int)*n;
    74 
    75         /*Recover some parameters*/
    76         int num_responses,num_controls,numberofvertices,solution_type;
    77         IssmDouble* scaling_factors = NULL;
    78         femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    79         femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    80         femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
    81         femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    82         numberofvertices=femmodel->vertices->NumberOfVertices();
    83 
    84         /*Constrain input vector and update controls*/
    85         double  *XL = NULL;
    86         double  *XU = NULL;
    87         GetPassiveVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
    88         GetPassiveVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
    89         for(int i=0;i<numberofvertices;i++){
    90                 for(int c=0;c<num_controls;c++){
    91                         int index = num_controls*i+c;
    92                         X[index] = X[index]*reCast<double>(scaling_factors[c]);
    93                         if(X[index]>XU[index]) X[index]=XU[index];
    94                         if(X[index]<XL[index]) X[index]=XL[index];
    95                 }
    96         }
    97 
    98         /*Start Tracing*/
    99         simul_starttrace(femmodel);
    100 
    101         /*Set X as our new control input abd as INDEPENDENT!!*/
    102 #ifdef _HAVE_AD_
    103         IssmDouble* aX=xNew<IssmDouble>(num_controls*numberofvertices,"t");
    104 #else
    105         IssmDouble* aX=xNew<IssmDouble>(num_controls*numberofvertices);
    106 #endif
    107         if(my_rank==0){
    108                 for(int i=0;i<intn;i++){
    109                         aX[i]<<=X[i];
    110                 }
    111         }
    112         ISSM_MPI_Bcast(aX,intn,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    113         SetControlInputsFromVectorx(femmodel,aX);
    114         xDelete<IssmDouble>(aX);
    115 
    116         /*Compute solution (forward)*/
    117         void (*solutioncore)(FemModel*)=NULL;
    118         CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
    119         solutioncore(femmodel);
    120 
    121         /*Get Dependents*/
    122         IssmDouble  output_value;
    123         int         num_dependents;
    124         IssmPDouble *dependents;
    125         DataSet*    dependent_objects=NULL;
    126 
    127         femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
    128         femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
    129 
    130         /*Go through our dependent variables, and compute the response:*/
    131         dependents=xNew<IssmPDouble>(num_dependents);
    132         for(int i=0;i<dependent_objects->Size();i++){
    133                 DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
    134                 dep->Responsex(&output_value,femmodel);
    135                 if (my_rank==0) {
    136                         output_value>>=dependents[i];
    137                 }
    138         }
    139 
    140         /*Turning off trace tape*/
    141         trace_off();
    142 
    143         /*Print tape statistics so that user can kill this run if something is off already:*/
    144         if(VerboseAutodiff()){ /*{{{*/
    145                 size_t  tape_stats[15];
    146                 tapestats(my_rank,tape_stats); //reading of tape statistics
    147                 int commSize=IssmComm::GetSize();
    148                 int *sstats=new int[7];
    149                 sstats[0]=tape_stats[NUM_OPERATIONS];
    150                 sstats[1]=tape_stats[OP_FILE_ACCESS];
    151                 sstats[2]=tape_stats[NUM_LOCATIONS];
    152                 sstats[3]=tape_stats[LOC_FILE_ACCESS];
    153                 sstats[4]=tape_stats[NUM_VALUES];
    154                 sstats[5]=tape_stats[VAL_FILE_ACCESS];
    155                 sstats[6]=tape_stats[TAY_STACK_SIZE];
    156                 int *rstats=NULL;
    157                 if (my_rank==0) rstats=new int[commSize*7];
    158                 ISSM_MPI_Gather(sstats,7,ISSM_MPI_INT,rstats,7,ISSM_MPI_INT,0,IssmComm::GetComm());
    159                 if (my_rank==0) {
    160                         int offset=50;
    161                         int rOffset=(commSize/10)+1;
    162                         _printf_("   ADOLC statistics: \n");
    163                         _printf_("     "<<setw(offset)<<left<<"#independents: " <<setw(12)<<right<<tape_stats[NUM_INDEPENDENTS] << "\n");
    164                         _printf_("     "<<setw(offset)<<left<<"#dependents: " <<setw(12)<<right<<tape_stats[NUM_DEPENDENTS] << "\n");
    165                         _printf_("     "<<setw(offset)<<left<<"max #live active variables: " <<setw(12)<<right<<tape_stats[NUM_MAX_LIVES] << "\n");
    166                         _printf_("     operations: entry size "<< sizeof(unsigned char) << " Bytes \n");
    167                         _printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffObufsizeEnum) " <<setw(12)<<right<<tape_stats[OP_BUFFER_SIZE] << "\n");
    168                         for (int r=0;r<commSize;++r)
    169                          _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+0] << (rstats[r*7+1]?" ->file":"") << "\n");
    170                         _printf_("     locations: entry size " << sizeof(locint) << " Bytes\n");
    171                         _printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffLbufsizeEnum) " <<setw(12)<<right<<tape_stats[LOC_BUFFER_SIZE] << "\n");
    172                         for (int r=0;r<commSize;++r)
    173                          _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+2] << (rstats[r*7+3]?" ->file":"") << "\n");
    174                         _printf_("     constant values: entry size " << sizeof(double) << " Bytes\n");
    175                         _printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffCbufsizeEnum) " <<setw(12)<<right<<tape_stats[VAL_BUFFER_SIZE] << "\n");
    176                         for (int r=0;r<commSize;++r)
    177                          _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+4] << (rstats[r*7+5]?" ->file":"") << "\n");
    178                         _printf_("     Taylor stack: entry size " << sizeof(revreal) << " Bytes\n");
    179                         _printf_("     "<<setw(offset)<<left<<"  #entries in buffer (AutodiffTbufsizeEnum) " <<setw(12)<<right<<tape_stats[TAY_BUFFER_SIZE] << "\n");
    180                         for (int r=0;r<commSize;++r)
    181                          _printf_("       ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+6] << (rstats[r*7+6]>tape_stats[TAY_BUFFER_SIZE]?" ->file":"") << "\n");
    182                         delete []rstats;
    183                 }
    184                 delete [] sstats;
    185         } /*}}}*/
    186 
    187         /*diverse: */
    188         int  dummy;
    189         int  num_independents=0;
    190 
    191         /*intermediary: */
    192         IssmPDouble *aWeightVector=NULL;
    193         IssmPDouble *weightVectorTimesJac=NULL;
    194 
    195         /*output: */
    196         IssmPDouble *totalgradient=NULL;
    197 
    198         /*retrieve parameters: */
    199         num_independents = intn;
    200 
    201         /*if no dependents, no point in running a driver: */
    202         if(!(num_dependents*num_independents)) _error_("this is not allowed");
    203 
    204         /*for adolc to run in parallel, we 0 out on rank~=0. But we still keep track of num_dependents:*/
    205         int num_dependents_old   = num_dependents;
    206         int num_independents_old = num_independents;
    207         if(my_rank!=0){
    208                 num_dependents   = 0;
    209                 num_independents = 0;
    210         }
    211 
    212         /*get the EDF pointer:*/
    213         ext_diff_fct *anEDF_for_solverx_p=xDynamicCast<GenericParam<Adolc_edf> * >(femmodel->parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p;
    214 
    215         /* these are always needed regardless of the interpreter */
    216         anEDF_for_solverx_p->dp_x=xNew<double>(anEDF_for_solverx_p->max_n);
    217         anEDF_for_solverx_p->dp_y=xNew<double>(anEDF_for_solverx_p->max_m);
    218 
    219         /* Ok, now we are going to call the fos_reverse in a loop on the index, from 0 to num_dependents, so
    220          * as to generate num_dependents gradients: */
    221         totalgradient=xNewZeroInit<IssmPDouble>(num_independents);
    222 
    223         for(int aDepIndex=0;aDepIndex<num_dependents_old;aDepIndex++){
    224 
    225                 /*initialize direction index in the weights vector: */
    226                 aWeightVector=xNewZeroInit<IssmPDouble>(num_dependents);
    227                 if (my_rank==0) aWeightVector[aDepIndex]=1.;
    228 
    229                 /*initialize output gradient: */
    230                 weightVectorTimesJac=xNew<IssmPDouble>(num_independents);
    231 
    232                 /*set the forward method function pointer: */
    233                 #ifdef _HAVE_GSL_
    234                 anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx;
    235                 #endif
    236                 #ifdef _HAVE_MUMPS_
    237                 anEDF_for_solverx_p->fos_reverse_iArr=fos_reverse_mumpsSolveEDF;
    238                 #endif
    239 
    240                 anEDF_for_solverx_p->dp_U=xNew<IssmPDouble>(anEDF_for_solverx_p->max_m);
    241                 anEDF_for_solverx_p->dp_Z=xNew<IssmPDouble>(anEDF_for_solverx_p->max_n);
    242 
    243                 /*call driver: */
    244                 fos_reverse(my_rank,num_dependents,num_independents, aWeightVector, weightVectorTimesJac );
    245 
    246                 /*Add to totalgradient: */
    247                 if(my_rank==0) for(int i=0;i<num_independents;i++) totalgradient[i]+=weightVectorTimesJac[i];
    248 
    249                 /*free resources :*/
    250                 xDelete(weightVectorTimesJac);
    251                 xDelete(aWeightVector);
    252         }
    253 
    254         /*Broadcast gradient to other ranks*/
    255         ISSM_MPI_Bcast(totalgradient,num_independents_old,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
    256 
    257         /*Check size of Jlist to avoid crashes*/
    258         _assert_((*Jlisti)<JlistM);
    259         _assert_(JlistN==num_responses+1);
    260 
    261         /*Compute objective function*/
    262         IssmDouble* Jtemp = NULL;
    263         IssmDouble J;
    264         femmodel->CostFunctionx(&J,&Jtemp,NULL);
    265         *pf = reCast<double>(J);
    266         _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
    267 
    268         /*Record cost function values and delete Jtemp*/
    269         for(int i=0;i<num_responses;i++) Jlist[(*Jlisti)*JlistN+i] = reCast<IssmPDouble>(Jtemp[i]);
    270         Jlist[(*Jlisti)*JlistN+num_responses] = reCast<IssmPDouble>(J);
    271         xDelete<IssmDouble>(Jtemp);
    272 
    273         if(*indic==0){
    274                 /*dry run, no gradient required*/
    275 
    276                 /*Retrieve objective functions independently*/
    277                 _printf0_("            N/A |\n");
    278                 for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
    279                 _printf0_("\n");
    280 
    281                 *Jlisti = (*Jlisti) +1;
    282                 xDelete<double>(XU);
    283                 xDelete<double>(XL);
    284                 return;
    285         }
    286 
    287         /*Compute gradient*/
    288         for(long i=0;i<num_independents_old;i++) G[i] = totalgradient[i];
    289 
    290         /*Constrain Gradient*/
    291         IssmDouble  Gnorm = 0.;
    292         for(int i=0;i<numberofvertices;i++){
    293                 for(int c=0;c<num_controls;c++){
    294                         int index = num_controls*i+c;
    295                         if(X[index]>=XU[index]) G[index]=0.;
    296                         if(X[index]<=XL[index]) G[index]=0.;
    297                         G[index] = G[index]*reCast<double>(scaling_factors[c]);
    298                         X[index] = X[index]/reCast<double>(scaling_factors[c]);
    299                         Gnorm += G[index]*G[index];
    300                 }
    301         }
    302         Gnorm = sqrt(Gnorm);
    303 
    304         /*Print info*/
    305         _printf0_("       "<<setw(12)<<setprecision(7)<<Gnorm<<" |");
    306         for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
    307         _printf0_("\n");
    308 
    309         /*Clean-up and return*/
    310         *Jlisti = (*Jlisti) +1;
    311         xDelete<double>(XU);
    312         xDelete<double>(XL);
    313         xDelete<IssmDouble>(scaling_factors);
    314         xDelete<IssmPDouble>(totalgradient);
    315 
    316 }/*}}}*/
    317 void controlm1qn3_core(FemModel* femmodel){/*{{{*/
     37void controlm1qn3_core(FemModel* femmodel){
    31838
    31939        /*Intermediaries*/
    32040        long         omode;
    321         double           f;
    322         double           dxmin,gttol;
     41        double       f,dxmin,gttol;
    32342        int          maxsteps,maxiter;
    32443        int          intn,numberofvertices,num_controls,num_cost_functions,solution_type;
    325         IssmDouble      *scaling_factors = NULL;
    326         double      *X  = NULL;
    327         double      *G  = NULL;
     44        IssmDouble  *scaling_factors = NULL;
     45        IssmDouble  *X  = NULL;
     46        IssmDouble  *G  = NULL;
    32847
    32948        /*Recover some parameters*/
    … …  
    33352        femmodel->parameters->FindParam(&maxsteps,InversionMaxstepsEnum);
    33453        femmodel->parameters->FindParam(&maxiter,InversionMaxiterEnum);
    335         femmodel->parameters->FindParamAndMakePassive(&dxmin,InversionDxminEnum);
    336         femmodel->parameters->FindParamAndMakePassive(&gttol,InversionGttolEnum);
     54        femmodel->parameters->FindParam(&dxmin,InversionDxminEnum);
     55        femmodel->parameters->FindParam(&gttol,InversionGttolEnum);
    33756        femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
    33857        femmodel->parameters->SetParam(false,SaveResultsEnum);
    … …  
    34160        /*Initialize M1QN3 parameters*/
    34261        if(VerboseControl())_printf0_("   Initialize M1QN3 parameters\n");
    343         SimulFunc simul_ptr    = &simul_ad; /*Cost function address*/
     62        SimulFunc costfuncion  = &simul;    /*Cost function address*/
    34463        void**    prosca       = &euclid_;  /*Dot product function (euclid is the default)*/
    34564        char      normtype[]   = "dfn";     /*Norm type: dfn = scalar product defined by prosca*/
    … …  
    35877
    35978        /*Get initial guess*/
    360         Vector<double> *Xpetsc = NULL;
    361         GetPassiveVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
     79        Vector<IssmDouble> *Xpetsc = NULL;
     80        GetVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
    36281        X = Xpetsc->ToMPISerial();
    36382        Xpetsc->GetSize(&intn);
    … …  
    37392                for(int c=0;c<num_controls;c++){
    37493                        int index = num_controls*i+c;
    375                         X[index] = X[index]/reCast<IssmPDouble>(scaling_factors[c]);
     94                        X[index] = X[index]/scaling_factors[c];
    37695                }
    37796        }
    … …  
    392111        mystruct.M        = maxiter;
    393112        mystruct.N        = num_cost_functions+1;
    394         mystruct.Jlist    = xNewZeroInit<IssmPDouble>(mystruct.M*mystruct.N);
     113        mystruct.Jlist    = xNewZeroInit<IssmDouble>(mystruct.M*mystruct.N);
    395114        mystruct.i        = xNewZeroInit<int>(1);
     115
    396116        /*Initialize Gradient and cost function of M1QN3*/
    397117        indic = 4; /*gradient required*/
    398         simul_ad(&indic,&n,X,&f,G,izs,rzs,(void*)&mystruct);
     118        simul(&indic,&n,X,&f,G,izs,rzs,(void*)&mystruct);
     119
    399120        /*Estimation of the expected decrease in f during the first iteration*/
    400121        double df1=f;
    401122
    402123        /*Call M1QN3 solver*/
    403         m1qn3_(simul_ptr,prosca,&ctonbe_,&ctcabe_,
     124        m1qn3_(costfuncion,prosca,&ctonbe_,&ctcabe_,
    404125                                &n,X,&f,G,&dxmin,&df1,
    405126                                &gttol,normtype,&impres,&io,imode,&omode,&niter,&nsim,iz,dz,&ndz,
    406127                                &reverse,&indic,izs,rzs,(void*)&mystruct);
     128
    407129        switch(int(omode)){
    408130                case 0:  _printf0_("   Stop requested (indic = 0)\n"); break;
    … …  
    416138                default: _printf0_("   Unknown end condition\n");
    417139        }
     140
    418141        /*Constrain solution vector*/
    419         double  *XL = NULL;
    420         double  *XU = NULL;
    421         GetPassiveVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
    422         GetPassiveVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
    423 
    424         for(int i=0;i<numberofvertices;i++){
    425                 for(int c=0;c<num_controls;c++){
    426                         int index = num_controls*i+c;
    427                         X[index] = X[index]*reCast<double>(scaling_factors[c]);
     142        IssmDouble  *XL = NULL;
     143        IssmDouble  *XU = NULL;
     144        GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
     145        GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
     146        for(int i=0;i<numberofvertices;i++){
     147                for(int c=0;c<num_controls;c++){
     148                        int index = num_controls*i+c;
     149                        X[index] = X[index]*scaling_factors[c];
    428150                        if(X[index]>XU[index]) X[index]=XU[index];
    429151                        if(X[index]<XL[index]) X[index]=XL[index];
    430152                }
    431153        }
    432                
    433         /*Set X as our new control*/
    434         IssmDouble* aX=xNew<IssmDouble>(intn);
    435         for(int i=0;i<intn;i++) aX[i] = reCast<IssmDouble>(X[i]);
    436         SetControlInputsFromVectorx(femmodel,aX);
    437         xDelete(aX);
    438 
    439         /*Set final gradient in inputs*/
    440         IssmDouble* aG=xNew<IssmDouble>(intn);
    441         for(int i=0;i<intn;i++) aG[i] = reCast<IssmDouble>(G[i]);
    442         ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG);
    443         xDelete(aG);
    444 
    445         /*Add last cost function to results*/
     154        SetControlInputsFromVectorx(femmodel,X);
     155        ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
    446156        femmodel->OutputControlsx(&femmodel->results);
    447         femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,mystruct.Jlist,(*mystruct.i),mystruct.N,0,0));
     157        femmodel->results->AddObject(new GenericExternalResult<double*>(femmodel->results->Size()+1,JEnum,mystruct.Jlist,(*mystruct.i),mystruct.N,0,0));
    448158
    449159        /*Finalize*/
    … …  
    454164        solutioncore(femmodel);
    455165
    456 
    457166        /*Clean-up and return*/
    458167        xDelete<double>(G);
    … …  
    461170        xDelete<double>(XU);
    462171        xDelete<double>(XL);
     172        xDelete<double>(scaling_factors);
     173        xDelete<double>(mystruct.Jlist);
     174        xDelete<int>(mystruct.i);
     175}
     176
     177/*Cost function definition*/
     178void simul(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs){
     179
     180        /*Recover Arguments*/
     181        m1qn3_struct *input_struct = (m1qn3_struct*)dzs;
     182        FemModel     *femmodel     = input_struct->femmodel;
     183        IssmDouble   *Jlist        = input_struct->Jlist;
     184        int           JlistM       = input_struct->M;
     185        int           JlistN       = input_struct->N;
     186        int          *Jlisti       = input_struct->i;
     187
     188        /*Recover some parameters*/
     189        int num_responses,num_controls,numberofvertices,solution_type;
     190        IssmDouble* scaling_factors = NULL;
     191        femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
     192        femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
     193        femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
     194        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
     195        numberofvertices=femmodel->vertices->NumberOfVertices();
     196
     197        /*Constrain input vector and update controls*/
     198        IssmDouble  *XL = NULL;
     199        IssmDouble  *XU = NULL;
     200        GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
     201        GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
     202        for(int i=0;i<numberofvertices;i++){
     203                for(int c=0;c<num_controls;c++){
     204                        int index = num_controls*i+c;
     205                        X[index] = X[index]*scaling_factors[c];
     206                        if(X[index]>XU[index]) X[index]=XU[index];
     207                        if(X[index]<XL[index]) X[index]=XL[index];
     208                }
     209        }
     210        SetControlInputsFromVectorx(femmodel,X);
     211
     212        /*Compute solution and adjoint*/
     213        void (*solutioncore)(FemModel*)=NULL;
     214        void (*adjointcore)(FemModel*)=NULL;
     215        CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
     216        solutioncore(femmodel);
     217
     218        /*Check size of Jlist to avoid crashes*/
     219        _assert_((*Jlisti)<JlistM);
     220        _assert_(JlistN==num_responses+1);
     221
     222        /*Compute objective function*/
     223        IssmDouble* Jtemp = NULL;
     224        femmodel->CostFunctionx(pf,&Jtemp,NULL);
     225        _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
     226
     227        /*Record cost function values and delete Jtemp*/
     228        for(int i=0;i<num_responses;i++) Jlist[(*Jlisti)*JlistN+i] = Jtemp[i];
     229        Jlist[(*Jlisti)*JlistN+num_responses] = *pf;
     230        xDelete<IssmDouble>(Jtemp);
     231
     232        if(*indic==0){
     233                /*dry run, no gradient required*/
     234
     235                /*Retrieve objective functions independently*/
     236                _printf0_("            N/A |\n");
     237                for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
     238                _printf0_("\n");
     239
     240                *Jlisti = (*Jlisti) +1;
     241                xDelete<IssmDouble>(XU);
     242                xDelete<IssmDouble>(XL);
     243                return;
     244        }
     245
     246        /*Compute Adjoint*/
     247        AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type);
     248        adjointcore(femmodel);
     249
     250        /*Compute gradient*/
     251        IssmDouble* G2 = NULL;
     252        Gradjx(&G2,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     253        for(long i=0;i<*n;i++) G[i] = -G2[i];
     254        xDelete<IssmDouble>(G2);
     255
     256        /*Constrain Gradient*/
     257        IssmDouble  Gnorm = 0.;
     258        for(int i=0;i<numberofvertices;i++){
     259                for(int c=0;c<num_controls;c++){
     260                        int index = num_controls*i+c;
     261                        if(X[index]>=XU[index]) G[index]=0.;
     262                        if(X[index]<=XL[index]) G[index]=0.;
     263                        G[index] = G[index]*scaling_factors[c];
     264                        X[index] = X[index]/scaling_factors[c];
     265                        Gnorm += G[index]*G[index];
     266                }
     267        }
     268        Gnorm = sqrt(Gnorm);
     269
     270        /*Print info*/
     271        _printf0_("       "<<setw(12)<<setprecision(7)<<Gnorm<<" |");
     272        for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[(*Jlisti)*JlistN+i]);
     273        _printf0_("\n");
     274
     275        /*Clean-up and return*/
     276        *Jlisti = (*Jlisti) +1;
     277        xDelete<IssmDouble>(XU);
     278        xDelete<IssmDouble>(XL);
    463279        xDelete<IssmDouble>(scaling_factors);
    464         xDelete<IssmPDouble>(mystruct.Jlist);
    465         xDelete<int>(mystruct.i);
    466 }/*}}}*/
     280}
    467281
    468282#else
    469 void controlm1qn3_core(FemModel* femmodel){_error_("M1QN3 not installed");}
     283void controlm1qn3_core(FemModel* femmodel){
     284        _error_("M1QN3 not installed");
     285}
    470286#endif //_HAVE_M1QN3_

  • issm/trunk-jpl/src/c/cores/stressbalance_core.cpp

    r22437 r22438  
    7979                femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
    8080        }
    81        
    82 //      if(solution_type==StressbalanceSolutionEnum)femmodel->RequestedDependentsx();
     81
     82        if(solution_type==StressbalanceSolutionEnum)femmodel->RequestedDependentsx();
    8383
    8484        /*Free ressources:*/   

  • issm/trunk-jpl/src/c/modules/SurfaceAbsVelMisfitx/SurfaceAbsVelMisfitx.cpp

    r22437 r22438  
    5959
    6060        /*Retrieve all inputs we will be needing: */
    61 //      Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     61        Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    6262        Input* vx_input     =topelement->GetInput(VxEnum);                                 _assert_(vx_input);
    6363        Input* vxobs_input  =topelement->GetInput(InversionVxObsEnum);                     _assert_(vxobs_input);
    … …  
    7979
    8080                /*Get all parameters at gaussian point*/
    81                 //weights_input->GetInputValue(&weight,gauss,SurfaceAbsVelMisfitEnum);
     81                weights_input->GetInputValue(&weight,gauss,SurfaceAbsVelMisfitEnum);
    8282                vx_input->GetInputValue(&vx,gauss);
    8383                vxobs_input->GetInputValue(&vxobs,gauss);
    … …  
    9898
    9999                /*Add to cost function*/
    100                 //Jelem+=misfit*weight*Jdet*gauss->weight;
    101                 Jelem+=misfit*Jdet*gauss->weight;
    102 
     100                Jelem+=misfit*weight*Jdet*gauss->weight;
    103101        }
    104102

  • issm/trunk-jpl/src/c/modules/SurfaceLogVelMisfitx/SurfaceLogVelMisfitx.cpp

    r22437 r22438  
    6262
    6363        /*Retrieve all inputs we will be needing: */
    64 //      Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     64        Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    6565        Input* vx_input     =topelement->GetInput(VxEnum);                                 _assert_(vx_input);
    6666        Input* vxobs_input  =topelement->GetInput(InversionVxObsEnum);                     _assert_(vxobs_input);
    … …  
    8282
    8383                /*Get all parameters at gaussian point*/
    84                 //weights_input->GetInputValue(&weight,gauss,SurfaceLogVelMisfitEnum);
     84                weights_input->GetInputValue(&weight,gauss,SurfaceLogVelMisfitEnum);
    8585                vx_input->GetInputValue(&vx,gauss);
    8686                vxobs_input->GetInputValue(&vxobs,gauss);
    … …  
    106106
    107107                misfit=4*pow(meanvel,2)*pow(log(velocity_mag/obs_velocity_mag),2);
    108                 weight = 1.;
     108
    109109                /*Add to cost function*/
    110110                Jelem+=misfit*weight*Jdet*gauss->weight;

  • issm/trunk-jpl/src/c/modules/SurfaceRelVelMisfitx/SurfaceRelVelMisfitx.cpp

    r22437 r22438  
    6161
    6262        /*Retrieve all inputs we will be needing: */
    63 //      Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
     63        Input* weights_input=topelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    6464        Input* vx_input     =topelement->GetInput(VxEnum);                                 _assert_(vx_input);
    6565        Input* vxobs_input  =topelement->GetInput(InversionVxObsEnum);                     _assert_(vxobs_input);
    … …  
    8181
    8282                /*Get all parameters at gaussian point*/
    83                 //weights_input->GetInputValue(&weight,gauss,SurfaceRelVelMisfitEnum);
     83                weights_input->GetInputValue(&weight,gauss,SurfaceRelVelMisfitEnum);
    8484                vx_input->GetInputValue(&vx,gauss);
    8585                vxobs_input->GetInputValue(&vxobs,gauss);
    … …  
    106106                        misfit=0.5*(scalex*pow((vx-vxobs),2));
    107107                }
    108                 weight = 1.;
     108
    109109                /*Add to cost function*/
    110110                Jelem+=misfit*weight*Jdet*gauss->weight;
Note: See TracChangeset for help on using the changeset viewer.