Changeset 24950


Ignore:
Timestamp:
06/01/20 17:42:54 (5 years ago)
Author:
Eric.Larour
Message:

CHG: fixing some bugs introduced in the ESA core.

Location:
issm/trunk-jpl/src/c
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/EsaAnalysis.cpp

    r24335 r24950  
    5959        /*some constant parameters: */
    6060        parameters->AddObject(iomodel->CopyConstantObject("md.esa.hemisphere",EsaHemisphereEnum));
    61 
     61        parameters->AddObject(iomodel->CopyConstantObject("md.slr.eartharea",SealevelEarthAreaEnum));
     62        parameters->AddObject(iomodel->CopyConstantObject("md.slr.rigid",SealevelriseRigidEnum));
     63        parameters->AddObject(iomodel->CopyConstantObject("md.slr.horiz",SealevelriseHorizEnum));
     64        parameters->AddObject(iomodel->CopyConstantObject("md.slr.elastic",SealevelriseElasticEnum));
     65        parameters->AddObject(iomodel->CopyConstantObject("md.slr.rotation",SealevelriseRotationEnum));
     66       
    6267        /*love numbers: */
    6368        iomodel->FetchData(&love_h,&nl,NULL,"md.esa.love_h");

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r24949 r24950  
    53365336
    53375337        /*Get area of element: precomputed in the sealevelrise_core_geometry:*/
    5338         this->GetInput2Value(&area,AreaEnum);
     5338        area=GetAreaSpherical();
    53395339
    53405340        /*element centroid (spherical): */
    … …  
    60846084        S=0; for(int i=0;i<NUMVERTICES;i++) S+=Sg_old[this->vertices[i]->Sid()]/NUMVERTICES;
    60856085
    6086         /*convert to SI: */
     6086        /*convert to kg/m^2: */
    60876087        S=S*rho_water;
    60886088
    60896089        for(int i=0;i<gsize;i++) Sgo[i]+=G[i]*S;
     6090        //cblas_daxpy(gsize,S,&G[0],1,&Sgo[0],1); //for later.
    60906091
    60916092        return;

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

    r23587 r24950  
    7272        if(isesa)femmodel->SetCurrentConfiguration(EsaAnalysisEnum);
    7373
     74        if(VerboseSolution()) _printf0_("   computing elastic geodetic core\n");
    7475        if(isesa){
    7576

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

    r24949 r24950  
    232232}
    233233/*}}}*/
    234 void steric_core(FemModel* femmodel){ /*{{{*/
    235 
    236         /*variables:*/
    237         Vector<IssmDouble> *bedrock  = NULL;
    238         Vector<IssmDouble> *SL  = NULL;
    239         Vector<IssmDouble> *steric_rate_g  = NULL;
    240         Vector<IssmDouble> *dynamic_rate_g = NULL;
    241         Vector<IssmDouble> *hydro_rate_g  = NULL;
    242         Vector<IssmDouble> *U_esa_rate= NULL;
    243         Vector<IssmDouble> *N_esa_rate= NULL;
    244         Vector<IssmDouble> *U_gia_rate= NULL;
    245         Vector<IssmDouble> *N_gia_rate= NULL;
    246 
    247         /*parameters: */
    248         bool isslr=0;
    249         int  solution_type;
    250         IssmDouble          dt;
    251         int  geodetic=0;
    252 
    253         /*Retrieve parameters:*/
    254         femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
    255         femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    256         femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    257         femmodel->parameters->FindParam(&geodetic,SealevelriseGeodeticEnum);
    258 
    259         /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/
    260         if(solution_type==SealevelriseSolutionEnum)isslr=1;
    261 
    262         /*Should we be here?:*/
    263         if(!isslr)return;
    264 
    265         /*Verbose: */
    266         if(VerboseSolution()) _printf0_("         computing steric sea level rise\n");
    267 
    268         /*Retrieve geoid viscous and elastic rates, bedrock uplift viscous and elastic rates + steric rate, as vectors:*/
    269         GetVectorFromInputsx(&bedrock,femmodel,BedEnum,VertexSIdEnum);
    270         GetVectorFromInputsx(&SL,femmodel,SealevelEnum,VertexSIdEnum);
    271         GetStericRate(&steric_rate_g,femmodel);
    272         GetDynamicRate(&dynamic_rate_g,femmodel);
    273         GetVectorFromInputsx(&hydro_rate_g,femmodel,SealevelriseHydroRateEnum,VertexSIdEnum);
    274         if(geodetic){
    275                 GetVectorFromInputsx(&U_esa_rate,femmodel,SealevelUEsaRateEnum,VertexSIdEnum);
    276                 GetVectorFromInputsx(&U_gia_rate,femmodel,SealevelUGiaRateEnum,VertexSIdEnum);
    277                 GetVectorFromInputsx(&N_esa_rate,femmodel,SealevelNEsaRateEnum,VertexSIdEnum);
    278                 GetVectorFromInputsx(&N_gia_rate,femmodel,SealevelNGiaRateEnum,VertexSIdEnum);
    279         }
    280 
    281         /*compute: sea level change = initial sea level + (N_gia_rate+N_esa_rate)  * dt + steric_rate + dynamic_rate + hydro_rate* dt*/
    282         if(geodetic){
    283                 SL->AXPY(N_gia_rate,dt);
    284                 SL->AXPY(N_esa_rate,dt);
    285         }
    286         SL->AXPY(steric_rate_g,dt);
    287         SL->AXPY(dynamic_rate_g,dt);
    288         SL->AXPY(hydro_rate_g,dt);
    289 
    290         /*compute new bedrock position: */
    291         if(geodetic){
    292                 bedrock->AXPY(U_esa_rate,dt);
    293                 bedrock->AXPY(U_gia_rate,dt);
    294         }
    295 
    296         /*update element inputs:*/
    297         InputUpdateFromVectorx(femmodel,bedrock,BedEnum,VertexSIdEnum);
    298         InputUpdateFromVectorx(femmodel,SL,SealevelEnum,VertexSIdEnum);
    299 
    300         /*Free ressources:*/   
    301         delete bedrock;
    302         delete SL;
    303         delete steric_rate_g;
    304         delete dynamic_rate_g;
    305         delete hydro_rate_g;
    306         if(geodetic){
    307                 delete U_esa_rate;
    308                 delete U_gia_rate;
    309                 delete N_esa_rate;
    310                 delete N_gia_rate;
    311         }
    312 }
    313 /*}}}*/
    314234SealevelMasks* sealevelrise_core_masks(FemModel* femmodel) {  /*{{{*/
    315235
    … …  
    321241        /*go through elements and fill the masks: */
    322242        for (int i=0;i<femmodel->elements->Size();i++){
    323 
    324243                Element*   element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    325244                element->SetSealevelMasks(masks);
    … …  
    699618}
    700619/*}}}*/
     620void steric_core(FemModel* femmodel){ /*{{{*/
     621
     622        /*variables:*/
     623        Vector<IssmDouble> *bedrock  = NULL;
     624        Vector<IssmDouble> *SL  = NULL;
     625        Vector<IssmDouble> *steric_rate_g  = NULL;
     626        Vector<IssmDouble> *dynamic_rate_g = NULL;
     627        Vector<IssmDouble> *hydro_rate_g  = NULL;
     628        Vector<IssmDouble> *U_esa_rate= NULL;
     629        Vector<IssmDouble> *N_esa_rate= NULL;
     630        Vector<IssmDouble> *U_gia_rate= NULL;
     631        Vector<IssmDouble> *N_gia_rate= NULL;
     632
     633        /*parameters: */
     634        bool isslr=0;
     635        int  solution_type;
     636        IssmDouble          dt;
     637        int  geodetic=0;
     638
     639        /*Retrieve parameters:*/
     640        femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
     641        femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     642        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
     643        femmodel->parameters->FindParam(&geodetic,SealevelriseGeodeticEnum);
     644
     645        /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/
     646        if(solution_type==SealevelriseSolutionEnum)isslr=1;
     647
     648        /*Should we be here?:*/
     649        if(!isslr)return;
     650
     651        /*Verbose: */
     652        if(VerboseSolution()) _printf0_("         computing steric sea level rise\n");
     653
     654        /*Retrieve geoid viscous and elastic rates, bedrock uplift viscous and elastic rates + steric rate, as vectors:*/
     655        GetVectorFromInputsx(&bedrock,femmodel,BedEnum,VertexSIdEnum);
     656        GetVectorFromInputsx(&SL,femmodel,SealevelEnum,VertexSIdEnum);
     657        GetStericRate(&steric_rate_g,femmodel);
     658        GetDynamicRate(&dynamic_rate_g,femmodel);
     659        GetVectorFromInputsx(&hydro_rate_g,femmodel,SealevelriseHydroRateEnum,VertexSIdEnum);
     660        if(geodetic){
     661                GetVectorFromInputsx(&U_esa_rate,femmodel,SealevelUEsaRateEnum,VertexSIdEnum);
     662                GetVectorFromInputsx(&U_gia_rate,femmodel,SealevelUGiaRateEnum,VertexSIdEnum);
     663                GetVectorFromInputsx(&N_esa_rate,femmodel,SealevelNEsaRateEnum,VertexSIdEnum);
     664                GetVectorFromInputsx(&N_gia_rate,femmodel,SealevelNGiaRateEnum,VertexSIdEnum);
     665        }
     666
     667        /*compute: sea level change = initial sea level + (N_gia_rate+N_esa_rate)  * dt + steric_rate + dynamic_rate + hydro_rate* dt*/
     668        if(geodetic){
     669                SL->AXPY(N_gia_rate,dt);
     670                SL->AXPY(N_esa_rate,dt);
     671        }
     672        SL->AXPY(steric_rate_g,dt);
     673        SL->AXPY(dynamic_rate_g,dt);
     674        SL->AXPY(hydro_rate_g,dt);
     675
     676        /*compute new bedrock position: */
     677        if(geodetic){
     678                bedrock->AXPY(U_esa_rate,dt);
     679                bedrock->AXPY(U_gia_rate,dt);
     680        }
     681
     682        /*update element inputs:*/
     683        InputUpdateFromVectorx(femmodel,bedrock,BedEnum,VertexSIdEnum);
     684        InputUpdateFromVectorx(femmodel,SL,SealevelEnum,VertexSIdEnum);
     685
     686        /*Free ressources:*/   
     687        delete bedrock;
     688        delete SL;
     689        delete steric_rate_g;
     690        delete dynamic_rate_g;
     691        delete hydro_rate_g;
     692        if(geodetic){
     693                delete U_esa_rate;
     694                delete U_gia_rate;
     695                delete N_esa_rate;
     696                delete N_gia_rate;
     697        }
     698}
     699/*}}}*/
    701700
    702701/*support routines:*/
Note: See TracChangeset for help on using the changeset viewer.