Changeset 22437


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

ADD: allows AD to be used with m1qn3 core, not compliled in now (dead code)

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

Legend:

Unmodified
Added
Removed

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

    r22427 r22437  
    1212
    1313#include "./classes.h"
     14#include "./FemModel.h"
    1415#include "./ExternalResults/ExternalResult.h"
    1516#include "./ExternalResults/Results.h"
    … …  
    1718#include "./Elements/Element.h"
    1819#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

    r22427 r22437  
    88/*Headers:*/
    99#include "./Definition.h"
    10 #include "../datastructures/datastructures.h"
    11 #include "./Elements/Element.h"
    12 #include "./Elements/Elements.h"
    1310#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"
    1811
    1912IssmDouble OutputDefinitionsResponsex(FemModel* femmodel,int output_enum);

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

    r22427 r22437  
    99/*{{{*/
    1010#include "./Definition.h"
    11 #include "../datastructures/datastructures.h"
    12 #include "./Elements/Element.h"
    13 #include "./Elements/Elements.h"
    1411#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"
    1912/*}}}*/
    2013

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

    r22422 r22437  
    88/*Headers:*/
    99#include "./Definition.h"
    10 #include "../datastructures/datastructures.h"
    11 #include "./Elements/Element.h"
    12 #include "./Elements/Elements.h"
    1310#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"
    2511
    2612
    … …  
    2915class Numberedcostfunction: public Object, public Definition{
    3016
    31         public:
     17public:
    3218
    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();
     19int   definitionenum;
     20char* name;
     21int   number_cost_functions;
     22int*  cost_functions_list;
    4223
    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);
     24/*Numberedcostfunction constructors, destructors :*/
     25Numberedcostfunction();
     26Numberedcostfunction(char* in_name, int in_definitionenum,int number_cost_functions_in,int* cost_functions_list_in);
     27~Numberedcostfunction();
    5028
    51                 /*Definition virtual function resolutoin: */
    52                 int             DefinitionEnum();
    53                 char*           Name();
    54                 IssmDouble Response(FemModel* femmodel);
     29/*Object virtual function resolutoin: */
     30Object* copy();
     31void            DeepEcho(void);
     32void            Echo(void);
     33int             Id(void);
     34void            Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction);
     35int             ObjectEnum(void);
     36
     37/*Definition virtual function resolutoin: */
     38int             DefinitionEnum();
     39char*           Name();
     40IssmDouble Response(FemModel* femmodel);
    5541};
    5642

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

    r22427 r22437  
    1515#include "./Elements/Element.h"
    1616#include "./Elements/Elements.h"
    17 #include "./FemModel.h"
    1817#include "../classes/Params/Parameters.h"
    1918

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

    r22427 r22437  
    99/*{{{*/
    1010#include "./Definition.h"
    11 #include "../datastructures/datastructures.h"
    12 #include "./Elements/Element.h"
    13 #include "./Elements/Elements.h"
    1411#include "./FemModel.h"
    15 #include "../classes/Params/Parameters.h"
    1612
    1713/*}}}*/

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

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

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

    r21450 r22437  
    1111#include "../solutionsequences/solutionsequences.h"
    1212
    13 #if defined (_HAVE_M1QN3_) & !defined(_HAVE_ADOLC_)
     13#if defined (_HAVE_M1QN3_)
     14//& !defined(_HAVE_ADOLC_)
    1415/*m1qn3 prototypes*/
    1516extern "C" void *ctonbe_; // DIS mode : Conversion
    1617extern "C" void *ctcabe_; // DIS mode : Conversion
    1718extern "C" void *euclid_; // Scalar product
    18 typedef void (*SimulFunc) (long* indic,long* n, double* x, double* pf,double* g,long [],float [],void* dzs);
    19 extern "C" void m1qn3_ (void f(long* indic,long* n, double* x, double* pf,double* g,long [],float [],void* dzs),
     19typedef void (*SimulFunc) (long* indic,long* n, double* x,double* pf,double* g,long [],float [],void* dzs);
     20extern "C" void m1qn3_ (void f(long* indic,long* n, double* x,double* pf,double* g,long [],float [],void* dzs),
    2021                        void **, void **, void **,
    21                         long *, double [], double *, double [], double*, double *,
    22                         double *, char [], long *, long *, long *, long *, long *, long *, long [], double [], long *,
     22                        long *, double [],double *, double[], double*, double *,
     23                        double *, char [], long *, long *, long *, long *, long *, long *, long [],double [], long *,
    2324                        long *, long *, long [], float [],void* );
    24 
    25 /*Cost function prototype*/
    26 void simul(long* indic,long* n,double* X,double* pf,double* G,long izs[1],float rzs[1],void* dzs);
    2725
    2826/*Use struct to provide arguments*/
    2927typedef struct{
    30         FemModel  * femmodel;
    31         IssmDouble* Jlist;
    32         int         M;
    33         int         N;
    34         int*        i;
     28        FemModel   * femmodel;
     29        IssmPDouble* Jlist;
     30        int          M;
     31        int          N;
     32        int*         i;
    3533} m1qn3_struct;
    3634
    37 void controlm1qn3_core(FemModel* femmodel){
     35/*m1qm3 functions*/
     36void 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}/*}}}*/
     61void 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}/*}}}*/
     317void controlm1qn3_core(FemModel* femmodel){/*{{{*/
    38318
    39319        /*Intermediaries*/
    40320        long         omode;
    41         double       f,dxmin,gttol;
     321        double           f;
     322        double           dxmin,gttol;
    42323        int          maxsteps,maxiter;
    43324        int          intn,numberofvertices,num_controls,num_cost_functions,solution_type;
    44         IssmDouble  *scaling_factors = NULL;
    45         IssmDouble  *X  = NULL;
    46         IssmDouble  *G  = NULL;
     325        IssmDouble      *scaling_factors = NULL;
     326        double      *X  = NULL;
     327        double      *G  = NULL;
    47328
    48329        /*Recover some parameters*/
    … …  
    52333        femmodel->parameters->FindParam(&maxsteps,InversionMaxstepsEnum);
    53334        femmodel->parameters->FindParam(&maxiter,InversionMaxiterEnum);
    54         femmodel->parameters->FindParam(&dxmin,InversionDxminEnum);
    55         femmodel->parameters->FindParam(&gttol,InversionGttolEnum);
     335        femmodel->parameters->FindParamAndMakePassive(&dxmin,InversionDxminEnum);
     336        femmodel->parameters->FindParamAndMakePassive(&gttol,InversionGttolEnum);
    56337        femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
    57338        femmodel->parameters->SetParam(false,SaveResultsEnum);
    … …  
    60341        /*Initialize M1QN3 parameters*/
    61342        if(VerboseControl())_printf0_("   Initialize M1QN3 parameters\n");
    62         SimulFunc costfuncion  = &simul;    /*Cost function address*/
     343        SimulFunc simul_ptr    = &simul_ad; /*Cost function address*/
    63344        void**    prosca       = &euclid_;  /*Dot product function (euclid is the default)*/
    64345        char      normtype[]   = "dfn";     /*Norm type: dfn = scalar product defined by prosca*/
    … …  
    77358
    78359        /*Get initial guess*/
    79         Vector<IssmDouble> *Xpetsc = NULL;
    80         GetVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
     360        Vector<double> *Xpetsc = NULL;
     361        GetPassiveVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
    81362        X = Xpetsc->ToMPISerial();
    82363        Xpetsc->GetSize(&intn);
    … …  
    92373                for(int c=0;c<num_controls;c++){
    93374                        int index = num_controls*i+c;
    94                         X[index] = X[index]/scaling_factors[c];
     375                        X[index] = X[index]/reCast<IssmPDouble>(scaling_factors[c]);
    95376                }
    96377        }
    … …  
    111392        mystruct.M        = maxiter;
    112393        mystruct.N        = num_cost_functions+1;
    113         mystruct.Jlist    = xNewZeroInit<IssmDouble>(mystruct.M*mystruct.N);
     394        mystruct.Jlist    = xNewZeroInit<IssmPDouble>(mystruct.M*mystruct.N);
    114395        mystruct.i        = xNewZeroInit<int>(1);
    115 
    116396        /*Initialize Gradient and cost function of M1QN3*/
    117397        indic = 4; /*gradient required*/
    118         simul(&indic,&n,X,&f,G,izs,rzs,(void*)&mystruct);
    119 
     398        simul_ad(&indic,&n,X,&f,G,izs,rzs,(void*)&mystruct);
    120399        /*Estimation of the expected decrease in f during the first iteration*/
    121400        double df1=f;
    122401
    123402        /*Call M1QN3 solver*/
    124         m1qn3_(costfuncion,prosca,&ctonbe_,&ctcabe_,
     403        m1qn3_(simul_ptr,prosca,&ctonbe_,&ctcabe_,
    125404                                &n,X,&f,G,&dxmin,&df1,
    126405                                &gttol,normtype,&impres,&io,imode,&omode,&niter,&nsim,iz,dz,&ndz,
    127406                                &reverse,&indic,izs,rzs,(void*)&mystruct);
    128 
    129407        switch(int(omode)){
    130408                case 0:  _printf0_("   Stop requested (indic = 0)\n"); break;
    … …  
    138416                default: _printf0_("   Unknown end condition\n");
    139417        }
    140 
    141418        /*Constrain solution vector*/
    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");
     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
    146424        for(int i=0;i<numberofvertices;i++){
    147425                for(int c=0;c<num_controls;c++){
    148426                        int index = num_controls*i+c;
    149                         X[index] = X[index]*scaling_factors[c];
     427                        X[index] = X[index]*reCast<double>(scaling_factors[c]);
    150428                        if(X[index]>XU[index]) X[index]=XU[index];
    151429                        if(X[index]<XL[index]) X[index]=XL[index];
    152430                }
    153431        }
    154         SetControlInputsFromVectorx(femmodel,X);
    155         ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
     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*/
    156446        femmodel->OutputControlsx(&femmodel->results);
    157         femmodel->results->AddObject(new GenericExternalResult<double*>(femmodel->results->Size()+1,JEnum,mystruct.Jlist,(*mystruct.i),mystruct.N,0,0));
     447        femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,mystruct.Jlist,(*mystruct.i),mystruct.N,0,0));
    158448
    159449        /*Finalize*/
    … …  
    164454        solutioncore(femmodel);
    165455
     456
    166457        /*Clean-up and return*/
    167458        xDelete<double>(G);
    … …  
    170461        xDelete<double>(XU);
    171462        xDelete<double>(XL);
    172         xDelete<double>(scaling_factors);
    173         xDelete<double>(mystruct.Jlist);
     463        xDelete<IssmDouble>(scaling_factors);
     464        xDelete<IssmPDouble>(mystruct.Jlist);
    174465        xDelete<int>(mystruct.i);
    175 }
    176 
    177 /*Cost function definition*/
    178 void 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);
    279         xDelete<IssmDouble>(scaling_factors);
    280 }
     466}/*}}}*/
    281467
    282468#else
    283 void controlm1qn3_core(FemModel* femmodel){
    284         _error_("M1QN3 not installed");
    285 }
     469void controlm1qn3_core(FemModel* femmodel){_error_("M1QN3 not installed");}
    286470#endif //_HAVE_M1QN3_

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

    r19527 r22437  
    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

    r18521 r22437  
    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;
     100                //Jelem+=misfit*weight*Jdet*gauss->weight;
     101                Jelem+=misfit*Jdet*gauss->weight;
     102
    101103        }
    102104

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

    r18521 r22437  
    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 
     108                weight = 1.;
    109109                /*Add to cost function*/
    110110                Jelem+=misfit*weight*Jdet*gauss->weight;

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

    r18521 r22437  
    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 
     108                weight = 1.;
    109109                /*Add to cost function*/
    110110                Jelem+=misfit*weight*Jdet*gauss->weight;
Note: See TracChangeset for help on using the changeset viewer.