Index: /issm/trunk-jpl/src/c/analyses/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/hydrology_core.cpp (revision 15190) +++ /issm/trunk-jpl/src/c/analyses/hydrology_core.cpp (revision 15191) @@ -93,4 +93,5 @@ if(isefficientlayer){ InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,EplHeadEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,HydrologydcMaskEplactiveEnum); } /*unload results*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15190) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15191) @@ -6395,4 +6395,6 @@ int *doflist = NULL; bool converged; + bool isefficientlayer; + IssmDouble activeEpl[numdof],activate[numdof]={0.0}; IssmDouble values[numdof]; IssmDouble residual[numdof]; @@ -6403,4 +6405,7 @@ GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + /*Get the flag to know if the efficient layer is present*/ + this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); + /*Use the dof list to index into the solution vector: */ for(int i=0;iparameters->FindParam(&kmax,HydrologySedimentKmaxEnum); this->parameters->FindParam(&penalty_factor,HydrologydcPenaltyFactorEnum); + + if(isefficientlayer)GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum); + kappa=kmax*pow(10.,penalty_factor); @@ -6422,12 +6430,16 @@ if(values[i]>h_max){ residual[i]=kappa*(values[i]-h_max); - if(reCast(dt))residual[i]=residual[i]; + if(isefficientlayer) activate[i]=1.0; } - else + else{ + if(isefficientlayer){ + if(activeEpl[i]==1.0)activate[i]=1.0; + } residual[i]=0.0; + } } } /*Add input to the element: */ - + if(isefficientlayer) this->inputs->AddInput(new TriaP1Input(HydrologydcMaskEplactiveEnum,activate)); this->inputs->AddInput(new TriaP1Input(SedimentHeadEnum,values)); this->inputs->AddInput(new TriaP1Input(SedimentHeadResidualEnum,residual)); @@ -6509,4 +6521,5 @@ IssmDouble sed_trans,sed_thick; IssmDouble leakage; + IssmDouble activeEpl[numdof]; IssmDouble wh_trans[numdof]={0.0}; IssmDouble storing[numdof]; @@ -6516,42 +6529,44 @@ /*Get the flag to know if the efficient layer is present*/ - this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); - - if(!isefficientlayer){ - transfer->SetValues(numdof,doflist,&wh_trans[0],ADD_VAL); - return; - } - /*Also get the flag to the transfer method (TO DO)*/ - this->parameters->FindParam(&transfermethod,HydrologydcTransferFlagEnum); - - /*Switch between the different transfer methods cases(TO DO)*/ - switch(transfermethod){ - case 0: - /*Just kipping the transfer to zero*/ - break; - case 1: - this->parameters->FindParam(&leakage,HydrologydcLeakageFactorEnum); - - sed_trans = matpar->GetSedimentTransmitivity(); - sed_thick = matpar->GetSedimentThickness(); - - GetInputListOnVertices(&sed_head[0],SedimentHeadEnum); - GetInputListOnVertices(&epl_head[0],EplHeadEnum); - - for(int i=0;iepl_head[i]){ - storing[i]=matpar->GetSedimentStoring(); - wh_trans[i]=sed_trans*storing[i]*(sed_head[i]-epl_head[i])/(leakage*sed_thick); + this->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); + + if(isefficientlayer){ + /*Also get the flag to the transfer method*/ + this->parameters->FindParam(&transfermethod,HydrologydcTransferFlagEnum); + + /*Switch between the different transfer methods cases*/ + switch(transfermethod){ + case 0: + /*Just kipping the transfer to zero*/ + break; + case 1: + + GetInputListOnVertices(&activeEpl[0],HydrologydcMaskEplactiveEnum); + GetInputListOnVertices(&sed_head[0],SedimentHeadEnum); + GetInputListOnVertices(&epl_head[0],EplHeadEnum); + + this->parameters->FindParam(&leakage,HydrologydcLeakageFactorEnum); + + sed_trans = matpar->GetSedimentTransmitivity(); + sed_thick = matpar->GetSedimentThickness(); + + for(int i=0;iepl_head[i]){ + storing[i]=matpar->GetSedimentStoring(); + wh_trans[i]=sed_trans*storing[i]*(sed_head[i]-epl_head[i])/(leakage*sed_thick); + } + else{ + storing[i]=matpar->GetEplStoring(); + wh_trans[i]=sed_trans*storing[i]*(sed_head[i]-epl_head[i])/(leakage*sed_thick); + } } - else{ - storing[i]=matpar->GetEplStoring(); - wh_trans[i]=sed_trans*storing[i]*(sed_head[i]-epl_head[i])/(leakage*sed_thick); - } + break; + default: + _error_("no case higher than 1 for the Transfer method"); } - break; - default: - _error_("no case higher than 1 for the Transfer method"); - } - + } /*Assign output pointer*/ transfer->SetValues(numdof,doflist,&wh_trans[0],ADD_VAL); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 15190) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 15191) @@ -92,5 +92,5 @@ bool IsNodeOnShelf(); bool IsNodeOnShelfFromFlags(IssmDouble* flags); - bool IsOnWater(); + bool IsOnWater(); void GetSolutionFromInputs(Vector* solution); void GetVectorFromInputs(Vector* vector, int name_enum); Index: /issm/trunk-jpl/src/m/classes/hydrologydc.m =================================================================== --- /issm/trunk-jpl/src/m/classes/hydrologydc.m (revision 15190) +++ /issm/trunk-jpl/src/m/classes/hydrologydc.m (revision 15191) @@ -22,4 +22,5 @@ spcepl_head = NaN; + mask_eplactive = NaN; epl_compressibility = 0; epl_porosity = 0; @@ -50,5 +51,5 @@ obj.sediment_compressibility = 1.0e-08; - obj.sediment_porosity = .4; + obj.sediment_porosity = 0.4; obj.sediment_thickness = 20.0; obj.sediment_transmitivity = 8.0e-04; @@ -88,4 +89,5 @@ if obj.isefficientlayer==1, md = checkfield(md,'hydrology.spcepl_head','forcing',1); + md = checkfield(md,'hydrology.mask_eplactive','size',[md.mesh.numberofvertices 1],'values',[0 1]); md = checkfield(md,'hydrology.epl_compressibility','>',0,'numel',1); md = checkfield(md,'hydrology.epl_porosity','>',0,'numel',1); @@ -125,4 +127,5 @@ disp(sprintf(' - for the epl layer')); fielddisplay(obj,'spcepl_head','epl water head constraints (NaN means no constraint) [m above MSL]'); + fielddisplay(obj,'mask_eplactive','active (1) or not (0) EPL'); fielddisplay(obj,'epl_compressibility','epl compressibility [Pa^-1]'); fielddisplay(obj,'epl_porosity','epl [dimensionless]'); @@ -133,5 +136,5 @@ end % }}} function marshall(obj,md,fid) % {{{ - WriteData(fid,'enum',HydrologyModelEnum(),'data',HydrologydcEnum(),'format','Integer'); + WriteData(fid,'enum',HydrologyModelEnum(),'data',HydrologydcEnum(),'format','Integer'); WriteData(fid,'object',obj,'fieldname','water_compressibility','format','Double'); WriteData(fid,'object',obj,'fieldname','isefficientlayer','format','Boolean'); @@ -154,5 +157,6 @@ if obj.isefficientlayer==1, - WriteData(fid,'object',obj,'fieldname','spcepl_head','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1); + WriteData(fid,'object',obj,'fieldname','spcepl_head','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1); + WriteData(fid,'object',obj,'fieldname','mask_eplactive','format','DoubleMat','mattype',1); WriteData(fid,'object',obj,'fieldname','epl_compressibility','format','Double'); WriteData(fid,'object',obj,'fieldname','epl_porosity','format','Double'); Index: /issm/trunk-jpl/test/NightlyRun/test332.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test332.m (revision 15190) +++ /issm/trunk-jpl/test/NightlyRun/test332.m (revision 15191) @@ -3,5 +3,5 @@ md=parameterize(md,'../Par/SquareSheetConstrained.par'); md=setflowequation(md,'macayeal','all'); -md.cluster=generic('name',oshostname(),'np',2); +md.cluster=generic('name',oshostname(),'np',1); md.hydrology=(hydrologydc); md.hydrology.isefficientlayer=0; @@ -16,5 +16,5 @@ md.timestepping.time_step=0; md.timestepping.final_time=1.0; -%md.verbose.convergence=1; +%md.verbose=verbose('1111111'); md=solve(md,HydrologySolutionEnum()); Index: /issm/trunk-jpl/test/NightlyRun/test333.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test333.m (revision 15190) +++ /issm/trunk-jpl/test/NightlyRun/test333.m (revision 15191) @@ -15,4 +15,5 @@ md.initialization.epl_head=0.0*ones(md.mesh.numberofvertices,1); md.hydrology.spcepl_head=NaN*ones(md.mesh.numberofvertices,1); +md.hydrology.mask_eplactive=0*ones(md.mesh.numberofvertices,1); md.basalforcings.melting_rate = 2.0*ones(md.mesh.numberofvertices,1); md.hydrology.epl_transmitivity=30; @@ -21,5 +22,21 @@ md.timestepping.final_time=2.0; +%md.verbose.solution=1; + md=solve(md,HydrologySolutionEnum()); + +%store=md.constants.g*md.hydrology.sediment_porosity* ... +% md.materials.rho_freshwater*((md.hydrology.sediment_compressibility/md.hydrology.sediment_porosity)+md.hydrology.water_compressibility) + +%sed=ones(1,size(md.results.HydrologySolution,2)); +%epl=ones(1,size(md.results.HydrologySolution,2)); +%res=ones(1,size(md.results.HydrologySolution,2)); +%input=ones(1,size(md.results.HydrologySolution,2)); +%for i= 1:size(md.results.HydrologySolution,2) +% sed(i)=mean(md.results.HydrologySolution(i).SedimentHead); +% res(i)=mean(md.results.HydrologySolution(i).SedimentHeadResidual); +% epl(i)=mean(md.results.HydrologySolution(i).EplHead); +% input(i)=2.0*(i*0.2); +%end %Fields and tolerances to track changes