Changeset 15062


Ignore:
Timestamp:
05/21/13 17:45:31 (12 years ago)
Author:
bdef
Message:

CHG: adding the convergence test in the outer loop of the hydrological system

Location:
issm/trunk-jpl/src
Files:
1 added
11 edited

Legend:

Unmodified
Added
Removed

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

    r15060 r15062  
    14851485        case HydrologyShreveAnalysisEnum:
    14861486                GetSolutionFromInputsHydrologyShreve(solution);
     1487                break;
     1488        case HydrologyDCInefficientAnalysisEnum:
     1489                GetSolutionFromInputsHydrologyDCInefficient(solution);
     1490                break;
     1491        case HydrologyDCEfficientAnalysisEnum:
     1492                GetSolutionFromInputsHydrologyDCEfficient(solution);
    14871493                break;
    14881494        #endif
    … …  
    63816387}
    63826388/*}}}*/
     6389/*FUNCTION Tria::GetSolutionFromInputsHydrologyDCInefficient{{{*/
     6390void  Tria::GetSolutionFromInputsHydrologyDCInefficient(Vector<IssmDouble>* solution){
     6391
     6392        const int    numdof=NDOF1*NUMVERTICES;
     6393
     6394        int         i;
     6395        int        *doflist = NULL;
     6396        IssmDouble  sedimenthead;
     6397        IssmDouble  values[numdof];
     6398        GaussTria  *gauss   = NULL;
     6399
     6400        /*Get dof list: */
     6401        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     6402
     6403        /*Get inputs*/
     6404        Input* sedimenthead_input=inputs->GetInput(SedimentHeadEnum); _assert_(sedimenthead_input);
     6405
     6406        /*P1 element only for now*/
     6407        gauss=new GaussTria();
     6408        for(i=0;i<NUMVERTICES;i++){
     6409
     6410                gauss->GaussVertex(i);
     6411                sedimenthead_input->GetInputValue(&sedimenthead,gauss);
     6412                values[i]=sedimenthead;
     6413        }
     6414
     6415        solution->SetValues(numdof,doflist,values,INS_VAL);
     6416
     6417        /*Free ressources:*/
     6418        delete gauss;
     6419        xDelete<int>(doflist);
     6420}
     6421/*}}}*/
     6422/*FUNCTION Tria::GetSolutionFromInputsHydrologyDCEfficient{{{*/
     6423void  Tria::GetSolutionFromInputsHydrologyDCEfficient(Vector<IssmDouble>* solution){
     6424
     6425        const int    numdof=NDOF1*NUMVERTICES;
     6426
     6427        int         i;
     6428        int        *doflist = NULL;
     6429        IssmDouble  eplhead;
     6430        IssmDouble  values[numdof];
     6431        GaussTria  *gauss   = NULL;
     6432
     6433        /*Get dof list: */
     6434        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     6435
     6436        /*Get inputs*/
     6437        Input* eplhead_input=inputs->GetInput(EplHeadEnum); _assert_(eplhead_input);
     6438
     6439        /*P1 element only for now*/
     6440        gauss=new GaussTria();
     6441        for(i=0;i<NUMVERTICES;i++){
     6442                gauss->GaussVertex(i);
     6443                eplhead_input->GetInputValue(&eplhead,gauss);
     6444                values[i]=eplhead;
     6445        }
     6446
     6447        solution->SetValues(numdof,doflist,values,INS_VAL);
     6448
     6449        /*Free ressources:*/
     6450        delete gauss;
     6451        xDelete<int>(doflist);
     6452}
     6453/*}}}*/
    63836454/*FUNCTION Tria::InputUpdateFromSolutionHydrologyShreve{{{*/
    63846455void  Tria::InputUpdateFromSolutionHydrologyShreve(IssmDouble* solution){

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

    r15060 r15062  
    249249                ElementVector* CreatePVectorHydrologyDCEfficient(void);
    250250                void      GetSolutionFromInputsHydrologyShreve(Vector<IssmDouble>* solution);
     251                void      GetSolutionFromInputsHydrologyDCInefficient(Vector<IssmDouble>* solution);
     252                void      GetSolutionFromInputsHydrologyDCEfficient(Vector<IssmDouble>* solution);
    251253                void    CreateHydrologyWaterVelocityInput(void);
    252254                void      InputUpdateFromSolutionHydrology(IssmDouble* solution);

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateParametersHydrologyDCInefficient.cpp

    r15020 r15062  
    1818        IssmDouble  penalty_factor;
    1919        IssmDouble  leakagefactor;
     20        IssmDouble  rel_tol;
    2021
    2122        /*Get parameters: */
    … …  
    3536        iomodel->FetchData(&transfer_flag,HydrologydcTransferFlagEnum);
    3637        iomodel->FetchData(&penalty_factor,HydrologydcPenaltyFactorEnum);
     38        iomodel->FetchData(&rel_tol,HydrologydcRelTolEnum);
    3739
    3840        if(sedimentlimit_flag==1){
    … …  
    5153        parameters->AddObject(new IntParam(HydrologydcSedimentlimitFlagEnum,sedimentlimit_flag));
    5254        parameters->AddObject(new IntParam(HydrologydcTransferFlagEnum,transfer_flag));
     55        parameters->AddObject(new DoubleParam(HydrologydcRelTolEnum,rel_tol));
    5356
    5457        /*Assign output pointer: */

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

    r15049 r15062  
    9090        SedimentHeadResidualEnum,
    9191        EplHeadEnum,
     92  HydrologydcRelTolEnum,
    9293        HydrologydcSpcsedimentHeadEnum,
    9394        HydrologydcSedimentCompressibilityEnum,

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

    r15049 r15062  
    9898                case SedimentHeadResidualEnum : return "SedimentHeadResidual";
    9999                case EplHeadEnum : return "EplHead";
     100                case HydrologydcRelTolEnum : return "HydrologydcRelTol";
    100101                case HydrologydcSpcsedimentHeadEnum : return "HydrologydcSpcsedimentHead";
    101102                case HydrologydcSedimentCompressibilityEnum : return "HydrologydcSedimentCompressibility";

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

    r15049 r15062  
    9898              else if (strcmp(name,"SedimentHeadResidual")==0) return SedimentHeadResidualEnum;
    9999              else if (strcmp(name,"EplHead")==0) return EplHeadEnum;
     100              else if (strcmp(name,"HydrologydcRelTol")==0) return HydrologydcRelTolEnum;
    100101              else if (strcmp(name,"HydrologydcSpcsedimentHead")==0) return HydrologydcSpcsedimentHeadEnum;
    101102              else if (strcmp(name,"HydrologydcSedimentCompressibility")==0) return HydrologydcSedimentCompressibilityEnum;
    … …  
    136137              else if (strcmp(name,"InversionNsteps")==0) return InversionNstepsEnum;
    137138              else if (strcmp(name,"InversionNumControlParameters")==0) return InversionNumControlParametersEnum;
    138               else if (strcmp(name,"InversionNumCostFunctions")==0) return InversionNumCostFunctionsEnum;
    139139         else stage=2;
    140140   }
    141141   if(stage==2){
    142               if (strcmp(name,"InversionStepThreshold")==0) return InversionStepThresholdEnum;
     142              if (strcmp(name,"InversionNumCostFunctions")==0) return InversionNumCostFunctionsEnum;
     143              else if (strcmp(name,"InversionStepThreshold")==0) return InversionStepThresholdEnum;
    143144              else if (strcmp(name,"InversionThicknessObs")==0) return InversionThicknessObsEnum;
    144145              else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum;
    … …  
    259260              else if (strcmp(name,"TimesteppingFinalTime")==0) return TimesteppingFinalTimeEnum;
    260261              else if (strcmp(name,"TimesteppingCflCoefficient")==0) return TimesteppingCflCoefficientEnum;
    261               else if (strcmp(name,"TimesteppingTimeAdapt")==0) return TimesteppingTimeAdaptEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"TimesteppingTimeStep")==0) return TimesteppingTimeStepEnum;
     265              if (strcmp(name,"TimesteppingTimeAdapt")==0) return TimesteppingTimeAdaptEnum;
     266              else if (strcmp(name,"TimesteppingTimeStep")==0) return TimesteppingTimeStepEnum;
    266267              else if (strcmp(name,"TransientIsdiagnostic")==0) return TransientIsdiagnosticEnum;
    267268              else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum;
    … …  
    382383              else if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum;
    383384              else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum;
    384               else if (strcmp(name,"StokesIceFront")==0) return StokesIceFrontEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
     388              if (strcmp(name,"StokesIceFront")==0) return StokesIceFrontEnum;
     389              else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
    389390              else if (strcmp(name,"StringParam")==0) return StringParamEnum;
    390391              else if (strcmp(name,"Tria")==0) return TriaEnum;
    … …  
    505506              else if (strcmp(name,"PatchNodes")==0) return PatchNodesEnum;
    506507              else if (strcmp(name,"PatchVertices")==0) return PatchVerticesEnum;
    507               else if (strcmp(name,"PentaP1ElementResult")==0) return PentaP1ElementResultEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum;
     511              if (strcmp(name,"PentaP1ElementResult")==0) return PentaP1ElementResultEnum;
     512              else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum;
    512513              else if (strcmp(name,"Step")==0) return StepEnum;
    513514              else if (strcmp(name,"Time")==0) return TimeEnum;

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

    r15020 r15062  
    99
    1010void solutionsequence_hydro_nonlinear(FemModel* femmodel){
    11 
     11       
    1212        /*solution : */
    1313        Vector<IssmDouble>* ug=NULL;
    14         Vector<IssmDouble>* uf=NULL;
    15         Vector<IssmDouble>* uf_old=NULL;
     14        Vector<IssmDouble>* uf_sed=NULL;
     15        Vector<IssmDouble>* uf_epl=NULL;
     16        Vector<IssmDouble>* old_uf=NULL;
     17        Vector<IssmDouble>* aged_uf_sed=NULL;
     18        Vector<IssmDouble>* aged_uf_epl=NULL;
    1619        Vector<IssmDouble>* ys=NULL;
    17         IssmDouble sediment_kmax,time;
    18 
    19         /*intermediary: */
     20        Vector<IssmDouble>* duf=NULL;
     21       
    2022        Matrix<IssmDouble>* Kff=NULL;
    2123        Matrix<IssmDouble>* Kfs=NULL;
    2224        Vector<IssmDouble>* pf=NULL;
    2325        Vector<IssmDouble>* df=NULL;
    24 
    25         bool converged,isefficientlayer;
    26         int  constraints_converged;
    27         int  num_unstable_constraints;
    28         int  count;
    29         int  hydro_maxiter;
     26       
     27        bool       sedconverged,eplconverged,hydroconverged;
     28        bool       isefficientlayer;
     29        int        constraints_converged;
     30        int        num_unstable_constraints;
     31        int        count;
     32        int        hydro_maxiter;
     33        IssmDouble sediment_kmax,time;
     34        IssmDouble eps_hyd;
     35        IssmDouble ndu_sed,nu_sed;
     36        IssmDouble ndu_epl,nu_epl;     
    3037
    3138        /*Recover parameters: */
     39        femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);//FIXME
    3240        femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
    33         femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);//FIXME
     41        femmodel->parameters->FindParam(&eps_hyd,HydrologydcRelTolEnum);
     42        femmodel->parameters->FindParam(&time,TimeEnum);
    3443        femmodel->BasisIntegralsx();
    3544        hydro_maxiter=100;
    3645        count=1;
    37         converged=false;
     46        sedconverged=false;
     47        eplconverged=false;
     48        hydroconverged=false;
    3849
    39         femmodel->parameters->FindParam(&time,TimeEnum);
     50        /*Iteration on the two layers + transfer*/
     51        femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);
     52        GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
     53        Reducevectorgtofx(&uf_sed, ug, femmodel->nodes,femmodel->parameters); _assert_(uf_sed);
     54        if(isefficientlayer) {
     55                femmodel->SetCurrentConfiguration(HydrologyDCEfficientAnalysisEnum);
     56                GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
     57                Reducevectorgtofx(&uf_epl, ug, femmodel->nodes,femmodel->parameters);
     58        }
     59
     60
    4061        for(;;){
    41                 /*First layer*/
     62                //save pointer to old velocity
     63                delete aged_uf_sed;aged_uf_sed=uf_sed;
     64                if(isefficientlayer){
     65                        delete aged_uf_epl;
     66                        aged_uf_epl=uf_epl;
     67                }
     68
    4269                femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);
    4370                InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,true,ResetPenaltiesEnum);
    … …  
    4673                femmodel->parameters->SetParam(HydrologySedimentEnum,HydrologyLayerEnum);
    4774                femmodel->HydrologyTransferx();
     75               
     76                /*Iteration on the sediment layer*/
    4877                for(;;){
    49 
    5078                        femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df, &sediment_kmax);
    5179                        CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCInefficientAnalysisEnum);
    52                         Reduceloadx(pf,Kfs,ys); delete Kfs; delete uf;
    53                         Solverx(&uf, Kff, pf,uf_old, df, femmodel->parameters);
    54                         delete uf_old; uf_old=uf->Duplicate();
     80                        Reduceloadx(pf,Kfs,ys); delete Kfs; delete uf_sed;
     81                        Solverx(&uf_sed, Kff, pf,old_uf, df, femmodel->parameters);
     82                        delete old_uf; old_uf=uf_sed->Duplicate();
    5583                        delete Kff; delete pf; delete ug; delete df;
    56                         Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters); delete ys;
     84                       
     85                        Mergesolutionfromftogx(&ug,uf_sed,ys,femmodel->nodes,femmodel->parameters); delete ys;
    5786                        InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
    58 
    5987                        ConstraintsStatex(&constraints_converged,&num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    6088
    61                         if (!converged){
     89                        if (!sedconverged){
    6290                                if(VerboseConvergence()) _pprintLine_("   #unstable constraints = " << num_unstable_constraints);
    63                                 if(num_unstable_constraints==0) converged = true;
     91                                if(num_unstable_constraints==0) sedconverged = true;
    6492                                if (count>=hydro_maxiter){
    6593                                        _error_("   maximum number of iterations (" << hydro_maxiter << ") exceeded");
    … …  
    6896                        count++;
    6997
    70                         if(converged){
     98                        if(sedconverged){
    7199                                femmodel->parameters->SetParam(sediment_kmax,HydrologySedimentKmaxEnum);
    72                                 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,converged,ConvergedEnum);
     100                                InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sedconverged,ConvergedEnum);
    73101                                InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
    74102                                InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,HydrologySedimentKmaxEnum);
    … …  
    84112                femmodel->UpdateConstraintsx();
    85113                femmodel->parameters->SetParam(HydrologyEfficientEnum,HydrologyLayerEnum);
    86                 converged = false;
     114                eplconverged = false;
     115                /*Iteration on the EPL layer*/
    87116                for(;;){
    88117                        femmodel->SystemMatricesx(&Kff, &Kfs, &pf,&df,NULL);
    89118                        CreateNodalConstraintsx(&ys,femmodel->nodes,HydrologyDCEfficientAnalysisEnum);
    90                         Reduceloadx(pf,Kfs,ys); delete Kfs; delete uf;
    91                         Solverx(&uf, Kff, pf,uf_old, df, femmodel->parameters);
    92                         delete uf_old; uf_old=uf->Duplicate();
     119                        Reduceloadx(pf,Kfs,ys); delete Kfs; delete uf_epl;
     120                        Solverx(&uf_epl, Kff, pf,old_uf, df, femmodel->parameters);
     121                        delete old_uf; old_uf=uf_epl->Duplicate();
    93122                        delete Kff;delete pf; delete ug; delete df;
    94                         Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters); delete ys;
     123                        Mergesolutionfromftogx(&ug,uf_epl,ys,femmodel->nodes,femmodel->parameters); delete ys;
    95124                        InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
    96125
    97126                        ConstraintsStatex(&constraints_converged,&num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    98127
    99                         if (!converged){
     128                        if (!eplconverged){
    100129                                if(VerboseConvergence()) _pprintLine_("   #unstable constraints = " << num_unstable_constraints);
    101                                 if(num_unstable_constraints==0) converged = true;
     130                                if(num_unstable_constraints==0) eplconverged = true;
    102131                                if (count>=hydro_maxiter){
    103132                                        _error_("   maximum number of iterations (" << hydro_maxiter << ") exceeded");
    … …  
    106135                        count++;
    107136
    108                         if(converged){
    109                                 InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,converged,ConvergedEnum);
     137                        if(eplconverged){
     138                                InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,eplconverged,ConvergedEnum);
    110139                                InputUpdateFromConstantx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,sediment_kmax,MeltingOffsetEnum);
    111140                                InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
    … …  
    115144
    116145                /*System convergence check*/
    117                 break;
     146                if(!hydroconverged){
     147                        //compute norm(du)/norm(u)
     148                        _assert_(aged_uf_sed); _assert_(uf_sed);
     149                        duf=uf_sed->Duplicate(); _assert_(duf);
     150                        aged_uf_sed->Copy(duf); duf->AYPX(uf_sed,-1.0);
     151                        ndu_sed=duf->Norm(NORM_TWO); nu_sed=aged_uf_sed->Norm(NORM_TWO);
     152                        if (xIsNan<IssmDouble>(ndu_sed) || xIsNan<IssmDouble>(nu_sed)) _error_("convergence criterion is NaN!");
     153                       
     154                        if(isefficientlayer){
     155                                duf=aged_uf_epl->Duplicate(); aged_uf_epl->Copy(duf); duf->AYPX(uf_epl,-1.0);
     156                                ndu_epl=duf->Norm(NORM_TWO); nu_epl=aged_uf_epl->Norm(NORM_TWO);
     157                                if (xIsNan<IssmDouble>(ndu_epl) || xIsNan<IssmDouble>(nu_epl)) _error_("convergence criterion is NaN!");
     158                        }
     159                        //print
     160                        if (!xIsNan<IssmDouble>(eps_hyd)){
     161                                if((ndu_sed/nu_sed)<eps_hyd){
     162                                        if(VerboseConvergence()) _pprintLine_(setw(50) << left << "   Sediment Convergence criterion: norm(du)/norm(u)" << ndu_sed/nu_sed*100 << " < " << eps_hyd*100 << " %");
     163                                        hydroconverged=true;
     164                                }
     165                                else{
     166                                        if(VerboseConvergence()) _pprintLine_(setw(50) << left << "   Sediment Convergence criterion: norm(du)/norm(u)" << ndu_sed/nu_sed*100 << " > " << eps_hyd*100 << " %");
     167                                        hydroconverged=false;
     168                                }
     169                                if(isefficientlayer){
     170                                        if((ndu_epl/nu_epl)<eps_hyd){
     171                                                if(VerboseConvergence()) _pprintLine_(setw(50) << left << "   EPL Convergence criterion: norm(du)/norm(u)" << ndu_epl/nu_epl*100 << " < " << eps_hyd*100 << " %");
     172                                        }
     173                                        else{
     174                                                if(VerboseConvergence()) _pprintLine_(setw(50) << left << "   EPL Convergence criterion: norm(du)/norm(u)" << ndu_epl/nu_epl*100 << " > " << eps_hyd*100 << " %");
     175                                                hydroconverged=false;
     176                                        }
     177                                }
     178                        }
     179                        else _pprintLine_(setw(50) << left << "   Convergence criterion: norm(du)/norm(u)" << ndu_sed/nu_sed*100 << " %");
     180                }
     181       
     182                if(hydroconverged)break;
    118183        }
    119184
    120185        InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug);
    121 
     186       
    122187        /*Free ressources: */
    123188        delete ug;
    124         delete uf;
    125         delete uf_old;
     189        delete uf_sed;
     190        delete uf_epl;
     191        delete old_uf;
     192        delete aged_uf_sed;
     193        delete aged_uf_epl;
     194        delete duf;
     195       
    126196}

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

    r15015 r15062  
    66classdef hydrologydc
    77        properties (SetAccess=public)
     8                water_compressibility    = 0;
     9                isefficientlayer         = 0;
     10                penalty_factor           = 0;
     11                rel_tol                  = 0;
     12                sedimentlimit_flag       = 0;
     13                sedimentlimit            = 0;
     14                transfer_flag            = 0;
     15                leakage_factor           = 0;
     16
    817                spcsediment_head         = NaN;
    918                sediment_compressibility = 0;
    … …  
    1120                sediment_thickness       = 0;
    1221                sediment_transmitivity   = 0;
    13                 water_compressibility    = 0;
    14                 isefficientlayer         = 0;
    15                 penalty_factor           = 0;
    16                 sedimentlimit_flag       = 0;
    17                 sedimentlimit            = 0;
    18                 transfer_flag            = 0;
    19                 leakage_factor           = 0;
    2022               
    21                 spcepl_head         = NaN;
    22                 epl_compressibility = 0;
    23                 epl_porosity        = 0;
    24                 epl_thickness       = 0;
    25                 epl_transmitivity   = 0;
     23                spcepl_head              = NaN;
     24                epl_compressibility      = 0;
     25                epl_porosity             = 0;
     26                epl_thickness            = 0;
     27                epl_transmitivity        = 0;
    2628               
    2729  end
    … …  
    3840                       
    3941                %Parameters from de Fleurian 2013
    40                         obj.sediment_compressibility = 1.0e-08;
    41                         obj.sediment_porosity        = .4;
    42                         obj.sediment_thickness       = 20.0;
    43                         obj.sediment_transmitivity   = 8.0e-04;
    4442                        obj.water_compressibility    = 5.04e-10;
    4543                        obj.isefficientlayer         = 1;
    4644                        obj.penalty_factor           = 3;
     45                        obj.rel_tol                  = 1.0e-06;
    4746                        obj.sedimentlimit_flag       = 0;
    4847                        obj.sedimentlimit            = 0;
    … …  
    5049                        obj.leakage_factor           = 10.0;
    5150
    52                         obj.epl_compressibility = 1.0e-08;
    53                         obj.epl_porosity        = 0.4;
    54                         obj.epl_thickness       = 1.0;
    55                         obj.epl_transmitivity   = 8.0e-02;
     51                        obj.sediment_compressibility = 1.0e-08;
     52                        obj.sediment_porosity        = .4;
     53                        obj.sediment_thickness       = 20.0;
     54                        obj.sediment_transmitivity   = 8.0e-04;
     55
     56                        obj.epl_compressibility      = 1.0e-08;
     57                        obj.epl_porosity             = 0.4;
     58                        obj.epl_thickness            = 1.0;
     59                        obj.epl_transmitivity        = 8.0e-02;
    5660                       
    5761                end % }}}
    … …  
    6367                        end
    6468                       
     69                        md = checkfield(md,'hydrology.water_compressibility','>',0,'numel',1);
     70                        md = checkfield(md,'hydrology.isefficientlayer','numel',[1],'values',[0 1]);
     71                        md = checkfield(md,'hydrology.penalty_factor','>',0,'numel',1);
     72                        md = checkfield(md,'hydrology.rel_tol','>',0,'numel',1);
     73                        md = checkfield(md,'hydrology.sedimentlimit_flag','numel',[1],'values',[0 1 2 3]);
     74                        md = checkfield(md,'hydrology.transfer_flag','numel',[1],'values',[0 1]);
     75                        if obj.sedimentlimit_flag==1,
     76                                md = checkfield(md,'hydrology.sedimentlimit','>',0,'numel',1);
     77            end
     78                        if obj.transfer_flag==1,
     79                                md = checkfield(md,'hydrology.leakage_factor','>',0,'numel',1);
     80            end
     81                       
    6582                        md = checkfield(md,'hydrology.spcsediment_head','forcing',1);
    6683                        md = checkfield(md,'hydrology.sediment_compressibility','>',0,'numel',1);
    … …  
    6885                        md = checkfield(md,'hydrology.sediment_thickness','>',0,'numel',1);
    6986                        md = checkfield(md,'hydrology.sediment_transmitivity','>',0,'numel',1);
    70                         md = checkfield(md,'hydrology.water_compressibility','>',0,'numel',1);
    71                         md = checkfield(md,'hydrology.isefficientlayer','numel',[1],'values',[0 1]);
    72                         md = checkfield(md,'hydrology.penalty_factor','>',0,'numel',1);
    73                         md = checkfield(md,'hydrology.sedimentlimit_flag','numel',[1],'values',[0 1 2 3]);
    74                         md = checkfield(md,'hydrology.transfer_flag','numel',[1],'values',[0 1]);
    75                        
    76                         if obj.sedimentlimit_flag==1,
    77                                 md = checkfield(md,'hydrology.sedimentlimit','>',0,'numel',1);
    78             end
    79                        
    80                         if obj.transfer_flag==1,
    81                                 md = checkfield(md,'hydrology.leakage_factor','>',0,'numel',1);
    82             end
    8387                       
    8488                        if obj.isefficientlayer==1,
    85                                
    8689                                md = checkfield(md,'hydrology.spcepl_head','forcing',1);
    8790                                md = checkfield(md,'hydrology.epl_compressibility','>',0,'numel',1);
    … …  
    9396                function disp(obj) % {{{
    9497                        disp(sprintf('   hydrology Dual Porous Continuum Equivalent parameters:'));
    95                         disp(sprintf('   - for the sediment layer'));
    96                        
    97                         fielddisplay(obj,'spcsediment_head','sediment water head constraints (NaN means no constraint) [m above MSL]');
    98                         fielddisplay(obj,'sediment_compressibility','sediment compressibility [Pa^-1]');
    99                         fielddisplay(obj,'sediment_porosity','sediment [dimensionless]');
    100                         fielddisplay(obj,'sediment_thickness','sediment thickness [m]');
    101                         fielddisplay(obj,'sediment_transmitivity','sediment transmitivity [m^2/s]');
     98                        disp(sprintf('   - general parameters'));
    10299                        fielddisplay(obj,'water_compressibility','compressibility of water [Pa^-1]');
    103100                        fielddisplay(obj,'isefficientlayer','do we use an efficient drainage system [1: true; 0: false]');
    104                         fielddisplay(obj,'penalty_factor','exponent of the value used in the penalisation method');
     101                        fielddisplay(obj,'penalty_factor','exponent of the value used in the penalisation method [dimensionless]');
     102                        fielddisplay(obj,'rel_tol','tolerance of the nonlinear iteration for the transfer between layers [dimensionless]');
    105103                        fielddisplay(obj,'sedimentlimit_flag','what kind of upper limit is applied for the inefficient layer');
    106104                        disp(sprintf('%55s  0: no limit',' '));
    … …  
    117115                                fielddisplay(obj,'leakage_factor','user defined leakage factor [m]');
    118116            end
    119 
     117                        disp(sprintf('   - for the sediment layer'));
     118                        fielddisplay(obj,'spcsediment_head','sediment water head constraints (NaN means no constraint) [m above MSL]');
     119                        fielddisplay(obj,'sediment_compressibility','sediment compressibility [Pa^-1]');
     120                        fielddisplay(obj,'sediment_porosity','sediment [dimensionless]');
     121                        fielddisplay(obj,'sediment_thickness','sediment thickness [m]');
     122                        fielddisplay(obj,'sediment_transmitivity','sediment transmitivity [m^2/s]');
     123                       
    120124                        if obj.isefficientlayer==1,
    121125                                disp(sprintf('   - for the epl layer'));
    122                                
    123126                                fielddisplay(obj,'spcepl_head','epl water head constraints (NaN means no constraint) [m above MSL]');
    124127                                fielddisplay(obj,'epl_compressibility','epl compressibility [Pa^-1]');
    … …  
    130133                end % }}}
    131134                function marshall(obj,fid) % {{{
    132                         WriteData(fid,'enum',HydrologyModelEnum(),'data',HydrologydcEnum(),'format','Integer');
     135                        WriteData(fid,'enum',HydrologyModelEnum(),'data',HydrologydcEnum(),'format','Integer');
     136                        WriteData(fid,'object',obj,'fieldname','water_compressibility','format','Double');
     137                        WriteData(fid,'object',obj,'fieldname','isefficientlayer','format','Boolean');
     138                        WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double');
     139                        WriteData(fid,'object',obj,'fieldname','rel_tol','format','Double');
     140                        WriteData(fid,'object',obj,'fieldname','sedimentlimit_flag','format','Integer');
     141                        WriteData(fid,'object',obj,'fieldname','transfer_flag','format','Integer');
     142                        if obj.sedimentlimit_flag==1,
     143                                WriteData(fid,'object',obj,'fieldname','sedimentlimit','format','Double');
     144            end
     145                        if obj.transfer_flag==1,
     146                                WriteData(fid,'object',obj,'fieldname','leakage_factor','format','Double');
     147            end
     148
    133149                        WriteData(fid,'object',obj,'fieldname','spcsediment_head','format','DoubleMat','mattype',1);
    134150                        WriteData(fid,'object',obj,'fieldname','sediment_compressibility','format','Double');                   
    135151                        WriteData(fid,'object',obj,'fieldname','sediment_porosity','format','Double');                 
    136152                        WriteData(fid,'object',obj,'fieldname','sediment_thickness','format','Double');
    137                         WriteData(fid,'object',obj,'fieldname','sediment_transmitivity','format','Double');                     
    138                         WriteData(fid,'object',obj,'fieldname','water_compressibility','format','Double');
    139                         WriteData(fid,'object',obj,'fieldname','isefficientlayer','format','Boolean');
    140                         WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double');
    141                         WriteData(fid,'object',obj,'fieldname','sedimentlimit_flag','format','Integer');
    142                         WriteData(fid,'object',obj,'fieldname','transfer_flag','format','Integer');
    143 
    144                         if obj.sedimentlimit_flag==1,
    145                                 WriteData(fid,'object',obj,'fieldname','sedimentlimit','format','Double');
    146             end
    147 
    148                         if obj.transfer_flag==1,
    149                                 WriteData(fid,'object',obj,'fieldname','leakage_factor','format','Double');
    150             end
    151 
    152                         if obj.isefficientlayer==1,
    153                                
     153                        WriteData(fid,'object',obj,'fieldname','sediment_transmitivity','format','Double');             
     154                       
     155                        if obj.isefficientlayer==1,     
    154156                                WriteData(fid,'object',obj,'fieldname','spcepl_head','format','DoubleMat','mattype',1);
    155157                                WriteData(fid,'object',obj,'fieldname','epl_compressibility','format','Double');                       

  • issm/trunk-jpl/src/m/contrib/paraview/writeVTKcell.m

    r15011 r15062  
    11function writeVTKcell(filename,model,Solution)
    22% vtk export
     3% function writeVTKcell(filename,model,Solution)
    34% creates a vtk-file filename.vtk containing simplicial mesh data
    45% (only work for triangle now)

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

    r15049 r15062  
    11431143        return StringToEnum('EplHead')[0]
    11441144
     1145def HydrologydcRelTolEnum():
     1146        """
     1147        HYDROLOGYDCRELTOLENUM - Enum of HydrologydcRelTol
     1148
     1149        WARNING: DO NOT MODIFY THIS FILE
     1150                                this file has been automatically generated by src/c/shared/Enum/Synchronize.sh
     1151                                Please read src/c/shared/Enum/README for more information
     1152
     1153           Usage:
     1154              macro=HydrologydcRelTolEnum()
     1155        """
     1156
     1157        return StringToEnum('HydrologydcRelTol')[0]
     1158
    11451159def HydrologydcSpcsedimentHeadEnum():
    11461160        """
    … …  
    77777791        """
    77787792
    7779         return 554
    7780 
     7793        return 555
     7794

  • issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m

    r15049 r15062  
    99%      macro=MaximumNumberOfEnums()
    1010
    11 macro=554;
     11macro=555;
Note: See TracChangeset for help on using the changeset viewer.