Changeset 16853


Ignore:
Timestamp:
11/21/13 11:31:43 (11 years ago)
Author:
Mathieu Morlighem
Message:

NEW: PVector DC

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

Legend:

Unmodified
Added
Removed

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

    r16782 r16853  
    100100}/*}}}*/
    101101ElementVector* HydrologyDCEfficientAnalysis::CreatePVector(Element* element){/*{{{*/
    102 _error_("not implemented yet");
     102
     103        /*Intermediaries*/
     104        int      meshtype;
     105        Element* basalelement;
     106
     107        /*Get basal element*/
     108        element->FindParam(&meshtype,MeshTypeEnum);
     109        switch(meshtype){
     110                case Mesh2DhorizontalEnum:
     111                        basalelement = element;
     112                        break;
     113                case Mesh3DEnum:
     114                        if(!element->IsOnBed()) return NULL;
     115                        basalelement = element->SpawnBasalElement();
     116                        break;
     117                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     118        }
     119
     120        /*Check that all nodes are active, else return empty matrix*/
     121        if(element->AllActive()) return NULL;
     122
     123        /*Intermediaries */
     124        IssmDouble dt,scalar,water_head;
     125        IssmDouble transfer,residual,epl_thickness;
     126        IssmDouble Jdet;
     127        IssmDouble *xyz_list     = NULL;
     128        Input*      old_wh_input = NULL;
     129
     130        /*Fetch number of nodes and dof for this finite element*/
     131        int numnodes = basalelement->GetNumberOfNodes();
     132
     133        /*Initialize Element vector*/
     134        ElementVector* pe    = basalelement->NewElementVector();
     135        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     136
     137        /*Retrieve all inputs and parameters*/
     138        basalelement->GetVerticesCoordinates(&xyz_list);
     139        IssmDouble epl_specificstoring = EplSpecificStoring(basalelement);
     140        basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
     141        Input* residual_input    = element->GetInput(SedimentHeadResidualEnum);    _assert_(residual_input);
     142        Input* transfer_input    = element->GetInput(WaterTransferEnum);           _assert_(transfer_input);
     143        Input* thickness_input   = element->GetInput(HydrologydcEplThicknessEnum); _assert_(thickness_input);
     144        if(dt!= 0.){old_wh_input = element->GetInput(EplHeadOldEnum);              _assert_(old_wh_input);}
     145
     146        /* Start  looping on the number of gaussian points: */
     147        Gauss* gauss=basalelement->NewGauss(2);
     148        for(int ig=gauss->begin();ig<gauss->end();ig++){
     149                gauss->GaussPoint(ig);
     150
     151                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     152                basalelement->NodalFunctions(basis,gauss);
     153
     154                /*Loading term*/
     155                transfer_input->GetInputValue(&transfer,gauss);
     156                scalar = Jdet*gauss->weight*(-transfer);
     157                if(dt!=0.) scalar = scalar*dt;
     158                for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
     159
     160                /*Transient term*/
     161                if(dt!=0.){
     162                        thickness_input->GetInputValue(&epl_thickness,gauss);
     163                        old_wh_input->GetInputValue(&water_head,gauss);
     164                        scalar = Jdet*gauss->weight*water_head*epl_specificstoring*epl_thickness;
     165                        for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
     166                }
     167        }
     168
     169        /*Clean up and return*/
     170        xDelete<IssmDouble>(xyz_list);
     171        xDelete<IssmDouble>(basis);
     172        delete gauss;
     173        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     174        return pe;
    103175}/*}}}*/
    104176void HydrologyDCEfficientAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    … …  
    118190        }
    119191}/*}}}*/
     192
     193/*Intermediaries*/
     194IssmDouble HydrologyDCEfficientAnalysis::EplSpecificStoring(Element* element){/*{{{*/
     195        IssmDouble rho_freshwater        = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
     196        IssmDouble g                     = element->GetMaterialParameter(ConstantsGEnum);
     197        IssmDouble epl_porosity          = element->GetMaterialParameter(HydrologydcEplPorosityEnum);
     198        IssmDouble epl_compressibility   = element->GetMaterialParameter(HydrologydcEplCompressibilityEnum);
     199        IssmDouble water_compressibility = element->GetMaterialParameter(HydrologydcWaterCompressibilityEnum);
     200        return rho_freshwater*g*epl_porosity*(water_compressibility+(epl_compressibility/epl_porosity));                 
     201}/*}}}*/

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

    r16782 r16853  
    2525                void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    2626                void InputUpdateFromSolution(IssmDouble* solution,Element* element);
     27
     28                /*Intermediaries*/
     29                IssmDouble EplSpecificStoring(Element* element);
    2730};
    2831#endif

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

    r16782 r16853  
    141141}/*}}}*/
    142142ElementVector* HydrologyDCInefficientAnalysis::CreatePVector(Element* element){/*{{{*/
    143 _error_("not implemented yet");
     143
     144        /*Intermediaries*/
     145        int      meshtype;
     146        Element* basalelement;
     147
     148        /*Get basal element*/
     149        element->FindParam(&meshtype,MeshTypeEnum);
     150        switch(meshtype){
     151                case Mesh2DhorizontalEnum:
     152                        basalelement = element;
     153                        break;
     154                case Mesh3DEnum:
     155                        if(!element->IsOnBed()) return NULL;
     156                        basalelement = element->SpawnBasalElement();
     157                        break;
     158                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     159        }
     160
     161        /*Intermediaries */
     162        IssmDouble dt,scalar,water_head;
     163        IssmDouble water_load,transfer;
     164        IssmDouble Jdet;
     165        IssmDouble *xyz_list  = NULL;
     166        Input*      old_wh_input      = NULL;
     167
     168        /*Fetch number of nodes and dof for this finite element*/
     169        int numnodes = basalelement->GetNumberOfNodes();
     170
     171        /*Initialize Element vector*/
     172        ElementVector* pe    = basalelement->NewElementVector();
     173        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     174
     175        /*Retrieve all inputs and parameters*/
     176        basalelement->GetVerticesCoordinates(&xyz_list);
     177        IssmDouble sediment_storing = SedimentStoring(basalelement);
     178        basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
     179        Input* water_input       = element->GetInput(BasalforcingsMeltingRateEnum); _assert_(water_input);
     180        Input* transfer_input    = element->GetInput(WaterTransferEnum);            _assert_(transfer_input);
     181        if(dt!= 0.){old_wh_input = element->GetInput(SedimentHeadOldEnum);          _assert_(old_wh_input);}
     182
     183        /* Start  looping on the number of gaussian points: */
     184        Gauss* gauss=basalelement->NewGauss(2);
     185        for(int ig=gauss->begin();ig<gauss->end();ig++){
     186                gauss->GaussPoint(ig);
     187
     188                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     189                basalelement->NodalFunctions(basis,gauss);
     190
     191                /*Loading term*/
     192                water_input->GetInputValue(&water_load,gauss);
     193                transfer_input->GetInputValue(&transfer,gauss);
     194                scalar = Jdet*gauss->weight*(water_load+transfer);
     195                if(dt!=0.) scalar = scalar*dt;
     196                for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
     197
     198                /*Transient term*/
     199                if(dt!=0.){
     200                        old_wh_input->GetInputValue(&water_head,gauss);
     201                        scalar = Jdet*gauss->weight*water_head*sediment_storing;
     202                        for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
     203                }
     204        }
     205
     206        /*Clean up and return*/
     207        xDelete<IssmDouble>(xyz_list);
     208        xDelete<IssmDouble>(basis);
     209        delete gauss;
     210        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     211        return pe;
    144212}/*}}}*/
    145213void HydrologyDCInefficientAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    … …  
    204272        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    205273}/*}}}*/
     274
     275/*Intermediaries*/
     276IssmDouble HydrologyDCInefficientAnalysis::SedimentStoring(Element* element){/*{{{*/
     277        IssmDouble rho_freshwater           = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
     278        IssmDouble g                        = element->GetMaterialParameter(ConstantsGEnum);
     279        IssmDouble sediment_porosity        = element->GetMaterialParameter(HydrologydcSedimentPorosityEnum);
     280        IssmDouble sediment_thickness       = element->GetMaterialParameter(HydrologydcSedimentThicknessEnum);
     281        IssmDouble sediment_compressibility = element->GetMaterialParameter(HydrologydcSedimentCompressibilityEnum);
     282        IssmDouble water_compressibility    = element->GetMaterialParameter(HydrologydcWaterCompressibilityEnum);
     283        return rho_freshwater*g*sediment_porosity*sediment_thickness*(water_compressibility+(sediment_compressibility/sediment_porosity));               
     284}/*}}}*/

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

    r16782 r16853  
    2525                void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    2626                void InputUpdateFromSolution(IssmDouble* solution,Element* element);
     27
     28                /*Intermediaries*/
     29                IssmDouble SedimentStoring(Element* element);
    2730};
    2831#endif

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

    r16782 r16853  
    8686}/*}}}*/
    8787ElementVector* HydrologyShreveAnalysis::CreatePVector(Element* element){/*{{{*/
    88 _error_("not implemented yet");
     88
     89        /*Skip if water or ice shelf element*/
     90        if(element->IsFloating()) return NULL;
     91
     92        /*Intermediaries */
     93        IssmDouble  Jdet,dt;
     94        IssmDouble  mb,oldw;
     95        IssmDouble* xyz_list = NULL;
     96
     97        /*Fetch number of nodes and dof for this finite element*/
     98        int numnodes = element->GetNumberOfNodes();
     99
     100        /*Initialize Element vector and other vectors*/
     101        ElementVector* pe    = element->NewElementVector();
     102        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     103
     104        /*Retrieve all inputs and parameters*/
     105        element->GetVerticesCoordinates(&xyz_list);
     106        element->FindParam(&dt,TimesteppingTimeStepEnum);
     107        Input* mb_input   = element->GetInput(BasalforcingsMeltingRateEnum); _assert_(mb_input);
     108        Input* oldw_input = element->GetInput(WaterColumnOldEnum);           _assert_(oldw);
     109
     110        /*Initialize mb_correction to 0, do not forget!:*/
     111        /* Start  looping on the number of gaussian points: */
     112        Gauss* gauss=element->NewGauss(2);
     113        for(int ig=gauss->begin();ig<gauss->end();ig++){
     114                gauss->GaussPoint(ig);
     115
     116                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     117                element->NodalFunctions(basis,gauss);
     118
     119                mb_input->GetInputValue(&mb,gauss);
     120                oldw_input->GetInputValue(&oldw,gauss);
     121
     122                if(dt!=0.){
     123                        for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(oldw+dt*mb)*basis[i];
     124                }
     125                else{
     126                        for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*mb*basis[i];
     127                }
     128        }
     129
     130        /*Clean up and return*/
     131        xDelete<IssmDouble>(xyz_list);
     132        xDelete<IssmDouble>(basis);
     133        delete gauss;
     134        return pe;
    89135}/*}}}*/
    90136void HydrologyShreveAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/

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

    r16782 r16853  
    7474}/*}}}*/
    7575ElementVector* L2ProjectionEPLAnalysis::CreatePVector(Element* element){/*{{{*/
    76 _error_("not implemented yet");
     76
     77        /*Intermediaries*/
     78        int      meshtype;
     79        Element* basalelement;
     80
     81        /*Get basal element*/
     82        element->FindParam(&meshtype,MeshTypeEnum);
     83        switch(meshtype){
     84                case Mesh2DhorizontalEnum:
     85                        basalelement = element;
     86                        break;
     87                case Mesh3DEnum:
     88                        if(!element->IsOnBed()) return NULL;
     89                        basalelement = element->SpawnBasalElement();
     90                        break;
     91                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     92        }
     93
     94        /*Intermediaries */
     95        int         input_enum,index;
     96        IssmDouble  Jdet,slopes[2];
     97        Input      *input     = NULL;
     98        IssmDouble *xyz_list  = NULL;
     99
     100        /*Fetch number of nodes and dof for this finite element*/
     101        int numnodes = basalelement->GetNumberOfNodes();
     102
     103        /*Initialize Element vector*/
     104        ElementVector* pe    = basalelement->NewElementVector();
     105        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     106
     107        /*Retrieve all inputs and parameters*/
     108        basalelement->GetVerticesCoordinates(&xyz_list);
     109        basalelement->FindParam(&input_enum,InputToL2ProjectEnum);
     110        switch(input_enum){
     111                case EplHeadSlopeXEnum: input = element->GetInput(EplHeadEnum); index = 0; _assert_(input); break;
     112                case EplHeadSlopeYEnum: input = element->GetInput(EplHeadEnum); index = 1; _assert_(input); break;
     113                default: _error_("not implemented");
     114        }
     115
     116        /* Start  looping on the number of gaussian points: */
     117        Gauss* gauss=basalelement->NewGauss(2);
     118        for(int ig=gauss->begin();ig<gauss->end();ig++){
     119                gauss->GaussPoint(ig);
     120
     121                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     122                basalelement->NodalFunctions(basis,gauss);
     123
     124                input->GetInputDerivativeValue(&slopes[0],xyz_list,gauss);
     125                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*slopes[index]*basis[i];
     126        }
     127
     128        /*Clean up and return*/
     129        xDelete<IssmDouble>(xyz_list);
     130        xDelete<IssmDouble>(basis);
     131        delete gauss;
     132        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
     133        return pe;
    77134}/*}}}*/
    78135void L2ProjectionEPLAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/

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

    r16850 r16853  
    3737                virtual        ~Element(){};
    3838
     39                virtual bool        AllActive(void)=0;
     40                virtual bool        AnyActive(void)=0;
    3941                virtual IssmDouble CharacteristicLength(void)=0;
    4042                virtual void   Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0;

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

    r16850 r16853  
    6868                /*}}}*/
    6969                /*Element virtual functions definitions: {{{*/
     70                bool        AllActive(void){_error_("not implemented yet");};
     71                bool        AnyActive(void){_error_("not implemented yet");};
    7072                void   BasalFrictionCreateInput(void);
    7173                IssmDouble CharacteristicLength(void){_error_("not implemented yet");};

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

    r16850 r16853  
    6767                /*}}}*/
    6868                /*Element virtual functions definitions: {{{*/
     69                bool        AllActive(void){_error_("not implemented yet");};
     70                bool        AnyActive(void){_error_("not implemented yet");};
    6971                void        AddBasalInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");};
    7072                void        AddInput(int input_enum, IssmDouble* values, int interpolation_enum){_error_("not implemented yet");};

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

    r16850 r16853  
    6767                /*}}}*/
    6868                /*Element virtual functions definitions: {{{*/
     69                bool        AllActive(void);
     70                bool        AnyActive(void);
    6971                IssmDouble  CharacteristicLength(void);
    7072                void        ComputeBasalStress(Vector<IssmDouble>* sigma_b);
    … …  
    375377                void           HydrologyEPLGetMask(Vector<IssmDouble>* vec_mask);
    376378                void           ComputeEPLThickness(void);
    377                 bool           AllActive(void);
    378                 bool           AnyActive(void);
     379
    379380                #endif
    380381

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

    r16812 r16853  
    254254                case MaterialsMixedLayerCapacityEnum:       return this->mixed_layer_capacity;
    255255                case MaterialsThermalExchangeVelocityEnum:  return this->thermal_exchange_velocity;
     256                case HydrologydcSedimentPorosityEnum:       return this->sediment_porosity;
     257                case HydrologydcSedimentThicknessEnum:      return this->sediment_thickness;
     258                case HydrologydcSedimentCompressibilityEnum:return this->sediment_compressibility;
     259                case HydrologydcEplPorosityEnum:            return this->epl_porosity;
     260                case HydrologydcEplCompressibilityEnum:     return this->epl_compressibility;
     261                case HydrologydcWaterCompressibilityEnum:   return this->water_compressibility;
    256262                case ConstantsGEnum:                        return this->g;
    257263                default: _error_("Enum "<<EnumToStringx(enum_in)<<" not supported yet");
Note: See TracChangeset for help on using the changeset viewer.