Changeset 25247


Ignore:
Timestamp:
07/10/20 04:52:59 (5 years ago)
Author:
bdef
Message:

CHG:clean up and parameter fetch optimisation in Hydrosrc/c/analyses/HydrologyDCEfficientAnalysis.cpp

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

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.h

    r25227 r25247  
    3535
    3636                /*Intermediaries*/
    37                 IssmDouble EplStoring(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input);
    38                 IssmDouble EplTransmitivity(Element* element,Gauss* gauss, Input2* epl_thick_input, Input2* epl_head_input, Input2* base_input);
     37                IssmDouble EplStoring(Element* element,Gauss* gauss, Input2* epl_thick_input);
     38                IssmDouble EplTransmitivity(Element* element,Gauss* gauss, Input2* epl_thick_input);
    3939                void GetHydrologyDCInefficientHmax(IssmDouble* ph_max,Element* element, Node* innode);
    4040                IssmDouble GetHydrologyKMatrixTransfer(Element* element);

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

    r25227 r25247  
    238238        Input2* sed_head_input = basalelement->GetInput2(SedimentHeadSubstepEnum);
    239239        Input2* base_input     = basalelement->GetInput2(BaseEnum);
    240         Input2* old_wh_input   = basalelement->GetInput2(SedimentHeadOldEnum);                  _assert_(old_wh_input);
    241240
    242241        /*Transfer related Inputs*/
    … …  
    244243                basalelement->GetInput2Value(&active_element,HydrologydcMaskEplactiveEltEnum);
    245244        }
    246 
    247245        /* Start  looping on the number of gaussian points: */
    248246        Gauss* gauss=basalelement->NewGauss(2);
    249247
    250248        for(int ig=gauss -> begin();ig<gauss->end();ig++){
    251                 gauss          -> GaussPoint(ig);
    252                 basalelement   -> JacobianDeterminant(&Jdet,xyz_list,gauss);
     249                gauss       ->GaussPoint(ig);
     250                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
    253251                basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    254252                basalelement->NodalFunctions(basis,gauss);
    … …  
    259257                /*Diffusivity*/
    260258                D_scalar=sediment_transmitivity*gauss->weight*Jdet;
    261                 //D_scalar=gauss->weight*Jdet;
    262259                if(dt!=0.) D_scalar=D_scalar*dt;
    263260                for(int i=0;i<numnodes;i++){
    … …  
    266263                        }
    267264                }
    268 
    269265                /*Transient*/
    270266                if(dt!=0.){
    271267                        D_scalar=sediment_storing*gauss->weight*Jdet;
    272                         //D_scalar=(sediment_storing/sediment_transmitivity)*gauss->weight*Jdet;
    273268                        for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
    274 
    275269                        /*Transfer EPL part*/
    276270                        if(isefficientlayer){
    … …  
    278272                                        transfer=GetHydrologyKMatrixTransfer(basalelement);
    279273                                        D_scalar=dt*transfer*gauss->weight*Jdet;
    280                                         //D_scalar=dt*(transfer/sediment_transmitivity)*gauss->weight*Jdet;
    281274                                        for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
    282275                                }
    … …  
    298291
    299292        /*Intermediaries*/
    300         bool             thawed_element;
    301         int                      domaintype;
     293        bool            thawed_element;
     294        int             domaintype;
    302295        Element* basalelement;
    303296
    … …  
    328321        /*Intermediaries */
    329322        bool       active_element,isefficientlayer;
    330         int        smb_model;
    331         int        smbsubstepping;
    332         int        hydrologysubstepping;
    333         int        smb_averaging;
     323        int        smb_model,smbsubstepping;
     324        int        hydrologysubstepping,smb_averaging;
    334325        IssmDouble dt,scalar,sediment_storing;
    335326        IssmDouble water_head,sediment_transmitivity;
    336327        IssmDouble water_load,runoff_value,transfer;
    337         IssmDouble Jdet;
     328        IssmDouble Jdet,time;
    338329
    339330        IssmDouble *xyz_list             = NULL;
    … …  
    355346        basalelement->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
    356347        basalelement->FindParam(&smb_model,SmbEnum);
    357         basalelement->FindParam(&smb_averaging,SmbAveragingEnum);
    358348
    359349        Input2* sed_head_input   = basalelement->GetInput2(SedimentHeadSubstepEnum);
    … …  
    363353        Input2* SedTrans_input   = basalelement->GetInput2(HydrologydcSedimentTransmitivityEnum); _assert_(SedTrans_input);
    364354
    365         IssmDouble time;
    366         basalelement->FindParam(&time,TimeEnum);
    367 
    368355        if(dt!= 0.){
    369356                old_wh_input = basalelement->GetInput2(SedimentHeadOldEnum); _assert_(old_wh_input);
    370357        }
    371358        if(smb_model==SMBgradientscomponentsEnum){
     359                basalelement->FindParam(&time,TimeEnum);
    372360                basalelement->FindParam(&smbsubstepping,SmbStepsPerStepEnum);
    373361                basalelement->FindParam(&hydrologysubstepping,HydrologyStepsPerStepEnum);
    … …  
    383371                else{
    384372                        //finer stepping in smb, we average the runoff from transient input
     373                        basalelement->FindParam(&smb_averaging,SmbAveragingEnum);
    385374                        dummy_input = basalelement->GetInput2(SmbRunoffTransientEnum,time-dt,time,smb_averaging); _assert_(dummy_input);
    386375                }
    … …  
    395384        /* Start  looping on the number of gaussian points: */
    396385        Gauss* gauss=basalelement->NewGauss(2);
    397 
    398         IssmDouble yts;
    399         basalelement->FindParam(&yts,ConstantsYtsEnum);
    400386
    401387        for(int ig=gauss->begin();ig<gauss->end();ig++){
    … …  
    411397                        else                     runoff_value = 0.;
    412398                        scalar = Jdet*gauss->weight*(water_load+runoff_value);
    413                         //scalar = Jdet*gauss->weight*(water_load)/sediment_transmitivity;
    414399                        if(dt!=0.) scalar = scalar*dt;
    415400                        for(int i=0;i<numnodes;i++){
    … …  
    424409                                else runoff_value = 0.;
    425410                                scalar = Jdet*gauss->weight*(water_load+runoff_value);
    426                                 //scalar = Jdet*gauss->weight*(water_load)/sediment_transmitivity;
    427411                                if(dt!=0.) scalar = scalar*dt;
    428412                                for(int i=0;i<numnodes;i++){
    … …  
    445429                                }
    446430                                scalar = Jdet*gauss->weight*((water_head*sediment_storing)+(dt*transfer));
    447                                 //scalar = Jdet*gauss->weight*((water_head*sediment_storing)+(dt*transfer))/sediment_transmitivity;
    448431                                for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
    449432                        }
    450433                        else{
    451434                                scalar = Jdet*gauss->weight*(water_head*sediment_storing);
    452                                 //scalar = Jdet*gauss->weight*(water_head*sediment_storing)/sediment_transmitivity;
    453435                                for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
    454436                        }
    … …  
    571553        IssmDouble storing,yield;
    572554        IssmDouble base_elev,prestep_head,water_sheet;
    573         IssmDouble rho_freshwater           = element->FindParam(MaterialsRhoFreshwaterEnum);
    574         IssmDouble g                        = element->FindParam(ConstantsGEnum);
    575         IssmDouble sediment_porosity        = element->FindParam(HydrologydcSedimentPorosityEnum);
     555        IssmDouble porewater_mass           = element->FindParam(HydrologydcSedimentPoreWaterMassEnum);
     556        IssmDouble layer_compressibility    = element->FindParam(HydrologydcSedimentLayerCompressibilityEnum);
    576557        IssmDouble sediment_thickness       = element->FindParam(HydrologydcSedimentThicknessEnum);
    577         IssmDouble sediment_compressibility = element->FindParam(HydrologydcSedimentCompressibilityEnum);
    578         IssmDouble water_compressibility    = element->FindParam(HydrologydcWaterCompressibilityEnum);
    579558        element->FindParam(&unconf_scheme,HydrologydcUnconfinedFlagEnum);
    580559        switch(unconf_scheme){
    581560        case 0:
    582                 sediment_storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
     561                sediment_storing=porewater_mass*sediment_thickness*layer_compressibility;
    583562                break;
    584563        case 1:
     564                yield = element->FindParam(HydrologydcSedimentPorosityEnum);
    585565                base_input->GetInputValue(&base_elev,gauss);
    586566                sed_head_input->GetInputValue(&prestep_head,gauss);
     567
    587568                water_sheet=max(0.0,(prestep_head-(base_elev-sediment_thickness)));
    588 
    589                 /* if (water_sheet<sediment_thickness){ */
    590                 /*      sediment_storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity)); */
    591                 /* } */
    592                 /* else{ */
    593                 /*      sediment_storing=sediment_porosity; */
    594                 /* } */
    595                 storing=rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));
     569                storing=porewater_mass*sediment_thickness*layer_compressibility;
    596570                //using logistic function for heavyside approximation
    597571                expfac=10.;
    598                 yield=sediment_porosity;
    599572                sediment_storing=yield+(storing-yield)/(1+exp(-2*expfac*(water_sheet-0.99*sediment_thickness)));
    600573                break;
    … …  
    606579IssmDouble HydrologyDCInefficientAnalysis::SedimentTransmitivity(Element* element,Gauss* gauss,Input2* sed_head_input, Input2* base_input,Input2* SedTrans_input){/*{{{*/
    607580        int unconf_scheme;
    608         IssmDouble ratio,expfac;
    609581        IssmDouble sediment_transmitivity;
    610582        IssmDouble FullLayer_transmitivity;
    611         IssmDouble meltingrate;
    612         IssmDouble groundedice;
    613583        IssmDouble base_elev,prestep_head,water_sheet;
    614584        IssmDouble sediment_thickness       = element->FindParam(HydrologydcSedimentThicknessEnum);

  • issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp

    r25030 r25247  
    420420                /*If already processed return*/
    421421                if(fabs(this->current_step-this_step)<1.e-5) return;
    422                 //              if(this->current_step>this_step-1.e-5 && this->current_step<this_step+1.e-5) return;
    423422
    424423                /*Prepare input*/

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r25150 r25247  
    11/*!\file: CreateParameters.cpp
    22 * \brief general driver for creating parameters dataset
    3  */ 
     3 */
    44
    55#ifdef HAVE_CONFIG_H
    … …  
    6363        parameters->AddObject(iomodel->CopyConstantObject("md.calving.law",CalvingLawEnum));
    6464        parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.parameterization",FrontalForcingsParamEnum));
    65         parameters->AddObject(new IntParam(SealevelriseRunCountEnum,1)); 
     65        parameters->AddObject(new IntParam(SealevelriseRunCountEnum,1));
    6666
    6767          {/*This is specific to ice...*/
    … …  
    9898          }
    9999
    100         /*amr properties*/     
     100        /*amr properties*/
    101101        int amrtype,amr_frequency;
    102102        iomodel->FindConstant(&amr_frequency,"md.transient.amr_frequency");
    … …  
    156156                case LinearFloatingMeltRateEnum:
    157157                        iomodel->FindConstant(&interp,"md.timestepping.interp_forcings");
    158                         iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.deepwater_melting_rate"); 
     158                        iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.deepwater_melting_rate");
    159159                        if(N==1){
    160160                                _assert_(M==1);
    … …  
    166166                        }
    167167                        xDelete<IssmDouble>(transparam);
    168                         iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.upperwater_melting_rate"); 
     168                        iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.upperwater_melting_rate");
    169169                        if(N==1){
    170170                                _assert_(M==1);
    … …  
    176176                        }
    177177                        xDelete<IssmDouble>(transparam);
    178                         iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.deepwater_elevation"); 
     178                        iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.deepwater_elevation");
    179179                        if(N==1){
    180180                                _assert_(M==1);
    … …  
    183183                        else{
    184184                                _assert_(N==2);
    185                                 parameters->AddObject(new TransientParam(BasalforcingsDeepwaterElevationEnum,&transparam[0],&transparam[M],interp,M)); 
    186                         }
    187                         xDelete<IssmDouble>(transparam);
    188                         iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.upperwater_elevation"); 
     185                                parameters->AddObject(new TransientParam(BasalforcingsDeepwaterElevationEnum,&transparam[0],&transparam[M],interp,M));
     186                        }
     187                        xDelete<IssmDouble>(transparam);
     188                        iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.upperwater_elevation");
    189189                        if(N==1){
    190190                                _assert_(M==1);
    … …  
    226226                        parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.isplume",BasalforcingsPicoIsplumeEnum));
    227227                        iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.farocean_temperature");
    228                         _assert_(M>=1 && N>=1); 
     228                        _assert_(M>=1 && N>=1);
    229229                        parameters->AddObject(new TransientArrayParam(BasalforcingsPicoFarOceantemperatureEnum,transparam,&transparam[N*(M-1)],interp,N,M));
    230230                        xDelete<IssmDouble>(transparam);
    231231                        iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.farocean_salinity");
    232                         _assert_(M>=1 && N>=1); 
     232                        _assert_(M>=1 && N>=1);
    233233                        parameters->AddObject(new TransientArrayParam(BasalforcingsPicoFarOceansalinityEnum,transparam,&transparam[N*(M-1)],interp,N,M));
    234234                        xDelete<IssmDouble>(transparam);
    … …  
    237237                        parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.num_basins",BasalforcingsIsmip6NumBasinsEnum));
    238238                        parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.gamma_0",BasalforcingsIsmip6Gamma0Enum));
    239                         parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.islocal",BasalforcingsIsmip6IsLocalEnum)); 
     239                        parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.islocal",BasalforcingsIsmip6IsLocalEnum));
    240240                        iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.delta_t");
    241241                        parameters->AddObject(new DoubleVecParam(BasalforcingsIsmip6DeltaTEnum,transparam,N));
    … …  
    279279        }
    280280        iomodel->FindConstant(&time,"md.timestepping.start_time");
    281         parameters->AddObject(new DoubleParam(TimeEnum,time)); 
    282         parameters->AddObject(new IntParam(StepEnum,0)); 
     281        parameters->AddObject(new DoubleParam(TimeEnum,time));
     282        parameters->AddObject(new IntParam(StepEnum,0));
    283283
    284284        /*By default, save all results*/
    … …  
    421421        iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
    422422        if(hydrology_model==HydrologydcEnum){
     423                IssmDouble sedcomp, sedporo, watcomp, rhofresh, g;
     424
    423425                /*FIXME: this cshould go to Analysis!!!*/
    424                 parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.sediment_compressibility",HydrologydcSedimentCompressibilityEnum));
     426                /* parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.sediment_compressibility",HydrologydcSedimentCompressibilityEnum)); */
     427                /* parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.water_compressibility",HydrologydcWaterCompressibilityEnum)); */
    425428                parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.sediment_porosity",HydrologydcSedimentPorosityEnum));
    426429                parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.sediment_thickness",HydrologydcSedimentThicknessEnum));
    427                 parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.water_compressibility",HydrologydcWaterCompressibilityEnum));
    428430                parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.isefficientlayer",HydrologydcIsefficientlayerEnum));
     431
     432                iomodel->FindConstant(&sedcomp,"md.hydrology.sediment_compressibility");
     433                iomodel->FindConstant(&sedporo,"md.hydrology.sediment_porosity");
     434                iomodel->FindConstant(&watcomp,"md.hydrology.water_compressibility");
     435                iomodel->FindConstant(&rhofresh,"md.materials.rho_freshwater");
     436                iomodel->FindConstant(&g,"md.constants.g");
     437
     438                parameters->AddObject(new DoubleParam(HydrologydcSedimentLayerCompressibilityEnum,(watcomp + sedcomp/sedporo)));
     439                parameters->AddObject(new DoubleParam(HydrologydcSedimentPoreWaterMassEnum,(rhofresh*g*sedporo)));
     440
    429441
    430442                bool isefficientlayer;
    431443                iomodel->FindConstant(&isefficientlayer,"md.hydrology.isefficientlayer");
    432444                if(isefficientlayer){
    433                         parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_compressibility",HydrologydcEplCompressibilityEnum));
    434                         parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_porosity",HydrologydcEplPorosityEnum));
     445                        IssmDouble eplcomp, eplporo;
     446                        /* parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_compressibility",HydrologydcEplCompressibilityEnum)); */
     447                        /* parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_porosity",HydrologydcEplPorosityEnum)); */
    435448                        parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_initial_thickness",HydrologydcEplInitialThicknessEnum));
    436449                        parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_colapse_thickness",HydrologydcEplColapseThicknessEnum));
    437450                        parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_max_thickness",HydrologydcEplMaxThicknessEnum));
    438451                        parameters->AddObject(iomodel->CopyConstantObject("md.hydrology.epl_conductivity",HydrologydcEplConductivityEnum));
     452
     453                        iomodel->FindConstant(&eplcomp,"md.hydrology.epl_compressibility");
     454                        iomodel->FindConstant(&eplporo,"md.hydrology.epl_porosity");
     455                        parameters->AddObject(new DoubleParam(HydrologydcEplLayerCompressibilityEnum,(watcomp + eplcomp/eplporo)));
     456                        parameters->AddObject(new DoubleParam(HydrologydcEplPoreWaterMassEnum,(rhofresh*g*eplporo)));
    439457                }
    440458        }
    … …  
    472490        if(mass_flux_present){
    473491
    474                 /*Fetch the mass flux segments necessary to compute the mass fluxes.  Build a DoubleMatArrayParam object out of them: */ 
     492                /*Fetch the mass flux segments necessary to compute the mass fluxes.  Build a DoubleMatArrayParam object out of them: */
    475493                iomodel->FetchData(&array,&mdims_array,&ndims_array,&mass_flux_num_profiles,"md.qmu.mass_flux_segments");
    476494                if(mass_flux_num_profiles==0)_error_("mass_flux_num_profiles is 0, when MassFlux computations were requested!");
    … …  
    518536                        xDelete<IssmDouble>(matrix);
    519537                }
    520                 xDelete<int>(mdims_array); 
     538                xDelete<int>(mdims_array);
    521539                xDelete<int>(ndims_array);
    522540                xDelete<IssmDouble*>(array);

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r25139 r25247  
    176176        HydrologyStorageEnum,
    177177        HydrologydcEplColapseThicknessEnum,
    178         HydrologydcEplCompressibilityEnum,
    179178        HydrologydcEplConductivityEnum,
    180179        HydrologydcEplInitialThicknessEnum,
     180        HydrologydcEplLayerCompressibilityEnum,
    181181        HydrologydcEplMaxThicknessEnum,
    182         HydrologydcEplPorosityEnum,
     182        HydrologydcEplPoreWaterMassEnum,
    183183        HydrologydcEplThickCompEnum,
    184184        HydrologydcEplflipLockEnum,
    … …  
    189189        HydrologydcPenaltyLockEnum,
    190190        HydrologydcRelTolEnum,
    191         HydrologydcSedimentCompressibilityEnum,
    192191        HydrologydcSedimentlimitEnum,
    193192        HydrologydcSedimentlimitFlagEnum,
     193        HydrologydcSedimentLayerCompressibilityEnum,
     194        HydrologydcSedimentPoreWaterMassEnum,
    194195        HydrologydcSedimentPorosityEnum,
    195196        HydrologydcSedimentThicknessEnum,
    196197        HydrologydcTransferFlagEnum,
    197198        HydrologydcUnconfinedFlagEnum,
    198         HydrologydcWaterCompressibilityEnum,
    199199        HydrologyshreveStabilizationEnum,
    200200        IcecapToEarthCommEnum,
Note: See TracChangeset for help on using the changeset viewer.