Changeset 17021


Ignore:
Timestamp:
12/16/13 10:31:04 (11 years ago)
Author:
bdef
Message:

CHG:the sediment transmitivity is now a mat and no more a constant

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

Legend:

Unmodified
Added
Removed

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

    r17012 r17021  
    247247                scalar = Jdet*gauss->weight*(-transfer);
    248248                if(dt!=0.) scalar = scalar*dt;
    249                 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
    250 
     249               
     250                for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
     251               
    251252                /*Transient term*/
    252253                if(dt!=0.){
    … …  
    254255                        old_wh_input->GetInputValue(&water_head,gauss);
    255256                        scalar = Jdet*gauss->weight*water_head*epl_specificstoring*epl_thickness;
    256                         for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
     257                       
     258                        for(int i=0;i<numnodes;i++)pe->values[i]+=scalar*basis[i];
    257259                }
    258260        }

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

    r17005 r17021  
    8585        iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
    8686        iomodel->FetchDataToInput(elements,SedimentHeadEnum);
     87        iomodel->FetchDataToInput(elements,HydrologydcSedimentTransmitivityEnum);
     88
    8789        if(isefficientlayer)iomodel->FetchDataToInput(elements,HydrologydcMaskEplactiveNodeEnum);
    8890}/*}}}*/
    … …  
    168170        /*Intermediaries */
    169171        IssmDouble  D_scalar,Jdet,dt;
     172        IssmDouble  sediment_transmitivity;
    170173        IssmDouble *xyz_list  = NULL;
    171174
    … …  
    177180        IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
    178181        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    179         IssmDouble     D[2][2]={0.};
     182        IssmDouble     D[2][2]= {0.};
    180183
    181184        /*Retrieve all inputs and parameters*/
    182         basalelement->GetVerticesCoordinates(&xyz_list);
    183         IssmDouble sediment_storing       = SedimentStoring(basalelement);
    184         IssmDouble sediment_transmitivity = basalelement->GetMaterialParameter(HydrologydcSedimentTransmitivityEnum);
    185         basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
     185        basalelement                       ->GetVerticesCoordinates(&xyz_list);
     186        IssmDouble sediment_storing        = SedimentStoring(basalelement);
     187        Input* SedTrans_input=basalelement ->GetInput(HydrologydcSedimentTransmitivityEnum); _assert_(SedTrans_input);
     188        basalelement                       ->FindParam(&dt,TimesteppingTimeStepEnum);
    186189
    187190        /* Start  looping on the number of gaussian points: */
    188191        Gauss* gauss=basalelement->NewGauss(2);
    189         for(int ig=gauss->begin();ig<gauss->end();ig++){
    190                 gauss->GaussPoint(ig);
    191 
    192                 basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     192        for(int ig=gauss -> begin();ig<gauss->end();ig++){
     193                gauss          -> GaussPoint(ig);
     194                basalelement   -> JacobianDeterminant(&Jdet,xyz_list,gauss);
     195                SedTrans_input -> GetInputValue(&sediment_transmitivity,gauss);
    193196
    194197                /*Diffusivity*/

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

    r17012 r17021  
    46604660        bool       active_element;
    46614661        int        transfermethod,step;
    4662         IssmDouble sed_trans,sed_thick;
    46634662        IssmDouble leakage,h_max;
    4664         IssmDouble wh_trans;
     4663        IssmDouble wh_trans,sed_thick;
    46654664        IssmDouble activeEpl[numdof],epl_thickness[numdof];
    46664665        IssmDouble epl_specificstoring[numdof],sedstoring[numdof];
    46674666        IssmDouble epl_head[numdof],sed_head[numdof];
    4668         IssmDouble old_transfer[numdof];
     4667        IssmDouble old_transfer[numdof],sed_trans[numdof];
    46694668
    46704669        Input* active_element_input=NULL;
    … …  
    46944693
    46954694                        GetInputListOnVertices(&sed_head[0],SedimentHeadEnum);
     4695                        GetInputListOnVertices(&sed_trans[0],HydrologydcSedimentTransmitivityEnum);
    46964696                        GetInputListOnVertices(&epl_head[0],EplHeadEnum);
    46974697                        GetInputListOnVertices(&epl_thickness[0],HydrologydcEplThicknessEnum);
    … …  
    46994699                        this->parameters->FindParam(&leakage,HydrologydcLeakageFactorEnum);
    47004700
    4701                         sed_trans = matpar->GetSedimentTransmitivity();
    47024701                        sed_thick = matpar->GetSedimentThickness();
    47034702
    47044703                        if(!active_element){
    47054704                                /*No transfer if the EPL is not active*/
    4706                                 /* for(int i=0;i<numdof;i++){ */
    4707                                 /*      wh_trans=0.0; */
    4708                                 /*      /\*Assign output pointer*\/ */
    4709                                 /*      transfer->SetValue(doflist[i],wh_trans,INS_VAL); */
    4710                                 /* } */
    47114705                        }
    47124706                        else{
    47134707                       
    4714                                 //GetInputListOnVertices(&old_transfer[0],WaterTransferEnum);
     4708                                GetInputListOnVertices(&old_transfer[0],WaterTransferEnum);
    47154709                       
    47164710                                for(int i=0;i<numdof;i++){
    … …  
    47184712                                        sedstoring[i]=matpar->GetSedimentStoring();
    47194713                                        this->GetHydrologyDCInefficientHmax(&h_max,i);
    4720                                                
     4714                                                                               
     4715                                        //avoiding transfer at first time
     4716                                        if(epl_head[i]==0.0)epl_head[i]=sed_head[i];
    47214717                                        /*EPL head higher than sediment head, transfer from the epl to the sediment*/
    47224718                                        if(epl_head[i]>sed_head[i]){
    4723                                                 wh_trans=epl_specificstoring[i]*epl_thickness[i]*sed_trans*(epl_head[i]-sed_head[i])/(leakage*sed_thick);                               
    4724 
     4719                                                wh_trans=epl_specificstoring[i]*epl_thickness[i]*sed_trans[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick);                           
    47254720                                                /*No transfer if the sediment head is allready at the maximum*/
    47264721                                                if(sed_head[i]>=h_max){
    … …  
    47304725                                        /*EPL head lower than sediment head, transfer from the sediment to the epl*/
    47314726                                        else if(epl_head[i]<=sed_head[i]){
    4732                                                 wh_trans=sedstoring[i]*sed_trans*(epl_head[i]-sed_head[i])/(leakage*sed_thick);
     4727                                                wh_trans=sedstoring[i]*sed_trans[i]*(epl_head[i]-sed_head[i])/(leakage*sed_thick);
    47334728                                        }
    47344729
    47354730                                        /*Introduce relaxation*/
    4736                                         //wh_trans=old_transfer[i]+0.9*(wh_trans-old_transfer[i]);
     4731                                        wh_trans=old_transfer[i]+0.66*(wh_trans-old_transfer[i]);
    47374732                               
    47384733                                        /*Assign output pointer*/
    47394734                                        transfer->SetValue(doflist[i],wh_trans,INS_VAL);
     4735
     4736                                       
    47404737                                }
    47414738                        }
    … …  
    48814878                                        /*And proceed to the real thing*/
    48824879                                        thickness[i] = old_thickness[i]*(1+((rho_water*gravity*dt)/(rho_ice*latentheat))*epl_conductivity*EPLgrad2-2.0*(A*dt/(pow(n,n)))*(pow(EPL_N,n)));
    4883 
    48844880                        }
    48854881                }

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

    r17006 r17021  
    5959                iomodel->Constant(&this->sediment_porosity,HydrologydcSedimentPorosityEnum);
    6060                iomodel->Constant(&this->sediment_thickness,HydrologydcSedimentThicknessEnum);
    61                 iomodel->Constant(&this->sediment_transmitivity,HydrologydcSedimentTransmitivityEnum);
    6261                iomodel->Constant(&this->water_compressibility,HydrologydcWaterCompressibilityEnum);
    6362                iomodel->Constant(&isefficientlayer,HydrologydcIsefficientlayerEnum);
    … …  
    262261                case HydrologydcEplInitialThicknessEnum:     return this->epl_init_thickness;
    263262                case HydrologydcWaterCompressibilityEnum:    return this->water_compressibility;
    264                 case HydrologydcSedimentTransmitivityEnum:   return this->sediment_transmitivity;
    265263                case HydrologyshreveCREnum:                  return this->hydro_CR;
    266264                case HydrologyshreveKnEnum:                  return this->hydro_kn;
    … …  
    390388        return this->rho_freshwater * this->g * this->epl_porosity *
    391389    (this->water_compressibility + (this->epl_compressibility / this->epl_porosity));           
    392 }               
    393 /*}}}*/
    394 /*FUNCTION Matpar::GetSedimentTransitivity {{{*/
    395 IssmDouble Matpar::GetSedimentTransmitivity(){
    396         return sediment_transmitivity;           
    397 }               
     390}                       
    398391/*}}}*/
    399392/*FUNCTION Matpar::GetSedimentThickness {{{*/

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

    r16812 r17021  
    4444                IssmDouble  sediment_porosity;   
    4545                IssmDouble  sediment_thickness;
    46                 IssmDouble  sediment_transmitivity;     
    4746                IssmDouble  water_compressibility;
    4847
    … …  
    122121                IssmDouble GetSedimentStoring();
    123122                IssmDouble GetEplSpecificStoring();
    124                 IssmDouble GetSedimentTransmitivity();
    125123                IssmDouble GetSedimentThickness();
    126124                IssmDouble GetEplConductivity();

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

    r17012 r17021  
    9797                sedconverged=false;
    9898                for(;;){
    99 
     99                        /* if(isefficientlayer){  */
     100                        /*      /\*Updating Elemental Mask*\/ */
     101                        /*      HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis(); */
     102                        /*      analysis->ElementizeEplMask(femmodel); */
     103                        /*      delete analysis; */
     104                        /*      femmodel->HydrologyTransferx(); */
     105                        /* } */
    100106                        uf_sed_sub_iter=uf_sed->Duplicate();_assert_(uf_sed_sub_iter);
    101107                        uf_sed->Copy(uf_sed_sub_iter);
    … …  
    123129
    124130                        if(sedconverged){
    125                                 if(isefficientlayer){
    126                                         /*Updating Elemental Mask*/
    127                                         HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis();
    128                                         analysis->ElementizeEplMask(femmodel);
    129                                         delete analysis;
    130                                         femmodel->HydrologyTransferx();
    131                                 }
    132131                                sedconverged=false;     
    133132                                /*Checking convegence on the value of the sediment head*/
    … …  
    149148
    150149                        if(sedconverged){
     150                                if(isefficientlayer){
     151                                        /*Updating Elemental Mask*/
     152                                        HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis();
     153                                        analysis->ElementizeEplMask(femmodel);
     154                                        delete analysis;
     155                                        femmodel->HydrologyTransferx();
     156                                }
    151157                                femmodel->parameters->SetParam(sediment_kmax,HydrologySedimentKmaxEnum);
    152158                                InputUpdateFromConstantx(femmodel,sedconverged,ConvergedEnum);
    … …  
    159165                /*Second layer*/
    160166                if(isefficientlayer){
    161 
     167                       
    162168                        femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum);
    163169                        InputUpdateFromConstantx(femmodel,true,ResetPenaltiesEnum);
    164170                        InputUpdateFromConstantx(femmodel,false,ConvergedEnum);
    165171                        femmodel->parameters->SetParam(HydrologyEfficientEnum,HydrologyLayerEnum);
     172                       
     173                               
    166174
    167175                        /*Iteration on the EPL layer*/
    168176                        eplconverged = false;
    169177                        for(;;){
    170 
    171                         /*Start by retrieving the EPL head slopes*/
     178                               
     179                                //updating mask after the computation of the epl thickness (Allow to close too thin EPL)
     180                                femmodel->HydrologyEPLupdateDomainx();
     181                                /*Start by retrieving the EPL head slopes and compute EPL Thickness*/
    172182                                if(VerboseSolution()) _printf0_("computing EPL Head slope...\n");
    173183                                femmodel->SetCurrentConfiguration(L2ProjectionEPLAnalysisEnum);
    … …  
    177187                                femmodel->parameters->SetParam(EplHeadSlopeYEnum,InputToL2ProjectEnum);
    178188                                solutionsequence_linear(femmodel);
    179 
    180                                
    181189                                femmodel->HydrologyEPLThicknessx();
    182 
     190                               
    183191                                femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum);
    184 
    185                                 //updating mask after the computation of the epl thickness (Allow to close too thin EPL)
    186                                 femmodel->HydrologyEPLupdateDomainx();
     192                       
     193                                /*Updating Elemental Mask*/
     194                                HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis();
     195                                analysis->ElementizeEplMask(femmodel);
     196                                delete analysis;
     197                                femmodel->HydrologyTransferx();
    187198                                       
    188199                                ug_epl_sub_iter=ug_epl->Duplicate();_assert_(ug_epl_sub_iter);
    189200                                ug_epl->Copy(ug_epl_sub_iter);
     201                               
    190202                                SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
    191203                                CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum);
    … …  
    213225                                eplcount++;
    214226
    215 
    216227                                if(eplconverged){
    217228                                        eplconverged=false;
    … …  
    232243
    233244                                if(eplconverged){
    234                                         /*Updating Elemental Mask*/
    235                                         HydrologyDCInefficientAnalysis* analysis = new HydrologyDCInefficientAnalysis();
    236                                         analysis->ElementizeEplMask(femmodel);
    237                                         delete analysis;
    238                                         femmodel->HydrologyTransferx();
    239                                        
     245                               
    240246                                        InputUpdateFromConstantx(femmodel,eplconverged,ConvergedEnum);
    241247                                        InputUpdateFromConstantx(femmodel,sediment_kmax,MeltingOffsetEnum);
Note: See TracChangeset for help on using the changeset viewer.