Changeset 25330


Ignore:
Timestamp:
07/31/20 21:56:02 (5 years ago)
Author:
Mathieu Morlighem
Message:

CHG: allow AD to deal with P0 inputs

Location:
issm/trunk-jpl/src/c/cores
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r24099 r25330  
    197197
    198198        /*Recover some parameters*/
    199         IssmDouble *scaling_factors = NULL;
    200         int        *N               = NULL;
    201         int        *control_enum    = NULL;
     199        double *scaling_factors = NULL;
     200        int    *M = NULL;
     201        int    *N = NULL;
     202        int    *control_enum    = NULL;
    202203        femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    203204        femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    204         femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
     205        femmodel->parameters->FindParamAndMakePassive(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
    205206        femmodel->parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
     207        femmodel->parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
    206208        femmodel->parameters->FindParam(&control_enum,NULL,InversionControlParametersEnum);
    207         int numberofvertices=femmodel->vertices->NumberOfVertices();
    208209
    209210        /*Constrain input vector and update controls*/
    … …  
    213214        GetPassiveVectorFromControlInputsx(&XU,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
    214215
    215         int N_add = 0;
     216        int offset = 0;
    216217        for (int c=0;c<num_controls;c++){
    217                 for(int i=0;i<numberofvertices*N[c];i++){
    218                         int index = N_add*numberofvertices+i;
    219                         X[index] = X[index]*reCast<double>(scaling_factors[c]);
     218                for(int i=0;i<M[c]*N[c];i++){
     219                        int index = offset+i;
     220                        X[index] = X[index]*scaling_factors[c];
    220221                        if(X[index]>XU[index]) X[index]=XU[index];
    221222                        if(X[index]<XL[index]) X[index]=XL[index];
    222223                }
    223                 N_add+=N[c];
     224                offset += M[c]*N[c];
    224225        }
    225226
    … …  
    439440        /*Constrain Gradient*/
    440441        IssmDouble  Gnorm = 0.;
    441         N_add = 0;
    442         for (int c=0;c<num_controls;c++){
    443                 for(int i=0;i<numberofvertices*N[c];i++){
    444                         int index = N_add*numberofvertices+i;
     442        offset = 0;
     443        for(int c=0;c<num_controls;c++){
     444                for(int i=0;i<M[c]*N[c];i++){
     445                        int index = offset+i;
    445446                        if(X[index]>=XU[index]) G[index]=0.;
    446447                        if(X[index]<=XL[index]) G[index]=0.;
    447                         G[index] = G[index]*reCast<double>(scaling_factors[c]);
    448                         X[index] = X[index]/reCast<double>(scaling_factors[c]);
     448                        G[index] = G[index]*scaling_factors[c];
     449                        X[index] = X[index]/scaling_factors[c];
    449450                        Gnorm += G[index]*G[index];
    450451                }
    451                 N_add+=N[c];
     452                offset += M[c]*N[c];
    452453        }
    453454        Gnorm = sqrt(Gnorm);
    … …  
    464465        xDelete<int>(control_enum);
    465466        xDelete<int>(N);
    466         xDelete<IssmDouble>(scaling_factors);
     467        xDelete<double>(scaling_factors);
    467468        xDelete<IssmPDouble>(totalgradient);
    468469}/*}}}*/
    … …  
    474475        double           dxmin,gttol;
    475476        int          maxsteps,maxiter;
    476         int          intn,numberofvertices,num_controls,num_cost_functions,solution_type;
    477         IssmDouble      *scaling_factors = NULL;
     477        int          intn,num_controls,num_cost_functions,solution_type;
     478        double*      scaling_factors = NULL;
    478479        double      *X  = NULL;
    479480        double      *G  = NULL;
    480481        int*                    N       = NULL;
    481         int offset = 0;
    482         int N_add;
     482   int*                 M       = NULL;
    483483        int* control_enum;
    484484
    … …  
    491491        femmodel->parameters->FindParamAndMakePassive(&dxmin,InversionDxminEnum);
    492492        femmodel->parameters->FindParamAndMakePassive(&gttol,InversionGttolEnum);
    493         femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
     493        femmodel->parameters->FindParamAndMakePassive(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
    494494        femmodel->parameters->FindParam(&control_enum,NULL,InversionControlParametersEnum);
    495495        femmodel->parameters->SetParam(false,SaveResultsEnum);
    496496        femmodel->parameters->FindParam(&N,NULL,ControlInputSizeNEnum);
    497         numberofvertices=femmodel->vertices->NumberOfVertices();
     497   femmodel->parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
    498498
    499499        /*Initialize M1QN3 parameters*/
    … …  
    517517        /*Get initial guess*/
    518518        GetPassiveVectorFromControlInputsx(&X,&intn,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
    519         //_assert_(intn==numberofvertices*num_controls);
    520519
    521520        /*Get problem dimension and initialize gradient and initial guess*/
    … …  
    524523
    525524        /*Scale control for M1QN3*/
    526         N_add = 0;
    527         for (int c=0;c<num_controls;c++){
    528                 for(int i=0;i<numberofvertices*N[c];i++){
    529                         int index = N_add*numberofvertices+i;
    530                         X[index] = X[index]/reCast<IssmPDouble>(scaling_factors[c]);
    531                 }
    532                 N_add+=N[c];
    533         }
     525   int offset = 0;
     526   for(int c=0;c<num_controls;c++){
     527      for(int i=0;i<M[c]*N[c];i++){
     528         int index = offset+i;
     529         X[index] = X[index]/scaling_factors[c];
     530      }
     531      offset += M[c]*N[c];
     532   }
    534533
    535534        /*Allocate m1qn3 working arrays (see documentation)*/
    … …  
    578577        GetPassiveVectorFromControlInputsx(&XU,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
    579578
    580         N_add = 0;
    581         for (int c=0;c<num_controls;c++){
    582                 for(int i=0;i<numberofvertices*N[c];i++){
    583                         int index = N_add*numberofvertices+i;
    584                         X[index] = X[index]*reCast<double>(scaling_factors[c]);
    585                         if(X[index]>XU[index]) X[index]=XU[index];
    586                         if(X[index]<XL[index]) X[index]=XL[index];
    587                 }
    588                 N_add +=N[c];
    589         }
     579   offset = 0;
     580   for (int c=0;c<num_controls;c++){
     581      for(int i=0;i<M[c]*N[c];i++){
     582         int index = offset+i;
     583         X[index] = X[index]*scaling_factors[c];
     584         if(X[index]>XU[index]) X[index]=XU[index];
     585         if(X[index]<XL[index]) X[index]=XL[index];
     586      }
     587      offset += M[c]*N[c];
     588   }
    590589
    591590        /*Set X as our new control*/
    … …  
    612611
    613612                        /*Disect results*/
    614                         GenericExternalResult<IssmPDouble*>* G_output = new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,Gradient1Enum+i,&G[offset],N[i],numberofvertices,1,0.);
    615                         GenericExternalResult<IssmPDouble*>* X_output = new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,control_enum[i],&X[offset],N[i],numberofvertices,1,0.);
     613                        GenericExternalResult<IssmPDouble*>* G_output = new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,Gradient1Enum+i,&G[offset],N[i],M[i],1,0.);
     614                        GenericExternalResult<IssmPDouble*>* X_output = new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,control_enum[i],&X[offset],N[i],M[i],1,0.);
    616615
    617616                        /*transpose for consistency with MATLAB's formating*/
    … …  
    623622                        femmodel->results->AddObject(X_output);
    624623
    625                         offset += N[i]*numberofvertices;
     624                        offset += N[i]*M[i];
    626625                }
    627626        }
    … …  
    649648        xDelete<double>(XU);
    650649        xDelete<double>(XL);
    651         xDelete<IssmDouble>(scaling_factors);
     650        xDelete<double>(scaling_factors);
    652651        xDelete<IssmPDouble>(mystruct.Jlist);
    653652        xDelete<int>(mystruct.i);

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

    r25317 r25330  
    4242        double  f,dxmin,gttol;
    4343        int     maxsteps,maxiter;
    44         int     intn,numberofvertices,num_controls,num_cost_functions,solution_type;
     44        int     intn,num_controls,num_cost_functions,solution_type;
    4545        double *scaling_factors = NULL;
    4646        double *X  = NULL;
    … …  
    6363        femmodel->parameters->FindParamAndMakePassive(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
    6464        femmodel->parameters->SetParam(false,SaveResultsEnum);
    65         numberofvertices=femmodel->vertices->NumberOfVertices();
    6665
    6766        /*Initialize M1QN3 parameters*/
    … …  
    211210
    212211        /*Recover some parameters*/
    213         int num_responses,num_controls,numberofvertices,solution_type;
     212        int num_responses,num_controls,solution_type;
    214213        double* scaling_factors = NULL;
    215214        int* M = NULL;
    … …  
    221220        femmodel->parameters->FindParamAndMakePassive(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
    222221        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    223         numberofvertices=femmodel->vertices->NumberOfVertices();
    224222
    225223        /*Constrain input vector and update controls*/
Note: See TracChangeset for help on using the changeset viewer.