Changeset 15191


Ignore:
Timestamp:
06/05/13 13:40:31 (12 years ago)
Author:
bdef
Message:

NEW:Adding a mask for the active EPL nodes

Location:
issm/trunk-jpl
Files:
6 edited

Legend:

Unmodified
Added
Removed

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

    r15161 r15191  
    9393                                if(isefficientlayer){
    9494                                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,EplHeadEnum);
     95                                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,HydrologydcMaskEplactiveEnum);
    9596                                }
    9697                                /*unload results*/

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

    r15182 r15191  
    63956395        int        *doflist        = NULL;
    63966396        bool        converged;
     6397        bool        isefficientlayer;
     6398        IssmDouble  activeEpl[numdof],activate[numdof]={0.0};
    63976399        IssmDouble  values[numdof];
    63986400        IssmDouble  residual[numdof];
    … …  
    64036405        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    64046406
     6407        /*Get the flag to know if the efficient layer is present*/
     6408                this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
     6409
    64056410        /*Use the dof list to index into the solution vector: */
    64066411        for(int i=0;i<numdof;i++){
    … …  
    64166421                this->parameters->FindParam(&kmax,HydrologySedimentKmaxEnum);
    64176422                this->parameters->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum);
     6423               
     6424                if(isefficientlayer)GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum);
     6425               
    64186426                kappa=kmax*pow(10.,penalty_factor);
    64196427               
    … …  
    64226430                        if(values[i]>h_max){
    64236431                                residual[i]=kappa*(values[i]-h_max);
    6424                                 if(reCast<bool,IssmDouble>(dt))residual[i]=residual[i];
     6432                                if(isefficientlayer) activate[i]=1.0;
    64256433                        }
    6426                         else
     6434                        else{
     6435                                if(isefficientlayer){
     6436                                        if(activeEpl[i]==1.0)activate[i]=1.0;
     6437                                }                               
    64276438                                residual[i]=0.0;
     6439                        }
    64286440                }
    64296441        }
    64306442        /*Add input to the element: */
    6431 
     6443        if(isefficientlayer) this->inputs->AddInput(new TriaP1Input(HydrologydcMaskEplactiveEnum,activate));
    64326444        this->inputs->AddInput(new TriaP1Input(SedimentHeadEnum,values));
    64336445        this->inputs->AddInput(new TriaP1Input(SedimentHeadResidualEnum,residual));
    … …  
    65096521        IssmDouble sed_trans,sed_thick;
    65106522        IssmDouble leakage;
     6523        IssmDouble activeEpl[numdof];
    65116524        IssmDouble wh_trans[numdof]={0.0};
    65126525        IssmDouble storing[numdof];
    … …  
    65166529
    65176530        /*Get the flag to know if the efficient layer is present*/
    6518         this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
    6519 
    6520         if(!isefficientlayer){
    6521                 transfer->SetValues(numdof,doflist,&wh_trans[0],ADD_VAL);
    6522                 return;
    6523         }       
    6524         /*Also get the flag to the transfer method (TO DO)*/
    6525         this->parameters->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
    6526 
    6527         /*Switch between the different transfer methods cases(TO DO)*/
    6528         switch(transfermethod){
    6529         case 0:
    6530                 /*Just kipping the transfer to zero*/
    6531                 break;
    6532         case 1:
    6533                 this->parameters->FindParam(&leakage,HydrologydcLeakageFactorEnum);
    6534 
    6535                 sed_trans = matpar->GetSedimentTransmitivity();
    6536                 sed_thick = matpar->GetSedimentThickness();
    6537 
    6538                 GetInputListOnVertices(&sed_head[0],SedimentHeadEnum);
    6539                 GetInputListOnVertices(&epl_head[0],EplHeadEnum);
    6540 
    6541                 for(int i=0;i<numdof;i++){
    6542                         if(sed_head[i]>epl_head[i]){
    6543                                 storing[i]=matpar->GetSedimentStoring();
    6544                                 wh_trans[i]=sed_trans*storing[i]*(sed_head[i]-epl_head[i])/(leakage*sed_thick);
     6531                this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
     6532
     6533        if(isefficientlayer){
     6534                /*Also get the flag to the transfer method*/
     6535                this->parameters->FindParam(&transfermethod,HydrologydcTransferFlagEnum);
     6536               
     6537                /*Switch between the different transfer methods cases*/
     6538                switch(transfermethod){
     6539                case 0:
     6540                        /*Just kipping the transfer to zero*/
     6541                        break;
     6542                case 1:
     6543                       
     6544                        GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum);
     6545                        GetInputListOnVertices(&sed_head[0],SedimentHeadEnum);
     6546                        GetInputListOnVertices(&epl_head[0],EplHeadEnum);
     6547
     6548                        this->parameters->FindParam(&leakage,HydrologydcLeakageFactorEnum);
     6549                       
     6550                        sed_trans = matpar->GetSedimentTransmitivity();
     6551                        sed_thick = matpar->GetSedimentThickness();
     6552                       
     6553                        for(int i=0;i<numdof;i++){
     6554                       
     6555                                if(activeEpl[i]==0.0)break;
     6556                               
     6557                                if(sed_head[i]>epl_head[i]){
     6558                                        storing[i]=matpar->GetSedimentStoring();
     6559                                        wh_trans[i]=sed_trans*storing[i]*(sed_head[i]-epl_head[i])/(leakage*sed_thick);
     6560                                }
     6561                                else{
     6562                                        storing[i]=matpar->GetEplStoring();
     6563                                        wh_trans[i]=sed_trans*storing[i]*(sed_head[i]-epl_head[i])/(leakage*sed_thick);
     6564                                }
    65456565                        }
    6546                         else{
    6547                                 storing[i]=matpar->GetEplStoring();
    6548                                 wh_trans[i]=sed_trans*storing[i]*(sed_head[i]-epl_head[i])/(leakage*sed_thick);
    6549                         }
     6566                        break;
     6567                default:
     6568                        _error_("no case higher than 1 for the Transfer method");
    65506569                }
    6551                 break;
    6552         default:
    6553                 _error_("no case higher than 1 for the Transfer method");
    6554         }
    6555 
     6570        }
    65566571        /*Assign output pointer*/
    65576572        transfer->SetValues(numdof,doflist,&wh_trans[0],ADD_VAL);

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

    r15130 r15191  
    9292                bool   IsNodeOnShelf();
    9393                bool   IsNodeOnShelfFromFlags(IssmDouble* flags);
    94                 bool   IsOnWater(); 
     94                bool   IsOnWater();
    9595                void   GetSolutionFromInputs(Vector<IssmDouble>* solution);
    9696                void   GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum);

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

    r15133 r15191  
    2222
    2323                spcepl_head              = NaN;
     24                mask_eplactive           = NaN;
    2425                epl_compressibility      = 0;
    2526                epl_porosity             = 0;
    … …  
    5051
    5152                        obj.sediment_compressibility = 1.0e-08;
    52                         obj.sediment_porosity        = .4;
     53                        obj.sediment_porosity        = 0.4;
    5354                        obj.sediment_thickness       = 20.0;
    5455                        obj.sediment_transmitivity   = 8.0e-04;
    … …  
    8889                        if obj.isefficientlayer==1,
    8990                                md = checkfield(md,'hydrology.spcepl_head','forcing',1);
     91                                md = checkfield(md,'hydrology.mask_eplactive','size',[md.mesh.numberofvertices 1],'values',[0 1]);
    9092                                md = checkfield(md,'hydrology.epl_compressibility','>',0,'numel',1);
    9193                                md = checkfield(md,'hydrology.epl_porosity','>',0,'numel',1);
    … …  
    125127                                disp(sprintf('   - for the epl layer'));
    126128                                fielddisplay(obj,'spcepl_head','epl water head constraints (NaN means no constraint) [m above MSL]');
     129                                fielddisplay(obj,'mask_eplactive','active (1) or not (0) EPL');
    127130                                fielddisplay(obj,'epl_compressibility','epl compressibility [Pa^-1]');
    128131                                fielddisplay(obj,'epl_porosity','epl [dimensionless]');
    … …  
    133136                end % }}}
    134137                function marshall(obj,md,fid) % {{{
    135                         WriteData(fid,'enum',HydrologyModelEnum(),'data',HydrologydcEnum(),'format','Integer'); 
     138                        WriteData(fid,'enum',HydrologyModelEnum(),'data',HydrologydcEnum(),'format','Integer');
    136139                        WriteData(fid,'object',obj,'fieldname','water_compressibility','format','Double');
    137140                        WriteData(fid,'object',obj,'fieldname','isefficientlayer','format','Boolean');
    … …  
    154157
    155158                        if obj.isefficientlayer==1,     
    156                                 WriteData(fid,'object',obj,'fieldname','spcepl_head','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);
     159                                WriteData(fid,'object',obj,'fieldname','spcepl_head','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);     
     160                                WriteData(fid,'object',obj,'fieldname','mask_eplactive','format','DoubleMat','mattype',1);
    157161                                WriteData(fid,'object',obj,'fieldname','epl_compressibility','format','Double');                       
    158162                                WriteData(fid,'object',obj,'fieldname','epl_porosity','format','Double');                       

  • issm/trunk-jpl/test/NightlyRun/test332.m

    r15078 r15191  
    33md=parameterize(md,'../Par/SquareSheetConstrained.par');
    44md=setflowequation(md,'macayeal','all');
    5 md.cluster=generic('name',oshostname(),'np',2);
     5md.cluster=generic('name',oshostname(),'np',1);
    66md.hydrology=(hydrologydc);
    77md.hydrology.isefficientlayer=0;
    … …  
    1616md.timestepping.time_step=0;
    1717md.timestepping.final_time=1.0;
    18 %md.verbose.convergence=1;
     18%md.verbose=verbose('1111111');
    1919md=solve(md,HydrologySolutionEnum());
    2020

  • issm/trunk-jpl/test/NightlyRun/test333.m

    r15162 r15191  
    1515md.initialization.epl_head=0.0*ones(md.mesh.numberofvertices,1);
    1616md.hydrology.spcepl_head=NaN*ones(md.mesh.numberofvertices,1);
     17md.hydrology.mask_eplactive=0*ones(md.mesh.numberofvertices,1);
    1718md.basalforcings.melting_rate = 2.0*ones(md.mesh.numberofvertices,1);
    1819md.hydrology.epl_transmitivity=30;
    … …  
    2122md.timestepping.final_time=2.0;
    2223
     24%md.verbose.solution=1;
     25
    2326md=solve(md,HydrologySolutionEnum());
     27
     28%store=md.constants.g*md.hydrology.sediment_porosity* ...
     29%                       md.materials.rho_freshwater*((md.hydrology.sediment_compressibility/md.hydrology.sediment_porosity)+md.hydrology.water_compressibility)
     30
     31%sed=ones(1,size(md.results.HydrologySolution,2));
     32%epl=ones(1,size(md.results.HydrologySolution,2));
     33%res=ones(1,size(md.results.HydrologySolution,2));
     34%input=ones(1,size(md.results.HydrologySolution,2));
     35%for i= 1:size(md.results.HydrologySolution,2)
     36%       sed(i)=mean(md.results.HydrologySolution(i).SedimentHead);
     37%       res(i)=mean(md.results.HydrologySolution(i).SedimentHeadResidual);
     38%       epl(i)=mean(md.results.HydrologySolution(i).EplHead);
     39%       input(i)=2.0*(i*0.2);
     40%end
    2441
    2542%Fields and tolerances to track changes
Note: See TracChangeset for help on using the changeset viewer.