source: issm/oecreview/Archive/16554-17801/ISSM-16637-16638.diff@ 17802

Last change on this file since 17802 was 17802, checked in by Mathieu Morlighem, 11 years ago

Added archives
File size: 34.1 KB

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

     
    5656        iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
    5757        iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
    5858        iomodel->FetchDataToInput(elements,EplHeadEnum);
     59        iomodel->FetchDataToInput(elements,HydrologydcEplThicknessEnum);
     60       
    5961}/*}}}*/
    6062void HydrologyDCEfficientAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
    6163

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

     
    113113        HydrologydcEplCompressibilityEnum,
    114114        HydrologydcEplPorosityEnum,
    115115        HydrologydcEplThicknessEnum,
    116         HydrologydcEplTransmitivityEnum,
     116        HydrologydcEplThicknessOldEnum,
     117        HydrologydcEplConductivityEnum,
    117118        HydrologydcIsefficientlayerEnum,
    118119        HydrologydcSedimentlimitFlagEnum,
    119120        HydrologydcSedimentlimitEnum,

  • ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

     
    121121                case HydrologydcEplCompressibilityEnum : return "HydrologydcEplCompressibility";
    122122                case HydrologydcEplPorosityEnum : return "HydrologydcEplPorosity";
    123123                case HydrologydcEplThicknessEnum : return "HydrologydcEplThickness";
    124                 case HydrologydcEplTransmitivityEnum : return "HydrologydcEplTransmitivity";
     124                case HydrologydcEplThicknessOldEnum : return "HydrologydcEplThicknessOld";
     125                case HydrologydcEplConductivityEnum : return "HydrologydcEplConductivity";
    125126                case HydrologydcIsefficientlayerEnum : return "HydrologydcIsefficientlayer";
    126127                case HydrologydcSedimentlimitFlagEnum : return "HydrologydcSedimentlimitFlag";
    127128                case HydrologydcSedimentlimitEnum : return "HydrologydcSedimentlimit";

  • ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

     
    121121              else if (strcmp(name,"HydrologydcEplCompressibility")==0) return HydrologydcEplCompressibilityEnum;
    122122              else if (strcmp(name,"HydrologydcEplPorosity")==0) return HydrologydcEplPorosityEnum;
    123123              else if (strcmp(name,"HydrologydcEplThickness")==0) return HydrologydcEplThicknessEnum;
    124               else if (strcmp(name,"HydrologydcEplTransmitivity")==0) return HydrologydcEplTransmitivityEnum;
     124              else if (strcmp(name,"HydrologydcEplThicknessOld")==0) return HydrologydcEplThicknessOldEnum;
     125              else if (strcmp(name,"HydrologydcEplConductivity")==0) return HydrologydcEplConductivityEnum;
    125126              else if (strcmp(name,"HydrologydcIsefficientlayer")==0) return HydrologydcIsefficientlayerEnum;
    126127              else if (strcmp(name,"HydrologydcSedimentlimitFlag")==0) return HydrologydcSedimentlimitFlagEnum;
    127128              else if (strcmp(name,"HydrologydcSedimentlimit")==0) return HydrologydcSedimentlimitEnum;
    … …  
    135136              else if (strcmp(name,"HydrologySedimentKmax")==0) return HydrologySedimentKmaxEnum;
    136137              else if (strcmp(name,"WaterTransfer")==0) return WaterTransferEnum;
    137138              else if (strcmp(name,"IndependentObject")==0) return IndependentObjectEnum;
    138               else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum;
    139139         else stage=2;
    140140   }
    141141   if(stage==2){
    142               if (strcmp(name,"InversionCostFunctionThreshold")==0) return InversionCostFunctionThresholdEnum;
     142              if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum;
     143              else if (strcmp(name,"InversionCostFunctionThreshold")==0) return InversionCostFunctionThresholdEnum;
    143144              else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum;
    144145              else if (strcmp(name,"InversionCostFunctions")==0) return InversionCostFunctionsEnum;
    145146              else if (strcmp(name,"InversionGradientScaling")==0) return InversionGradientScalingEnum;
    … …  
    258259              else if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum;
    259260              else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum;
    260261              else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum;
    261               else if (strcmp(name,"Surface")==0) return SurfaceEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"SurfaceforcingsPrecipitation")==0) return SurfaceforcingsPrecipitationEnum;
     265              if (strcmp(name,"Surface")==0) return SurfaceEnum;
     266              else if (strcmp(name,"SurfaceforcingsPrecipitation")==0) return SurfaceforcingsPrecipitationEnum;
    266267              else if (strcmp(name,"SurfaceforcingsMassBalance")==0) return SurfaceforcingsMassBalanceEnum;
    267268              else if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum;
    268269              else if (strcmp(name,"SurfaceforcingsDesfac")==0) return SurfaceforcingsDesfacEnum;
    … …  
    381382              else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum;
    382383              else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum;
    383384              else if (strcmp(name,"Element")==0) return ElementEnum;
    384               else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"FileParam")==0) return FileParamEnum;
     388              if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum;
     389              else if (strcmp(name,"FileParam")==0) return FileParamEnum;
    389390              else if (strcmp(name,"Input")==0) return InputEnum;
    390391              else if (strcmp(name,"IntInput")==0) return IntInputEnum;
    391392              else if (strcmp(name,"InputToExtrude")==0) return InputToExtrudeEnum;
    … …  
    504505              else if (strcmp(name,"BasalFriction")==0) return BasalFrictionEnum;
    505506              else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum;
    506507              else if (strcmp(name,"QmuTemperature")==0) return QmuTemperatureEnum;
    507               else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum;
     511              if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum;
     512              else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum;
    512513              else if (strcmp(name,"StressTensor")==0) return StressTensorEnum;
    513514              else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum;
    514515              else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;

  • ../trunk-jpl/src/c/modules/InputDuplicatex/InputDuplicatex.cpp

     
    88#include "../../toolkits/toolkits.h"
    99
    1010void InputDuplicatex(FemModel* femmodel,int original_enum, int new_enum){
    11 
    1211        /*Go through elemnets, and ask to reinitialie the input: */
    1312        for(int i=0;i<femmodel->elements->Size();i++){
    1413                Element* element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));

  • ../trunk-jpl/src/c/solutionsequences/solutionsequence_hydro_nonlinear.cpp

     
    128128                        femmodel->HydrologyEPLupdateDomainx();
    129129                        femmodel->parameters->SetParam(HydrologyEfficientEnum,HydrologyLayerEnum);
    130130
     131                        /*Reset constraint on the ZigZag Lock*/
     132                        ResetConstraintsx(femmodel);
    131133                        /*Iteration on the EPL layer*/
    132134                        eplconverged = false;
    133135                        for(;;){
     136                                femmodel->HydrologyEPLThicknessx();
    134137                                femmodel->HydrologyTransferx();
    135138                                SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
    136139                                CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum);

  • ../trunk-jpl/src/c/cores/hydrology_core.cpp

     
    8383                        femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
    8484                        if (isefficientlayer){
    8585                                InputDuplicatex(femmodel,EplHeadEnum,EplHeadOldEnum);
     86                                InputDuplicatex(femmodel,HydrologydcEplThicknessEnum,HydrologydcEplThicknessOldEnum);
    8687                        }
    8788                       
    8889                        /*Proceed now to heads computations*/
    … …  
    9192                        if(save_results && ((i+1)%output_frequency==0 || (i+1)==nsteps)){
    9293                                if(VerboseSolution()) _printf0_("   saving results \n");
    9394                                if(isefficientlayer){
    94                                         int outputs[6] = {SedimentHeadEnum,SedimentHeadResidualEnum,EplHeadEnum,HydrologydcMaskEplactiveEnum,EplHeadSlopeXEnum,EplHeadSlopeYEnum};
    95                                         femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],6);
     95                                        int outputs[7] = {SedimentHeadEnum,SedimentHeadResidualEnum,EplHeadEnum,HydrologydcMaskEplactiveEnum,EplHeadSlopeXEnum,EplHeadSlopeYEnum,HydrologydcEplThicknessEnum};
     96                                        femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],7);
    9697                                }
    9798                                else{
    9899                                        int outputs[2] = {SedimentHeadEnum,SedimentHeadResidualEnum};

  • ../trunk-jpl/src/c/classes/Materials/Matice.cpp

     
    144144        return pow(B,-n);
    145145}
    146146/*}}}*/
     147/*FUNCTION Matice::GetAbar {{{*/
     148IssmDouble Matice::GetAbar(){
     149        /*
     150         * A = 1/B^n
     151         */
     152
     153        IssmDouble B,n;
     154
     155        inputs->GetInputAverage(&B,MaterialsRheologyBbarEnum);
     156        n=this->GetN();
     157
     158        return pow(B,-n);
     159}
     160/*}}}*/
    147161/*FUNCTION Matice::GetB {{{*/
    148162IssmDouble Matice::GetB(){
    149163

  • ../trunk-jpl/src/c/classes/Materials/Material.h

     
    3535                virtual void       GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* pepsilon)=0;
    3636                virtual void       GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon)=0;
    3737                virtual IssmDouble GetA()=0;
     38                virtual IssmDouble GetAbar()=0;
    3839                virtual IssmDouble GetB()=0;
    3940                virtual IssmDouble GetBbar()=0;
    4041                virtual IssmDouble GetN()=0;

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

     
    6565                if(isefficientlayer){
    6666                                iomodel->Constant(&this->epl_compressibility,HydrologydcEplCompressibilityEnum);
    6767                                iomodel->Constant(&this->epl_porosity,HydrologydcEplPorosityEnum);
    68                                 iomodel->Constant(&this->epl_thickness,HydrologydcEplThicknessEnum);
    69                                 iomodel->Constant(&this->epl_transmitivity,HydrologydcEplTransmitivityEnum);
     68                                iomodel->Constant(&this->epl_conductivity,HydrologydcEplConductivityEnum);
    7069                }
    7170        }
    7271        else{
    … …  
    353352    ( this->water_compressibility+( this->sediment_compressibility/ this->sediment_porosity));           
    354353}               
    355354/*}}}*/
    356 /*FUNCTION Matpar::GetEplStoring {{{*/
    357 IssmDouble Matpar::GetEplStoring(){
    358         return this->rho_freshwater* this->g* this->epl_porosity* this->epl_thickness*
     355/*FUNCTION Matpar::GetEplSpecificStoring {{{*/
     356IssmDouble Matpar::GetEplSpecificStoring(){
     357        return this->rho_freshwater* this->g* this->epl_porosity*
    359358    ( this->water_compressibility+( this->epl_compressibility/ this->epl_porosity));             
    360359}               
    361360/*}}}*/
    … …  
    368367IssmDouble Matpar::GetSedimentThickness(){
    369368        return sediment_thickness;               
    370369}               
    371 /*}}}*/ 
    372 /*FUNCTION Matpar::GetEplTransitivity {{{*/
    373 IssmDouble Matpar::GetEplTransmitivity(){
    374         return epl_transmitivity;               
     370/*}}}*/ 
     371/*FUNCTION Matpar::GetEplConductivity {{{*/
     372IssmDouble Matpar::GetEplConductivity(){
     373        return epl_conductivity;                 
    375374}               
    376375/*}}}*/                 
    377376/*FUNCTION Matpar::TMeltingPoint {{{*/

  • ../trunk-jpl/src/c/classes/Materials/Matice.h

     
    6262                void       GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* pepsilon);
    6363                void       GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
    6464                IssmDouble GetA();
     65                IssmDouble GetAbar();
    6566                IssmDouble GetB();
    6667                IssmDouble GetBbar();
    6768                IssmDouble GetD();

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

     
    4848
    4949                IssmDouble  epl_compressibility;
    5050                IssmDouble  epl_porosity;       
    51                 IssmDouble  epl_thickness;
    52                 IssmDouble  epl_transmitivity;   
     51                IssmDouble  epl_conductivity;   
    5352
    5453                /*gia: */
    5554                IssmDouble lithosphere_shear_modulus;
    … …  
    9392                void       GetViscosityDerivativeEpsSquareFS(IssmDouble* pmu_prime, IssmDouble* pepsilon){_error_("not supported");};
    9493                void       GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon){_error_("not supported");};
    9594                IssmDouble GetA(){_error_("not supported");};
     95                IssmDouble GetAbar(){_error_("not supported");};
    9696                IssmDouble GetB(){_error_("not supported");};
    9797                IssmDouble GetBbar(){_error_("not supported");};
    9898                IssmDouble GetN(){_error_("not supported");};
    … …  
    120120                IssmDouble GetHydrologyCR();
    121121                IssmDouble GetHydrologyN();
    122122                IssmDouble GetSedimentStoring();
    123                 IssmDouble GetEplStoring();
     123                IssmDouble GetEplSpecificStoring();
    124124                IssmDouble GetSedimentTransmitivity();
    125125                IssmDouble GetSedimentThickness();
    126                 IssmDouble GetEplTransmitivity();
     126                IssmDouble GetEplConductivity();
    127127                IssmDouble TMeltingPoint(IssmDouble pressure);
    128128                IssmDouble PureIceEnthalpy(IssmDouble pressure);
    129129                IssmDouble GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure);

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

     
    15771577void  Tria::InputDuplicate(int original_enum,int new_enum){
    15781578
    15791579        /*Call inputs method*/
    1580         if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum);
    1581 
     1580        if (IsInput(original_enum)) {
     1581                inputs->DuplicateInput(original_enum,new_enum);
     1582        }
    15821583}
    15831584/*}}}*/
    15841585/*FUNCTION Tria::InputScale{{{*/
    … …  
    20962097                                name==SedimentHeadEnum ||
    20972098                                name==EplHeadOldEnum ||
    20982099                                name==EplHeadEnum ||
     2100                                name==HydrologydcEplThicknessOldEnum ||
     2101                                name==HydrologydcEplThicknessEnum ||
    20992102                                name==HydrologydcMaskEplactiveEnum ||
    21002103                                name==MeshVertexonbedEnum ||
    21012104                                name==WaterTransferEnum ||
    … …  
    24542457                                this->ComputeStressTensor();
    24552458                                input=this->inputs->GetInput(output_enum);
    24562459                                break;
    2457                 case SurfaceNormalVelocityEnum:
    2458                               this->ComputeSurfaceNormalVelocity();
    2459                                 input=this->inputs->GetInput(output_enum);
    2460                                 break;
    24612460                        default:
    24622461                                _error_("input "<<EnumToStringx(output_enum)<<" not found in element");
    24632462                }
    … …  
    24892488                        for(int i=0;i<NUMVERTICES;i++) values[i] = values[i]/reCast<IssmDouble>(connectivity[i]);
    24902489
    24912490                        vector->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL);
    2492                        
    2493                         /*FIXME: Inputs needs to be updated first in Tria::ResultInterpolation, this is a hack*/
    2494                         this->inputs->DeleteInput(SurfaceNormalVelocityEnum);
    24952491                        break;
    24962492                                        }
    24972493                default:
    … …  
    26502646        *(surface_normal+2) = normal[2]/normal_norm;
    26512647}
    26522648/*}}}*/
    2653 /*FUNCTION Tria::ComputeSurfaceNormalVelocity{{{*/
    2654 void Tria::ComputeSurfaceNormalVelocity(){
    2655 
    2656   IssmDouble      sum,tangential_vector[2],normal_vector[2],time,ux,uy;
    2657   IssmDouble      normal_velocity[NUMVERTICES],xyz_list[NUMVERTICES][3];
    2658   IssmDouble      value[NUMVERTICES],verticesonsurface[NUMVERTICES];
    2659 
    2660   for(int iv=0;iv<NUMVERTICES;iv++){
    2661     normal_velocity[iv]=0.;
    2662     value[iv]=0.;
    2663   }
    2664 
    2665   GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    2666 
    2667   GaussTria* gauss=new GaussTria();
    2668   Input* slope_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(slope_input);
    2669   //  Input* slope_input= inputs->GetInput(SurfaceEnum); _assert_(slope_input);
    2670   Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    2671   Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    2672 
    2673 
    2674   /*Get list of nodes on surface*/
    2675   GetInputListOnVertices(&verticesonsurface[0],MeshVertexonsurfaceEnum); 
    2676   sum = verticesonsurface[0]+verticesonsurface[1]+verticesonsurface[2];
    2677   _assert_(sum==0. || sum==1. || sum==2.);
    2678 
    2679   /*Compute normal velocity for surface nodes from L2 projected slope*/
    2680   for(int iv=0;iv<NUMVERTICES;iv++){
    2681     if(verticesonsurface[iv] == 1){
    2682         gauss->GaussNode(P1Enum,iv);   
    2683         slope_input->GetInputValue(&value[iv],gauss);
    2684         //      slope_input->GetInputDerivativeValue(&value[iv],&xyz_list[0][0],gauss);
    2685         vx_input->GetInputValue(&ux,gauss);
    2686         vy_input->GetInputValue(&uy,gauss);
    2687         tangential_vector[0]=sqrt(1./(pow(value[iv],2.)+1.));
    2688         tangential_vector[1]=value[iv]*tangential_vector[0];
    2689         normal_vector[0]=-1.*tangential_vector[1];
    2690         normal_vector[1]=tangential_vector[0];
    2691         normal_velocity[iv]=ux*normal_vector[0]+uy*normal_vector[1];
    2692     }
    2693   }
    2694 
    2695   delete gauss;
    2696   this->inputs->AddInput(new TriaInput(SurfaceNormalVelocityEnum,&normal_velocity[0],P1Enum));
    2697 
    2698 }
    2699 /*}}}*/
    27002649/*FUNCTION Tria::TimeAdapt{{{*/
    27012650IssmDouble  Tria::TimeAdapt(void){
    27022651
    … …  
    32773226        /*Now get the average SMB over the element*/
    32783227        Input* smb_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(smb_input);
    32793228        smb_input->GetInputAverage(&smb);                                                                                                                                                                                               // average smb on element in m ice s-1
    3280         Total_Smb=rho_ice*base*smb;                                                                                                                                                                                                                     // smb on element in kg s-1
     3229   Total_Smb=rho_ice*base*smb;                                                                                                                                                                                                                  // smb on element in kg s-1
    32813230
    32823231        /*Return: */
    32833232        return Total_Smb;
    … …  
    67336682ElementMatrix* Tria::CreateKMatrixHydrologyDCEfficient(void){
    67346683
    67356684        /* Intermediaries */
    6736         IssmDouble  D_scalar,Jdet;
    6737         IssmDouble      epl_transmitivity,dt;
    6738         IssmDouble  epl_storing;
     6685        IssmDouble  D_scalar,Jdet,dt;
     6686        IssmDouble  epl_thickness;
     6687        IssmDouble      epl_conductivity;
     6688        IssmDouble  epl_specificstoring;
    67396689        IssmDouble  xyz_list[NUMVERTICES][3];
    67406690
    67416691        /*Check that all nodes are active, else return empty matrix*/
    … …  
    67546704
    67556705        /*Retrieve all inputs and parameters*/
    67566706        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     6707        Input* thickness_input=inputs->GetInput(HydrologydcEplThicknessEnum); _assert_(thickness);
    67576708        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    6758         epl_storing       = matpar->GetEplStoring();
    6759         epl_transmitivity = matpar->GetEplTransmitivity();
     6709        epl_specificstoring = matpar->GetEplSpecificStoring();
     6710        epl_conductivity    = matpar->GetEplConductivity();
    67606711
     6712
    67616713        /* Start looping on the number of gaussian points: */
    67626714        GaussTria* gauss=new GaussTria(2);
    67636715        for(int ig=gauss->begin();ig<gauss->end();ig++){
    6764 
     6716               
     6717               
    67656718                gauss->GaussPoint(ig);
    67666719                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     6720                thickness_input->GetInputValue(&epl_thickness,gauss);
    67676721
    67686722                /*Diffusivity*/
    6769                 D_scalar=epl_transmitivity*gauss->weight*Jdet;
     6723                D_scalar=epl_conductivity*epl_thickness*gauss->weight*Jdet;
    67706724                if(reCast<bool,IssmDouble>(dt)) D_scalar=D_scalar*dt;
    67716725                D[0][0]=D_scalar; D[0][1]=0.;
    67726726                D[1][0]=0.;       D[1][1]=D_scalar;
    … …  
    67796733                /*Transient*/
    67806734                if(reCast<bool,IssmDouble>(dt)){
    67816735                        GetNodalFunctions(basis,gauss);
    6782                         D_scalar=epl_storing*gauss->weight*Jdet;
     6736                        D_scalar=epl_specificstoring*epl_thickness*gauss->weight*Jdet;
    67836737
    67846738                        TripleMultiply(basis,numnodes,1,0,
    67856739                                                &D_scalar,1,1,0,
    … …  
    69156869        IssmDouble xyz_list[NUMVERTICES][3];
    69166870        IssmDouble dt,scalar,water_head;
    69176871        IssmDouble transfer,residual;
    6918         IssmDouble epl_storing;
     6872        IssmDouble epl_thickness;
     6873        IssmDouble epl_specificstoring;
    69196874        GaussTria* gauss = NULL;
    69206875
    69216876        /*Check that all nodes are active, else return empty matrix*/
    … …  
    69316886        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    69326887
    69336888        /*Retrieve all inputs and parameters*/
    6934         epl_storing = matpar->GetEplStoring();
     6889        epl_specificstoring = matpar->GetEplSpecificStoring();
    69356890        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    69366891        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    6937         Input* residual_input=inputs->GetInput(SedimentHeadResidualEnum);  _assert_(residual_input);
    6938         Input* transfer_input=inputs->GetInput(WaterTransferEnum);  _assert_(transfer_input);
     6892        Input* residual_input=inputs->GetInput(SedimentHeadResidualEnum);     _assert_(residual_input);
     6893        Input* transfer_input=inputs->GetInput(WaterTransferEnum);            _assert_(transfer_input);
     6894        Input* thickness_input=inputs->GetInput(HydrologydcEplThicknessEnum); _assert_(thickness_input);
    69396895        Input* old_wh_input=NULL;
    69406896
    69416897        if(reCast<bool,IssmDouble>(dt)){
    … …  
    69586914
    69596915                /*Transient term*/
    69606916                if(reCast<bool,IssmDouble>(dt)){
     6917                        thickness_input->GetInputValue(&epl_thickness,gauss);
    69616918                        old_wh_input->GetInputValue(&water_head,gauss);
    6962                         scalar = Jdet*gauss->weight*water_head*epl_storing;
     6919                        scalar = Jdet*gauss->weight*water_head*epl_specificstoring*epl_thickness;
    69636920                        for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
    69646921                }
    69656922        }
    … …  
    71417098        IssmDouble sed_trans,sed_thick;
    71427099        IssmDouble leakage,h_max;
    71437100        IssmDouble wh_trans;
    7144         IssmDouble activeEpl[numdof];
    7145         IssmDouble eplstoring[numdof],sedstoring[numdof];
     7101        IssmDouble activeEpl[numdof],epl_thickness[numdof];
     7102        IssmDouble epl_specificstoring[numdof],sedstoring[numdof];
    71467103        IssmDouble epl_head[numdof],sed_head[numdof];
    71477104
    71487105        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    … …  
    71647121                        GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum);
    71657122                        GetInputListOnVertices(&sed_head[0],SedimentHeadEnum);
    71667123                        GetInputListOnVertices(&epl_head[0],EplHeadEnum);
     7124                        GetInputListOnVertices(&epl_thickness[0],HydrologydcEplThicknessEnum);
    71677125
    71687126                        this->parameters->FindParam(&leakage,HydrologydcLeakageFactorEnum);
    71697127
    … …  
    71767134                                        wh_trans=0.0;
    71777135                                }
    71787136                                else{
    7179                                         eplstoring[i]=matpar->GetEplStoring();         
     7137                                        epl_specificstoring[i]=matpar->GetEplSpecificStoring();         
    71807138                                        sedstoring[i]=matpar->GetSedimentStoring();
    71817139
    71827140                                        /*EPL head higher than sediment head, transfer from the epl to the sediment*/
    71837141                                        if(epl_head[i]>sed_head[i]){
    7184                                                 wh_trans=eplstoring[i]*sed_trans*(epl_head[i]-sed_head[i])/(leakage*sed_thick);                         
     7142                                                wh_trans=epl_specificstoring[i]*epl_thickness[i]*sed_trans*(epl_head[i]-sed_head[i])/(leakage*sed_thick);                               
    71857143
    71867144                                                /*No transfer if the sediment head is allready at the maximum*/
    71877145                                                this->GetHydrologyDCInefficientHmax(&h_max,nodes[i]);
    … …  
    72457203        GetInputListOnVertices(&eplhead[0],EplHeadEnum);
    72467204        GetInputListOnVertices(&residual[0],SedimentHeadResidualEnum);
    72477205
    7248         /*Get minimum sediment head*/
     7206        /*Get minimum sediment head of the element*/
    72497207        sedheadmin=sedhead[0];
    72507208        for(i=1;i<numdof;i++) if(sedhead[i]<=sedheadmin)sedheadmin=sedhead[i];
    72517209
    … …  
    72797237/*}}}*/
    72807238/*FUNCTION Tria::ComputeEPLThickness{{{*/
    72817239void  Tria::ComputeEPLThickness(void){
     7240
     7241        int         i;
     7242        const int   numdof         = NDOF1 *NUMVERTICES;
     7243        bool        isefficientlayer;
     7244        IssmDouble  n,A,dt;
     7245        IssmDouble  rho_water,rho_ice;
     7246        IssmDouble  gravity,latentheat,EPLgrad;
     7247        IssmDouble  EPL_N,epl_conductivity;
     7248        IssmDouble  activeEpl[numdof],thickness[numdof];
     7249        IssmDouble  eplhead[numdof], old_thickness[numdof];
     7250        IssmDouble  epl_slopeX[numdof],epl_slopeY[numdof];
     7251        IssmDouble  ice_thickness[numdof],bed[numdof];
     7252
     7253
     7254        /*Get the flag to know if the efficient layer is present*/
     7255        this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
     7256        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     7257
     7258        if(isefficientlayer){
     7259                /*For now, assuming just one way to compute EPL thickness*/
     7260                rho_water  = matpar->GetRhoWater();
     7261                rho_ice    = matpar->GetRhoIce();
     7262                gravity    = matpar->GetG();
     7263                latentheat = matpar->GetLatentHeat();
     7264                n          = material->GetN();
     7265                A          = material->GetAbar();
     7266               
     7267                GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum);
     7268                GetInputListOnVertices(&eplhead[0],EplHeadEnum);
     7269                GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum);
     7270                GetInputListOnVertices(&epl_slopeY[0],EplHeadSlopeYEnum);
     7271                GetInputListOnVertices(&old_thickness[0],HydrologydcEplThicknessOldEnum);
     7272                GetInputListOnVertices(&ice_thickness[0],ThicknessEnum);
     7273                GetInputListOnVertices(&bed[0],BedEnum);
     7274               
     7275                epl_conductivity = matpar->GetEplConductivity();
     7276               
     7277                for(int i=0;i<numdof;i++){
     7278                        /*Keeping thickness to 1 if EPL is not active*/
     7279                        if(activeEpl[i]==0.0){
     7280                                thickness[i]=1.0;
     7281                        }
     7282                        else{
     7283
     7284                                /*Compute first the effective pressure in the EPL*/
     7285                                EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i])));
     7286                               
     7287                                /*Get then the gradient of EPL heads*/
     7288                                EPLgrad = epl_slopeX[i]+epl_slopeY[i];
     7289                               
     7290                                /*And proceed to the real thing*/
     7291                                thickness[i] = old_thickness[i]*(1-((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*pow(EPLgrad,2.0)-2.0*(A*dt/(pow(n,n)))*EPL_N);
     7292
     7293                        }
     7294                }
     7295                this->inputs->AddInput(new TriaInput(HydrologydcEplThicknessEnum,thickness,P1Enum));
     7296        }
    72827297}
    72837298/*}}}*/
    72847299

  • ../trunk-jpl/src/c/classes/Elements/Penta.cpp

     
    26442644                                name==EplHeadEnum ||
    26452645                                name==SedimentHeadOldEnum ||
    26462646                                name==EplHeadOldEnum ||
     2647                                name==HydrologydcEplThicknessOldEnum ||
     2648                                name==HydrologydcEplThicknessEnum ||
    26472649                                name==HydrologydcMaskEplactiveEnum ||
    26482650                                name==WaterTransferEnum
    26492651
    … …  
    28702872
    28712873        switch(input->GetResultInterpolation()){
    28722874                case P0Enum:
    2873                         _error_("not implemented...");
     2875                        _error_("P0 not implemented yet for input "<<EnumToStringx(output_enum));
    28742876                        break;
    28752877                case P1Enum:{
    28762878                        IssmDouble  values[NUMVERTICES];
    … …  
    1084110843/*FUNCTION Penta::ComputeEPLThickness{{{*/
    1084210844void  Penta::ComputeEPLThickness(void){
    1084310845
     10846        int         i;
     10847        const int   numdof   = NDOF1 *NUMVERTICES;
     10848        const int   numdof2d = NDOF1 *NUMVERTICES2D;
     10849        bool        isefficientlayer;
     10850        IssmDouble  n,A,dt;
     10851        IssmDouble  rho_water,rho_ice;
     10852        IssmDouble  gravity,latentheat,EPLgrad;
     10853        IssmDouble  EPL_N,epl_conductivity;
     10854        IssmDouble  activeEpl[numdof],thickness[numdof];
     10855        IssmDouble  eplhead[numdof], old_thickness[numdof];
     10856        IssmDouble  epl_slopeX[numdof],epl_slopeY[numdof];
     10857        IssmDouble  ice_thickness[numdof],bed[numdof];
     10858        Penta       *penta = NULL;
     10859        /*If not on bed, return*/
    1084410860        if (!IsOnBed())return;
    1084510861
    10846         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    10847         tria->ComputeEPLThickness();
    10848         delete tria->material; delete tria;
     10862        /*Get the flag to know if the efficient layer is present*/
     10863        this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
     10864        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    1084910865
     10866        if(isefficientlayer){
     10867                /*For now, assuming just one way to compute EPL thickness*/
     10868                rho_water        = matpar->GetRhoWater();
     10869                rho_ice          = matpar->GetRhoIce();
     10870                gravity          = matpar->GetG();
     10871                latentheat       = matpar->GetLatentHeat();
     10872                epl_conductivity = matpar->GetEplConductivity();
     10873                n                = material->GetN();
     10874                A                = material->GetA();
     10875               
     10876               
     10877                GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum);
     10878                GetInputListOnVertices(&eplhead[0],EplHeadEnum);
     10879                GetInputListOnVertices(&epl_slopeX[0],EplHeadSlopeXEnum);
     10880                GetInputListOnVertices(&epl_slopeY[0],EplHeadSlopeYEnum);
     10881                GetInputListOnVertices(&old_thickness[0],HydrologydcEplThicknessOldEnum);
     10882                GetInputListOnVertices(&ice_thickness[0],ThicknessEnum);
     10883                GetInputListOnVertices(&bed[0],BedEnum);
     10884               
     10885                for(int i=0;i<numdof2d;i++){
     10886                        /*Keeping thickness to 1 if EPL is not active*/
     10887                        if(activeEpl[i]==0.0){
     10888                                thickness[i]=1.0;
     10889                                thickness[i+numdof2d]=thickness[i];
     10890                        }
     10891                        else{
     10892
     10893                                /*Compute first the effective pressure in the EPL*/
     10894                                EPL_N=gravity*((rho_ice*ice_thickness[i])-(rho_water*(eplhead[i]-bed[i])));
     10895                               
     10896                                /*Get then the gradient of EPL heads*/
     10897                                EPLgrad = epl_slopeX[i]+epl_slopeY[i];
     10898                               
     10899                                /*And proceed to the real thing*/
     10900                                thickness[i] = old_thickness[i]*(1-((rho_water*gravity*dt)/(rho_ice* latentheat))*epl_conductivity*pow(EPLgrad,2.0)-2.0*(A*dt/(pow(n,n)))*EPL_N);
     10901                                thickness[i+numdof2d]=thickness[i];
     10902                        }
     10903                }
     10904                penta=this;
     10905                for(;;){
     10906
     10907                        /*Add input to the element: */                 
     10908                        penta->inputs->AddInput(new PentaInput(HydrologydcEplThicknessEnum,thickness,P1Enum));
     10909                       
     10910                        /*Stop if we have reached the surface*/
     10911                        if (penta->IsOnSurface()) break;
     10912                       
     10913                        /* get upper Penta*/
     10914                        penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id);
     10915                }
     10916        }
    1085010917}
    1085110918/*}}}*/
    1085210919/*FUNCTION Penta::InputUpdateFromSolutionHydrologyDCInefficient{{{*/

  • ../trunk-jpl/src/m/classes/hydrologydc.m

     
    2626                epl_compressibility      = 0;
    2727                epl_porosity             = 0;
    2828                epl_thickness            = 0;
    29                 epl_transmitivity        = 0;
     29                epl_conductivity         = 0;
    3030  end
    3131        methods
    3232                % {{{ function obj = hydrologydc(varargin)
    … …  
    6060                        obj.epl_compressibility      = 1.0e-08;
    6161                        obj.epl_porosity             = 0.4;
    6262                        obj.epl_thickness            = 1.0;
    63                         obj.epl_transmitivity        = 8.0e-02;
     63                        obj.epl_conductivity         = 8.0e-02;
    6464
    6565                end
    6666                % }}}
    … …  
    9696                                md = checkfield(md,'hydrology.epl_compressibility','>',0,'numel',1);
    9797                                md = checkfield(md,'hydrology.epl_porosity','>',0,'numel',1);
    9898                                md = checkfield(md,'hydrology.epl_thickness','>',0,'numel',1);
    99                                 md = checkfield(md,'hydrology.epl_transmitivity','>',0,'numel',1);
     99                                md = checkfield(md,'hydrology.epl_conductivity','>',0,'numel',1);
    100100            end
    101101                end
    102102                % }}}
    … …  
    136136                                fielddisplay(obj,'mask_eplactive','active (1) or not (0) EPL');
    137137                                fielddisplay(obj,'epl_compressibility','epl compressibility [Pa^-1]');
    138138                                fielddisplay(obj,'epl_porosity','epl [dimensionless]');
    139                                 fielddisplay(obj,'epl_thickness','epl thickness [m]');
    140                                 fielddisplay(obj,'epl_transmitivity','epl transmitivity [m^2/s]');
     139                                fielddisplay(obj,'epl_thickness','epl initial thickness [m]');
     140                                fielddisplay(obj,'epl_conductivity','epl conductivity [m^2/s]');
    141141            end
    142142
    143143                end
    … …  
    171171                                WriteData(fid,'object',obj,'fieldname','epl_compressibility','format','Double');                       
    172172                                WriteData(fid,'object',obj,'fieldname','epl_porosity','format','Double');                       
    173173                                WriteData(fid,'object',obj,'fieldname','epl_thickness','format','Double');
    174                                 WriteData(fid,'object',obj,'fieldname','epl_transmitivity','format','Double');
     174                                %WriteData(fid,'object',obj,'fieldname','epl_transmitivity','format','Double');
     175                                WriteData(fid,'object',obj,'fieldname','epl_conductivity','format','Double');
    175176                        end
    176177                end
    177178% }}}

  • ../trunk-jpl/src/m/enum/HydrologydcEplConductivityEnum.m

     
     1function macro=HydrologydcEplConductivityEnum()
     2%HYDROLOGYDCEPLCONDUCTIVITYENUM - Enum of HydrologydcEplConductivity
     3%
     4%   WARNING: DO NOT MODIFY THIS FILE
     5%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
     6%            Please read src/c/shared/Enum/README for more information
     7%
     8%   Usage:
     9%      macro=HydrologydcEplConductivityEnum()
     10
     11macro=StringToEnum('HydrologydcEplConductivity');

  • ../trunk-jpl/src/m/enum/EnumDefinitions.py

     
    113113def HydrologydcEplCompressibilityEnum(): return StringToEnum("HydrologydcEplCompressibility")[0]
    114114def HydrologydcEplPorosityEnum(): return StringToEnum("HydrologydcEplPorosity")[0]
    115115def HydrologydcEplThicknessEnum(): return StringToEnum("HydrologydcEplThickness")[0]
    116 def HydrologydcEplTransmitivityEnum(): return StringToEnum("HydrologydcEplTransmitivity")[0]
     116def HydrologydcEplThicknessOldEnum(): return StringToEnum("HydrologydcEplThicknessOld")[0]
     117def HydrologydcEplConductivityEnum(): return StringToEnum("HydrologydcEplConductivity")[0]
    117118def HydrologydcIsefficientlayerEnum(): return StringToEnum("HydrologydcIsefficientlayer")[0]
    118119def HydrologydcSedimentlimitFlagEnum(): return StringToEnum("HydrologydcSedimentlimitFlag")[0]
    119120def HydrologydcSedimentlimitEnum(): return StringToEnum("HydrologydcSedimentlimit")[0]

  • ../trunk-jpl/src/m/enum/HydrologydcEplThicknessOldEnum.m

     
     1function macro=HydrologydcEplThicknessOldEnum()
     2%HYDROLOGYDCEPLTHICKNESSOLDENUM - Enum of HydrologydcEplThicknessOld
     3%
     4%   WARNING: DO NOT MODIFY THIS FILE
     5%            this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
     6%            Please read src/c/shared/Enum/README for more information
     7%
     8%   Usage:
     9%      macro=HydrologydcEplThicknessOldEnum()
     10
     11macro=StringToEnum('HydrologydcEplThicknessOld');
Note: See TracBrowser for help on using the repository browser.