Changeset 18757

Timestamp:

11/07/14 15:38:27 (10 years ago)

Author:

srebuffi

Message:

CHG: added calving class

Location:

issm/trunk-jpl/src/m/classes

Files:

• calving.m (added)
• calving.py (added)
• calvinglevermann.m (added)
• masstransport.m (modified) (6 diffs)
• model.m (modified) (4 diffs)
• model.py (modified) (5 diffs)
• transient.m (modified) (7 diffs)

issm/trunk-jpl/src/m/classes/masstransport.m

@@ -7 +7 @@
         properties (SetAccess=public) 
                  spcthickness           = NaN;
-                 iscalvingrate          = 0;
-                 calvingrate            = NaN;
-                 levermann_calving_coeff= 0;
                  isfreesurface          = 0;
                  min_thickness          = 0;
@@ -45 +42 @@
         end
         methods
-         function createxml(obj,fid) % {{{
-            fprintf(fid, '<!-- masstransport -->\n');            
-                    
-            % Masstransport solution parameters
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="Masstransport solution parameters">','<section name="masstransport" />');                    
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="spcthickness" type="',class(obj.spcthickness),'" default="',convert2str(obj.spcthickness),'">','     <section name="masstransport" />','     <help> thickness constraints (NaN means no constraint) [m] </help>','  </parameter>');
-                                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="iscalvingrate" type="',class(obj.iscalvingrate),'" default="',convert2str(obj.iscalvingrate),'">','     <section name="masstransport" />','     <help> do we use calving (loss of ice) </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="calvingrate" type="',class(obj.calvingrate),'" default="',convert2str(obj.calvingrate),'">','     <section name="masstransport" />','     <help> calving rate at given location [m/a] </help>','  </parameter>');
-                                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="levermann_calving_coeff" type="',class(obj.levermann_calving_coeff),'" default="',convert2str(obj.levermann_calving_coeff),'">','     <section name="masstransport" />','     <help> Proportionality coefficient in Levermann model  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="isfreesurface" type="',class(obj.isfreesurface),'" default="',convert2str(obj.isfreesurface),'">','     <section name="masstransport" />','     <help> do we use free surfaces (FS only) are mass conservation </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="min_thickness" type="',class(obj.min_thickness),'" default="',convert2str(obj.min_thickness),'">','     <section name="masstransport" />','     <help> minimum ice thickness allowed [m] </help>','  </parameter>');
-            
-            % hydrostatic_adjustment drop-down (incremental or absolute )
-            fprintf(fid,'%s\n%s\n%s\n%s\n','  <parameter key ="hydrostatic_adjustment" type="alternative" optional ="false">','     <section name="masstransport" />','     <help> adjustment of ice shelves surface and bed elevations: ''Incremental'' or ''Absolute''  </help>');
-            fprintf(fid,'%s\n','       <option value="Incremental" type="string" default="true"> </option>');
-            fprintf(fid,'%s\n%s\n','       <option value="Absolute" type="string" default="false"> </option>','</parameter>');
-            
-            %stabilization drop-down (0, 1, 2, or 3)
-            fprintf(fid,'%s\n%s\n%s\n%s\n','  <parameter key ="stabilization" type="alternative" optional="false">','     <section name="masstransport" />','     <help> 0: no, 1:  artificial_diffusivity, 2: streamline upwinding, 3: discontinuous Galerkin </help>');
-            fprintf(fid,'%s\n','       <option value="0" type="string" default="true"> </option>');
-            fprintf(fid,'%s\n','       <option value="1" type="string" default="false"> </option>');
-            fprintf(fid,'%s\n','       <option value="2" type="string" default="false"> </option>');
-            fprintf(fid,'%s\n%s\n','       <option value="3" type="string" default="false"> </option>','</parameter>');
-
-            fprintf(fid,'%s\n%s\n','</frame>');    
-            
-            %Penalty options
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="2" label="Penalty options">','<section name="masstransport" />'); 
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vertex_pairing" type="',class(obj.vertex_pairing),'" default="',convert2str(obj.vertex_pairing),'">','     <section name="masstransport" />','     <help> offset used by penalties: penalty = Kmax*10^offset </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="penalty_factor" type="',class(obj.penalty_factor),'" default="',convert2str(obj.penalty_factor),'">','     <section name="masstransport" />','     <help> pairs of vertices that are penalized </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="requested_outputs" type="',class(obj.requested_outputs),'" default="',convert2str(obj.requested_outputs),'">','     <section name="masstransport" />','     <help> additional outputs requested </help>','  </parameter>');
-            fprintf(fid,'%s\n%s\n','</frame>');
-        
-        end % }}}
                 function obj = masstransport(varargin) % {{{
                         switch nargin
@@ -104 +67 @@
                 function obj = setdefaultparameters(obj) % {{{
 
-                        %Proportionality coefficient in Levermann model
-                        obj.levermann_calving_coeff=8e21;
-                        
                         %Type of stabilization to use 0:nothing 1:artificial_diffusivity 3:Discontinuous Galerkin
                         obj.stabilization=1;
@@ -128 +88 @@
 
                         md = checkfield(md,'fieldname','masstransport.spcthickness','forcing',1);
-                        if(ismember(LevelsetAnalysisEnum(), analyses) & md.transient.islevelset)
-                                md = checkfield(md,'fieldname','masstransport.iscalvingrate','values',[0 1]);
-                                md = checkfield(md,'fieldname','masstransport.calvingrate','NaN',1,'size',[md.mesh.numberofvertices 1],'>=',0);
-                                md = checkfield(md,'fieldname','masstransport.levermann_calving_coeff','>',0);
-                        end
-
                         md = checkfield(md,'fieldname','masstransport.isfreesurface','values',[0 1]);
                         md = checkfield(md,'fieldname','masstransport.hydrostatic_adjustment','values',{'Absolute' 'Incremental'});
@@ -144 +98 @@
                         disp(sprintf('   Masstransport solution parameters:'));
                         fielddisplay(obj,'spcthickness','thickness constraints (NaN means no constraint) [m]');
-                        fielddisplay(obj,'iscalvingrate','do we use calving (loss of ice)');
-                        fielddisplay(obj,'calvingrate','calving rate at given location [m/a]');
-                        fielddisplay(obj,'levermann_calving_coeff','proportionality coefficient in Levermann model');
                         fielddisplay(obj,'isfreesurface','do we use free surfaces (FS only) are mass conservation');
                         fielddisplay(obj,'min_thickness','minimum ice thickness allowed [m]');
@@ -163 +114 @@
 
                         WriteData(fid,'object',obj,'fieldname','spcthickness','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);
-                        WriteData(fid,'object',obj,'fieldname','calvingrate','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1,'scale',1./yts);
-                        WriteData(fid,'object',obj,'fieldname','iscalvingrate','format','Boolean');
                         WriteData(fid,'object',obj,'fieldname','isfreesurface','format','Boolean');
-                        WriteData(fid,'object',obj,'fieldname','levermann_calving_coeff','format','Double');
                         WriteData(fid,'object',obj,'fieldname','min_thickness','format','Double');
                         WriteData(fid,'data',StringToEnum(obj.hydrostatic_adjustment),'format','Integer','enum',MasstransportHydrostaticAdjustmentEnum());

issm/trunk-jpl/src/m/classes/model.m

@@ -37 +37 @@
                 steadystate      = 0;
                 transient        = 0;
+                calving          = 0;
                 gia              = 0;
 
@@ -796 +797 @@
                         end
                         md.masstransport.spcthickness=project3d(md,'vector',md.masstransport.spcthickness,'type','node');
-                        md.masstransport.calvingrate=project3d(md,'vector',md.masstransport.calvingrate,'type','node');
                         md.balancethickness.spcthickness=project3d(md,'vector',md.balancethickness.spcthickness,'type','node');
                         md.damage.spcdamage=project3d(md,'vector',md.damage.spcdamage,'type','node');
                         md.stressbalance.referential=project3d(md,'vector',md.stressbalance.referential,'type','node');
                         md.stressbalance.loadingforce=project3d(md,'vector',md.stressbalance.loadingforce,'type','node');
+
+                        % Calving variables
+                        if isa(md.calving,'calving'),md.calving.calvingrate=project3d(md,'vector',md.calving.calvingrate,'type','node');end;
+                        if isa(md.calving,'calvinglevermann'),md.calving.coeff=project3d(md,'vector',md.calving.coeff,'type','node');end;
 
                         % Hydrologydc variables
@@ -1160 +1164 @@
                         md.steadystate      = steadystate();
                         md.transient        = transient();
+                        md.calving          = calving();
                         md.gia              = gia();
                         md.seaice           = seaice();
@@ -1320 +1325 @@
                         disp(sprintf('%19s: %-22s -- %s','verbose'         ,['[1x1 ' class(obj.verbose) ']'],'verbosity level in solve'));
                         disp(sprintf('%19s: %-22s -- %s','settings'        ,['[1x1 ' class(obj.settings) ']'],'settings properties'));
-                        disp(sprintf('%19s: %-22s -- %s','toolkits'          ,['[1x1 ' class(obj.toolkits) ']'],'PETSc options for each solution'));
+                        disp(sprintf('%19s: %-22s -- %s','toolkits'        ,['[1x1 ' class(obj.toolkits) ']'],'PETSc options for each solution'));
                         disp(sprintf('%19s: %-22s -- %s','cluster'         ,['[1x1 ' class(obj.cluster) ']'],'cluster parameters (number of cpus...)'));
                         disp(sprintf('%19s: %-22s -- %s','balancethickness',['[1x1 ' class(obj.balancethickness) ']'],'parameters for balancethickness solution'));
-                        disp(sprintf('%19s: %-22s -- %s','stressbalance'      ,['[1x1 ' class(obj.stressbalance) ']'],'parameters for stressbalance solution'));
+                        disp(sprintf('%19s: %-22s -- %s','stressbalance'   ,['[1x1 ' class(obj.stressbalance) ']'],'parameters for stressbalance solution'));
                         disp(sprintf('%19s: %-22s -- %s','groundingline'   ,['[1x1 ' class(obj.groundingline) ']'],'parameters for groundingline solution'));
                         disp(sprintf('%19s: %-22s -- %s','hydrology'       ,['[1x1 ' class(obj.hydrology) ']'],'parameters for hydrology solution'));
-                        disp(sprintf('%19s: %-22s -- %s','masstransport'      ,['[1x1 ' class(obj.masstransport) ']'],'parameters for masstransport solution'));
+                        disp(sprintf('%19s: %-22s -- %s','masstransport'   ,['[1x1 ' class(obj.masstransport) ']'],'parameters for masstransport solution'));
                         disp(sprintf('%19s: %-22s -- %s','thermal'         ,['[1x1 ' class(obj.thermal) ']'],'parameters for thermal solution'));
                         disp(sprintf('%19s: %-22s -- %s','steadystate'     ,['[1x1 ' class(obj.steadystate) ']'],'parameters for steadystate solution'));
                         disp(sprintf('%19s: %-22s -- %s','transient'       ,['[1x1 ' class(obj.transient) ']'],'parameters for transient solution'));
-                        disp(sprintf('%19s: %-22s -- %s','gia'       ,['[1x1 ' class(obj.gia) ']'],'parameters for gia solution'));
+                        disp(sprintf('%19s: %-22s -- %s','calving'         ,['[1x1 ' class(obj.calving) ']'],'parameters for calving'));
+                        disp(sprintf('%19s: %-22s -- %s','gia'             ,['[1x1 ' class(obj.gia) ']'],'parameters for gia solution'));
                         disp(sprintf('%19s: %-22s -- %s','autodiff'        ,['[1x1 ' class(obj.autodiff) ']'],'automatic differentiation parameters'));
                         disp(sprintf('%19s: %-22s -- %s','flaim'           ,['[1x1 ' class(obj.flaim) ']'],'flaim parameters'));

issm/trunk-jpl/src/m/classes/model.py

@@ -15 +15 @@
 from basalforcings import basalforcings
 from matice import matice
+from calving import calving
 from damage import damage
 from friction import friction
@@ -88 +89 @@
                 self.steadystate      = steadystate()
                 self.transient        = transient()
+                self.calving          = calving()
                 self.gia              = gia()
 
@@ -129 +131 @@
                         'steadystate',\
                         'transient',\
+                        'calving',\
                         'gia',\
                         'autodiff',\
@@ -168 +171 @@
                 string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("steadystate","[%s,%s]" % ("1x1",obj.steadystate.__class__.__name__),"parameters for steadystate solution"))
                 string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("transient","[%s,%s]" % ("1x1",obj.transient.__class__.__name__),"parameters for transient solution"))
+                string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("calving","[%s,%s]" % ("1x1",obj.calving.__class__.__name__),"parameters for calving"))
                 string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("autodiff","[%s,%s]" % ("1x1",obj.autodiff.__class__.__name__),"automatic differentiation parameters"))
                 string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("flaim","[%s,%s]" % ("1x1",obj.flaim.__class__.__name__),"flaim parameters"))
@@ -658 +662 @@
                 md.stressbalance.referential=project3d(md,'vector',md.stressbalance.referential,'type','node')
                 md.stressbalance.loadingforce=project3d(md,'vector',md.stressbalance.loadingforce,'type','node')
-                md.masstransport.calvingrate=project3d(md,'vector',md.masstransport.calvingrate,'type','node')
+
+                # Calving variables
+                if hasattr(md.calving,'calving'):
+                        md.calving.calvingrate=project3d(md,'vector',md.calving.calvingrate,'type','node')
+                if hasattr(md.calving,'calvinglevermann'):
+                        md.calving.coeff=project3d(md,'vector',md.calving.coeff,'type','node')
 
                 # Hydrologydc variables

issm/trunk-jpl/src/m/classes/transient.m

@@ -13 +13 @@
                 isdamageevolution = 0;
                 islevelset        = 0;
+                iscalving         = 0;
                 ishydrology       = 0;
                 requested_outputs = {};
@@ -29 +30 @@
                 fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="isdamageevolution" type="',class(obj.isdamageevolution),'" default="',convert2str(obj.isdamageevolution),'">','     <section name="transient" />','     <help> indicates whether damage evolution is used in the transient </help>','  </parameter>');
                 fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="islevelset" type="',class(obj.islevelset),'" default="',convert2str(obj.islevelset),'">','     <section name="transient" />','     <help> LEVEL SET DESCRIPTION...  </help>','  </parameter>');
+                fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="iscalving" type="',class(obj.iscalving),'" default="',convert2str(obj.iscalving),'">','     <section name="transient" />','     <help> indicates whether calving is used in the transient </help>','  </parameter>');
                 fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="requested_outputs" type="',class(obj.requested_outputs),'" default="',convert2str(obj.requested_outputs),'">','     <section name="transient" />','     <help> list of additional outputs requested </help>','  </parameter>');
             fprintf(fid,'%s\n%s\n','</frame>');    
@@ -51 +53 @@
                         obj.isdamageevolution = 0;
                         obj.islevelset      = 0;
+                        obj.iscalving       =0;
                         obj.ishydrology     = 0;
 
@@ -66 +69 @@
                         obj.isdamageevolution = 0;
                         obj.islevelset      = 0;
+                        obj.iscalving       = 0;
                         obj.ishydrology     = 0;
 
@@ -90 +94 @@
                         md = checkfield(md,'fieldname','transient.isdamageevolution','numel',[1],'values',[0 1]);
                         md = checkfield(md,'fieldname','transient.islevelset','numel',[1],'values',[0 1]);
+                        md = checkfield(md,'fieldname','transient.iscalving','numel',[1],'values',[0 1]);
                         md = checkfield(md,'fieldname','transient.ishydrology','numel',[1],'values',[0 1]);
                         md = checkfield(md,'fieldname','transient.requested_outputs','stringrow',1);
@@ -104 +109 @@
                         fielddisplay(obj,'isdamageevolution','indicates whether damage evolution is used in the transient');
                         fielddisplay(obj,'islevelset','LEVEL SET DESCRIPTION...');
+                        fielddisplay(obj,'iscalving','indicates whether calving is used in the transient');
                         fielddisplay(obj,'ishydrology','indicates whether an hydrology model is used');
                         fielddisplay(obj,'requested_outputs','list of additional outputs requested');
@@ -117 +123 @@
                         WriteData(fid,'object',obj,'fieldname','ishydrology','format','Boolean');
                         WriteData(fid,'object',obj,'fieldname','islevelset','format','Boolean');
+                        WriteData(fid,'object',obj,'fieldname','iscalving','format','Boolean');
 
                         %process requested outputs