Index: /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 21468) +++ /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 21469) @@ -105,4 +105,11 @@ iomodel->FetchDataToInput(elements,"md.smb.b_neg",SmbBNegEnum); break; + case SMBgradientselaEnum: + iomodel->FetchDataToInput(elements,"md.smb.ela",SmbElaEnum); + iomodel->FetchDataToInput(elements,"md.smb.b_pos",SmbBPosEnum); + iomodel->FetchDataToInput(elements,"md.smb.b_neg",SmbBNegEnum); + iomodel->FetchDataToInput(elements,"md.smb.b_max",SmbBMaxEnum); + iomodel->FetchDataToInput(elements,"md.smb.b_min",SmbBMinEnum); + break; case SMBhenningEnum: iomodel->FetchDataToInput(elements,"md.smb.smbref",SmbSmbrefEnum,0.); @@ -217,4 +224,7 @@ /*Nothing to add to parameters*/ break; + case SMBgradientselaEnum: + /*Nothing to add to parameters*/ + break; case SMBhenningEnum: /*Nothing to add to parameters*/ @@ -282,4 +292,8 @@ SmbGradientsx(femmodel); break; + case SMBgradientselaEnum: + if(VerboseSolution())_printf0_(" call smb gradients ela module\n"); + SmbGradientsElax(femmodel); + break; case SMBhenningEnum: if(VerboseSolution())_printf0_(" call smb Henning module\n"); Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 21468) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 21469) @@ -17,4 +17,5 @@ IssmDouble rho_water; // density of fresh water IssmDouble rho_ice; // density of ice + IssmDouble yts; // conversion factor year to second /*Loop over all the elements of this partition*/ @@ -44,4 +45,7 @@ rho_water=element->matpar->GetMaterialParameter(MaterialsRhoFreshwaterEnum); + /* Get constants */ + femmodel->parameters->FindParam(&yts,ConstantsYtsEnum); + // loop over all vertices for(v=0;v(s); xDelete(smb); + } + +}/*}}}*/ +void SmbGradientsElax(FemModel* femmodel){/*{{{*/ + + // void SurfaceMassBalancex(hd,agd,ni){ + // INPUT parameters: ni: working size of arrays + // INPUT: surface elevation (m): hd(NA) + // OUTPUT: surface mass-balance (m/yr ice): agd(NA) + int v; + + /*Loop over all the elements of this partition*/ + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + + /*Allocate all arrays*/ + int numvertices = element->GetNumberOfVertices(); + IssmDouble* ela = xNew(numvertices); // Equilibrium Line Altitude (m a.s.l) to which deviations are used to calculate the SMB + IssmDouble* b_pos = xNew(numvertices); // SMB gradient above ELA (m ice eq. per m elevation change) + IssmDouble* b_neg = xNew(numvertices); // SMB gradient below ELA (m ice eq. per m elevation change) + IssmDouble* b_max = xNew(numvertices); // Upper cap on SMB rate (m/y ice eq.) + IssmDouble* b_min = xNew(numvertices); // Lower cap on SMB rate (m/y ice eq.) + IssmDouble* s = xNew(numvertices); // Surface elevation (m a.s.l.) + IssmDouble* smb = xNew(numvertices); // SMB (m/y ice eq.) + + /*Recover ELA, SMB gradients, and caps*/ + element->GetInputListOnVertices(ela,SmbElaEnum); + element->GetInputListOnVertices(b_pos,SmbBPosEnum); + element->GetInputListOnVertices(b_neg,SmbBNegEnum); + element->GetInputListOnVertices(b_max,SmbBMaxEnum); + element->GetInputListOnVertices(b_min,SmbBMinEnum); + + /*Recover surface elevation at vertices: */ + element->GetInputListOnVertices(s,SurfaceEnum); + + /*Loop over all vertices, calculate SMB*/ + for(v=0;vela[v]){ + smb[v]=b_pos[v]*(s[v]-ela[v]); + } + // if surface is below or equal to the ELA + else{ + smb[v]=b_neg[v]*(s[v]-ela[v]); + } + + // if SMB is larger than upper cap, set SMB to upper cap + if(smb[v]>b_max[v]){ + smb[v]=b_max[v]; + } + // if SMB is smaller than lower cap, set SMB to lower cap + if(smb[v]AddInput(SmbMassBalanceEnum,smb,P1Enum); + xDelete(ela); + xDelete(b_pos); + xDelete(b_neg); + xDelete(b_max); + xDelete(b_min); + xDelete(s); + xDelete(smb); + } Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 21468) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 21469) @@ -11,4 +11,5 @@ void SurfaceMassBalancex(FemModel* femmodel); void SmbGradientsx(FemModel* femmodel); +void SmbGradientsElax(FemModel* femmodel); void Delta18oParameterizationx(FemModel* femmodel); void MungsmtpParameterizationx(FemModel* femmodel); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21468) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21469) @@ -464,4 +464,9 @@ SmbRefreezeEnum, SMBgcmEnum, + /*SMBgradientsela*/ + SMBgradientselaEnum, + SmbElaEnum, + SmbBMaxEnum, + SmbBMinEnum, /*}}}*/ /*Inputs {{{*/ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21468) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21469) @@ -465,4 +465,8 @@ case SmbRefreezeEnum : return "SmbRefreeze"; case SMBgcmEnum : return "SMBgcm"; + case SMBgradientselaEnum : return "SMBgradientsela"; + case SmbElaEnum : return "SmbEla"; + case SmbBMaxEnum : return "SmbBMax"; + case SmbBMinEnum : return "SmbBMin"; case AdjointpEnum : return "Adjointp"; case AdjointxEnum : return "Adjointx"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21468) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21469) @@ -474,4 +474,8 @@ else if (strcmp(name,"SmbRefreeze")==0) return SmbRefreezeEnum; else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum; + else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum; + else if (strcmp(name,"SmbEla")==0) return SmbElaEnum; + else if (strcmp(name,"SmbBMax")==0) return SmbBMaxEnum; + else if (strcmp(name,"SmbBMin")==0) return SmbBMinEnum; else if (strcmp(name,"Adjointp")==0) return AdjointpEnum; else if (strcmp(name,"Adjointx")==0) return AdjointxEnum; @@ -502,12 +506,12 @@ else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum; else if (strcmp(name,"TemperaturePDD")==0) return TemperaturePDDEnum; - else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum; + else stage=5; + } + if(stage==5){ + if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum; else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum; else if (strcmp(name,"Vel")==0) return VelEnum; else if (strcmp(name,"Velocity")==0) return VelocityEnum; - else stage=5; - } - if(stage==5){ - if (strcmp(name,"VxAverage")==0) return VxAverageEnum; + else if (strcmp(name,"VxAverage")==0) return VxAverageEnum; else if (strcmp(name,"Vx")==0) return VxEnum; else if (strcmp(name,"VxPicard")==0) return VxPicardEnum; @@ -625,12 +629,12 @@ else if (strcmp(name,"Outputdefinition40")==0) return Outputdefinition40Enum; else if (strcmp(name,"Outputdefinition41")==0) return Outputdefinition41Enum; - else if (strcmp(name,"Outputdefinition42")==0) return Outputdefinition42Enum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"Outputdefinition42")==0) return Outputdefinition42Enum; else if (strcmp(name,"Outputdefinition43")==0) return Outputdefinition43Enum; else if (strcmp(name,"Outputdefinition44")==0) return Outputdefinition44Enum; else if (strcmp(name,"Outputdefinition45")==0) return Outputdefinition45Enum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"Outputdefinition46")==0) return Outputdefinition46Enum; + else if (strcmp(name,"Outputdefinition46")==0) return Outputdefinition46Enum; else if (strcmp(name,"Outputdefinition47")==0) return Outputdefinition47Enum; else if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum; @@ -748,12 +752,12 @@ else if (strcmp(name,"Gsl")==0) return GslEnum; else if (strcmp(name,"Cuffey")==0) return CuffeyEnum; - else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum; else if (strcmp(name,"CuffeyTemperate")==0) return CuffeyTemperateEnum; else if (strcmp(name,"Paterson")==0) return PatersonEnum; else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum; - else stage=7; - } - if(stage==7){ - if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum; + else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum; else if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum; else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum; @@ -871,12 +875,12 @@ else if (strcmp(name,"BalancethicknessSolution")==0) return BalancethicknessSolutionEnum; else if (strcmp(name,"Balancethickness2Analysis")==0) return Balancethickness2AnalysisEnum; - else if (strcmp(name,"Balancethickness2Solution")==0) return Balancethickness2SolutionEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"Balancethickness2Solution")==0) return Balancethickness2SolutionEnum; else if (strcmp(name,"BalancethicknessSoftAnalysis")==0) return BalancethicknessSoftAnalysisEnum; else if (strcmp(name,"BalancethicknessSoftSolution")==0) return BalancethicknessSoftSolutionEnum; else if (strcmp(name,"BalancevelocityAnalysis")==0) return BalancevelocityAnalysisEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum; + else if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum; else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum; else if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum; @@ -994,7 +998,10 @@ else if (strcmp(name,"Parameters")==0) return ParametersEnum; else if (strcmp(name,"Vertices")==0) return VerticesEnum; - else if (strcmp(name,"Results")==0) return ResultsEnum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"Results")==0) return ResultsEnum; else if (strcmp(name,"MaximumNumberOfDefinitions")==0) return MaximumNumberOfDefinitionsEnum; - else stage=9; + else stage=10; } /*If we reach this point, the string provided has not been found*/ Index: /issm/trunk-jpl/src/m/classes/SMBgradientsela.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBgradientsela.m (revision 21469) +++ /issm/trunk-jpl/src/m/classes/SMBgradientsela.m (revision 21469) @@ -0,0 +1,88 @@ +%SMBgradientsela Class definition +% +% Usage: +% SMBgradientsela=SMBgradientsela(); + +classdef SMBgradientsela + properties (SetAccess=public) + ela = NaN; + b_pos = NaN; + b_neg = NaN; + b_max = NaN; + b_min = NaN; + requested_outputs = {}; + end + methods + function self = SMBgradientsela(varargin) % {{{ + switch nargin + case 0 + self=setdefaultparameters(self); + otherwise + error('constructor not supported'); + end + end % }}} + function self = extrude(self,md) % {{{ + + %Nothing for now + + end % }}} + function list = defaultoutputs(self,md) % {{{ + list = {''}; + end % }}} + function self = initialize(self,md) % {{{ + + %Nothing done for now + + end % }}} + function self = setdefaultparameters(self) % {{{ + + self.b_max=9999; + self.b_min=-9999; + + end % }}} + function md = checkconsistency(self,md,solution,analyses) % {{{ + + if ismember('MasstransportAnalysis',analyses), + md = checkfield(md,'fieldname','smb.ela','timeseries',1,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','smb.b_pos','timeseries',1,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','smb.b_neg','timeseries',1,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','smb.b_max','timeseries',1,'NaN',1,'Inf',1); + md = checkfield(md,'fieldname','smb.b_min','timeseries',1,'NaN',1,'Inf',1); + end + md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1); + end % }}} + function disp(self) % {{{ + disp(sprintf(' surface forcings parameters:')); + + disp(sprintf('\n SMB gradients ela parameters:')); + fielddisplay(self,'ela',' equilibrium line altitude from which deviation is used to calculate smb using the smb gradients ela method [m a.s.l.]'); + fielddisplay(self,'b_pos',' vertical smb gradient (dB/dz) above ela'); + fielddisplay(self,'b_neg',' vertical smb gradient (dB/dz) below ela'); + fielddisplay(self,'b_max',' upper cap on smb rate, default: 9999 (no cap) [m ice eq./yr] '); + fielddisplay(self,'b_min',' lower cap on smb rate, default: -9999 (no cap) [m ice eq./yr]'); + fielddisplay(self,' requested_outputs','additional outputs requested'); + + end % }}} + function marshall(self,prefix,md,fid) % {{{ + + yts=md.constants.yts; + + WriteData(fid,prefix,'name','md.smb.model','data',9,'format','Integer'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','ela','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','b_pos','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','b_neg','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','b_max','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','b_min','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + + %process requested outputs + outputs = self.requested_outputs; + pos = find(ismember(outputs,'default')); + if ~isempty(pos), + outputs(pos) = []; %remove 'default' from outputs + outputs = [outputs defaultoutputs(self,md)]; %add defaults + end + WriteData(fid,prefix,'data',outputs,'name','md.smb.requested_outputs','format','StringArray'); + + end % }}} + end +end Index: /issm/trunk-jpl/src/m/classes/SMBgradientsela.py =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBgradientsela.py (revision 21469) +++ /issm/trunk-jpl/src/m/classes/SMBgradientsela.py (revision 21469) @@ -0,0 +1,80 @@ +from fielddisplay import fielddisplay +from checkfield import checkfield +from WriteData import WriteData +from project3d import project3d + +class SMBgradientsela(object): + """ + SMBgradientsela Class definition + + Usage: + SMBgradientsela=SMBgradientsela(); + """ + + def __init__(self): # {{{ + self.ela = float('NaN') + self.b_pos = float('NaN') + self.b_neg = float('NaN') + self.b_max = float('NaN') + self.b_min = float('NaN') + self.requested_outputs = [] + #}}} + def __repr__(self): # {{{ + string=" surface forcings parameters:" + + string="%s\n%s"%(string,fielddisplay(self,'issmbgradientsela','is smb gradients ela method activated (0 or 1, default is 0)')) + string="%s\n%s"%(string,fielddisplay(self,'ela',' equilibrium line altitude from which deviation is used to calculate smb using the smb gradients ela method [m a.s.l.]')) + string="%s\n%s"%(string,fielddisplay(self,'b_pos',' vertical smb gradient (dB/dz) above ela')) + string="%s\n%s"%(string,fielddisplay(self,'b_neg',' vertical smb gradient (dB/dz) below ela')) + string="%s\n%s"%(string,fielddisplay(self,'b_max',' upper cap on smb rate, default: 9999 (no cap) [m ice eq./yr]')) + string="%s\n%s"%(string,fielddisplay(self,'b_min',' lower cap on smb rate, default: -9999 (no cap) [m ice eq./yr]')) + string="%s\n%s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested')) + + return string + #}}} + def extrude(self,md): # {{{ + + #Nothing for now + return self + #}}} + def defaultoutputs(self,md): # {{{ + return [] + #}}} + def initialize(self,md): # {{{ + + #Nothing for now + + return self + #}}} + def checkconsistency(self,md,solution,analyses): # {{{ + + if 'MasstransportAnalysis' in analyses: + md = checkfield(md,'fieldname','smb.ela','timeseries',1,'NaN',1,'Inf',1) + md = checkfield(md,'fieldname','smb.b_pos','timeseries',1,'NaN',1,'Inf',1) + md = checkfield(md,'fieldname','smb.b_neg','timeseries',1,'NaN',1,'Inf',1) + md = checkfield(md,'fieldname','smb.b_max','timeseries',1,'NaN',1,'Inf',1) + md = checkfield(md,'fieldname','smb.b_min','timeseries',1,'NaN',1,'Inf',1) + + md = checkfield(md,'fieldname','masstransport.requested_outputs','stringrow',1) + return md + # }}} + def marshall(self,prefix,md,fid): # {{{ + + yts=md.constants.yts + + WriteData(fid,prefix,'name','md.smb.model','data',9,'format','Integer'); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','ela','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) + WriteData(fid,prefix,'object',self,'class','smb','fieldname','b_pos','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) + WriteData(fid,prefix,'object',self,'class','smb','fieldname','b_neg','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) + WriteData(fid,prefix,'object',self,'class','smb','fieldname','b_max','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) + WriteData(fid,prefix,'object',self,'class','smb','fieldname','b_min','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts) + + #process requested outputs + outputs = self.requested_outputs + indices = [i for i, x in enumerate(outputs) if x == 'default'] + if len(indices) > 0: + outputscopy=outputs[0:max(0,indices[0]-1)]+self.defaultoutputs(md)+outputs[indices[0]+1:] + outputs =outputscopy + WriteData(fid,prefix,'data',outputs,'name','md.smb.requested_outputs','format','StringArray') + + # }}} Index: /issm/trunk-jpl/test/NightlyRun/test343.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test343.m (revision 21469) +++ /issm/trunk-jpl/test/NightlyRun/test343.m (revision 21469) @@ -0,0 +1,50 @@ +%Test Name: SquareSheetConstrainedSmbGradientsEla2d +md=triangle(model(),'../Exp/Square.exp',150000.); +md=setmask(md,'',''); +md=parameterize(md,'../Par/SquareSheetConstrained.par'); +md=setflowequation(md,'SSA','all'); +md.smb = SMBgradientsela(); +md.smb.ela=1500.*ones(md.mesh.numberofvertices+1,1); +md.smb.b_pos=0.002.*ones(md.mesh.numberofvertices+1,1); +md.smb.b_neg=0.005.*ones(md.mesh.numberofvertices+1,1); +md.smb.b_max=4.*(md.materials.rho_freshwater/md.materials.rho_ice).*ones(md.mesh.numberofvertices+1,1); +md.smb.b_min=-4.*(md.materials.rho_freshwater/md.materials.rho_ice).*ones(md.mesh.numberofvertices+1,1); + +%Change geometry +md.geometry.thickness=md.geometry.surface*30; +md.geometry.surface=md.geometry.base+md.geometry.thickness; + +%Transient options +md.transient.requested_outputs={'default','TotalSmb'}; +md.cluster=generic('name',oshostname(),'np',3); +md=solve(md,'Transient'); + +%Fields and tolerances to track changes +field_names ={'Vx1','Vy1','Vel1','Bed1','Surface1','Thickness1','SMB1','TotalSmb1','Vx2','Vy2','Vel2','Bed2','Surface2','Thickness2','SMB2','TotalSmb2','Vx3','Vy3','Vel3','Bed3','Surface3','Thickness3','SMB3','TotalSmb3'}; +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,1e-13,1e-13,1e-13,1e-13,1.5e-13,1e-13}; +field_values={... + (md.results.TransientSolution(1).Vx),... + (md.results.TransientSolution(1).Vy),... + (md.results.TransientSolution(1).Vel),... + (md.results.TransientSolution(1).Base),... + (md.results.TransientSolution(1).Surface),... + (md.results.TransientSolution(1).Thickness),... + (md.results.TransientSolution(1).SmbMassBalance),... + (md.results.TransientSolution(1).TotalSmb),... + (md.results.TransientSolution(2).Vx),... + (md.results.TransientSolution(2).Vy),... + (md.results.TransientSolution(2).Vel),... + (md.results.TransientSolution(2).Base),... + (md.results.TransientSolution(2).Surface),... + (md.results.TransientSolution(2).Thickness),... + (md.results.TransientSolution(2).TotalSmb),... + (md.results.TransientSolution(2).SmbMassBalance),... + (md.results.TransientSolution(3).Vx),... + (md.results.TransientSolution(3).Vy),... + (md.results.TransientSolution(3).Vel),... + (md.results.TransientSolution(3).Base),... + (md.results.TransientSolution(3).Surface),... + (md.results.TransientSolution(3).Thickness),... + (md.results.TransientSolution(3).SmbMassBalance),... + (md.results.TransientSolution(3).TotalSmb),... + }; Index: /issm/trunk-jpl/test/NightlyRun/test344.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test344.m (revision 21469) +++ /issm/trunk-jpl/test/NightlyRun/test344.m (revision 21469) @@ -0,0 +1,61 @@ +%Test Name: SquareSheetConstrainedSmbGradientsEla3d +md=triangle(model(),'../Exp/Square.exp',150000.); +md=setmask(md,'',''); +md=parameterize(md,'../Par/SquareSheetConstrained.par'); + +%Change geometry +md.geometry.thickness=md.geometry.surface*30; +md.geometry.surface=md.geometry.base+md.geometry.thickness; + +md=extrude(md,3,1.); +md=setflowequation(md,'HO','all'); +md.smb = SMBgradientsela(); +md.smb.ela=1500.*ones(md.mesh.numberofvertices+1,1); +md.smb.b_pos=0.002.*ones(md.mesh.numberofvertices+1,1); +md.smb.b_neg=0.005.*ones(md.mesh.numberofvertices+1,1); +md.smb.b_max=4.*(md.materials.rho_freshwater/md.materials.rho_ice).*ones(md.mesh.numberofvertices+1,1); +md.smb.b_min=-4.*(md.materials.rho_freshwater/md.materials.rho_ice).*ones(md.mesh.numberofvertices+1,1); + + +%Transient options +md.transient.requested_outputs={'default','TotalSmb'}; +md.cluster=generic('name',oshostname(),'np',3); +md=solve(md,'Transient'); + +%Fields and tolerances to track changes +field_names ={'Vx1','Vy1','Vz1','Vel1','Bed1','Surface1','Thickness1','Temperature1','SMB1','TotalSmb1','Vx2','Vy2','Vz2','Vel2','Bed2','Surface2','Thickness2','Temperature2','SMB2','TotalSmb2','Vx3','Vy3','Vz3','Vel3','Bed3','Surface3','Thickness3','Temperature3','SMB3','TotalSmb3'}; +field_tolerances={1e-09,1e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,1e-10,... + 1e-09,1e-09,1e-10,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,1e-10,... + 1e-09,5e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,1e-10}; +field_values={... + (md.results.TransientSolution(1).Vx),... + (md.results.TransientSolution(1).Vy),... + (md.results.TransientSolution(1).Vz),... + (md.results.TransientSolution(1).Vel),... + (md.results.TransientSolution(1).Base),... + (md.results.TransientSolution(1).Surface),... + (md.results.TransientSolution(1).Thickness),... + (md.results.TransientSolution(1).Temperature),... + (md.results.TransientSolution(1).SmbMassBalance),... + (md.results.TransientSolution(1).TotalSmb),... + (md.results.TransientSolution(2).Vx),... + (md.results.TransientSolution(2).Vy),... + (md.results.TransientSolution(2).Vz),... + (md.results.TransientSolution(2).Vel),... + (md.results.TransientSolution(2).Base),... + (md.results.TransientSolution(2).Surface),... + (md.results.TransientSolution(2).Thickness),... + (md.results.TransientSolution(2).Temperature),... + (md.results.TransientSolution(2).SmbMassBalance),... + (md.results.TransientSolution(2).TotalSmb),... + (md.results.TransientSolution(3).Vx),... + (md.results.TransientSolution(3).Vy),... + (md.results.TransientSolution(3).Vz),... + (md.results.TransientSolution(3).Vel),... + (md.results.TransientSolution(3).Base),... + (md.results.TransientSolution(3).Surface),... + (md.results.TransientSolution(3).Thickness),... + (md.results.TransientSolution(3).Temperature),... + (md.results.TransientSolution(3).SmbMassBalance),... + (md.results.TransientSolution(3).TotalSmb),... + };