Changeset 14695


Ignore:
Timestamp:
04/22/13 13:02:02 (12 years ago)
Author:
bdef
Message:

CHG: computing the real hydrologicalDC matrix and adding InputUpdateFormSolution

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

Legend:

Unmodified
Added
Removed

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

    r14679 r14695  
    58705870        /* Intermediaries */
    58715871        IssmDouble  D_scalar,Jdet;
    5872         IssmDouble  sediment_compressibility, sediment_porosity, sediment_thickness,
    5873                 sediment_transmitivity, water_compressibility, rho_freshwater, gravity, dt;
     5872        IssmDouble      sediment_transmitivity,dt;
    58745873        IssmDouble  sediment_storing;
    58755874        IssmDouble  xyz_list[NUMVERTICES][3];
    … …  
    58845883        /*Retrieve all inputs and parameters*/
    58855884        GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES);
    5886         sediment_compressibility                        = matpar->GetSedimentCompressibility();
    5887         sediment_porosity                                                       = matpar->GetSedimentPorosity();
    5888         sediment_thickness                                              = matpar->GetSedimentThickness();
    5889         sediment_transmitivity                          = matpar->GetSedimentTransmitivity();
    5890         water_compressibility                                   = matpar->GetWaterCompressibility();
    5891         rho_freshwater                                                          = matpar->GetRhoFreshwater();
    5892         gravity                                                                                         = matpar->GetG();
     5885
     5886        sediment_storing                         = matpar->GetSedimentStoring();
     5887        sediment_transmitivity = matpar->GetSedimentTransmitivity();
    58935888       
    58945889        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    5895 
    5896         /*Compute Storing coefficient from the above */
    5897 
    5898         sediment_storing = rho_freshwater*gravity*sediment_porosity*sediment_thickness*
    5899                 (water_compressibility+(sediment_compressibility/sediment_porosity));
    59005890
    59015891        /* Start looping on the number of gaussian points: */
    … …  
    60106000        IssmDouble Jdet;
    60116001        IssmDouble xyz_list[NUMVERTICES][3];
    6012         IssmDouble dt,gravity,scalar,water_head;
     6002        IssmDouble dt,scalar,water_head;
     6003        IssmDouble water_load;
     6004        IssmDouble sediment_storing;
    60136005        IssmDouble basis[numdof];
    60146006        GaussTria* gauss=NULL;
    … …  
    60186010
    60196011        /*Retrieve all inputs and parameters*/
     6012
     6013        sediment_storing = matpar->GetSedimentStoring();
     6014
    60206015        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    6021         gravity = matpar->GetG();
     6016
     6017        Input* water_input=inputs->GetInput(BasalforcingsMeltingRateEnum);  _assert_(water_input);
    60226018        Input* old_wh_input=NULL;
     6019       
    60236020        if(reCast<bool,IssmDouble>(dt)){
    60246021                old_wh_input=inputs->GetInput(SedimentHeadEnum); _assert_(old_wh_input);
    60256022        }
    6026 
     6023       
    60276024        /* Start  looping on the number of gaussian points: */
    60286025        gauss=new GaussTria(2);
    … …  
    60346031                GetNodalFunctions(basis, gauss);
    60356032
    6036                 /*Gravity term*/
    6037                 scalar = Jdet*gauss->weight*gravity;
     6033                /*Loading term*/
     6034                water_input->GetInputValue(&water_load,gauss);
     6035                scalar = Jdet*gauss->weight*water_load;
    60386036                if(reCast<bool,IssmDouble>(dt)) scalar = scalar*dt;
    60396037                for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
    … …  
    60426040                if(reCast<bool,IssmDouble>(dt)){
    60436041                        old_wh_input->GetInputValue(&water_head,gauss);
    6044                         scalar = Jdet*gauss->weight*water_head;
     6042                        scalar = Jdet*gauss->weight*water_head*sediment_storing;
    60456043                        for(int i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
    60466044                }
    … …  
    60916089void  Tria::InputUpdateFromSolutionHydrology(IssmDouble* solution){
    60926090
     6091        int hydrology_model;
     6092  parameters->FindParam(&hydrology_model,HydrologyEnum);
     6093       
     6094  switch(hydrology_model){
     6095        case HydrologyshreveEnum:
     6096                return InputUpdateFromSolutionHydrologyShreve(solution);
     6097        case HydrologydcEnum:
     6098                return InputUpdateFromSolutionHydrologyDC(solution);
     6099        default:
     6100                _error_("Approximation " << EnumToStringx(hydrology_model) << " not supported yet");
     6101  }
     6102}
     6103/*}}}*/
     6104/*FUNCTION Tria::InputUpdateFromSolutionHydrologyShreve{{{*/
     6105void  Tria::InputUpdateFromSolutionHydrologyShreve(IssmDouble* solution){
     6106
    60936107        /*Intermediaries*/
    60946108        const int   numdof         = NDOF1 *NUMVERTICES;
    … …  
    61086122        /*Add input to the element: */
    61096123        this->inputs->AddInput(new TriaP1Input(WatercolumnEnum,values));
     6124
     6125        /*Free ressources:*/
     6126        xDelete<int>(doflist);
     6127}
     6128/*}}}*/
     6129/*FUNCTION Tria::InputUpdateFromSolutionHydrologyDC{{{*/
     6130void  Tria::InputUpdateFromSolutionHydrologyDC(IssmDouble* solution){
     6131
     6132        /*Intermediaries*/
     6133        const int   numdof         = NDOF1 *NUMVERTICES;
     6134        int        *doflist        = NULL;
     6135        IssmDouble  values[numdof];
     6136
     6137        /*Get dof list: */
     6138        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     6139
     6140        /*Use the dof list to index into the solution vector: */
     6141        for(int i=0;i<numdof;i++){
     6142                values[i]=solution[doflist[i]];
     6143                if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
     6144        }
     6145
     6146        /*Add input to the element: */
     6147        this->inputs->AddInput(new TriaP1Input(SedimentHeadEnum,values));
    61106148
    61116149        /*Free ressources:*/

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

    r14650 r14695  
    244244                ElementVector* CreatePVectorHydrologyShreve(void);
    245245                ElementVector* CreatePVectorHydrologyDC(void);
    246                 void      CreateHydrologyWaterVelocityInput(void);
     246                void    CreateHydrologyWaterVelocityInput(void);
    247247                void      GetSolutionFromInputsHydrology(Vector<IssmDouble>* solution);
    248248                void      InputUpdateFromSolutionHydrology(IssmDouble* solution);
     249                void      InputUpdateFromSolutionHydrologyShreve(IssmDouble* solution);
     250                void      InputUpdateFromSolutionHydrologyDC(IssmDouble* solution);
    249251                #endif
    250252                #ifdef _HAVE_BALANCED_

  • issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.cpp

    r14655 r14695  
    315315}               
    316316/*}}}*/
    317 /*FUNCTION Matpar::GetSedimentCompressibility {{{*/
    318 IssmDouble Matpar::GetSedimentCompressibility(){
    319         return sediment_compressibility;                 
    320 }               
    321 /*}}}*/
    322 /*FUNCTION Matpar::GetSedimentPorosity {{{*/
    323 IssmDouble Matpar::GetSedimentPorosity(){
    324         return sediment_porosity;               
    325 }               
    326 /*}}}*/
    327 /*FUNCTION Matpar::GetSedimentThickness {{{*/
    328 IssmDouble Matpar::GetSedimentThickness(){
    329         return sediment_thickness;               
     317/*FUNCTION Matpar::GetSedimentStoring {{{*/
     318IssmDouble Matpar::GetSedimentStoring(){
     319        return this->rho_freshwater* this->g* this->sediment_porosity* this->sediment_thickness*
     320    ( this->water_compressibility+( this->sediment_compressibility/ this->sediment_porosity));           
    330321}               
    331322/*}}}*/
    … …  
    334325        return sediment_transmitivity;           
    335326}               
    336 /*}}}*/
    337 /*FUNCTION Matpar::GetWaterCompressibility {{{*/
    338 IssmDouble Matpar::GetWaterCompressibility(){
    339         return water_compressibility;           
    340 }               
    341327/*}}}*/                 
    342328/*FUNCTION Matpar::TMeltingPoint {{{*/

  • issm/trunk-jpl/src/c/classes/objects/Materials/Matpar.h

    r14660 r14695  
    106106                IssmDouble GetHydrologyCR();
    107107                IssmDouble GetHydrologyN();
    108                 IssmDouble GetSedimentCompressibility();
    109                 IssmDouble GetSedimentPorosity();
    110                 IssmDouble GetSedimentThickness();
     108                IssmDouble GetSedimentStoring();
    111109                IssmDouble GetSedimentTransmitivity();
    112                 IssmDouble GetWaterCompressibility();
    113110                IssmDouble TMeltingPoint(IssmDouble pressure);
    114111                IssmDouble PureIceEnthalpy(IssmDouble pressure);
Note: See TracChangeset for help on using the changeset viewer.