Index: /issm/trunk-jpl/externalpackages/matlab/install.sh =================================================================== --- /issm/trunk-jpl/externalpackages/matlab/install.sh (revision 15613) +++ /issm/trunk-jpl/externalpackages/matlab/install.sh (revision 15614) @@ -1,35 +1,3 @@ #!/bin/bash -set -eu - -if [ ! -h install ]; -then - echo " ****************************************************************" - echo " * *" - echo " * you need create a symbolic link from matlab's root directory *" - echo " * *" - echo " * First, you need to locate matlab's root directory *" - echo " * The easiest way to find it is to launch matlab and run the *" - echo " * following command: matlabroot *" - echo " * Then, create a symbolic link using the command 'ln -s ' the *" - echo " * target directory should be named 'install' *" - echo " * *" - echo " * Examples: *" - echo " * ln -s /usr/local/matlab712/ install *" - echo " * ln -s /Applications/MATLAB_R2012a.app/ install *" - echo " * *" - echo " * run ./install.sh again to check that the link is correct *" - echo " * *" - echo " ****************************************************************" -else - - #Read link - DIR=$(readlink install) - if [ -d $DIR ] - then - echo "Symbolic link correctly set" - else - rm -rf install - echo "Incorrect symbolic link:" - echo "provided target ($DIR) does not exist or is not a directory" - fi -fi +rm -rf install +ln -s /usr/local/matlab80/ install Index: /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m =================================================================== --- /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m (revision 15613) +++ /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.m (revision 15614) @@ -27,10 +27,5 @@ end -%segment on neumann (Ice Front) -> none -if (md.mesh.dimension==2) - md.diagnostic.icefront=zeros(0,4); -else - md.diagnostic.icefront=zeros(0,6); -end +%No ice front: do nothing %Create zeros basal melting rate and surface mass balance if not specified Index: /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py =================================================================== --- /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py (revision 15613) +++ /issm/trunk-jpl/src/m/boundaryconditions/SetIceSheetBC.py (revision 15614) @@ -32,9 +32,5 @@ print " boundary conditions for diagnostic model: spc set as zero" - #segment on neumann (Ice Front) -> none - if md.mesh.dimension==2: - md.diagnostic.icefront=numpy.zeros((0,4)) - else: - md.diagnostic.icefront=numpy.zeros((0,6)) + #No ice front -> do nothing #Create zeros basal melting rate and surface mass balance if not specified Index: /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.m =================================================================== --- /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.m (revision 15613) +++ /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.m (revision 15614) @@ -46,22 +46,5 @@ end -%segment on Ice Front -%segment on Neumann (Ice Front) -pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2))); -if (md.mesh.dimension==2) - pressureload=md.mesh.segments(pos,:); -elseif md.mesh.dimension==3 - pressureload_layer1=[md.mesh.segments(pos,1:2) md.mesh.segments(pos,2)+md.mesh.numberofvertices2d md.mesh.segments(pos,1)+md.mesh.numberofvertices2d md.mesh.segments(pos,3)]; - pressureload=[]; - for i=1:md.mesh.numberoflayers-1, - pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ]; - end -end - -%Add water or air enum depending on the element -pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))]; - -%plug onto model -md.diagnostic.icefront=pressureload; +%Ice front position md.mask.icelevelset(find(nodeonicefront))=0; Index: /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py =================================================================== --- /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py (revision 15613) +++ /issm/trunk-jpl/src/m/boundaryconditions/SetIceShelfBC.py (revision 15614) @@ -54,24 +54,5 @@ print " boundary conditions for diagnostic model: spc set as zero" - #segment on Ice Front - #segment on Neumann (Ice Front) -# pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2))); - pos=numpy.nonzero(numpy.logical_or(nodeonicefront[md.mesh.segments[:,0]-1],nodeonicefront[md.mesh.segments[:,1]-1]))[0] - if md.mesh.dimension==2: - pressureload=md.mesh.segments[pos,:] - elif md.mesh.dimension==3: -# pressureload_layer1=[md.mesh.segments(pos,1:2) md.mesh.segments(pos,2)+md.mesh.numberofvertices2d md.mesh.segments(pos,1)+md.mesh.numberofvertices2d md.mesh.segments(pos,3)]; - pressureload_layer1=numpy.hstack((md.mesh.segments[pos,0:2],md.mesh.segments[pos,1]+md.mesh.numberofvertices2d,md.mesh.segments[pos,0]+md.mesh.numberofvertices2d,md.mesh.segments[pos,2])) - pressureload=numpy.zeros((0,5),int) - for i in xrange(1,md.mesh.numberoflayers): -# pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ]; - pressureload=numpy.vstack((pressureload,numpy.hstack((pressureload_layer1[:,0:4]+(i-1)*md.mesh.numberofvertices2d,pressureload_layer1[:,4]+(i-1)*md.mesh.numberofelements2d)))) - - #Add water or air enum depending on the element -# pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))]; - pressureload=numpy.hstack((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1]-1].reshape(-1,1))) - - #plug onto model - md.diagnostic.icefront=pressureload + #Icefront position pos=numpy.nonzero(nodeonicefront)[0] md.mask.icelevelset[pos]=0 Index: /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.m =================================================================== --- /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.m (revision 15613) +++ /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.m (revision 15614) @@ -57,21 +57,5 @@ md.hydrology.spcwatercolumn(pos,1)=1; -%segment on Neumann (Ice Front) -pos=find(vertexonicefront(md.mesh.segments(:,1)) | vertexonicefront(md.mesh.segments(:,2))); -if (md.mesh.dimension==2) - pressureload=md.mesh.segments(pos,:); -elseif md.mesh.dimension==3 - pressureload_layer1=[md.mesh.segments(pos,1:2) md.mesh.segments(pos,2)+md.mesh.numberofvertices2d md.mesh.segments(pos,1)+md.mesh.numberofvertices2d md.mesh.segments(pos,3)]; - pressureload=[]; - for i=1:md.mesh.numberoflayers-1, - pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ]; - end -end - -%Add water or air enum depending on the element -pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))+ 0*md.mask.elementongroundedice(pressureload(:,end))]; - -%plug onto model -md.diagnostic.icefront=pressureload; +%Position of ice front md.mask.icelevelset(find(vertexonicefront))=0; Index: /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py =================================================================== --- /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py (revision 15613) +++ /issm/trunk-jpl/src/m/boundaryconditions/SetMarineIceSheetBC.py (revision 15614) @@ -61,23 +61,5 @@ md.hydrology.spcwatercolumn[pos,0]=1 - #segment on Neumann (Ice Front) -# pos=find(vertexonicefront(md.mesh.segments(:,1)) | vertexonicefront(md.mesh.segments(:,2))); - pos=numpy.nonzero(numpy.logical_or(vertexonicefront[md.mesh.segments[:,0]-1],vertexonicefront[md.mesh.segments[:,1]-1]))[0] - if md.mesh.dimension==2: - pressureload=md.mesh.segments[pos,:] - elif md.mesh.dimension==3: -# pressureload_layer1=[md.mesh.segments(pos,1:2) md.mesh.segments(pos,2)+md.mesh.numberofvertices2d md.mesh.segments(pos,1)+md.mesh.numberofvertices2d md.mesh.segments(pos,3)]; - pressureload_layer1=numpy.hstack((md.mesh.segments[pos,0:2],md.mesh.segments[pos,1]+md.mesh.numberofvertices2d,md.mesh.segments[pos,0]+md.mesh.numberofvertices2d,md.mesh.segments[pos,2])) - pressureload=numpy.zeros((0,5),int) - for i in xrange(1,md.mesh.numberoflayers): -# pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ]; - pressureload=numpy.vstack((pressureload,numpy.hstack((pressureload_layer1[:,0:4]+(i-1)*md.mesh.numberofvertices2d,pressureload_layer1[:,4]+(i-1)*md.mesh.numberofelements2d)))) - - #Add water or air enum depending on the element -# pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))+ 0*md.mask.elementongroundedice(pressureload(:,end))]; - pressureload=numpy.hstack((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1]-1].reshape(-1,1)+0*md.mask.elementongroundedice[pressureload[:,-1]-1].reshape(-1,1))) - - #plug onto model - md.diagnostic.icefront=pressureload + #Position of ice front pos=numpy.nonzero(vertexonicefront)[0] md.mask.icelevelset[pos]=0 Index: sm/trunk-jpl/src/m/classes/diagnostic.m =================================================================== --- /issm/trunk-jpl/src/m/classes/diagnostic.m (revision 15613) +++ (revision ) @@ -1,183 +1,0 @@ -%DIAGNOSTIC class definition -% -% Usage: -% diagnostic=diagnostic(); - -classdef diagnostic - properties (SetAccess=public) - spcvx = NaN; - spcvy = NaN; - spcvz = NaN; - restol = 0; - reltol = 0; - abstol = 0; - isnewton = 0; - FSreconditioning = 0; - viscosity_overshoot = 0; - icefront = NaN; - maxiter = 0; - shelf_dampening = 0; - vertex_pairing = NaN; - penalty_factor = NaN; - rift_penalty_lock = NaN; - rift_penalty_threshold = 0; - referential = NaN; - loadingforce = NaN; - requested_outputs = NaN; - end - methods - function obj = diagnostic(varargin) % {{{ - switch nargin - case 0 - obj=setdefaultparameters(obj); - otherwise - error('constructor not supported'); - end - end % }}} - function obj = setdefaultparameters(obj) % {{{ - - %maximum of non-linear iterations. - obj.maxiter=100; - - %Convergence criterion: absolute, relative and residual - obj.restol=10^-4; - obj.reltol=0.01; - obj.abstol=10; - - obj.FSreconditioning=10^13; - obj.shelf_dampening=0; - - %Penalty factor applied kappa=max(stiffness matrix)*10^penalty_factor - obj.penalty_factor=3; - - %coefficient to update the viscosity between each iteration of - %a diagnostic according to the following formula - %viscosity(n)=viscosity(n)+viscosity_overshoot(viscosity(n)-viscosity(n-1)) - obj.viscosity_overshoot=0; - - %Stop the iterations of rift if below a threshold - obj.rift_penalty_threshold=0; - - %in some solutions, it might be needed to stop a run when only - %a few constraints remain unstable. For thermal computation, this - %parameter is often used. - obj.rift_penalty_lock=10; - - end % }}} - function md = checkconsistency(obj,md,solution,analyses) % {{{ - - %Early return - if ~ismember(DiagnosticHorizAnalysisEnum(),analyses), return; end - %if ~ismember(DiagnosticHorizAnalysisEnum(),analyses) | (solution==TransientSolutionEnum() & md.transient.isdiagnostic==0), return; end - - md = checkfield(md,'diagnostic.spcvx','forcing',1); - md = checkfield(md,'diagnostic.spcvy','forcing',1); - if md.mesh.dimension==3, md = checkfield(md,'diagnostic.spcvz','forcing',1); end - md = checkfield(md,'diagnostic.restol','size',[1 1],'>',0,'NaN',1); - md = checkfield(md,'diagnostic.reltol','size',[1 1]); - md = checkfield(md,'diagnostic.abstol','size',[1 1]); - md = checkfield(md,'diagnostic.isnewton','numel',[1],'values',[0 1 2]); - md = checkfield(md,'diagnostic.FSreconditioning','size',[1 1],'NaN',1); - md = checkfield(md,'diagnostic.viscosity_overshoot','size',[1 1],'NaN',1); - if md.mesh.dimension==2, - md = checkfield(md,'diagnostic.icefront','size',[NaN 4],'NaN',1); - else - md = checkfield(md,'diagnostic.icefront','size',[NaN 6],'NaN',1); - end - md = checkfield(md,'diagnostic.icefront(:,end)','values',[0 1 2]); - md = checkfield(md,'diagnostic.maxiter','size',[1 1],'>=',1); - md = checkfield(md,'diagnostic.referential','size',[md.mesh.numberofvertices 6]); - md = checkfield(md,'diagnostic.loadingforce','size',[md.mesh.numberofvertices 3]); - if ~isempty(md.diagnostic.requested_outputs), - md = checkfield(md,'diagnostic.requested_outputs','size',[NaN 1]); - end - - %singular solution - if ~(any(~isnan(md.diagnostic.spcvx)) & any(~isnan(md.diagnostic.spcvy))), - md = checkmessage(md,['model is not well posed (singular). You need at least one node with fixed velocity!']); - end - %CHECK THAT EACH LINES CONTAINS ONLY NAN VALUES OR NO NAN VALUES - if any(sum(isnan(md.diagnostic.referential),2)~=0 & sum(isnan(md.diagnostic.referential),2)~=6), - md = checkmessage(md,['Each line of diagnostic.referential should contain either only NaN values or no NaN values']); - end - %CHECK THAT THE TWO VECTORS PROVIDED ARE ORTHOGONAL - if any(sum(isnan(md.diagnostic.referential),2)==0), - pos=find(sum(isnan(md.diagnostic.referential),2)==0); - if any(abs(dot(md.diagnostic.referential(pos,1:3),md.diagnostic.referential(pos,4:6),2))>eps), - md = checkmessage(md,['Vectors in diagnostic.referential (columns 1 to 3 and 4 to 6) must be orthogonal']); - end - end - %CHECK THAT NO rotation specified for FS Grounded ice at base - if md.mesh.dimension==3 & md.flowequation.isFS, - pos=find(md.mask.vertexongroundedice & md.mesh.vertexonbed); - if any(~isnan(md.diagnostic.referential(pos,:))), - md = checkmessage(md,['no referential should be specified for basal vertices of grounded ice']); - end - md = checkfield(md,'diagnostic.FSreconditioning','>',0); - end - end % }}} - function disp(obj) % {{{ - - disp(sprintf(' Diagnostic solution parameters:')); - - disp(sprintf('\n %s','Convergence criteria:')); - fielddisplay(obj,'restol','mechanical equilibrium residual convergence criterion'); - fielddisplay(obj,'reltol','velocity relative convergence criterion, NaN: not applied'); - fielddisplay(obj,'abstol','velocity absolute convergence criterion, NaN: not applied'); - fielddisplay(obj,'isnewton','0: Picard''s fixed point, 1: Newton''s method, 2: hybrid'); - fielddisplay(obj,'maxiter','maximum number of nonlinear iterations'); - fielddisplay(obj,'viscosity_overshoot','over-shooting constant new=new+C*(new-old)'); - - disp(sprintf('\n %s','boundary conditions:')); - fielddisplay(obj,'spcvx','x-axis velocity constraint (NaN means no constraint) [m/yr]'); - fielddisplay(obj,'spcvy','y-axis velocity constraint (NaN means no constraint) [m/yr]'); - fielddisplay(obj,'spcvz','z-axis velocity constraint (NaN means no constraint) [m/yr]'); - fielddisplay(obj,'icefront','segments on ice front list (last column 0: Air, 1: Water, 2: Ice)'); - - disp(sprintf('\n %s','Rift options:')); - fielddisplay(obj,'rift_penalty_threshold','threshold for instability of mechanical constraints'); - fielddisplay(obj,'rift_penalty_lock','number of iterations before rift penalties are locked'); - - disp(sprintf('\n %s','Penalty options:')); - fielddisplay(obj,'penalty_factor','offset used by penalties: penalty = Kmax*10^offset'); - fielddisplay(obj,'vertex_pairing','pairs of vertices that are penalized'); - - disp(sprintf('\n %s','Other:')); - fielddisplay(obj,'shelf_dampening','use dampening for floating ice ? Only for FS model'); - fielddisplay(obj,'FSreconditioning','multiplier for incompressibility equation. Only for FS model'); - fielddisplay(obj,'referential','local referential'); - fielddisplay(obj,'loadingforce','loading force applied on each point [N/m^3]'); - fielddisplay(obj,'requested_outputs','additional outputs requested'); - - end % }}} - function marshall(obj,md,fid) % {{{ - WriteData(fid,'object',obj,'fieldname','spcvx','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1); - WriteData(fid,'object',obj,'fieldname','spcvy','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1); - WriteData(fid,'object',obj,'fieldname','spcvz','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1); - WriteData(fid,'object',obj,'fieldname','restol','format','Double'); - WriteData(fid,'object',obj,'fieldname','reltol','format','Double'); - WriteData(fid,'object',obj,'fieldname','abstol','format','Double'); - WriteData(fid,'object',obj,'fieldname','isnewton','format','Integer'); - WriteData(fid,'object',obj,'fieldname','FSreconditioning','format','Double'); - WriteData(fid,'object',obj,'fieldname','viscosity_overshoot','format','Double'); - WriteData(fid,'object',obj,'fieldname','maxiter','format','Integer'); - WriteData(fid,'object',obj,'fieldname','shelf_dampening','format','Integer'); - WriteData(fid,'object',obj,'fieldname','vertex_pairing','format','DoubleMat','mattype',3); - WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double'); - WriteData(fid,'object',obj,'fieldname','rift_penalty_lock','format','Integer'); - WriteData(fid,'object',obj,'fieldname','rift_penalty_threshold','format','Integer'); - WriteData(fid,'object',obj,'fieldname','referential','format','DoubleMat','mattype',1); - WriteData(fid,'data',obj.loadingforce(:,1),'format','DoubleMat','mattype',1,'enum',LoadingforceXEnum); - WriteData(fid,'data',obj.loadingforce(:,2),'format','DoubleMat','mattype',1,'enum',LoadingforceYEnum); - WriteData(fid,'data',obj.loadingforce(:,3),'format','DoubleMat','mattype',1,'enum',LoadingforceZEnum); - WriteData(fid,'object',obj,'fieldname','requested_outputs','format','DoubleMat','mattype',3); - - %marshall ice front - data=obj.icefront; - pos=find(data(:,end)==0); data(pos,end)=AirEnum(); - pos=find(data(:,end)==1); data(pos,end)=WaterEnum(); - pos=find(data(:,end)==2); data(pos,end)=IceEnum(); - WriteData(fid,'data',data,'enum',DiagnosticIcefrontEnum(),'format','DoubleMat','mattype',3); - end % }}} - end -end Index: sm/trunk-jpl/src/m/classes/diagnostic.py =================================================================== --- /issm/trunk-jpl/src/m/classes/diagnostic.py (revision 15613) +++ (revision ) @@ -1,190 +1,0 @@ -import numpy -import sys -import copy -from fielddisplay import fielddisplay -from EnumDefinitions import * -from checkfield import * -from WriteData import * - -class diagnostic(object): - """ - DIAGNOSTIC class definition - - Usage: - diagnostic=diagnostic(); - """ - - def __init__(self): # {{{ - self.spcvx = float('NaN') - self.spcvy = float('NaN') - self.spcvz = float('NaN') - self.restol = 0 - self.reltol = 0 - self.abstol = 0 - self.isnewton = 0 - self.FSreconditioning = 0 - self.viscosity_overshoot = 0 - self.icefront = float('NaN') - self.maxiter = 0 - self.shelf_dampening = 0 - self.vertex_pairing = float('NaN') - self.penalty_factor = float('NaN') - self.rift_penalty_lock = float('NaN') - self.rift_penalty_threshold = 0 - self.referential = float('NaN') - self.loadingforce = float('NaN') - self.requested_outputs = float('NaN') - - #set defaults - self.setdefaultparameters() - - #}}} - def __repr__(self): # {{{ - - string=' Diagnostic solution parameters:' - string="%s\n%s"%(string,' Convergence criteria:') - string="%s\n%s"%(string,fielddisplay(self,'restol','mechanical equilibrium residual convergence criterion')) - string="%s\n%s"%(string,fielddisplay(self,'reltol','velocity relative convergence criterion, NaN: not applied')) - string="%s\n%s"%(string,fielddisplay(self,'abstol','velocity absolute convergence criterion, NaN: not applied')) - string="%s\n%s"%(string,fielddisplay(self,'isnewton',"0: Picard's fixed point, 1: Newton's method, 2: hybrid")) - string="%s\n%s"%(string,fielddisplay(self,'maxiter','maximum number of nonlinear iterations')) - string="%s\n%s"%(string,fielddisplay(self,'viscosity_overshoot','over-shooting constant new=new+C*(new-old)')) - - string="%s\n%s"%(string,'\n boundary conditions:') - string="%s\n%s"%(string,fielddisplay(self,'spcvx','x-axis velocity constraint (NaN means no constraint) [m/yr]')) - string="%s\n%s"%(string,fielddisplay(self,'spcvy','y-axis velocity constraint (NaN means no constraint) [m/yr]')) - string="%s\n%s"%(string,fielddisplay(self,'spcvz','z-axis velocity constraint (NaN means no constraint) [m/yr]')) - string="%s\n%s"%(string,fielddisplay(self,'icefront','segments on ice front list (last column 0: Air, 1: Water, 2: Ice')) - - string="%s\n%s"%(string,'\n Rift options:') - string="%s\n%s"%(string,fielddisplay(self,'rift_penalty_threshold','threshold for instability of mechanical constraints')) - string="%s\n%s"%(string,fielddisplay(self,'rift_penalty_lock','number of iterations before rift penalties are locked')) - - string="%s\n%s"%(string,'\n Penalty options:') - string="%s\n%s"%(string,fielddisplay(self,'penalty_factor','offset used by penalties: penalty = Kmax*10^offset')) - string="%s\n%s"%(string,fielddisplay(self,'vertex_pairing','pairs of vertices that are penalized')) - - string="%s\n%s"%(string,'\n Other:') - string="%s\n%s"%(string,fielddisplay(self,'shelf_dampening','use dampening for floating ice ? Only for FS model')) - string="%s\n%s"%(string,fielddisplay(self,'FSreconditioning','multiplier for incompressibility equation. Only for FS model')) - string="%s\n%s"%(string,fielddisplay(self,'referential','local referential')) - string="%s\n%s"%(string,fielddisplay(self,'loadingforce','loading force applied on each point [N/m^3]')) - string="%s\n%s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested')) - - return string - #}}} - def setdefaultparameters(self): # {{{ - #maximum of non-linear iterations. - self.maxiter=100 - - #Convergence criterion: absolute, relative and residual - self.restol=10**-4 - self.reltol=0.01 - self.abstol=10 - - self.FSreconditioning=10**13 - self.shelf_dampening=0 - - #Penalty factor applied kappa=max(stiffness matrix)*10^penalty_factor - self.penalty_factor=3 - - #coefficient to update the viscosity between each iteration of - #a diagnostic according to the following formula - #viscosity(n)=viscosity(n)+viscosity_overshoot(viscosity(n)-viscosity(n-1)) - self.viscosity_overshoot=0 - - #Stop the iterations of rift if below a threshold - self.rift_penalty_threshold=0 - - #in some solutions, it might be needed to stop a run when only - #a few constraints remain unstable. For thermal computation, this - #parameter is often used. - self.rift_penalty_lock=10 - - return self - #}}} - def checkconsistency(self,md,solution,analyses): # {{{ - - #Early return - if DiagnosticHorizAnalysisEnum() not in analyses: - return md - #if (DiagnosticHorizAnalysisEnum() not in analyses) | (solution==TransientSolutionEnum() and not md.transient.isdiagnostic): - # return md - - md = checkfield(md,'diagnostic.spcvx','forcing',1) - md = checkfield(md,'diagnostic.spcvy','forcing',1) - if md.mesh.dimension==3: - md = checkfield(md,'diagnostic.spcvz','forcing',1) - md = checkfield(md,'diagnostic.restol','size',[1],'>',0) - md = checkfield(md,'diagnostic.reltol','size',[1]) - md = checkfield(md,'diagnostic.abstol','size',[1]) - md = checkfield(md,'diagnostic.isnewton','numel',[1],'values',[0,1,2]) - md = checkfield(md,'diagnostic.FSreconditioning','size',[1],'NaN',1) - md = checkfield(md,'diagnostic.viscosity_overshoot','size',[1],'NaN',1) - if md.mesh.dimension==2: - md = checkfield(md,'diagnostic.icefront','size',[float('NaN'),4],'NaN',1) - else: - md = checkfield(md,'diagnostic.icefront','size',[float('NaN'),6],'NaN',1) - md = checkfield(md,'diagnostic.icefront[:,-1]','values',[0,1,2]) - md = checkfield(md,'diagnostic.maxiter','size',[1],'>=',1) - md = checkfield(md,'diagnostic.referential','size',[md.mesh.numberofvertices,6]) - md = checkfield(md,'diagnostic.loadingforce','size',[md.mesh.numberofvertices,3]) - if not md.diagnostic.requested_outputs: - md = checkfield(md,'diagnostic.requested_outputs','size',[float('NaN'),1]) - - #singular solution -# if ~any((~isnan(md.diagnostic.spcvx)+~isnan(md.diagnostic.spcvy))==2), - if not numpy.any(numpy.logical_and(numpy.logical_not(numpy.isnan(md.diagnostic.spcvx)),numpy.logical_not(numpy.isnan(md.diagnostic.spcvy)))): - md.checkmessage("model is not well posed (singular). You need at least one node with fixed velocity!") - #CHECK THAT EACH LINES CONTAINS ONLY NAN VALUES OR NO NAN VALUES -# if any(sum(isnan(md.diagnostic.referential),2)~=0 & sum(isnan(md.diagnostic.referential),2)~=6), - if numpy.any(numpy.logical_and(numpy.sum(numpy.isnan(md.diagnostic.referential),axis=1)!=0,numpy.sum(numpy.isnan(md.diagnostic.referential),axis=1)!=6)): - md.checkmessage("Each line of diagnostic.referential should contain either only NaN values or no NaN values") - #CHECK THAT THE TWO VECTORS PROVIDED ARE ORTHOGONAL -# if any(sum(isnan(md.diagnostic.referential),2)==0), - if numpy.any(numpy.sum(numpy.isnan(md.diagnostic.referential),axis=1)==0): - pos=[i for i,item in enumerate(numpy.sum(numpy.isnan(md.diagnostic.referential),axis=1)) if item==0] -# numpy.inner (and numpy.dot) calculate all the dot product permutations, resulting in a full matrix multiply -# if numpy.any(numpy.abs(numpy.inner(md.diagnostic.referential[pos,0:2],md.diagnostic.referential[pos,3:5]).diagonal())>sys.float_info.epsilon): -# md.checkmessage("Vectors in diagnostic.referential (columns 1 to 3 and 4 to 6) must be orthogonal") - for item in md.diagnostic.referential[pos,:]: - if numpy.abs(numpy.inner(item[0:2],item[3:5]))>sys.float_info.epsilon: - md.checkmessage("Vectors in diagnostic.referential (columns 1 to 3 and 4 to 6) must be orthogonal") - #CHECK THAT NO rotation specified for FS Grounded ice at base -# if md.mesh.dimension==3 & md.flowequation.isFS, - if md.mesh.dimension==3 and md.flowequation.isFS: - pos=numpy.nonzero(numpy.logical_and(md.mask.vertexongroundedice,md.mesh.vertexonbed)) - if numpy.any(numpy.logical_not(numpy.isnan(md.diagnostic.referential[pos,:]))): - md.checkmessage("no referential should be specified for basal vertices of grounded ice") - - return md - # }}} - def marshall(self,md,fid): # {{{ - WriteData(fid,'object',self,'fieldname','spcvx','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1) - WriteData(fid,'object',self,'fieldname','spcvy','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1) - WriteData(fid,'object',self,'fieldname','spcvz','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1) - WriteData(fid,'object',self,'fieldname','restol','format','Double') - WriteData(fid,'object',self,'fieldname','reltol','format','Double') - WriteData(fid,'object',self,'fieldname','abstol','format','Double') - WriteData(fid,'object',self,'fieldname','isnewton','format','Integer') - WriteData(fid,'object',self,'fieldname','FSreconditioning','format','Double') - WriteData(fid,'object',self,'fieldname','viscosity_overshoot','format','Double') - WriteData(fid,'object',self,'fieldname','maxiter','format','Integer') - WriteData(fid,'object',self,'fieldname','shelf_dampening','format','Integer') - WriteData(fid,'object',self,'fieldname','vertex_pairing','format','DoubleMat','mattype',3) - WriteData(fid,'object',self,'fieldname','penalty_factor','format','Double') - WriteData(fid,'object',self,'fieldname','rift_penalty_lock','format','Integer') - WriteData(fid,'object',self,'fieldname','rift_penalty_threshold','format','Integer') - WriteData(fid,'object',self,'fieldname','referential','format','DoubleMat','mattype',1) - WriteData(fid,'data',self.loadingforce[:,0],'format','DoubleMat','mattype',1,'enum',LoadingforceXEnum()) - WriteData(fid,'data',self.loadingforce[:,1],'format','DoubleMat','mattype',1,'enum',LoadingforceYEnum()) - WriteData(fid,'data',self.loadingforce[:,2],'format','DoubleMat','mattype',1,'enum',LoadingforceZEnum()) - WriteData(fid,'object',self,'fieldname','requested_outputs','format','DoubleMat','mattype',3) - - #marshall ice front - data=copy.deepcopy(self.icefront) - data[numpy.nonzero(data[:,-1]==0),-1]=AirEnum() - data[numpy.nonzero(data[:,-1]==1),-1]=WaterEnum() - data[numpy.nonzero(data[:,-1]==2),-1]=IceEnum() - WriteData(fid,'data',data,'enum',DiagnosticIcefrontEnum(),'format','DoubleMat','mattype',3) - # }}} Index: /issm/trunk-jpl/src/m/classes/model/model.m =================================================================== --- /issm/trunk-jpl/src/m/classes/model/model.m (revision 15613) +++ /issm/trunk-jpl/src/m/classes/model/model.m (revision 15614) @@ -91,4 +91,14 @@ disp('Recovering old hydrology class'); md.hydrology=hydrologyshreve(md.materials); + end + %2013 Juli 25th + if isa(md.diagnostic,'diagnostic'), + disp('Recovering old stressbalance class'); + icefront = md.diagnostic.icefront; + md.diagnostic=stressbalance(md.diagnostic); + + %Deal with front + md.mask.icelevelset=ones(md.mesh.numberofvertices,1); + md.mask.icelevelset(icefront(:,1:end-2))=0; end @@ -184,10 +194,4 @@ md.prognostic.spcthickness=project2d(md,md.prognostic.spcthickness,md.mesh.numberoflayers); md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers); - - %Extrusion of Neumann BC - if ~isnan(md.diagnostic.icefront), - numberofneumann2d=size(md.diagnostic.icefront,1)/(md.mesh.numberoflayers-1); - md.diagnostic.icefront=[md.diagnostic.icefront(1:numberofneumann2d,1:2) md.diagnostic.icefront(1:numberofneumann2d,5:6)]; %Add two columns on the first layer - end %materials @@ -499,16 +503,4 @@ end - %Diagnostic - if ~isnan(md2.diagnostic.icefront) - md2.diagnostic.icefront(:,1)=Pnode(md1.diagnostic.icefront(:,1)); - md2.diagnostic.icefront(:,2)=Pnode(md1.diagnostic.icefront(:,2)); - md2.diagnostic.icefront(:,end-1)=Pelem(md1.diagnostic.icefront(:,end-1)); - if md1.mesh.dimension==3 - md2.diagnostic.icefront(:,3)=Pnode(md1.diagnostic.icefront(:,3)); - md2.diagnostic.icefront(:,4)=Pnode(md1.diagnostic.icefront(:,4)); - end - md2.diagnostic.icefront=md2.diagnostic.icefront(find(md2.diagnostic.icefront(:,1) & md2.diagnostic.icefront(:,2) & md2.diagnostic.icefront(:,end)),:); - end - %Results fields if isstruct(md1.results), @@ -732,12 +724,4 @@ md.diagnostic.referential=project3d(md,'vector',md.diagnostic.referential,'type','node'); md.diagnostic.loadingforce=project3d(md,'vector',md.diagnostic.loadingforce,'type','node'); - - %in 3d, pressureload: [node1 node2 node3 node4 element] - pressureload_layer1=[md.diagnostic.icefront(:,1:2) md.diagnostic.icefront(:,2)+md.mesh.numberofvertices2d md.diagnostic.icefront(:,1)+md.mesh.numberofvertices2d md.diagnostic.icefront(:,3:4)]; %Add two columns on the first layer - pressureload=[]; - for i=1:numlayers-1, - pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d pressureload_layer1(:,6)]; - end - md.diagnostic.icefront=pressureload; %connectivity @@ -1099,5 +1083,5 @@ md.cluster = generic(); md.balancethickness = balancethickness(); - md.diagnostic = diagnostic(); + md.diagnostic = stressbalance(); md.hydrology = hydrologyshreve(); md.prognostic = prognostic(); Index: /issm/trunk-jpl/src/m/classes/model/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model/model.py (revision 15613) +++ /issm/trunk-jpl/src/m/classes/model/model.py (revision 15614) @@ -72,5 +72,5 @@ self.balancethickness = balancethickness() - self.diagnostic = diagnostic() + self.diagnostic = stressbalance() self.groundingline = groundingline() self.hydrology = hydrologyshreve() @@ -392,14 +392,4 @@ md2.thermal.spctemperature[nodestoflag2,0]=1 - #Diagnostic - if numpy.any(numpy.logical_not(numpy.isnan(md2.diagnostic.icefront))): - md2.diagnostic.icefront[:,0]=Pnode[md1.diagnostic.icefront[:,0].astype(int)-1] - md2.diagnostic.icefront[:,1]=Pnode[md1.diagnostic.icefront[:,1].astype(int)-1] - md2.diagnostic.icefront[:,-2]=Pelem[md1.diagnostic.icefront[:,-2].astype(int)-1] - if md1.mesh.dimension==3: - md2.diagnostic.icefront[:,2]=Pnode[md1.diagnostic.icefront[:,2].astype(int)-1] - md2.diagnostic.icefront[:,3]=Pnode[md1.diagnostic.icefront[:,3].astype(int)-1] - md2.diagnostic.icefront=md2.diagnostic.icefront[numpy.nonzero(logical_and_n(md2.diagnostic.icefront[:,0],md2.diagnostic.icefront[:,1],md2.diagnostic.icefront[:,-1]))[0],:] - #Results fields if md1.results: @@ -627,11 +617,4 @@ md.diagnostic.loadingforce=project3d(md,'vector',md.diagnostic.loadingforce,'type','node') - #in 3d, pressureload: [node1 node2 node3 node4 element] - pressureload_layer1=numpy.hstack((md.diagnostic.icefront[:,0:2],md.diagnostic.icefront[:,1:2]+md.mesh.numberofvertices2d,md.diagnostic.icefront[:,0:1]+md.mesh.numberofvertices2d,md.diagnostic.icefront[:,2:4])) #Add two columns on the first layer - pressureload=numpy.empty((0,6),int) - for i in xrange(numlayers-1): - pressureload=numpy.vstack((pressureload,numpy.hstack((pressureload_layer1[:,0:4]+i*md.mesh.numberofvertices2d,pressureload_layer1[:,4:5]+i*md.mesh.numberofelements2d,pressureload_layer1[:,5:6])))) - md.diagnostic.icefront=pressureload - #connectivity md.mesh.elementconnectivity=numpy.tile(md.mesh.elementconnectivity,(numlayers-1,1)) Index: /issm/trunk-jpl/src/m/classes/oldclasses/diagnostic.m =================================================================== --- /issm/trunk-jpl/src/m/classes/oldclasses/diagnostic.m (revision 15614) +++ /issm/trunk-jpl/src/m/classes/oldclasses/diagnostic.m (revision 15614) @@ -0,0 +1,183 @@ +%DIAGNOSTIC class definition +% +% Usage: +% diagnostic=diagnostic(); + +classdef diagnostic + properties (SetAccess=public) + spcvx = NaN; + spcvy = NaN; + spcvz = NaN; + restol = 0; + reltol = 0; + abstol = 0; + isnewton = 0; + FSreconditioning = 0; + viscosity_overshoot = 0; + icefront = NaN; + maxiter = 0; + shelf_dampening = 0; + vertex_pairing = NaN; + penalty_factor = NaN; + rift_penalty_lock = NaN; + rift_penalty_threshold = 0; + referential = NaN; + loadingforce = NaN; + requested_outputs = NaN; + end + methods + function obj = diagnostic(varargin) % {{{ + switch nargin + case 0 + obj=setdefaultparameters(obj); + otherwise + error('constructor not supported'); + end + end % }}} + function obj = setdefaultparameters(obj) % {{{ + + %maximum of non-linear iterations. + obj.maxiter=100; + + %Convergence criterion: absolute, relative and residual + obj.restol=10^-4; + obj.reltol=0.01; + obj.abstol=10; + + obj.FSreconditioning=10^13; + obj.shelf_dampening=0; + + %Penalty factor applied kappa=max(stiffness matrix)*10^penalty_factor + obj.penalty_factor=3; + + %coefficient to update the viscosity between each iteration of + %a diagnostic according to the following formula + %viscosity(n)=viscosity(n)+viscosity_overshoot(viscosity(n)-viscosity(n-1)) + obj.viscosity_overshoot=0; + + %Stop the iterations of rift if below a threshold + obj.rift_penalty_threshold=0; + + %in some solutions, it might be needed to stop a run when only + %a few constraints remain unstable. For thermal computation, this + %parameter is often used. + obj.rift_penalty_lock=10; + + end % }}} + function md = checkconsistency(obj,md,solution,analyses) % {{{ + + %Early return + if ~ismember(DiagnosticHorizAnalysisEnum(),analyses), return; end + %if ~ismember(DiagnosticHorizAnalysisEnum(),analyses) | (solution==TransientSolutionEnum() & md.transient.isdiagnostic==0), return; end + + md = checkfield(md,'diagnostic.spcvx','forcing',1); + md = checkfield(md,'diagnostic.spcvy','forcing',1); + if md.mesh.dimension==3, md = checkfield(md,'diagnostic.spcvz','forcing',1); end + md = checkfield(md,'diagnostic.restol','size',[1 1],'>',0,'NaN',1); + md = checkfield(md,'diagnostic.reltol','size',[1 1]); + md = checkfield(md,'diagnostic.abstol','size',[1 1]); + md = checkfield(md,'diagnostic.isnewton','numel',[1],'values',[0 1 2]); + md = checkfield(md,'diagnostic.FSreconditioning','size',[1 1],'NaN',1); + md = checkfield(md,'diagnostic.viscosity_overshoot','size',[1 1],'NaN',1); + if md.mesh.dimension==2, + md = checkfield(md,'diagnostic.icefront','size',[NaN 4],'NaN',1); + else + md = checkfield(md,'diagnostic.icefront','size',[NaN 6],'NaN',1); + end + md = checkfield(md,'diagnostic.icefront(:,end)','values',[0 1 2]); + md = checkfield(md,'diagnostic.maxiter','size',[1 1],'>=',1); + md = checkfield(md,'diagnostic.referential','size',[md.mesh.numberofvertices 6]); + md = checkfield(md,'diagnostic.loadingforce','size',[md.mesh.numberofvertices 3]); + if ~isempty(md.diagnostic.requested_outputs), + md = checkfield(md,'diagnostic.requested_outputs','size',[NaN 1]); + end + + %singular solution + if ~(any(~isnan(md.diagnostic.spcvx)) & any(~isnan(md.diagnostic.spcvy))), + md = checkmessage(md,['model is not well posed (singular). You need at least one node with fixed velocity!']); + end + %CHECK THAT EACH LINES CONTAINS ONLY NAN VALUES OR NO NAN VALUES + if any(sum(isnan(md.diagnostic.referential),2)~=0 & sum(isnan(md.diagnostic.referential),2)~=6), + md = checkmessage(md,['Each line of diagnostic.referential should contain either only NaN values or no NaN values']); + end + %CHECK THAT THE TWO VECTORS PROVIDED ARE ORTHOGONAL + if any(sum(isnan(md.diagnostic.referential),2)==0), + pos=find(sum(isnan(md.diagnostic.referential),2)==0); + if any(abs(dot(md.diagnostic.referential(pos,1:3),md.diagnostic.referential(pos,4:6),2))>eps), + md = checkmessage(md,['Vectors in diagnostic.referential (columns 1 to 3 and 4 to 6) must be orthogonal']); + end + end + %CHECK THAT NO rotation specified for FS Grounded ice at base + if md.mesh.dimension==3 & md.flowequation.isFS, + pos=find(md.mask.vertexongroundedice & md.mesh.vertexonbed); + if any(~isnan(md.diagnostic.referential(pos,:))), + md = checkmessage(md,['no referential should be specified for basal vertices of grounded ice']); + end + md = checkfield(md,'diagnostic.FSreconditioning','>',0); + end + end % }}} + function disp(obj) % {{{ + + disp(sprintf(' Diagnostic solution parameters:')); + + disp(sprintf('\n %s','Convergence criteria:')); + fielddisplay(obj,'restol','mechanical equilibrium residual convergence criterion'); + fielddisplay(obj,'reltol','velocity relative convergence criterion, NaN: not applied'); + fielddisplay(obj,'abstol','velocity absolute convergence criterion, NaN: not applied'); + fielddisplay(obj,'isnewton','0: Picard''s fixed point, 1: Newton''s method, 2: hybrid'); + fielddisplay(obj,'maxiter','maximum number of nonlinear iterations'); + fielddisplay(obj,'viscosity_overshoot','over-shooting constant new=new+C*(new-old)'); + + disp(sprintf('\n %s','boundary conditions:')); + fielddisplay(obj,'spcvx','x-axis velocity constraint (NaN means no constraint) [m/yr]'); + fielddisplay(obj,'spcvy','y-axis velocity constraint (NaN means no constraint) [m/yr]'); + fielddisplay(obj,'spcvz','z-axis velocity constraint (NaN means no constraint) [m/yr]'); + fielddisplay(obj,'icefront','segments on ice front list (last column 0: Air, 1: Water, 2: Ice)'); + + disp(sprintf('\n %s','Rift options:')); + fielddisplay(obj,'rift_penalty_threshold','threshold for instability of mechanical constraints'); + fielddisplay(obj,'rift_penalty_lock','number of iterations before rift penalties are locked'); + + disp(sprintf('\n %s','Penalty options:')); + fielddisplay(obj,'penalty_factor','offset used by penalties: penalty = Kmax*10^offset'); + fielddisplay(obj,'vertex_pairing','pairs of vertices that are penalized'); + + disp(sprintf('\n %s','Other:')); + fielddisplay(obj,'shelf_dampening','use dampening for floating ice ? Only for FS model'); + fielddisplay(obj,'FSreconditioning','multiplier for incompressibility equation. Only for FS model'); + fielddisplay(obj,'referential','local referential'); + fielddisplay(obj,'loadingforce','loading force applied on each point [N/m^3]'); + fielddisplay(obj,'requested_outputs','additional outputs requested'); + + end % }}} + function marshall(obj,md,fid) % {{{ + WriteData(fid,'object',obj,'fieldname','spcvx','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1); + WriteData(fid,'object',obj,'fieldname','spcvy','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1); + WriteData(fid,'object',obj,'fieldname','spcvz','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1); + WriteData(fid,'object',obj,'fieldname','restol','format','Double'); + WriteData(fid,'object',obj,'fieldname','reltol','format','Double'); + WriteData(fid,'object',obj,'fieldname','abstol','format','Double'); + WriteData(fid,'object',obj,'fieldname','isnewton','format','Integer'); + WriteData(fid,'object',obj,'fieldname','FSreconditioning','format','Double'); + WriteData(fid,'object',obj,'fieldname','viscosity_overshoot','format','Double'); + WriteData(fid,'object',obj,'fieldname','maxiter','format','Integer'); + WriteData(fid,'object',obj,'fieldname','shelf_dampening','format','Integer'); + WriteData(fid,'object',obj,'fieldname','vertex_pairing','format','DoubleMat','mattype',3); + WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double'); + WriteData(fid,'object',obj,'fieldname','rift_penalty_lock','format','Integer'); + WriteData(fid,'object',obj,'fieldname','rift_penalty_threshold','format','Integer'); + WriteData(fid,'object',obj,'fieldname','referential','format','DoubleMat','mattype',1); + WriteData(fid,'data',obj.loadingforce(:,1),'format','DoubleMat','mattype',1,'enum',LoadingforceXEnum); + WriteData(fid,'data',obj.loadingforce(:,2),'format','DoubleMat','mattype',1,'enum',LoadingforceYEnum); + WriteData(fid,'data',obj.loadingforce(:,3),'format','DoubleMat','mattype',1,'enum',LoadingforceZEnum); + WriteData(fid,'object',obj,'fieldname','requested_outputs','format','DoubleMat','mattype',3); + + %marshall ice front + data=obj.icefront; + pos=find(data(:,end)==0); data(pos,end)=AirEnum(); + pos=find(data(:,end)==1); data(pos,end)=WaterEnum(); + pos=find(data(:,end)==2); data(pos,end)=IceEnum(); + WriteData(fid,'data',data,'enum',DiagnosticIcefrontEnum(),'format','DoubleMat','mattype',3); + end % }}} + end +end Index: /issm/trunk-jpl/src/m/classes/stressbalance.m =================================================================== --- /issm/trunk-jpl/src/m/classes/stressbalance.m (revision 15614) +++ /issm/trunk-jpl/src/m/classes/stressbalance.m (revision 15614) @@ -0,0 +1,177 @@ +%STRESSBALANC Eclass definition +% +% Usage: +% stressbalance=stressbalance(); + +classdef stressbalance + properties (SetAccess=public) + spcvx = NaN; + spcvy = NaN; + spcvz = NaN; + restol = 0; + reltol = 0; + abstol = 0; + isnewton = 0; + FSreconditioning = 0; + viscosity_overshoot = 0; + maxiter = 0; + shelf_dampening = 0; + vertex_pairing = NaN; + penalty_factor = NaN; + rift_penalty_lock = NaN; + rift_penalty_threshold = 0; + referential = NaN; + loadingforce = NaN; + requested_outputs = NaN; + end + methods + function obj = stressbalance(varargin) % {{{ + switch nargin + case 0 + obj=setdefaultparameters(obj); + case 1 + inputstruct=varargin{1}; + list1 = properties('stressbalance'); + list2 = fieldnames(inputstruct); + for i=1:length(list1) + fieldname = list1{i}; + if ismember(fieldname,list2), + obj.(fieldname) = inputstruct.(fieldname); + end + end + otherwise + error('constructor not supported'); + end + end % }}} + function obj = setdefaultparameters(obj) % {{{ + + %maximum of non-linear iterations. + obj.maxiter=100; + + %Convergence criterion: absolute, relative and residual + obj.restol=10^-4; + obj.reltol=0.01; + obj.abstol=10; + + obj.FSreconditioning=10^13; + obj.shelf_dampening=0; + + %Penalty factor applied kappa=max(stiffness matrix)*10^penalty_factor + obj.penalty_factor=3; + + %coefficient to update the viscosity between each iteration of + %a stressbalance according to the following formula + %viscosity(n)=viscosity(n)+viscosity_overshoot(viscosity(n)-viscosity(n-1)) + obj.viscosity_overshoot=0; + + %Stop the iterations of rift if below a threshold + obj.rift_penalty_threshold=0; + + %in some solutions, it might be needed to stop a run when only + %a few constraints remain unstable. For thermal computation, this + %parameter is often used. + obj.rift_penalty_lock=10; + + end % }}} + function md = checkconsistency(obj,md,solution,analyses) % {{{ + + %Early return + if ~ismember(DiagnosticHorizAnalysisEnum(),analyses), return; end + + md = checkfield(md,'diagnostic.spcvx','forcing',1); + md = checkfield(md,'diagnostic.spcvy','forcing',1); + if md.mesh.dimension==3, md = checkfield(md,'diagnostic.spcvz','forcing',1); end + md = checkfield(md,'diagnostic.restol','size',[1 1],'>',0,'NaN',1); + md = checkfield(md,'diagnostic.reltol','size',[1 1]); + md = checkfield(md,'diagnostic.abstol','size',[1 1]); + md = checkfield(md,'diagnostic.isnewton','numel',[1],'values',[0 1 2]); + md = checkfield(md,'diagnostic.FSreconditioning','size',[1 1],'NaN',1); + md = checkfield(md,'diagnostic.viscosity_overshoot','size',[1 1],'NaN',1); + md = checkfield(md,'diagnostic.maxiter','size',[1 1],'>=',1); + md = checkfield(md,'diagnostic.referential','size',[md.mesh.numberofvertices 6]); + md = checkfield(md,'diagnostic.loadingforce','size',[md.mesh.numberofvertices 3]); + if ~isempty(md.diagnostic.requested_outputs), + md = checkfield(md,'diagnostic.requested_outputs','size',[NaN 1]); + end + + %singular solution + if ~(any(~isnan(md.diagnostic.spcvx)) & any(~isnan(md.diagnostic.spcvy))), + md = checkmessage(md,['model is not well posed (singular). You need at least one node with fixed velocity!']); + end + %CHECK THAT EACH LINES CONTAINS ONLY NAN VALUES OR NO NAN VALUES + if any(sum(isnan(md.diagnostic.referential),2)~=0 & sum(isnan(md.diagnostic.referential),2)~=6), + md = checkmessage(md,['Each line of diagnostic.referential should contain either only NaN values or no NaN values']); + end + %CHECK THAT THE TWO VECTORS PROVIDED ARE ORTHOGONAL + if any(sum(isnan(md.diagnostic.referential),2)==0), + pos=find(sum(isnan(md.diagnostic.referential),2)==0); + if any(abs(dot(md.diagnostic.referential(pos,1:3),md.diagnostic.referential(pos,4:6),2))>eps), + md = checkmessage(md,['Vectors in diagnostic.referential (columns 1 to 3 and 4 to 6) must be orthogonal']); + end + end + %CHECK THAT NO rotation specified for FS Grounded ice at base + if md.mesh.dimension==3 & md.flowequation.isFS, + pos=find(md.mask.vertexongroundedice & md.mesh.vertexonbed); + if any(~isnan(md.diagnostic.referential(pos,:))), + md = checkmessage(md,['no referential should be specified for basal vertices of grounded ice']); + end + md = checkfield(md,'diagnostic.FSreconditioning','>',0); + end + end % }}} + function disp(obj) % {{{ + + disp(sprintf(' StressBalance solution parameters:')); + + disp(sprintf('\n %s','Convergence criteria:')); + fielddisplay(obj,'restol','mechanical equilibrium residual convergence criterion'); + fielddisplay(obj,'reltol','velocity relative convergence criterion, NaN: not applied'); + fielddisplay(obj,'abstol','velocity absolute convergence criterion, NaN: not applied'); + fielddisplay(obj,'isnewton','0: Picard''s fixed point, 1: Newton''s method, 2: hybrid'); + fielddisplay(obj,'maxiter','maximum number of nonlinear iterations'); + fielddisplay(obj,'viscosity_overshoot','over-shooting constant new=new+C*(new-old)'); + + disp(sprintf('\n %s','boundary conditions:')); + fielddisplay(obj,'spcvx','x-axis velocity constraint (NaN means no constraint) [m/yr]'); + fielddisplay(obj,'spcvy','y-axis velocity constraint (NaN means no constraint) [m/yr]'); + fielddisplay(obj,'spcvz','z-axis velocity constraint (NaN means no constraint) [m/yr]'); + + disp(sprintf('\n %s','Rift options:')); + fielddisplay(obj,'rift_penalty_threshold','threshold for instability of mechanical constraints'); + fielddisplay(obj,'rift_penalty_lock','number of iterations before rift penalties are locked'); + + disp(sprintf('\n %s','Penalty options:')); + fielddisplay(obj,'penalty_factor','offset used by penalties: penalty = Kmax*10^offset'); + fielddisplay(obj,'vertex_pairing','pairs of vertices that are penalized'); + + disp(sprintf('\n %s','Other:')); + fielddisplay(obj,'shelf_dampening','use dampening for floating ice ? Only for FS model'); + fielddisplay(obj,'FSreconditioning','multiplier for incompressibility equation. Only for FS model'); + fielddisplay(obj,'referential','local referential'); + fielddisplay(obj,'loadingforce','loading force applied on each point [N/m^3]'); + fielddisplay(obj,'requested_outputs','additional outputs requested'); + + end % }}} + function marshall(obj,md,fid) % {{{ + WriteData(fid,'object',obj,'class','diagnostic','fieldname','spcvx','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1); + WriteData(fid,'object',obj,'class','diagnostic','fieldname','spcvy','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1); + WriteData(fid,'object',obj,'class','diagnostic','fieldname','spcvz','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1); + WriteData(fid,'object',obj,'class','diagnostic','fieldname','restol','format','Double'); + WriteData(fid,'object',obj,'class','diagnostic','fieldname','reltol','format','Double'); + WriteData(fid,'object',obj,'class','diagnostic','fieldname','abstol','format','Double'); + WriteData(fid,'object',obj,'class','diagnostic','fieldname','isnewton','format','Integer'); + WriteData(fid,'object',obj,'class','diagnostic','fieldname','FSreconditioning','format','Double'); + WriteData(fid,'object',obj,'class','diagnostic','fieldname','viscosity_overshoot','format','Double'); + WriteData(fid,'object',obj,'class','diagnostic','fieldname','maxiter','format','Integer'); + WriteData(fid,'object',obj,'class','diagnostic','fieldname','shelf_dampening','format','Integer'); + WriteData(fid,'object',obj,'class','diagnostic','fieldname','vertex_pairing','format','DoubleMat','mattype',3); + WriteData(fid,'object',obj,'class','diagnostic','fieldname','penalty_factor','format','Double'); + WriteData(fid,'object',obj,'class','diagnostic','fieldname','rift_penalty_lock','format','Integer'); + WriteData(fid,'object',obj,'class','diagnostic','fieldname','rift_penalty_threshold','format','Integer'); + WriteData(fid,'object',obj,'class','diagnostic','fieldname','referential','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'class','diagnostic','fieldname','requested_outputs','format','DoubleMat','mattype',3); + WriteData(fid,'data',obj.loadingforce(:,1),'format','DoubleMat','mattype',1,'enum',LoadingforceXEnum); + WriteData(fid,'data',obj.loadingforce(:,2),'format','DoubleMat','mattype',1,'enum',LoadingforceYEnum); + WriteData(fid,'data',obj.loadingforce(:,3),'format','DoubleMat','mattype',1,'enum',LoadingforceZEnum); + end % }}} + end +end Index: /issm/trunk-jpl/src/m/classes/stressbalance.py =================================================================== --- /issm/trunk-jpl/src/m/classes/stressbalance.py (revision 15614) +++ /issm/trunk-jpl/src/m/classes/stressbalance.py (revision 15614) @@ -0,0 +1,181 @@ +import numpy +import sys +import copy +from fielddisplay import fielddisplay +from EnumDefinitions import * +from checkfield import * +from WriteData import * + +class stressbalance(object): + """ + STRESSBALANCE class definition + + Usage: + stressbalance=stressbalance(); + """ + + def __init__(self): # {{{ + self.spcvx = float('NaN') + self.spcvy = float('NaN') + self.spcvz = float('NaN') + self.restol = 0 + self.reltol = 0 + self.abstol = 0 + self.isnewton = 0 + self.FSreconditioning = 0 + self.viscosity_overshoot = 0 + self.icefront = float('NaN') + self.maxiter = 0 + self.shelf_dampening = 0 + self.vertex_pairing = float('NaN') + self.penalty_factor = float('NaN') + self.rift_penalty_lock = float('NaN') + self.rift_penalty_threshold = 0 + self.referential = float('NaN') + self.loadingforce = float('NaN') + self.requested_outputs = float('NaN') + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(self): # {{{ + + string=' StressBalance solution parameters:' + string="%s\n%s"%(string,' Convergence criteria:') + string="%s\n%s"%(string,fielddisplay(self,'restol','mechanical equilibrium residual convergence criterion')) + string="%s\n%s"%(string,fielddisplay(self,'reltol','velocity relative convergence criterion, NaN: not applied')) + string="%s\n%s"%(string,fielddisplay(self,'abstol','velocity absolute convergence criterion, NaN: not applied')) + string="%s\n%s"%(string,fielddisplay(self,'isnewton',"0: Picard's fixed point, 1: Newton's method, 2: hybrid")) + string="%s\n%s"%(string,fielddisplay(self,'maxiter','maximum number of nonlinear iterations')) + string="%s\n%s"%(string,fielddisplay(self,'viscosity_overshoot','over-shooting constant new=new+C*(new-old)')) + + string="%s\n%s"%(string,'\n boundary conditions:') + string="%s\n%s"%(string,fielddisplay(self,'spcvx','x-axis velocity constraint (NaN means no constraint) [m/yr]')) + string="%s\n%s"%(string,fielddisplay(self,'spcvy','y-axis velocity constraint (NaN means no constraint) [m/yr]')) + string="%s\n%s"%(string,fielddisplay(self,'spcvz','z-axis velocity constraint (NaN means no constraint) [m/yr]')) + string="%s\n%s"%(string,fielddisplay(self,'icefront','segments on ice front list (last column 0: Air, 1: Water, 2: Ice')) + + string="%s\n%s"%(string,'\n Rift options:') + string="%s\n%s"%(string,fielddisplay(self,'rift_penalty_threshold','threshold for instability of mechanical constraints')) + string="%s\n%s"%(string,fielddisplay(self,'rift_penalty_lock','number of iterations before rift penalties are locked')) + + string="%s\n%s"%(string,'\n Penalty options:') + string="%s\n%s"%(string,fielddisplay(self,'penalty_factor','offset used by penalties: penalty = Kmax*10^offset')) + string="%s\n%s"%(string,fielddisplay(self,'vertex_pairing','pairs of vertices that are penalized')) + + string="%s\n%s"%(string,'\n Other:') + string="%s\n%s"%(string,fielddisplay(self,'shelf_dampening','use dampening for floating ice ? Only for FS model')) + string="%s\n%s"%(string,fielddisplay(self,'FSreconditioning','multiplier for incompressibility equation. Only for FS model')) + string="%s\n%s"%(string,fielddisplay(self,'referential','local referential')) + string="%s\n%s"%(string,fielddisplay(self,'loadingforce','loading force applied on each point [N/m^3]')) + string="%s\n%s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested')) + + return string + #}}} + def setdefaultparameters(self): # {{{ + #maximum of non-linear iterations. + self.maxiter=100 + + #Convergence criterion: absolute, relative and residual + self.restol=10**-4 + self.reltol=0.01 + self.abstol=10 + + self.FSreconditioning=10**13 + self.shelf_dampening=0 + + #Penalty factor applied kappa=max(stiffness matrix)*10^penalty_factor + self.penalty_factor=3 + + #coefficient to update the viscosity between each iteration of + #a stressbalance according to the following formula + #viscosity(n)=viscosity(n)+viscosity_overshoot(viscosity(n)-viscosity(n-1)) + self.viscosity_overshoot=0 + + #Stop the iterations of rift if below a threshold + self.rift_penalty_threshold=0 + + #in some solutions, it might be needed to stop a run when only + #a few constraints remain unstable. For thermal computation, this + #parameter is often used. + self.rift_penalty_lock=10 + + return self + #}}} + def checkconsistency(self,md,solution,analyses): # {{{ + + #Early return + if DiagnosticHorizAnalysisEnum() not in analyses: + return md + + md = checkfield(md,'diagnostic.spcvx','forcing',1) + md = checkfield(md,'diagnostic.spcvy','forcing',1) + if md.mesh.dimension==3: + md = checkfield(md,'diagnostic.spcvz','forcing',1) + md = checkfield(md,'diagnostic.restol','size',[1],'>',0) + md = checkfield(md,'diagnostic.reltol','size',[1]) + md = checkfield(md,'diagnostic.abstol','size',[1]) + md = checkfield(md,'diagnostic.isnewton','numel',[1],'values',[0,1,2]) + md = checkfield(md,'diagnostic.FSreconditioning','size',[1],'NaN',1) + md = checkfield(md,'diagnostic.viscosity_overshoot','size',[1],'NaN',1) + if md.mesh.dimension==2: + md = checkfield(md,'diagnostic.icefront','size',[float('NaN'),4],'NaN',1) + else: + md = checkfield(md,'diagnostic.icefront','size',[float('NaN'),6],'NaN',1) + md = checkfield(md,'diagnostic.icefront[:,-1]','values',[0,1,2]) + md = checkfield(md,'diagnostic.maxiter','size',[1],'>=',1) + md = checkfield(md,'diagnostic.referential','size',[md.mesh.numberofvertices,6]) + md = checkfield(md,'diagnostic.loadingforce','size',[md.mesh.numberofvertices,3]) + if not md.diagnostic.requested_outputs: + md = checkfield(md,'diagnostic.requested_outputs','size',[float('NaN'),1]) + + #singular solution +# if ~any((~isnan(md.diagnostic.spcvx)+~isnan(md.diagnostic.spcvy))==2), + if not numpy.any(numpy.logical_and(numpy.logical_not(numpy.isnan(md.diagnostic.spcvx)),numpy.logical_not(numpy.isnan(md.diagnostic.spcvy)))): + md.checkmessage("model is not well posed (singular). You need at least one node with fixed velocity!") + #CHECK THAT EACH LINES CONTAINS ONLY NAN VALUES OR NO NAN VALUES +# if any(sum(isnan(md.diagnostic.referential),2)~=0 & sum(isnan(md.diagnostic.referential),2)~=6), + if numpy.any(numpy.logical_and(numpy.sum(numpy.isnan(md.diagnostic.referential),axis=1)!=0,numpy.sum(numpy.isnan(md.diagnostic.referential),axis=1)!=6)): + md.checkmessage("Each line of diagnostic.referential should contain either only NaN values or no NaN values") + #CHECK THAT THE TWO VECTORS PROVIDED ARE ORTHOGONAL +# if any(sum(isnan(md.diagnostic.referential),2)==0), + if numpy.any(numpy.sum(numpy.isnan(md.diagnostic.referential),axis=1)==0): + pos=[i for i,item in enumerate(numpy.sum(numpy.isnan(md.diagnostic.referential),axis=1)) if item==0] +# numpy.inner (and numpy.dot) calculate all the dot product permutations, resulting in a full matrix multiply +# if numpy.any(numpy.abs(numpy.inner(md.diagnostic.referential[pos,0:2],md.diagnostic.referential[pos,3:5]).diagonal())>sys.float_info.epsilon): +# md.checkmessage("Vectors in diagnostic.referential (columns 1 to 3 and 4 to 6) must be orthogonal") + for item in md.diagnostic.referential[pos,:]: + if numpy.abs(numpy.inner(item[0:2],item[3:5]))>sys.float_info.epsilon: + md.checkmessage("Vectors in diagnostic.referential (columns 1 to 3 and 4 to 6) must be orthogonal") + #CHECK THAT NO rotation specified for FS Grounded ice at base +# if md.mesh.dimension==3 & md.flowequation.isFS, + if md.mesh.dimension==3 and md.flowequation.isFS: + pos=numpy.nonzero(numpy.logical_and(md.mask.vertexongroundedice,md.mesh.vertexonbed)) + if numpy.any(numpy.logical_not(numpy.isnan(md.diagnostic.referential[pos,:]))): + md.checkmessage("no referential should be specified for basal vertices of grounded ice") + + return md + # }}} + def marshall(self,md,fid): # {{{ + WriteData(fid,'object',self,'class','diagnostic','fieldname','spcvx','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1) + WriteData(fid,'object',self,'class','diagnostic','fieldname','spcvy','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1) + WriteData(fid,'object',self,'class','diagnostic','fieldname','spcvz','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1) + WriteData(fid,'object',self,'class','diagnostic','fieldname','restol','format','Double') + WriteData(fid,'object',self,'class','diagnostic','fieldname','reltol','format','Double') + WriteData(fid,'object',self,'class','diagnostic','fieldname','abstol','format','Double') + WriteData(fid,'object',self,'class','diagnostic','fieldname','isnewton','format','Integer') + WriteData(fid,'object',self,'class','diagnostic','fieldname','FSreconditioning','format','Double') + WriteData(fid,'object',self,'class','diagnostic','fieldname','viscosity_overshoot','format','Double') + WriteData(fid,'object',self,'class','diagnostic','fieldname','maxiter','format','Integer') + WriteData(fid,'object',self,'class','diagnostic','fieldname','shelf_dampening','format','Integer') + WriteData(fid,'object',self,'class','diagnostic','fieldname','vertex_pairing','format','DoubleMat','mattype',3) + WriteData(fid,'object',self,'class','diagnostic','fieldname','penalty_factor','format','Double') + WriteData(fid,'object',self,'class','diagnostic','fieldname','rift_penalty_lock','format','Integer') + WriteData(fid,'object',self,'class','diagnostic','fieldname','rift_penalty_threshold','format','Integer') + WriteData(fid,'object',self,'class','diagnostic','fieldname','referential','format','DoubleMat','mattype',1) + WriteData(fid,'object',self,'class','diagnostic','fieldname','requested_outputs','format','DoubleMat','mattype',3) + WriteData(fid,'data',self.loadingforce[:,0],'format','DoubleMat','mattype',1,'enum',LoadingforceXEnum()) + WriteData(fid,'data',self.loadingforce[:,1],'format','DoubleMat','mattype',1,'enum',LoadingforceYEnum()) + WriteData(fid,'data',self.loadingforce[:,2],'format','DoubleMat','mattype',1,'enum',LoadingforceZEnum()) + # }}} Index: /issm/trunk-jpl/src/m/consistency/checkfield.m =================================================================== --- /issm/trunk-jpl/src/m/consistency/checkfield.m (revision 15613) +++ /issm/trunk-jpl/src/m/consistency/checkfield.m (revision 15614) @@ -23,6 +23,4 @@ % Example: % md = checkfield(md,'mesh.elementonbed','size',[md.mesh.numberofelements 1],'values',[0 1]); -% md = checkfield(md,'diagnostic.icefront','size',[NaN 4],'NaN',1); -% md = checkfield(md,'diagnostic.icefront(:,end)','values',[0 1 2]); %get options Index: /issm/trunk-jpl/src/m/consistency/checkfield.py =================================================================== --- /issm/trunk-jpl/src/m/consistency/checkfield.py (revision 15613) +++ /issm/trunk-jpl/src/m/consistency/checkfield.py (revision 15614) @@ -29,6 +29,4 @@ Example: md = checkfield(md,'mesh.elementonbed','size',[md.mesh.numberofelements 1],'values',[0 1]); - md = checkfield(md,'diagnostic.icefront','size',[NaN 4],'NaN',1); - md = checkfield(md,'diagnostic.icefront(:,end)','values',[0 1 2]); """ Index: /issm/trunk-jpl/test/NightlyRun/test201.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test201.m (revision 15613) +++ /issm/trunk-jpl/test/NightlyRun/test201.m (revision 15614) @@ -5,4 +5,6 @@ md.cluster=generic('name',oshostname(),'np',3); md.verbose=verbose('all'); +md.mask.icelevelset(find(md.mask.icelevelset==0))=-1; + md=solve(md,DiagnosticSolutionEnum()); Index: /issm/trunk-jpl/test/NightlyRun/test217.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test217.m (revision 15613) +++ /issm/trunk-jpl/test/NightlyRun/test217.m (revision 15614) @@ -54,7 +54,5 @@ nodeonicefront=zeros(md.mesh.numberofvertices,1); pos=find(md.mesh.y==ymax); nodeonicefront(pos)=1; -pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2))); diagnostic.icefront=md.mesh.segments(pos,:); -diagnostic.icefront=[diagnostic.icefront 1*md.mask.elementonfloatingice(diagnostic.icefront(:,end))]; -md.diagnostic.icefront=diagnostic.icefront; +md.mask.icelevelset=1-nodeonicefront; md=solve(md,DiagnosticSolutionEnum()); Index: /issm/trunk-jpl/test/NightlyRun/test217.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test217.py (revision 15613) +++ /issm/trunk-jpl/test/NightlyRun/test217.py (revision 15614) @@ -66,8 +66,5 @@ pos=numpy.nonzero(md.mesh.y==ymax) nodeonicefront[pos]=1 -pos=numpy.nonzero(numpy.logical_or(nodeonicefront[md.mesh.segments[:,0].astype(int)-1],nodeonicefront[md.mesh.segments[:,1].astype(int)-1]))[0] -diagnostic.icefront=md.mesh.segments[pos,:] -diagnostic.icefront=numpy.hstack((diagnostic.icefront,1.*md.mask.elementonfloatingice[diagnostic.icefront[:,-1].astype(int)-1].reshape(-1,1))) -md.diagnostic.icefront=diagnostic.icefront +md.mask.icelevelset=1-nodeonicefront md=solve(md,DiagnosticSolutionEnum()) Index: /issm/trunk-jpl/test/Par/Pig.par =================================================================== --- /issm/trunk-jpl/test/Par/Pig.par (revision 15613) +++ /issm/trunk-jpl/test/Par/Pig.par (revision 15614) @@ -47,5 +47,4 @@ %Boundary conditions: md=SetMarineIceSheetBC(md); -md.diagnostic.icefront(:,end)=1; %Change name so that no test have the same name Index: /issm/trunk-jpl/test/Par/Pig.py =================================================================== --- /issm/trunk-jpl/test/Par/Pig.py (revision 15613) +++ /issm/trunk-jpl/test/Par/Pig.py (revision 15614) @@ -58,5 +58,4 @@ #Boundary conditions: md=SetMarineIceSheetBC(md) -md.diagnostic.icefront[:,-1]=1 #Change name so that no test have the same name Index: /issm/trunk-jpl/test/Par/RoundSheetShelf.par =================================================================== --- /issm/trunk-jpl/test/Par/RoundSheetShelf.par (revision 15613) +++ /issm/trunk-jpl/test/Par/RoundSheetShelf.par (revision 15614) @@ -101,7 +101,4 @@ md.diagnostic.spcvy(pos)=0; -pressureload=md.mesh.segments; -pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end)) + 0*md.mask.elementongroundedice(pressureload(:,end))]; -md.diagnostic.icefront=pressureload; md.mask.icelevelset(find(md.mesh.vertexonboundary))=0; md.balancethickness.spcthickness=NaN*ones(md.mesh.numberofvertices,1); Index: /issm/trunk-jpl/test/Par/RoundSheetShelf.py =================================================================== --- /issm/trunk-jpl/test/Par/RoundSheetShelf.py (revision 15613) +++ /issm/trunk-jpl/test/Par/RoundSheetShelf.py (revision 15614) @@ -110,7 +110,4 @@ md.diagnostic.spcvy[pos]=0 -pressureload=copy.deepcopy(md.mesh.segments) -pressureload=numpy.hstack((pressureload,(1*md.mask.elementonfloatingice[pressureload[:,-1].astype(int)-1] + 0*md.mask.elementongroundedice[pressureload[:,-1].astype(int)-1]).reshape(-1,1))) -md.diagnostic.icefront=pressureload pos=numpy.nonzero(md.mesh.vertexonboundary) md.mask.icelevelset[pos]=0