Changeset 23971


Ignore:
Timestamp:
05/31/19 13:02:22 (6 years ago)
Author:
Mathieu Morlighem
Message:

CHG: fixed a few things

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

Legend:

Unmodified
Added
Removed

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

    r23970 r23971  
    122122        iomodel->FetchDataToInput(elements,"md.hydrology.bump_height",HydrologyBumpHeightEnum);
    123123        iomodel->FetchDataToInput(elements,"md.hydrology.sheet_conductivity",HydrologySheetConductivityEnum);
    124         iomodel->FetchDataToInput(elements,"md.hydrology.neumannflux",HydrologyNeumannfluxEnum,0.);
     124        iomodel->FetchDataToInput(elements,"md.hydrology.neumannflux",HydrologyNeumannfluxEnum);
     125        iomodel->FetchDataToInput(elements,"md.hydrology.moulin_input",HydrologyMoulinInputEnum);
    125126        iomodel->FetchDataToInput(elements,"md.initialization.watercolumn",HydrologySheetThicknessEnum);
    126127        iomodel->FetchDataToInput(elements,"md.initialization.hydraulic_potential",HydraulicPotentialEnum);

  • issm/trunk-jpl/src/c/classes/Loads/Channel.cpp

    r23968 r23971  
    348348        /*Intermediaries */
    349349        IssmDouble  Jdet,v1,qc,fFactor,Afactor,Bfactor,Xifactor;
    350         IssmDouble  A,B,n,phi_old,phi,phi_0,dPw;
     350        IssmDouble  A,B,n,phi_old,phi,phi_0,dPw,ks,Ngrad;
    351351        IssmDouble  H,h,b,dphi[2],dphids,dphimds,db[2],dbds;
    352352        IssmDouble  xyz_list[NUMVERTICES][3];
    … …  
    369369        IssmDouble g         = element->FindParam(ConstantsGEnum);
    370370        IssmDouble kc        = element->FindParam(HydrologyChannelConductivityEnum);
    371         IssmDouble ks        = element->FindParam(HydrologySheetConductivityEnum);
    372371        IssmDouble lc        = element->FindParam(HydrologyChannelSheetWidthEnum);
    373372        IssmDouble c_t       = element->FindParam(HydrologyPressureMeltCoefficientEnum);
    … …  
    378377        Input* B_input      = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
    379378        Input* n_input      = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
     379        Input* ks_input     = element->GetInput(HydrologySheetConductivityEnum); _assert_(ks_input);
    380380        Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum);      _assert_(phiold_input);
    381381        Input* phi_input    = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
    … …  
    402402                phi_input->GetInputDerivativeValue(&dphi[0],&xyz_list_tria[0][0],gauss);
    403403                h_input->GetInputValue(&h,gauss);
     404                ks_input->GetInputValue(&ks,gauss);
    404405                B_input->GetInputValue(&B,gauss);
    405406                n_input->GetInputValue(&n,gauss);
    … …  
    412413                dphids  = dphi[0]*tx + dphi[1]*ty;
    413414                dphimds = rho_water*g*(db[0]*tx + db[1]*ty);
     415                Ngrad   = fabs(dphids);
     416                if(Ngrad<1.e-12) Ngrad = 1.e-12;
    414417
    415418                /*Compute the effective conductivity Kc = k h^alpha |grad Phi|^{beta-2} (same for sheet)*/
    416                 IssmDouble Kc = kc * pow(this->S,ALPHA_C) * pow(fabs(dphids),BETA_C-2.);
    417                 IssmDouble Ks = ks * pow(h      ,ALPHA_S) * pow(fabs(dphids),BETA_S-2.);
     419                IssmDouble Kc = kc * pow(this->S,ALPHA_C) * pow(Ngrad,BETA_C-2.);
     420                IssmDouble Ks = ks * pow(h      ,ALPHA_S) * pow(Ngrad,BETA_S-2.);
    418421
    419422                /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/
    … …  
    479482        IssmDouble  Jdet,v2,Afactor,Bfactor,fFactor;
    480483        IssmDouble  A,B,n,phi_old,phi,phi_0,dphimds;
    481         IssmDouble  H,h,b,db[2],dphids,qc,dPw;
     484        IssmDouble  H,h,b,db[2],dphids,qc,dPw,ks;
    482485        IssmDouble  xyz_list[NUMVERTICES][3];
    483486        IssmDouble  xyz_list_tria[3][3];
    … …  
    496499        IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
    497500        IssmDouble kc        = element->FindParam(HydrologyChannelConductivityEnum);
    498         IssmDouble ks        = element->FindParam(HydrologySheetConductivityEnum);
    499501        IssmDouble g         = element->FindParam(ConstantsGEnum);
    500502        IssmDouble lc        = element->FindParam(HydrologyChannelSheetWidthEnum);
    … …  
    506508        Input* B_input      = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
    507509        Input* n_input      = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
     510        Input* ks_input     = element->GetInput(HydrologySheetConductivityEnum); _assert_(ks_input);
    508511        Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum);      _assert_(phiold_input);
    509512        Input* phi_input    = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
    … …  
    526529                b_input->GetInputDerivativeValue(&db[0],&xyz_list_tria[0][0],gauss);
    527530                h_input->GetInputValue(&h,gauss);
     531                ks_input->GetInputValue(&ks,gauss);
    528532                B_input->GetInputValue(&B,gauss);
    529533                n_input->GetInputValue(&n,gauss);
    … …  
    585589
    586590        /*Intermediaries */
    587         IssmDouble  A,B,n,phi,phi_0;
     591        IssmDouble  A,B,n,phi,phi_0,ks,Ngrad;
    588592        IssmDouble  H,h,b,dphi[2],dphids,dphimds,db[2],dbds;
    589593        IssmDouble  xyz_list[NUMVERTICES][3];
    … …  
    599603        IssmDouble g         = element->FindParam(ConstantsGEnum);
    600604        IssmDouble kc        = element->FindParam(HydrologyChannelConductivityEnum);
    601         IssmDouble ks        = element->FindParam(HydrologySheetConductivityEnum);
    602605        IssmDouble lc        = element->FindParam(HydrologyChannelSheetWidthEnum);
    603606        IssmDouble c_t       = element->FindParam(HydrologyPressureMeltCoefficientEnum);
    … …  
    609612        Input* B_input      = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
    610613        Input* n_input      = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
     614        Input* ks_input     = element->GetInput(HydrologySheetConductivityEnum); _assert_(ks_input);
    611615        Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum);      _assert_(phiold_input);
    612616        Input* phi_input    = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
    … …  
    626630        phi_input->GetInputDerivativeValue(&dphi[0],&xyz_list_tria[0][0],gauss);
    627631        h_input->GetInputValue(&h,gauss);
     632        ks_input->GetInputValue(&ks,gauss);
    628633        B_input->GetInputValue(&B,gauss);
    629634        n_input->GetInputValue(&n,gauss);
    … …  
    636641        dphids  = dphi[0]*tx + dphi[1]*ty;
    637642        dphimds = rho_water*g*(db[0]*tx + db[1]*ty);
     643        Ngrad   = fabs(dphids);
     644        if(Ngrad<1.e-12) Ngrad = 1.e-12;
    638645
    639646        /*Compute the effective conductivity Kc = k h^alpha |grad Phi|^{beta-2} (same for sheet)*/
    640         IssmDouble Kc = kc * pow(this->S,ALPHA_C) * pow(fabs(dphids),BETA_C-2.);
    641         IssmDouble Ks = ks * pow(h      ,ALPHA_S) * pow(fabs(dphids),BETA_S-2.);
     647        IssmDouble Kc = kc * pow(this->S,ALPHA_C) * pow(Ngrad,BETA_C-2.);
     648        IssmDouble Ks = ks * pow(h      ,ALPHA_S) * pow(Ngrad,BETA_S-2.);
    642649
    643650        /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/

  • issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp

    r23960 r23971  
    159159
    160160        switch(analysis_type){
     161                case HydrologyGlaDSAnalysisEnum:
     162                        pe = this->CreatePVectorHydrologyGlaDS();
     163                        break;
    161164                case HydrologyShaktiAnalysisEnum:
    162165                        pe = this->CreatePVectorHydrologyShakti();
    … …  
    310313/*}}}*/
    311314
     315ElementVector* Moulin::CreatePVectorHydrologyGlaDS(void){/*{{{*/
     316
     317        /*If this node is not the master node (belongs to another partition of the
     318         * mesh), don't add the moulin input a second time*/
     319        if(node->IsClone()) return NULL;
     320
     321        IssmDouble moulin_load;
     322
     323        /*Initialize Element matrix*/
     324        ElementVector* pe=new ElementVector(&node,1,this->parameters);
     325
     326        this->element->GetInputValue(&moulin_load,node,HydrologyMoulinInputEnum);
     327        pe->values[0]=moulin_load;
     328
     329        /*Clean up and return*/
     330        return pe;
     331}
     332/*}}}*/
    312333ElementVector* Moulin::CreatePVectorHydrologyShakti(void){/*{{{*/
    313334

  • issm/trunk-jpl/src/c/classes/Loads/Moulin.h

    r23959 r23971  
    7777
    7878                ElementVector* CreatePVectorHydrologyShakti(void);
     79                ElementVector* CreatePVectorHydrologyGlaDS(void);
    7980                ElementVector* CreatePVectorHydrologyDCInefficient(void);
    8081                ElementVector* CreatePVectorHydrologyDCEfficient(void);
Note: See TracChangeset for help on using the changeset viewer.