Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 26486) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 26487) @@ -93,5 +93,5 @@ break; case FrontalForcingsRignotEnum: - iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.basin",FrontalForcingsBasinIdEnum); + iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.basin_id",FrontalForcingsBasinIdEnum); iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.subglacial_discharge",FrontalForcingsSubglacialDischargeEnum); iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.thermalforcing",FrontalForcingsThermalForcingEnum); @@ -141,5 +141,5 @@ break; case FrontalForcingsRignotEnum: - parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.numberofbasins",FrontalForcingsNumberofBasinsEnum)); + parameters->AddObject(iomodel->CopyConstantObject("md.frontalforcings.num_basins",FrontalForcingsNumberofBasinsEnum)); break; default: Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26486) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 26487) @@ -3378,4 +3378,59 @@ } /*}}}*/ +void Element::RignotMeltParameterization(){/*{{{*/ + + const int numvertices = this->GetNumberOfVertices(); + IssmDouble A, B, alpha, beta; + IssmDouble bed,qsg,qsg_basin,TF,yts; + int numbasins,basinid; + IssmDouble* basin_icefront_area=NULL; + + /* Coefficients */ + A = 3e-4; + B = 0.15; + alpha = 0.39; + beta = 1.18; + + /*Get inputs*/ + this->GetInputValue(&basinid,FrontalForcingsBasinIdEnum); + Input* bed_input = this->GetInput(BedEnum); _assert_(bed_input); + Input* qsg_input = this->GetInput(FrontalForcingsSubglacialDischargeEnum); _assert_(qsg_input); + Input* TF_input = this->GetInput(FrontalForcingsThermalForcingEnum); _assert_(TF_input); + + this->FindParam(&yts, ConstantsYtsEnum); + this->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum); + this->parameters->FindParam(&basin_icefront_area,&numbasins,FrontalForcingsBasinIcefrontAreaEnum); + IssmDouble meltrates[numvertices]; + + /* Start looping on the number of vertices: */ + GaussTria gauss; + for(int iv=0;ivGetInputValue(&bed,&gauss); + qsg_input->GetInputValue(&qsg,&gauss); + TF_input->GetInputValue(&TF,&gauss); + + if(basin_icefront_area[basinid]==0.) meltrates[iv]=0.; + else{ + /* change the unit of qsg (m^3/d -> m/d) with ice front area */ + qsg_basin=qsg/basin_icefront_area[basinid]; + + /* calculate melt rates */ + meltrates[iv]=((A*max(-bed,0.)*pow(max(qsg_basin,0.),alpha)+B)*pow(max(TF,0.),beta))/86400; //[m/s] + } + + if(xIsNan(meltrates[iv])) _error_("NaN found in vector"); + if(xIsInf(meltrates[iv])) _error_("Inf found in vector"); + } + + /*Add input*/ + this->AddInput(CalvingMeltingrateEnum,&meltrates[0],P1Enum); + + /*Cleanup and return*/ + xDelete(basin_icefront_area); +} +/*}}}*/ void Element::SetBoolInput(Inputs* inputs,int enum_in,bool value){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26486) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26487) @@ -167,4 +167,5 @@ void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac); void PositiveDegreeDaySicopolis(bool isfirnwarming); + void RignotMeltParameterization(); void ResultInterpolation(int* pinterpolation,int*nodesperelement,int* parray_size, int output_enum); void ResultToPatch(IssmDouble* values,int nodesperelement,int output_enum); @@ -345,5 +346,4 @@ virtual void ResetFSBasalBoundaryCondition()=0; virtual void ResetHooks()=0; - virtual void RignotMeltParameterization(void){_error_("not implemented yet");}; virtual void SetElementInput(int enum_in,IssmDouble value){_error_("not implemented yet");}; virtual void SetElementInput(int enum_in,IssmDouble value, int type)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 26486) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 26487) @@ -3402,64 +3402,4 @@ if(this->hneighbors) this->hneighbors->reset(); -} -/*}}}*/ -void Penta::RignotMeltParameterization(){/*{{{*/ - - if(!this->IsOnBase()) return; - - IssmDouble A, B, alpha, beta; - IssmDouble bed,qsg,qsg_basin,TF,yts; - int numbasins; - IssmDouble basinid[NUMVERTICES]; - IssmDouble* basin_icefront_area=NULL; - - /* Coefficients */ - A = 3e-4; - B = 0.15; - alpha = 0.39; - beta = 1.18; - - /*Get inputs*/ - Input* bed_input = this->GetInput(BedEnum); _assert_(bed_input); - Input* qsg_input = this->GetInput(FrontalForcingsSubglacialDischargeEnum); _assert_(qsg_input); - Input* TF_input = this->GetInput(FrontalForcingsThermalForcingEnum); _assert_(TF_input); - Element::GetInputListOnVertices(&basinid[0],FrontalForcingsBasinIdEnum); - - this->FindParam(&yts, ConstantsYtsEnum); - this->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum); - this->parameters->FindParam(&basin_icefront_area,&numbasins,FrontalForcingsBasinIcefrontAreaEnum); - - IssmDouble meltrates[NUMVERTICES2D]; //frontal melt-rate - - /* Start looping on the number of vertices: */ - GaussPenta gauss; - for(int iv=0;ivGetInputValue(&bed,&gauss); - qsg_input->GetInputValue(&qsg,&gauss); - TF_input->GetInputValue(&TF,&gauss); - - if(basinid[iv]==0 || basin_icefront_area[reCast(basinid[iv])-1]==0.) meltrates[iv]=0.; - else{ - /* change the unit of qsg (m^3/d -> m/d) with ice front area */ - qsg_basin=qsg/basin_icefront_area[reCast(basinid[iv])-1]; - - /* calculate melt rates */ - meltrates[iv]=((A*max(-bed,0.)*pow(max(qsg_basin,0.),alpha)+B)*pow(max(TF,0.),beta))/86400; //[m/s] - } - - if(xIsNan(meltrates[iv])) _error_("NaN found in vector"); - if(xIsInf(meltrates[iv])) _error_("Inf found in vector"); - } - - /*Add input*/ - this->AddInput(CalvingMeltingrateEnum,&meltrates[0],P1Enum); - - this->InputExtrude(CalvingMeltingrateEnum,-1); - - /*Cleanup and return*/ - xDelete(basin_icefront_area); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26486) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26487) @@ -171,5 +171,4 @@ void ResetFSBasalBoundaryCondition(void); void ResetHooks(); - void RignotMeltParameterization(); void SetElementInput(int enum_in,IssmDouble value); void SetElementInput(int enum_in,IssmDouble value,int type){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26486) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26487) @@ -4717,60 +4717,4 @@ if(this->hneighbors) this->hneighbors->reset(); -} -/*}}}*/ -void Tria::RignotMeltParameterization(){/*{{{*/ - - IssmDouble A, B, alpha, beta; - IssmDouble bed,qsg,qsg_basin,TF,yts; - int numbasins; - IssmDouble basinid[NUMVERTICES]; - IssmDouble* basin_icefront_area=NULL; - - /* Coefficients */ - A = 3e-4; - B = 0.15; - alpha = 0.39; - beta = 1.18; - - /*Get inputs*/ - Input* bed_input = this->GetInput(BedEnum); _assert_(bed_input); - Input* qsg_input = this->GetInput(FrontalForcingsSubglacialDischargeEnum); _assert_(qsg_input); - Input* TF_input = this->GetInput(FrontalForcingsThermalForcingEnum); _assert_(TF_input); - Element::GetInputListOnVertices(&basinid[0],FrontalForcingsBasinIdEnum); - - this->FindParam(&yts, ConstantsYtsEnum); - this->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum); - this->parameters->FindParam(&basin_icefront_area,&numbasins,FrontalForcingsBasinIcefrontAreaEnum); - - IssmDouble meltrates[NUMVERTICES]; //frontal melt-rate - - /* Start looping on the number of vertices: */ - GaussTria gauss; - for(int iv=0;ivGetInputValue(&bed,&gauss); - qsg_input->GetInputValue(&qsg,&gauss); - TF_input->GetInputValue(&TF,&gauss); - - if(basin_icefront_area[reCast(basinid[iv])-1]==0.) meltrates[iv]=0.; - else{ - /* change the unit of qsg (m^3/d -> m/d) with ice front area */ - qsg_basin=qsg/basin_icefront_area[reCast(basinid[iv])-1]; - - /* calculate melt rates */ - meltrates[iv]=((A*max(-bed,0.)*pow(max(qsg_basin,0.),alpha)+B)*pow(max(TF,0.),beta))/86400; //[m/s] - } - - if(xIsNan(meltrates[iv])) _error_("NaN found in vector"); - if(xIsInf(meltrates[iv])) _error_("Inf found in vector"); - } - - /*Add input*/ - this->AddInput(CalvingMeltingrateEnum,&meltrates[0],P1Enum); - - /*Cleanup and return*/ - xDelete(basin_icefront_area); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 26486) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 26487) @@ -132,5 +132,4 @@ void ResetHooks(); void ResetLevelsetFromSegmentlist(IssmDouble* segments,int numsegments); - void RignotMeltParameterization(); void SetElementInput(int enum_in,IssmDouble value); void SetElementInput(int enum_in,IssmDouble value,int type); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 26486) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 26487) @@ -1602,34 +1602,27 @@ void FemModel::IcefrontAreax(){/*{{{*/ - int numvertices = 6; - int numbasins; - IssmDouble* BasinId = xNew(numvertices); + int numbasins,BasinId; this->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum); - IssmDouble* basin_icefront_area = xNewZeroInit(numbasins); - - for(int basin=1;basin(numbasins); + + for(int basin=0;basinelements->objects){ Element* element= xDynamicCast(object); - element->GetInputListOnVertices(BasinId,FrontalForcingsBasinIdEnum); - for(int j=0;jGetIcefrontArea(); - break; - } - } + element->GetInputValue(&BasinId,FrontalForcingsBasinIdEnum); + if(BasinId==basin) local_icefront_area+=element->GetIcefrontArea(); } ISSM_MPI_Reduce(&local_icefront_area,&total_icefront_area,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm()); ISSM_MPI_Bcast(&total_icefront_area,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - basin_icefront_area[basin-1]=total_icefront_area; - } - + basin_icefront_area[basin]=total_icefront_area; + } + this->parameters->AddObject(new DoubleVecParam(FrontalForcingsBasinIcefrontAreaEnum,basin_icefront_area,numbasins)); xDelete(basin_icefront_area); - xDelete(BasinId); }/*}}}*/ void FemModel::IcefrontMassFluxx(IssmDouble* pM, bool scaled){/*{{{*/ Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26486) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 26487) @@ -243,5 +243,4 @@ bool randomflag{}; femmodel->parameters->FindParam(&randomflag,SmbRandomflagEnum); - _printf_("randomflag: "<',0,'nan',1,'Inf',1); - md = checkfield(md,'fieldname','frontalforcings.numberofbasins','numel',[1]); + md = checkfield(md,'fieldname','frontalforcings.num_basins','numel',1,'NaN',1,'Inf',1,'>',0); + md = checkfield(md,'fieldname','frontalforcings.basin_id','Inf',1,'>=',0,'<=',md.frontalforcings.num_basins,'size',[md.mesh.numberofelements 1]); md = checkfield(md,'fieldname','frontalforcings.subglacial_discharge','>=',0,'nan',1,'Inf',1,'timeseries',1); md = checkfield(md,'fieldname','frontalforcings.thermalforcing','nan',1,'Inf',1,'timeseries',1); @@ -52,6 +52,6 @@ function disp(self) % {{{ disp(sprintf(' Frontalforcings parameters:')); - fielddisplay(self,'basin','basin ID for vertices'); - fielddisplay(self,'numberofbasins', 'number of basins'); + fielddisplay(self,'basin_id','basin ID for elements'); + fielddisplay(self,'num_basins', 'number of basins'); fielddisplay(self,'subglacial_discharge','sum of subglacial discharge for each basin [m/d]'); fielddisplay(self,'thermalforcing','thermal forcing [∘C]'); @@ -59,6 +59,6 @@ function marshall(self,prefix,md,fid) % {{{ WriteData(fid,prefix,'name','md.frontalforcings.parameterization','data',2,'format','Integer'); - WriteData(fid,prefix,'object',self,'fieldname','basin','format','DoubleMat','mattype',1); - WriteData(fid,prefix,'object',self,'fieldname','numberofbasins','format','Integer'); + WriteData(fid,prefix,'object',self,'fieldname','basin_id','data',self.basin_id-1,'name','md.frontalforcings.basin_id','format','IntMat','mattype',2); %0-indexed + WriteData(fid,prefix,'object',self,'fieldname','num_basins','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','subglacial_discharge','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'fieldname','thermalforcing','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); Index: /issm/trunk-jpl/test/NightlyRun/test257.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test257.m (revision 26487) +++ /issm/trunk-jpl/test/NightlyRun/test257.m (revision 26487) @@ -0,0 +1,69 @@ +%Test Name: SquareShelfSMBautoregression +md=triangle(model(),'../Exp/Square.exp',80000.); +md=setmask(md,'all',''); +md=parameterize(md,'../Par/SquareShelfConstrained.par'); +md=setflowequation(md,'SSA','all'); +md.cluster=generic('name',oshostname(),'np',3); +md.transient.requested_outputs={'default','IceVolume','SmbMassBalance',}; + +ymax = max(md.mesh.y); +xmax = max(md.mesh.x); +%Generate basin IDs for 3 basins +idbasin = zeros(md.mesh.numberofelements,1); +iid1 = find(md.mesh.y>=2/3*ymax); +iid2 = intersect(find(md.mesh.y<2/3*ymax),find(md.mesh.x>=1/3*xmax)); +iid3 = intersect(find(md.mesh.y<2/3*ymax),find(md.mesh.x<1/3*xmax)); +for ii=1:md.mesh.numberofelements + for vertex=1:3 + if any(iid1==md.mesh.elements(ii,vertex)) %one vertex in basin 1 + idbasin(ii) = 1; + end + end + if idbasin(ii)==0 %no vertex was found in basin 1 + for vertex=1:3 + if any(iid2==md.mesh.elements(ii,vertex)) %one vertex in basin 2 + idbasin(ii) = 2; + end + end + end + if idbasin(ii)==0 %no vertex was found in basin 1 and 2 + idbasin(ii) = 3; + end +end + +md=solve(md,'Transient'); + +%Fields and tolerances to track changes +field_names ={'Vx1','Vy1','Vel1','Thickness1','Volume1','SmbMassBalance1','Vx2','Vy2','Vel2','Thickness2','Volume2','SmbMassBalance2','Vx3','Vy3','Vel3','Thickness3','Volume3','SmbMassBalance3'}; +field_tolerances={1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,... + 1e-13,1e-13,1e-13,1e-13,1e-13,1e-13... + 1e-13,1e-13,1e-13,1e-13,1e-13,1e-13}; +field_values={... + (md.results.TransientSolution(1).Vx),... + (md.results.TransientSolution(1).Vy),... + (md.results.TransientSolution(1).Vel),... + (md.results.TransientSolution(1).Pressure),... + (md.results.TransientSolution(1).Base),... + (md.results.TransientSolution(1).Surface),... + (md.results.TransientSolution(1).Thickness),... + (md.results.TransientSolution(1).IceVolume),... + (md.results.TransientSolution(1).SmbMassBalance),... + (md.results.TransientSolution(2).Vx),... + (md.results.TransientSolution(2).Vy),... + (md.results.TransientSolution(2).Vel),... + (md.results.TransientSolution(2).Pressure),... + (md.results.TransientSolution(2).Base),... + (md.results.TransientSolution(2).Surface),... + (md.results.TransientSolution(2).Thickness),... + (md.results.TransientSolution(2).IceVolume),... + (md.results.TransientSolution(2).SmbMassBalance),... + (md.results.TransientSolution(3).Vx),... + (md.results.TransientSolution(3).Vy),... + (md.results.TransientSolution(3).Vel),... + (md.results.TransientSolution(3).Pressure),... + (md.results.TransientSolution(3).Base),... + (md.results.TransientSolution(3).Surface),... + (md.results.TransientSolution(3).Thickness),... + (md.results.TransientSolution(3).IceVolume),... + (md.results.TransientSolution(3).SmbMassBalance),... + }; Index: /issm/trunk-jpl/test/NightlyRun/test542.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test542.m (revision 26487) +++ /issm/trunk-jpl/test/NightlyRun/test542.m (revision 26487) @@ -0,0 +1,83 @@ +%Test Name: PigTranFrontalforcingsrignot +md=triangle(model(),'../Exp/Pig.exp',10000.); +md=setmask(md,'../Exp/PigShelves.exp','../Exp/PigIslands.exp'); +md=parameterize(md,'../Par/Pig.par'); +md=setflowequation(md,'SSA','all'); +md.timestepping.time_step=2; +md.timestepping.final_time=50; + +%separate domain in 2 basins +idbasin = zeros(md.mesh.numberofelements,1); +iid1 = find(md.mesh.y>=mean(md.mesh.y)); +for ii=1:md.mesh.numberofelements + for vertex=1:3 + if any(iid1==md.mesh.elements(ii,vertex)) %one vertex in basin 1 + idbasin(ii) = 1; + end + end + if idbasin(ii)==0 %no vertex was found in basin 1 + idbasin(ii) = 2; + end +end + +%calving and frontalforcings parameters +md.mask.ice_levelset = 1e4*(md.mask.ice_levelset + 0.5); +md.calving=calvingvonmises(); +md.frontalforcings=frontalforcingsrignot(); +md.frontalforcings.num_basins = 2; +md.frontalforcings.basin_id = idbasin; +md.frontalforcings.subglacial_discharge = 0.1*ones(md.mesh.numberofvertices,1); +md.frontalforcings.thermalforcing = 0.5*ones(md.mesh.numberofvertices,1); +for elem=1:md.mesh.numberofelements + if idbasin(elem)==2 + md.frontalforcings.thermalforcing(md.mesh.elements(elem,1:3)) = 1.5; + end +end +md.transient.ismovingfront = 1; +md.levelset.spclevelset = NaN(md.mesh.numberofvertices,1); +md.levelset.migration_max = 1e10; + +md.transient.requested_outputs = {'default','CalvingMeltingrate'}; +md.cluster=generic('name',oshostname(),'np',2); +md=solve(md,'Transient'); + +%Fields and tolerances to track changes +field_names ={... + 'Vx1' ,'Vy1' ,'Vel1' ,'Pressure1' ,'Bed1' ,'Surface1' ,'Thickness1' ,'MaskIceLevelset1', 'CalvingMeltingrate1',... + 'Vx2' ,'Vy2' ,'Vel2' ,'Pressure2' ,'Bed2' ,'Surface2' ,'Thickness2' ,'MaskIceLevelset2', 'CalvingMeltingrate2',... + 'Vx10','Vy10','Vel10','Pressure10','Bed10','Surface10','Thickness10','MaskIceLevelset10', 'CalvingMeltingrate10',... + }; +field_tolerances={... + 1e-11,2e-11,2e-11,1e-12,2e-11,6e-12,9e-12,1e-12,1e-9,... + 2e-11,1e-11,1e-11,9e-12,2e-1,2e-11,2e-11,1e-11,1e-9,... + 2e-10,1e-10,1e-10,1e-10,1e-10,1e-10,1e-10,1e-9,1e-9,... + }; +field_values={... + (md.results.TransientSolution(1).Vx),... + (md.results.TransientSolution(1).Vy),... + (md.results.TransientSolution(1).Vel),... + (md.results.TransientSolution(1).Pressure),... + (md.results.TransientSolution(1).Base),... + (md.results.TransientSolution(1).Surface),... + (md.results.TransientSolution(1).Thickness),... + (md.results.TransientSolution(1).CalvingMeltingrate),... + (md.results.TransientSolution(1).MaskIceLevelset),... + (md.results.TransientSolution(2).Vx),... + (md.results.TransientSolution(2).Vy),... + (md.results.TransientSolution(2).Vel),... + (md.results.TransientSolution(2).Pressure),... + (md.results.TransientSolution(2).Base),... + (md.results.TransientSolution(2).Surface),... + (md.results.TransientSolution(2).Thickness),... + (md.results.TransientSolution(2).MaskIceLevelset),... + (md.results.TransientSolution(2).CalvingMeltingrate),... + (md.results.TransientSolution(10).Vx),... + (md.results.TransientSolution(10).Vy),... + (md.results.TransientSolution(10).Vel),... + (md.results.TransientSolution(10).Pressure),... + (md.results.TransientSolution(10).Base),... + (md.results.TransientSolution(10).Surface),... + (md.results.TransientSolution(10).Thickness),... + (md.results.TransientSolution(10).MaskIceLevelset),... + (md.results.TransientSolution(10).CalvingMeltingrate),... + };