Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 26554) +++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 26555) @@ -20,7 +20,6 @@ //#define FSANALYTICAL 10 //#define LATERALFRICTION 1 //#define DISCSLOPE 1 //testing for SSA -#define NOSPCSHEARVEL 1 //MLHO /*Model processing*/ void StressbalanceAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ @@ -195,16 +194,10 @@ } } else{//MLHO - #ifdef NOSPCSHEARVEL - IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,0); - IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,2); - #else - /*Default: apply spcs to shear vx and shear vy*/ - IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,0); - IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx",StressbalanceAnalysisEnum,finiteelement,1); - IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,2); - IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy",StressbalanceAnalysisEnum,finiteelement,3); - #endif + IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx_base",StressbalanceAnalysisEnum,finiteelement,0); + IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvx_shear",StressbalanceAnalysisEnum,finiteelement,1); + IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy_base",StressbalanceAnalysisEnum,finiteelement,2); + IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvy_shear",StressbalanceAnalysisEnum,finiteelement,3); } } Index: ../trunk-jpl/src/m/boundaryconditions/SetMLHOBC.m =================================================================== --- ../trunk-jpl/src/m/boundaryconditions/SetMLHOBC.m (nonexistent) +++ ../trunk-jpl/src/m/boundaryconditions/SetMLHOBC.m (revision 26555) @@ -0,0 +1,16 @@ +function md=SetMLHOBC(md) +%SETMLHOBC - Create the boundary conditions for stressbalance for MLHO: VxBase, VyBase, VxShear, VyShear +% +% Usage: +% md=SetMLHOBC(md) +% + + +%node on Dirichlet +if md.flowequation.isMLHO + md.stressbalance.spcvx_base=md.stressbalance.spcvx; + md.stressbalance.spcvy_base=md.stressbalance.spcvy; + + md.stressbalance.spcvx_shear=NaN*ones(size(md.stressbalance.spcvx_base)); + md.stressbalance.spcvy_shear=NaN*ones(size(md.stressbalance.spcvy_base)); +end Index: ../trunk-jpl/src/m/classes/flowequation.m =================================================================== --- ../trunk-jpl/src/m/classes/flowequation.m (revision 26554) +++ ../trunk-jpl/src/m/classes/flowequation.m (revision 26555) @@ -61,10 +61,9 @@ if isfield(objstruct,'borderstokes'), self.borderFS = objstruct.borderstokes; end; end - %Oct 8 2020 - %Warning if MLHO + %Nov 6 2021 if any(self.vertex_equation==4) - warning(['Monolayer Higher-Order (MLHO) detected in md.flowequation, this is probably a mistake, you need to run setflowequation again']); + disp(['Monolayer Higher-Order (MLHO) detected in md.flowequation, this is still under development. Please double check your settings.']); end end% }}} Index: ../trunk-jpl/src/m/classes/stressbalance.m =================================================================== --- ../trunk-jpl/src/m/classes/stressbalance.m (revision 26554) +++ ../trunk-jpl/src/m/classes/stressbalance.m (revision 26555) @@ -7,6 +7,10 @@ properties (SetAccess=public) spcvx = NaN; spcvy = NaN; + spcvx_base = NaN; + spcvy_base = NaN; + spcvx_shear = NaN; + spcvy_shear = NaN; spcvz = NaN; restol = 0; reltol = 0; @@ -31,6 +35,13 @@ self.referential=project3d(md,'vector',self.referential,'type','node'); self.loadingforce=project3d(md,'vector',self.loadingforce,'type','node'); + % for MLHO + if md.flowequation.isMLHO + self.spcvx_base=project3d(md,'vector',self.spcvx_base,'type','node'); + self.spcvy_base=project3d(md,'vector',self.spcvy_base,'type','node'); + self.spcvx_shear=project3d(md,'vector',self.spcvx_shear,'type','poly','degree',4); + self.spcvy_shear=project3d(md,'vector',self.spcvy_shear,'type','poly','degree',4); + end end % }}} function self = stressbalance(varargin) % {{{ switch nargin @@ -122,6 +133,13 @@ end md = checkfield(md,'fieldname','stressbalance.FSreconditioning','>',0); end + % CHECK THIS ONLY WORKS FOR MLHO + if md.flowequation.isMLHO + md = checkfield(md,'fieldname','stressbalance.spcvx_base','Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','stressbalance.spcvy_base','Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','stressbalance.spcvx_shear','Inf',1,'timeseries',1); + md = checkfield(md,'fieldname','stressbalance.spcvy_shear','Inf',1,'timeseries',1); + end end % }}} function list=defaultoutputs(self,md) % {{{ @@ -150,6 +168,12 @@ fielddisplay(self,'spcvy','y-axis velocity constraint (NaN means no constraint) [m/yr]'); fielddisplay(self,'spcvz','z-axis velocity constraint (NaN means no constraint) [m/yr]'); + disp(sprintf('\n %s','MLHO boundary conditions:')); + fielddisplay(self,'spcvx_base','x-axis basal velocity constraint (NaN means no constraint) [m/yr]'); + fielddisplay(self,'spcvy_base','y-axis basal velocity constraint (NaN means no constraint) [m/yr]'); + fielddisplay(self,'spcvx_shear','x-axis shear velocity constraint (NaN means no constraint) [m/yr]'); + fielddisplay(self,'spcvy_shear','y-axis shear velocity constraint (NaN means no constraint) [m/yr]'); + disp(sprintf('\n %s','Rift options:')); fielddisplay(self,'rift_penalty_threshold','threshold for instability of mechanical constraints'); fielddisplay(self,'rift_penalty_lock','number of iterations before rift penalties are locked'); @@ -201,6 +225,13 @@ outputs = [outputs defaultoutputs(self,md)]; %add defaults end WriteData(fid,prefix,'data',outputs,'name','md.stressbalance.requested_outputs','format','StringArray'); + % for MLHO + if (md.flowequation.isMLHO) + WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','spcvx_base','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','spcvy_base','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','spcvx_shear','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','spcvy_shear','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); + end end % }}} function savemodeljs(self,fid,modelname) % {{{ @@ -222,6 +253,10 @@ writejs2Darray(fid,[modelname '.stressbalance.loadingforce'],self.loadingforce); writejscellstring(fid,[modelname '.stressbalance.requested_outputs'],self.requested_outputs); + writejs1Darray(fid,[modelname '.stressbalance.spcvx_base'],self.spcvx_shear); + writejs1Darray(fid,[modelname '.stressbalance.spcvy_base'],self.spcvy_shear); + writejs1Darray(fid,[modelname '.stressbalance.spcvx_shear'],self.spcvx_shear); + writejs1Darray(fid,[modelname '.stressbalance.spcvy_shear'],self.spcvy_shear); end % }}} end end Index: ../trunk-jpl/src/m/classes/stressbalance.py =================================================================== --- ../trunk-jpl/src/m/classes/stressbalance.py (revision 26554) +++ ../trunk-jpl/src/m/classes/stressbalance.py (revision 26555) @@ -19,6 +19,10 @@ def __init__(self): # {{{ self.spcvx = np.nan self.spcvy = np.nan + self.spcvx_base = np.nan + self.spcvy_base = np.nan + self.spcvx_shear = np.nan + self.spcvy_shear = np.nan self.spcvz = np.nan self.restol = 0 self.reltol = 0 @@ -54,6 +58,11 @@ s += '{}\n'.format(fielddisplay(self, 'spcvy', 'y-axis velocity constraint (NaN means no constraint) [m / yr]')) s += '{}\n'.format(fielddisplay(self, 'spcvz', 'z-axis velocity constraint (NaN means no constraint) [m / yr]')) s += '{}\n'.format(fielddisplay(self, 'icefront', 'segments on ice front list (last column 0: Air, 1: Water, 2: Ice')) + s += ' MLHO boundary conditions:\n' + s += '{}\n'.format(fielddisplay(self, 'spcvx_base', 'x-axis basal velocity constraint (NaN means no constraint) [m / yr]')) + s += '{}\n'.format(fielddisplay(self, 'spcvy_base', 'y-axis basal velocity constraint (NaN means no constraint) [m / yr]')) + s += '{}\n'.format(fielddisplay(self, 'spcvx_shear', 'x-axis shear velocity constraint (NaN means no constraint) [m / yr]')) + s += '{}\n'.format(fielddisplay(self, 'spcvy_shear', 'y-axis shear velocity constraint (NaN means no constraint) [m / yr]')) s += ' Rift options:\n' s += '{}\n'.format(fielddisplay(self, 'rift_penalty_threshold', 'threshold for instability of mechanical constraints')) s += '{}\n'.format(fielddisplay(self, 'rift_penalty_lock', 'number of iterations before rift penalties are locked')) @@ -76,6 +85,12 @@ self.referential = project3d(md, 'vector', self.referential, 'type', 'node') self.loadingforce = project3d(md, 'vector', self.loadingforce, 'type', 'node') + if md.flowequation.isMLHO: + self.spcvx_base = project3d(md, 'vector', self.spcvx_base, 'type', 'node') + self.spcvy_base = project3d(md, 'vector', self.spcvy_base, 'type', 'node') + self.spcvx_shear = project3d(md, 'vector', self.spcvx_shear, 'type', 'poly', 'degree', 4) + self.spcvy_shear = project3d(md, 'vector', self.spcvy_shear, 'type', 'poly', 'degree', 4) + return self # }}} @@ -160,6 +175,11 @@ pos = np.nonzero(np.logical_and(md.mask.ocean_levelset, md.mesh.vertexonbase)) if np.any(np.logical_not(np.isnan(md.stressbalance.referential[pos, :]))): md.checkmessage("no referential should be specified for basal vertices of grounded ice") + if md.flowequation.isMLHO: + md = checkfield(md, 'fieldname', 'stressbalance.spcvx_base', 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'fieldname', 'stressbalance.spcvy_base', 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'fieldname', 'stressbalance.spcvx_shear', 'Inf', 1, 'timeseries', 1) + md = checkfield(md, 'fieldname', 'stressbalance.spcvy_shear', 'Inf', 1, 'timeseries', 1) return md # }}} @@ -194,4 +214,10 @@ outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:] outputs = outputscopy WriteData(fid, prefix, 'data', outputs, 'name', 'md.stressbalance.requested_outputs', 'format', 'StringArray') + # MLHO + if md.flowequation.isMLHO: + WriteData(fid, prefix, 'object', self, 'class', 'stressbalance', 'fieldname', 'spcvx_base', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) + WriteData(fid, prefix, 'object', self, 'class', 'stressbalance', 'fieldname', 'spcvy_base', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) + WriteData(fid, prefix, 'object', self, 'class', 'stressbalance', 'fieldname', 'spcvx_shear', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) + WriteData(fid, prefix, 'object', self, 'class', 'stressbalance', 'fieldname', 'spcvy_shear', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) # }}} Index: ../trunk-jpl/src/m/extrusion/project3d.m =================================================================== --- ../trunk-jpl/src/m/extrusion/project3d.m (revision 26554) +++ ../trunk-jpl/src/m/extrusion/project3d.m (revision 26555) @@ -6,11 +6,12 @@ % element vector of size (md.mesh.numberofelements2d,N/A). % arguments: % 'vector': 2d vector -% 'type': 'element' or 'node'. +% 'type': 'element' or 'node' or 'poly' % options: % 'layer' a layer number where vector should keep its values. If not specified, all layers adopt the % value of the 2d vector. % 'padding': default to 0 (value adopted by other 3d layers not being projected0 +% 'degree': degree of polynomials when extrude from bottom to the top % % Egs: % extruded_vector=project3d(md,'vector',vector2d,'type','node','layer',1,'padding',NaN); @@ -32,6 +33,7 @@ type = getfieldvalue(options,'type'); %mandatory layer = getfieldvalue(options,'layer',0); %optional (do all layers otherwise) paddingvalue = getfieldvalue(options,'padding',0); %0 by default +polyexponent = getfieldvalue(options,'degree',0); %0 by default, 0-degree polynomial if length(vector2d)==1, projected_vector=vector2d; @@ -79,6 +81,29 @@ projected_vector( ((layer-1)*md.mesh.numberofelements2d+1):(layer*md.mesh.numberofelements2d),:)=vector2d; end +elseif strcmpi(type,'poly'), % interpolate values from 0 to 1 with a polynomial degree n + %Initialize 3d vector + if size(vector2d,1)==md.mesh.numberofvertices2d + projected_vector=paddingvalue*ones(md.mesh.numberofvertices, size(vector2d,2)); + elseif size(vector2d,1)==md.mesh.numberofvertices2d+1 + projected_vector=paddingvalue*ones(md.mesh.numberofvertices+1,size(vector2d,2)); + projected_vector(end,:)=vector2d(end,:); + vector2d=vector2d(1:end-1,:); + else + error('vector length not supported') + end + + polycoeff = [0:1./(md.mesh.numberoflayers-1):1]; + + %Fill in + if layer==0, + for i=1:md.mesh.numberoflayers, + projected_vector(((i-1)*md.mesh.numberofvertices2d+1):(i*md.mesh.numberofvertices2d),:)=vector2d*(1-(1-polycoeff(i)).^polyexponent); + end + else + projected_vector(((layer-1)*md.mesh.numberofvertices2d+1):(layer*md.mesh.numberofvertices2d),:)=vector2d*(1-(1-polycoeff(layer)).^polyexponent); + end + else error('project3d error message: unknown projection type'); end Index: ../trunk-jpl/test/NightlyRun/test127.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test127.m (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test127.m (revision 26555) @@ -17,7 +17,7 @@ massfluxatgate('name','MassFlux5','profilename',['../Exp/MassFlux5.exp'],'definitionstring','Outputdefinition5'),... massfluxatgate('name','MassFlux6','profilename',['../Exp/MassFlux6.exp'],'definitionstring','Outputdefinition6')... }; - +md=SetMLHOBC(md); md=solve(md,'Stressbalance'); %Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test128.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test128.m (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test128.m (revision 26555) @@ -6,6 +6,7 @@ md.cluster=generic('name',oshostname(),'np',3); md.transient.requested_outputs={'IceVolume','VxShear','VyShear','VxBase','VyBase','VxSurface','VySurface'}; +md=SetMLHOBC(md); md=solve(md,'Transient'); %Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test128.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test128.py (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test128.py (revision 26555) @@ -15,6 +15,7 @@ md.cluster = generic('name', gethostname(), 'np', 3) md.transient.requested_outputs = ['IceVolume','VxSurface','VySurface','VxShear','VyShear','VxBase','VyBase'] +md = SetMLHOBC(md); md = solve(md, 'Transient') #Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test129.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test129.m (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test129.m (revision 26555) @@ -13,6 +13,7 @@ md.timestepping.final_time=19; md.settings.output_frequency=2; +md=SetMLHOBC(md); md=solve(md,'Transient'); md2=solve(md,'Transient','restart',1); Index: ../trunk-jpl/test/NightlyRun/test129.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test129.py (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test129.py (revision 26555) @@ -22,6 +22,7 @@ # time steps and resolution md.timestepping.final_time = 19 md.settings.output_frequency = 2 +md = SetMLHOBC(md); md = solve(md, 'Transient') md2 = copy.deepcopy(md) Index: ../trunk-jpl/test/NightlyRun/test248.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test248.m (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test248.m (revision 26555) @@ -5,6 +5,7 @@ md=setflowequation(md,'MLHO','all'); md.cluster=generic('name',oshostname(),'np',3); md.stressbalance.requested_outputs={'default','VxSurface','VySurface','VxShear','VyShear','VxBase','VyBase'}; +md=SetMLHOBC(md); md=solve(md,'Stressbalance'); %Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test248.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test248.py (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test248.py (revision 26555) @@ -15,6 +15,7 @@ md = setflowequation(md, 'MLHO', 'all') md.cluster = generic('name', gethostname(), 'np', 3) md.stressbalance.requested_outputs = ['default', 'VxSurface', 'VySurface', 'VxShear', 'VyShear', 'VxBase', 'VyBase'] +md = SetMLHOBC(md); md = solve(md, 'Stressbalance') #Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test249.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test249.m (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test249.m (revision 26555) @@ -6,6 +6,7 @@ md.cluster=generic('name',oshostname(),'np',3); md.transient.requested_outputs={'default','FloatingArea','GroundedArea','TotalGroundedBmb','TotalFloatingBmb'}; md.basalforcings.floatingice_melting_rate(:)=1; +md=SetMLHOBC(md); md=solve(md,'Transient'); %Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test249.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test249.py (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test249.py (revision 26555) @@ -16,6 +16,7 @@ md = setflowequation(md, 'MLHO', 'all') md.cluster = generic('name', gethostname(), 'np', 3) md.transient.requested_outputs = ['default', 'FloatingArea', 'GroundedArea', 'TotalFloatingBmb', 'TotalGroundedBmb'] +md = SetMLHOBC(md); md = solve(md, 'Transient') Index: ../trunk-jpl/test/NightlyRun/test254.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test254.m (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test254.m (revision 26555) @@ -57,6 +57,7 @@ md.mask.ice_levelset=-1+nodeonicefront; md.stressbalance.requested_outputs={'default','VySurface','VyShear','VyBase'}; +md=SetMLHOBC(md); md=solve(md,'Stressbalance'); %create analytical solution: strain rate is constant = ((rho_ice*g*h)/4B)^3 (Paterson, 4th Edition, page 292. Index: ../trunk-jpl/test/NightlyRun/test254.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test254.py (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test254.py (revision 26555) @@ -69,6 +69,7 @@ md.mask.ice_levelset = -1 + nodeonicefront md.stressbalance.requested_outputs = ['default', 'VySurface', 'VyShear', 'VyBase'] +md = SetMLHOBC(md); md = solve(md, 'Stressbalance') # create analytical solution: strain rate is constant = ((rho_ice * g * h) / 4B)^3 (Paterson, 4th Edition, page 292. Index: ../trunk-jpl/test/NightlyRun/test255.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test255.m (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test255.m (revision 26555) @@ -5,6 +5,7 @@ md=setflowequation(md,'MLHO','all'); md.cluster=generic('name',oshostname(),'np',3); md.masstransport.hydrostatic_adjustment='Incremental'; +md=SetMLHOBC(md); md=solve(md,'Transient'); %Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test255.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test255.py (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test255.py (revision 26555) @@ -14,6 +14,7 @@ md = setflowequation(md, 'MLHO', 'all') md.cluster = generic('name', gethostname(), 'np', 3) md.masstransport.hydrostatic_adjustment = 'Incremental' +md = SetMLHOBC(md); md = solve(md, 'Transient') #Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test256.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test256.m (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test256.m (revision 26555) @@ -6,6 +6,7 @@ md.geometry.base=md.geometry.base+50.; md.geometry.surface=md.geometry.surface+50.; md.cluster=generic('name',oshostname(),'np',1); md.masstransport.hydrostatic_adjustment='Incremental'; +md=SetMLHOBC(md); md=solve(md,'Transient'); %Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test256.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test256.py (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test256.py (revision 26555) @@ -16,6 +16,7 @@ md.geometry.surface = md.geometry.surface + 50. md.cluster = generic('name', gethostname(), 'np', 1) md.masstransport.hydrostatic_adjustment = 'Incremental' +md = SetMLHOBC(md); md = solve(md, 'Transient') #Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test330.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test330.m (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test330.m (revision 26555) @@ -5,6 +5,7 @@ md=setflowequation(md,'MLHO','all'); md.cluster=generic('name',oshostname(),'np',3); md.stressbalance.requested_outputs={'default','VxSurface','VySurface','VxShear','VyShear','VxBase','VyBase'}; +md=SetMLHOBC(md); md=solve(md,'Stressbalance'); %Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test330.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test330.py (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test330.py (revision 26555) @@ -14,6 +14,7 @@ md = setflowequation(md, 'MLHO', 'all') md.cluster = generic('name', gethostname(), 'np', 3) md.stressbalance.requested_outputs = ['default', 'VxSurface', 'VySurface', 'VxShear', 'VyShear', 'VxBase', 'VyBase'] +md = SetMLHOBC(md); md = solve(md, 'Stressbalance') #Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test332.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test332.m (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test332.m (revision 26555) @@ -4,6 +4,7 @@ md=parameterize(md,'../Par/SquareSheetConstrained.par'); md=setflowequation(md,'MLHO','all'); md.cluster=generic('name',oshostname(),'np',3); +md=SetMLHOBC(md); md=solve(md,'Transient'); %Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test332.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test332.py (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test332.py (revision 26555) @@ -13,6 +13,7 @@ md = parameterize(md, '../Par/SquareSheetConstrained.py') md = setflowequation(md, 'MLHO', 'all') md.cluster = generic('name', gethostname(), 'np', 3) +md = SetMLHOBC(md); md = solve(md, 'Transient') #Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test333.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test333.m (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test333.m (revision 26555) @@ -18,6 +18,7 @@ md.inversion.vx_obs=md.initialization.vx; md.inversion.vy_obs=md.initialization.vy; md.cluster=generic('name',oshostname(),'np',3); +md=SetMLHOBC(md); md=solve(md,'Stressbalance'); %Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test333.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test333.py (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test333.py (revision 26555) @@ -31,6 +31,7 @@ md.inversion.vy_obs = md.initialization.vy md.cluster = generic('name', gethostname(), 'np', 3) +md = SetMLHOBC(md); md = solve(md, 'Stressbalance') #Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test334.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test334.m (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test334.m (revision 26555) @@ -19,6 +19,7 @@ md.verbose.control=true; md.cluster=generic('name',oshostname(),'np',3); +md=SetMLHOBC(md); md=solve(md,'Stressbalance'); %Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test334.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test334.py (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test334.py (revision 26555) @@ -31,6 +31,7 @@ md.cluster = generic('name', gethostname(), 'np', 3) +md = SetMLHOBC(md); md = solve(md, 'Stressbalance') Index: ../trunk-jpl/test/NightlyRun/test335.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test335.m (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test335.m (revision 26555) @@ -6,6 +6,7 @@ md = extrude(md, 5, 1); md.cluster=generic('name',oshostname(),'np',3); md.stressbalance.requested_outputs={'default','VxSurface','VySurface','VxShear','VyShear','VxBase','VyBase'}; +md=SetMLHOBC(md); md=solve(md,'Stressbalance'); %Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test335.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test335.py (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test335.py (revision 26555) @@ -12,6 +12,7 @@ md = setmask(md, '', '') md = parameterize(md, '../Par/SquareSheetConstrained.py') md = setflowequation(md, 'MLHO', 'all') +md = SetMLHOBC(md); md.extrude(5, 1.) md.cluster = generic('name', gethostname(), 'np', 3) md.stressbalance.requested_outputs = ['default', 'VxSurface', 'VySurface', 'VxShear', 'VyShear', 'VxBase', 'VyBase'] Index: ../trunk-jpl/test/NightlyRun/test446.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test446.m (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test446.m (revision 26555) @@ -5,6 +5,7 @@ md=setflowequation(md,'MLHO','all'); md.cluster=generic('name',oshostname(),'np',3); md.stressbalance.requested_outputs={'default','VxSurface','VySurface','VxShear','VyShear','VxBase','VyBase'}; +md=SetMLHOBC(md); md=solve(md,'Stressbalance'); %Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test446.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test446.py (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test446.py (revision 26555) @@ -14,6 +14,7 @@ md = setflowequation(md, 'MLHO', 'all') md.cluster = generic('name', gethostname(), 'np', 3) md.stressbalance.requested_outputs = ['default', 'VxSurface', 'VySurface', 'VxShear', 'VyShear', 'VxBase', 'VyBase'] +md = SetMLHOBC(md); md = solve(md, 'Stressbalance') #Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test447.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test447.m (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test447.m (revision 26555) @@ -15,6 +15,7 @@ md=setflowequation(md,'MLHO','all'); md.cluster=generic('name',oshostname(),'np',3); md.transient.requested_outputs={'default','GroundedArea','FloatingArea','TotalFloatingBmb','TotalGroundedBmb','TotalSmb'}; +md=SetMLHOBC(md); md=solve(md,'Transient'); %Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test447.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test447.py (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test447.py (revision 26555) @@ -26,6 +26,7 @@ md = setflowequation(md, 'MLHO', 'all') md.cluster = generic('name', gethostname(), 'np', 3) md.transient.requested_outputs = ['default', 'GroundedArea', 'FloatingArea', 'TotalFloatingBmb', 'TotalGroundedBmb', 'TotalSmb'] +md = SetMLHOBC(md); md = solve(md, 'Transient') #Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test448.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test448.m (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test448.m (revision 26555) @@ -24,6 +24,7 @@ md.transient.isstressbalance=1; md.transient.isgroundingline=1; +md=SetMLHOBC(md); %test different grounding line dynamics. md.groundingline.migration='AggressiveMigration'; md=solve(md,'Transient'); Index: ../trunk-jpl/test/NightlyRun/test448.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test448.py (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test448.py (revision 26555) @@ -37,6 +37,7 @@ #test different grounding line dynamics. md.groundingline.migration = 'AggressiveMigration' +md = SetMLHOBC(md); md = solve(md, 'Transient') element_on_iceshelf_agressive = md.results.TransientSolution[0].MaskOceanLevelset vel_agressive = md.results.TransientSolution[0].Vel Index: ../trunk-jpl/test/NightlyRun/test449.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test449.m (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test449.m (revision 26555) @@ -44,6 +44,7 @@ md.timestepping.time_step=10; md.cluster=generic('name',oshostname(),'np',3); +md=SetMLHOBC(md); md=solve(md,'Transient'); %Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test449.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test449.py (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test449.py (revision 26555) @@ -55,6 +55,7 @@ md.timestepping.time_step = 10 md.cluster = generic('name', gethostname(), 'np', 3) +md = SetMLHOBC(md); md = solve(md, 'Transient') #print md.results.TransientSolution[0].BasalforcingsFloatingiceMeltingRate #print md.results.TransientSolution[1].BasalforcingsFloatingiceMeltingRate Index: ../trunk-jpl/test/NightlyRun/test518.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test518.m (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test518.m (revision 26555) @@ -5,6 +5,7 @@ md=setflowequation(md,'MLHO','all'); md.cluster=generic('name',oshostname(),'np',3); md.stressbalance.requested_outputs={'default','VxSurface','VySurface','VxShear','VyShear','VxBase','VyBase'}; +md=SetMLHOBC(md); md=solve(md,'Stressbalance'); %Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test518.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test518.py (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test518.py (revision 26555) @@ -14,6 +14,7 @@ md = setflowequation(md, 'MLHO', 'all') md.cluster = generic('name', gethostname(), 'np', 3) md.stressbalance.requested_outputs = ['default', 'VxSurface', 'VySurface', 'VxShear', 'VyShear', 'VxBase', 'VyBase'] +md = SetMLHOBC(md); md = solve(md, 'Stressbalance') #Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test519.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test519.m (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test519.m (revision 26555) @@ -8,6 +8,7 @@ [md.mesh.lat,md.mesh.long] = xy2ll(md.mesh.x,md.mesh.y,-1); md.mesh.scale_factor=0.9*ones(md.mesh.numberofvertices,1); md.transient.requested_outputs={'default','IceVolume','IceVolumeScaled','GroundedArea','GroundedAreaScaled','FloatingArea','FloatingAreaScaled','TotalSmb','TotalSmbScaled','TotalFloatingBmb','TotalFloatingBmbScaled'}; +md=SetMLHOBC(md); md=solve(md,'Transient'); %Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test519.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test519.py (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test519.py (revision 26555) @@ -15,6 +15,7 @@ md.mesh.scale_factor = 0.9 * np.ones((md.mesh.numberofvertices)) md.transient.requested_outputs = ['default', 'IceVolume', 'IceVolumeScaled', 'GroundedArea', 'GroundedAreaScaled', 'FloatingArea', 'FloatingAreaScaled', 'TotalSmb', 'TotalSmbScaled', 'TotalFloatingBmb', 'TotalFloatingBmbScaled'] md.cluster = generic('name', gethostname(), 'np', 3) +md = SetMLHOBC(md); md = solve(md, 'Transient') # Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test810.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test810.m (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test810.m (revision 26555) @@ -14,6 +14,7 @@ md.transient.isthermal=0; md.transient.isgroundingline=1; +md=SetMLHOBC(md); md=solve(md,'Transient'); %Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test810.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test810.py (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test810.py (revision 26555) @@ -21,6 +21,7 @@ md.transient.issmb = True md.transient.isthermal = False md.transient.isgroundingline = True +md = SetMLHOBC(md); md = solve(md, 'Transient') Index: ../trunk-jpl/test/NightlyRun/test811.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test811.m (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test811.m (revision 26555) @@ -17,6 +17,7 @@ md.frontalforcings.meltingrate=zeros(md.mesh.numberofvertices,1); md.levelset.migration_max = 1e10; +md=SetMLHOBC(md); md=solve(md,'Transient'); %Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test811.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test811.py (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test811.py (revision 26555) @@ -27,6 +27,7 @@ md.frontalforcings.meltingrate = np.zeros((md.mesh.numberofvertices)) md.levelset.migration_max = 1e10 +md = SetMLHOBC(md); md = solve(md, 'Transient') #Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test812.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test812.m (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test812.m (revision 26555) @@ -34,6 +34,7 @@ md.transient.requested_outputs={'default','StrainRateparallel','StrainRateperpendicular','Calvingratex','Calvingratey','CalvingCalvingrate'}; +md=SetMLHOBC(md); md=solve(md,'Transient'); %Fields and tolerances to track changes Index: ../trunk-jpl/test/NightlyRun/test812.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test812.py (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test812.py (revision 26555) @@ -43,6 +43,7 @@ md.transient.requested_outputs = ['default', 'StrainRateparallel', 'StrainRateperpendicular', 'Calvingratex', 'Calvingratey', 'CalvingCalvingrate'] +md = SetMLHOBC(md); md = solve(md, 'Transient') #Fields and tolerances to track changes Index: ../trunk-jpl/src/m/boundaryconditions/SetMLHOBC.py =================================================================== --- ../trunk-jpl/src/m/boundaryconditions/SetMLHOBC.py (nonexistent) +++ ../trunk-jpl/src/m/boundaryconditions/SetMLHOBC.py (revision 26555) @@ -0,0 +1,21 @@ +import os +import MatlabFuncs as m + + +def SetMLHOBC(md): + """ + SETMLHOBC - Create the boundary conditions for stressbalance for MLHO: VxBase, VyBase, VxShear, VyShear + + Usage: + md = SetIceShelfBC(md, varargin) + + Example: + md = SetIceShelfBC(md) + + """ + + #node on Dirichlet (boundary and ~icefront) + md.stressbalance.spcvx_base = md.stressbalance.spcvx + md.stressbalance.spcvy_base = md.stressbalance.spcvy + md.stressbalance.spcvx_shear = float('nan') * md.stressbalance.spcvx + md.stressbalance.spcvy_shear = float('nan') * md.stressbalance.spcvy Index: ../trunk-jpl/src/m/extrusion/project3d.py =================================================================== --- ../trunk-jpl/src/m/extrusion/project3d.py (revision 26554) +++ ../trunk-jpl/src/m/extrusion/project3d.py (revision 26555) @@ -14,7 +14,7 @@ arguments: 'vector': 2d vector - 'type': 'element' or 'node' + 'type': 'element' or 'node' or 'poly' options: 'layer' a layer number where vector should keep its values. If @@ -22,6 +22,7 @@ vector. 'padding': default to 0 (value adopted by other 3d layers not being projected. + 'degree': degree of polynomials when extrude from bottom to the top Examples: extruded_vector = project3d(md, 'vector', vector2d, 'type', 'node', 'layer', 1, 'padding', NaN) @@ -41,6 +42,7 @@ vectype = options.getfieldvalue('type') #mandatory layer = options.getfieldvalue('layer', 0) #optional (do all layers otherwise) paddingvalue = options.getfieldvalue('padding', 0) #0 by default + polyexponent = options.getfieldvalue('degree', 0) #0 by default, 0-degree polynomial #Handle special case where vector2d is single element (differs from representation in MATLAB) if isinstance(vector2d, (bool, int, float)): @@ -114,6 +116,37 @@ projected_vector[(i * md.mesh.numberofelements2d):((i + 1) * md.mesh.numberofelements2d), :] = vector2d else: projected_vector[((layer - 1) * md.mesh.numberofelements2d):(layer * md.mesh.numberofelements2d), :] = vector2d + elif vectype.lower() == 'poly': + #Initialize 3d vector + if np.ndim(vector2d) == 1: + if vector2d.shape[0] == md.mesh.numberofelements2d: + projected_vector = (paddingvalue * np.ones((md.mesh.numberofelements))).astype(vector2d.dtype) + elif vector2d.shape[0] == md.mesh.numberofelements2d + 1: + projected_vector = (paddingvalue * np.ones((md.mesh.numberofelements + 1))).astype(vector2d.dtype) + projected_vector[-1] = vector2d[-1] + vector2d = vector2d[:-1] + else: + raise TypeError("vector length not supported") + #Fill in + if layer == 0: + for i in range(md.mesh.numberoflayers - 1): + projected_vector[(i * md.mesh.numberofelements2d):((i + 1) * md.mesh.numberofelements2d)] = vector2d*(1.0-(1.0-i/(md.mesh.numberoflayers - 1.0))**polyexponent) + else: + projected_vector[((layer - 1) * md.mesh.numberofelements2d):(layer * md.mesh.numberofelements2d)] = vector2d*(1.0-(1.0-layer/(md.mesh.numberoflayers - 1.0))**polyexponent) + else: + if vector2d.shape[0] == md.mesh.numberofelements2d: + projected_vector = (paddingvalue * np.ones((md.mesh.numberofelements, np.size(vector2d, axis=1)))).astype(vector2d.dtype) + elif vector2d.shape[0] == md.mesh.numberofelements2d + 1: + projected_vector = (paddingvalue * np.ones((md.mesh.numberofelements + 1, np.size(vector2d, axis=1)))).astype(vector2d.dtype) + projected_vector[-1, :] = vector2d[-1, :] + vector2d = vector2d[:-1, :] + else: + raise TypeError("vector length not supported") + #Fill in + if layer == 0: + for i in range(md.mesh.numberoflayers - 1): + projected_vector[(i * md.mesh.numberofelements2d):((i + 1) * md.mesh.numberofelements2d), :] = vector2d*(1.0-(1.0-i/(md.mesh.numberoflayers - 1.0))**polyexponent) + else: else: raise TypeError("project3d error message: unknown projection type") Index: ../trunk-jpl/test/NightlyRun/test127.py =================================================================== --- ../trunk-jpl/test/NightlyRun/test127.py (revision 26554) +++ ../trunk-jpl/test/NightlyRun/test127.py (revision 26555) @@ -24,7 +24,7 @@ massfluxatgate('name', 'MassFlux4', 'profilename', '../Exp/MassFlux4.exp', 'definitionstring', 'Outputdefinition4'), massfluxatgate('name', 'MassFlux5', 'profilename', '../Exp/MassFlux5.exp', 'definitionstring', 'Outputdefinition5'), massfluxatgate('name', 'MassFlux6', 'profilename', '../Exp/MassFlux6.exp', 'definitionstring', 'Outputdefinition6')] - +md = SetMLHOBC(md); md = solve(md, 'Stressbalance') #Fields and tolerances to track changes