Index: /issm/trunk/src/m/Makefile.am =================================================================== --- /issm/trunk/src/m/Makefile.am (revision 12328) +++ /issm/trunk/src/m/Makefile.am (revision 12329) @@ -13,9 +13,5 @@ ./planet/*.m \ ./qmu/*.m \ - ./shared/*.m \ - ./solutions/*.m \ - ./solvers/*.m \ ./utils/*.m \ - ./utils/Analysis/*.m \ ./utils/Array/*.m \ ./utils/BC/*.m \ Index: /issm/trunk/src/m/classes/autodiff.py =================================================================== --- /issm/trunk/src/m/classes/autodiff.py (revision 12329) +++ /issm/trunk/src/m/classes/autodiff.py (revision 12329) @@ -0,0 +1,28 @@ +#module imports +from fielddisplay import fielddisplay + +class autodiff: + #properties + def __init__(self): + # {{{ Properties + self.isautodiff = False + self.forward = True + self.reverse = False + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + string=' automatic differentiation parameters:' + string="%s\n%s"%(string,fielddisplay(obj,'isautodiff','indicates if the automatic differentiation is activated')) + string="%s\n%s"%(string,fielddisplay(obj,'forward','forward differentiation')) + string="%s\n%s"%(string,fielddisplay(obj,'reverse','backward differentiation')) + return string + #}}} + def setdefaultparameters(obj): + # {{{setdefaultparameters + return obj + #}}} + Index: /issm/trunk/src/m/classes/balancethickness.py =================================================================== --- /issm/trunk/src/m/classes/balancethickness.py (revision 12329) +++ /issm/trunk/src/m/classes/balancethickness.py (revision 12329) @@ -0,0 +1,33 @@ +#module imports +from fielddisplay import fielddisplay +class balancethickness: + #properties + def __init__(self): + # {{{ Properties + self.spcthickness = float('NaN') + self.thickening_rate = float('NaN') + self.stabilization = 0 + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + + string=' balance thickness solution parameters:' + + string="%s\n\n%s"%(string,fielddisplay(obj,'spcthickness','thickness constraints (NaN means no constraint)')) + string="%s\n%s"%(string,fielddisplay(obj,'thickening_rate','ice thickening rate used in the mass conservation (dh/dt)')) + string="%s\n%s"%(string,fielddisplay(obj,'stabilization','0: None, 1: SU, 2: MacAyeal''s artificial diffusivity, 3:DG')) + return string + #}}} + def setdefaultparameters(obj): + # {{{setdefaultparameters + + #Type of stabilization used + obj.stabilization=1 + + return obj + #}}} + Index: /issm/trunk/src/m/classes/basalforcings.py =================================================================== --- /issm/trunk/src/m/classes/basalforcings.py (revision 12329) +++ /issm/trunk/src/m/classes/basalforcings.py (revision 12329) @@ -0,0 +1,29 @@ +#module imports +from fielddisplay import fielddisplay + +class basalforcings: + #properties + def __init__(self): + # {{{ Properties + self.melting_rate = float('NaN') + self.melting_rate_correction = float('NaN') + self.geothermalflux = float('NaN') + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + string=" basal forcings parameters:" + + string="%s\n\n%s"%(string,fielddisplay(obj,"melting_rate","basal melting rate (positive if melting)")) + string="%s\n%s"%(string,fielddisplay(obj,"melting_rate_correction","additional melting applied when the grounding line retreats")) + string="%s\n%s"%(string,fielddisplay(obj,"geothermalflux","geothermal heat flux [W/m^2]")) + return string + #}}} + def setdefaultparameters(obj): + # {{{setdefaultparameters + return obj + #}}} + Index: /issm/trunk/src/m/classes/clusters/generic.m =================================================================== --- /issm/trunk/src/m/classes/clusters/generic.m (revision 12328) +++ /issm/trunk/src/m/classes/clusters/generic.m (revision 12329) @@ -2,5 +2,4 @@ % % Usage: -% cluster=generic('name','astrid',); % cluster=generic('name','astrid','np',3); % cluster=generic('name',oshostname(),'np',3,'login','username'); @@ -14,9 +13,9 @@ port=0; interactive=1; - codepath=[issmtier() '/bin']; - executionpath=[issmtier() '/execution']; - valgrind=[issmtier() '/externalpackages/valgrind/install/bin/valgrind']; - valgrindlib=[issmtier() '/externalpackages/valgrind/install/lib/libmpidebug.so']; - valgrindsup=[issmtier() '/externalpackages/valgrind/issm.supp']; + codepath=[issmdir() '/bin']; + executionpath=[issmdir() '/execution']; + valgrind=[issmdir() '/externalpackages/valgrind/install/bin/valgrind']; + valgrindlib=[issmdir() '/externalpackages/valgrind/install/lib/libmpidebug.so']; + valgrindsup=[issmdir() '/externalpackages/valgrind/issm.supp']; %}}} end @@ -28,6 +27,5 @@ %get name - if ~exist(options,'name'), error('option ''name'' has not been provided'); end - cluster.name=getfieldvalue(options,'name'); + cluster.name=getfieldvalue(options,'name',oshostname()); %initialize cluster using user settings if provided Index: /issm/trunk/src/m/classes/clusters/generic.py =================================================================== --- /issm/trunk/src/m/classes/clusters/generic.py (revision 12329) +++ /issm/trunk/src/m/classes/clusters/generic.py (revision 12329) @@ -0,0 +1,217 @@ +#GENERIC cluster class definition +# +# Usage: +# cluster=generic('name','astrid',); +# cluster=generic('name','astrid','np',3); +# cluster=generic('name',oshostname(),'np',3,'login','username'); + + +class generic: + #properties + def __init__(self): + # {{{ Properties + self.name='' + self.login='' + self.np=1 + self.port=0 + self.interactive=1 + self.codepath=issmdir() + '/bin' + self.executionpath=issmdir() + '/execution' + self.valgrind=issmdir() + '/externalpackages/valgrind/install/bin/valgrind' + self.valgrindlib=issmdir() + '/externalpackages/valgrind/install/lib/libmpidebug.so' + self.valgrindsup=issmdir() + '/externalpackages/valgrind/issm.supp' + #}}} + def __repr__(obj): + # {{{ Display + + string="class 'generic' object:" + string="%s\n\n%s"%(string,"%s%s"%(' name: ',obj.name)) + string="%s\n%s"%(string,"%s%i"%(' np: ',obj.np)) + string="%s\n%s"%(string,"%s%i"%(' port: ',obj.port)) + string="%s\n%s"%(string,"%s%s"%(' codepath: ',obj.codepath)) + string="%s\n%s"%(string,"%s%s"%(' executionpath: ',obj.executionpath)) + string="%s\n%s"%(string,"%s%s"%(' valgrind: ',obj.valgrind)) + string="%s\n%s"%(string,"%s%s"%(' valgrindlib: ',obj.valgrindlib)) + string="%s\n%s"%(string,"%s%s"%(' valgrindsup: ',obj.valgrindsup)) + return string + #}}} + + +#old matlab +# function cluster=generic(varargin) % {{{1 +# +# %use provided options to change fields +# options=pairoptions(varargin{:}); +# +# %get name +# if ~exist(options,'name'), error('option ''name'' has not been provided'); end +# cluster.name=getfieldvalue(options,'name'); +# +# %initialize cluster using user settings if provided +# if (exist([cluster.name '_settings'])==2), eval([cluster.name '_settings']); end +# +# %OK get other fields +# for i=1:size(options.list,1), +# fieldname=options.list{i,1}; +# fieldvalue=options.list{i,2}; +# if ismember(fieldname,properties('generic')), +# cluster.(fieldname)=fieldvalue; +# else +# disp(['''' fieldname ''' is not a property of cluster generic']); +# end +# end +# end +# %}}} +# function checkconsistency(cluster,md,solution,analyses) % {{{1 +# if cluster.np<1 +# checkmessage(['number of processors should be at least 1']); +# end +# if isnan(cluster.np), +# checkessage('number of processors should not be NaN!'); +# end +# end +# %}}} +# function BuildQueueScript(cluster,md) % {{{1 +# +# %retrieve parameters +# modelname=md.miscellaneous.name; +# solution=md.private.solution; +# isvalgrind=md.debug.valgrind; +# isgprof=md.debug.gprof; +# +# %open file for writing: +# if ~ispc, +# fid=fopen([modelname '.queue'],'w'); +# else +# fid=fopen([modelname '.bat'],'w'); +# end +# +# %write instructions for launching a job on the cluster +# if ~ispc, +# fprintf(fid,'#!/bin/sh\n'); +# else +# fprintf(fid,'@echo off\n'); +# end +# +# if ~isvalgrind, +# if cluster.interactive +# if ~ispc, +# fprintf(fid,'mpiexec -np %i %s/issm.exe %s %s %s ',cluster.np,cluster.codepath,EnumToString(solution),cluster.executionpath,modelname); +# else +# fprintf(fid,'"%s/issm.exe" %s "%s" %s ',cluster.codepath,EnumToString(solution),cluster.executionpath,modelname); +# end +# else +# if ~ispc, +# fprintf(fid,'mpiexec -np %i %s/issm.exe %s %s %s 2> %s.errlog >%s.outlog ',cluster.np,cluster.codepath,EnumToString(solution),cluster.executionpath,modelname,modelname,modelname); +# else +# fprintf(fid,'"%s/issm.exe" %s "%s" %s 2> %s.errlog >%s.outlog ',cluster.codepath,EnumToString(solution),cluster.executionpath,modelname,modelname,modelname); +# end +# end +# else +# if ~ispc, +# %Add --gen-suppressions=all to get suppression lines +# fprintf(fid,'LD_PRELOAD=%s \\\n',cluster.valgrindlib); +# fprintf(fid,'mpiexec -np %i %s --leak-check=full --suppressions=%s %s/issm.exe %s %s %s 2> %s.errlog >%s.outlog ',... +# cluster.np,cluster.valgrind,cluster.valgrindsup, cluster.codepath,EnumToString(solution),cluster.executionpath,modelname,modelname,modelname); +# else +# error('valgrind not supported on windows platforms'); +# end +# end +# +# if isgprof, +# if ~ispc, +# fprintf(fid,'\n gprof %s/issm.exe gmon.out > %s.performance',cluster.codepath,modelname); +# else +# error('gprof not supported on windows platforms'); +# end +# +# end +# +# if ~md.settings.io_gather, +# if ~ispc, +# %concatenate the output files: +# fprintf(fid,'\ncat %s.outbin.* > %s.outbin',modelname,modelname); +# else +# error('iogather not supported on windows platforms'); +# end +# +# end +# +# %close file: +# fclose(fid); +# +# %in interactive mode, create a run file, and errlog and outlog file +# if cluster.interactive, +# fid=fopen([modelname '.errlog'],'w'); fclose(fid); +# fid=fopen([modelname '.outlog'],'w'); fclose(fid); +# end +# +# +# end +# %}}} +# function LaunchQueueJob(cluster,md,options)% {{{1 +# +# if ~ispc, +# %lauch command, to be executed via ssh +# launchcommand=['cd ' cluster.executionpath ' && rm -rf ./' md.private.runtimename ' && mkdir ' md.private.runtimename ... +# ' && cd ' md.private.runtimename ' && mv ../' md.private.runtimename '.tar.gz ./ && tar -zxf ' md.private.runtimename '.tar.gz && source ' md.miscellaneous.name '.queue ']; +# +# if ~strcmpi(options.batch,'yes'), +# +# %compress the files into one zip. +# compressstring=['tar -zcf ' md.private.runtimename '.tar.gz ' md.miscellaneous.name '.bin ' md.miscellaneous.name '.queue ' md.miscellaneous.name '.petsc ']; +# if md.qmu.isdakota, +# compressstring=[compressstring md.miscellaneous.name '.qmu.in']; +# end +# if cluster.interactive, +# compressstring=[compressstring ' ' md.miscellaneous.name '.errlog ' md.miscellaneous.name '.outlog ']; +# end +# system(compressstring); +# +# disp('uploading input file and queueing script'); +# issmscpout(cluster.name,cluster.executionpath,cluster.login,cluster.port,{[md.private.runtimename '.tar.gz']}); +# +# disp('launching solution sequence on remote cluster'); +# issmssh(cluster.name,cluster.login,cluster.port,launchcommand); +# else +# disp('batch mode requested: not launching job interactively'); +# disp('launch solution sequence on remote cluster by hand'); +# end +# else +# %launch right here, do not compress or archive. +# system([md.miscellaneous.name '.bat']); +# end +# +# end %}}} +# function Download(cluster,md)% {{{1 +# +# if ~ispc, +# %some check +# if isempty(md.private.runtimename), +# error('supply runtime name for results to be loaded!'); +# end +# +# %Figure out the directory where all the files are in: +# directory=[cluster.executionpath '/' md.private.runtimename '/']; +# +# %What packages are we picking up from remote cluster +# packages={[md.miscellaneous.name '.outlog'],[md.miscellaneous.name '.errlog']}; +# if md.qmu.isdakota, +# packages{end+1}=[md.miscellaneous.name '.qmu.err']; +# packages{end+1}=[md.miscellaneous.name '.qmu.out']; +# if isfield(md.qmu.params,'tabular_graphics_data'), +# if md.qmu.params.tabular_graphics_data==true, +# packages{end+1}='dakota_tabular.dat'; +# end +# end +# else +# packages{end+1}=[md.miscellaneous.name '.outbin']; +# end +# +# %copy files from cluster to present directory +# issmscpin(cluster.name, cluster.login, cluster.port, directory, packages); +# else +# %do nothing! +# end +# end %}}} +# Index: /issm/trunk/src/m/classes/clusters/none.m =================================================================== --- /issm/trunk/src/m/classes/clusters/none.m (revision 12328) +++ /issm/trunk/src/m/classes/clusters/none.m (revision 12329) @@ -12,18 +12,5 @@ methods function cluster=none(varargin) % {{{1 - cluster=AssignObjectFields(pairoptions(varargin{:}),cluster); - end - %}}} - function disp(cluster) % {{{1 - % display the object - disp(sprintf('cluster class ''%s'' object ''%s'' = ',class(cluster),inputname(1))); - disp(sprintf(' name: %s',cluster.name)); - end - %}}} - function checkconsistency(cluster,md,solution,analyses) % {{{1 - end - %}}} - function BuildQueueScript(cluster,md) % {{{1 - error('none.BuildQueueScript error message: serial cluster cannot build queue script'); + error('Cannot assign md.cluster to ''none'': ISSM is not available in serial model anymore'); end %}}} Index: /issm/trunk/src/m/classes/clusters/none.py =================================================================== --- /issm/trunk/src/m/classes/clusters/none.py (revision 12329) +++ /issm/trunk/src/m/classes/clusters/none.py (revision 12329) @@ -0,0 +1,24 @@ +#NONE cluster class definition +# +# Usage: +# cluster=none(); +# cluster=none('np',3); +# cluster=none('np',3,'login','username'); + +class none: + #properties + def __init__(self): + # {{{ Properties + self.name='none' + #}}} + def __repr__(obj): + # {{{ Display + + string="class 'none' object:" + string="%s\n\n%s"%(string,"%s%s"%(' name: ',obj.name)) + return string + #}}} + def checkconsistency(cluster,md,solution,analyses): + pass + def BuildQueueScript(cluster,md): + raise RuntimeError('none.BuildQueueScript error message: serial cluster cannot build queue script') Index: /issm/trunk/src/m/classes/constants.py =================================================================== --- /issm/trunk/src/m/classes/constants.py (revision 12329) +++ /issm/trunk/src/m/classes/constants.py (revision 12329) @@ -0,0 +1,41 @@ +#module imports +from fielddisplay import fielddisplay + +class constants: + #properties + def __init__(self): + # {{{ Properties + self.g = 0 + self.yts = 0 + self.referencetemperature = 0 + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + string=" constants parameters:" + string="%s\n\n%s"%(string,fielddisplay(obj,"g","gravitational acceleration")) + string="%s\n%s"%(string,fielddisplay(obj,"yts","number of seconds in a year")) + string="%s\n%s"%(string,fielddisplay(obj,"referencetemperature","reference temperature used in the enthalpy model")) + + + return string + #}}} + + def setdefaultparameters(obj): + # {{{setdefaultparameters + + #acceleration due to gravity (m/s^2) + obj.g=9.81 + + #converstion from year to seconds + obj.yts=365*24*3600 + + #the reference temperature for enthalpy model (cf Aschwanden) + obj.referencetemperature=223.15 + + return obj + #}}} + Index: /issm/trunk/src/m/classes/debug.py =================================================================== --- /issm/trunk/src/m/classes/debug.py (revision 12329) +++ /issm/trunk/src/m/classes/debug.py (revision 12329) @@ -0,0 +1,28 @@ +#module imports +from fielddisplay import fielddisplay + +class debug: + #properties + def __init__(self): + # {{{ Properties + self.valgrind=False + self.gprof = False + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + string=" debug parameters:" + + string="%s\n\n%s"%(string,fielddisplay(obj,"valgrind","use Valgrind to debug (0 or 1)")) + string="%s\n%s"%(string,fielddisplay(obj,"gprof","use gnu-profiler to find out where the time is spent")) + return string + #}}} + + def setdefaultparameters(obj): + # {{{setdefaultparameters + return obj + #}}} + Index: /issm/trunk/src/m/classes/diagnostic.py =================================================================== --- /issm/trunk/src/m/classes/diagnostic.py (revision 12329) +++ /issm/trunk/src/m/classes/diagnostic.py (revision 12329) @@ -0,0 +1,99 @@ +#module imports +from fielddisplay import fielddisplay + +class diagnostic: + #properties + def __init__(self): + # {{{ Properties + 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.stokesreconditioning = 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.requested_outputs = float('NaN') + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + + + string='\n Diagnostic solution parameters:' + string="%s\n\n%s"%(string,' Convergence criteria:') + + string="%s\n%s"%(string,fielddisplay(obj,'restol','mechanical equilibrium residual convergence criterion')) + string="%s\n%s"%(string,fielddisplay(obj,'reltol','velocity relative convergence criterion, NaN -> not applied')) + string="%s\n%s"%(string,fielddisplay(obj,'abstol','velocity absolute convergence criterion, NaN -> not applied')) + string="%s\n%s"%(string,fielddisplay(obj,'isnewton','Apply Newton''s method instead of a Picard fixed point method')) + string="%s\n%s"%(string,fielddisplay(obj,'maxiter','maximum number of nonlinear iterations')) + string="%s\n%s"%(string,fielddisplay(obj,'viscosity_overshoot','over-shooting constant new=new+C*(new-old)')) + + string="%s\n%s"%(string,' boundary conditions:') + + string="%s\n%s"%(string,fielddisplay(obj,'spcvx','x-axis velocity constraint (NaN means no constraint)')) + string="%s\n%s"%(string,fielddisplay(obj,'spcvy','y-axis velocity constraint (NaN means no constraint)')) + string="%s\n%s"%(string,fielddisplay(obj,'spcvz','z-axis velocity constraint (NaN means no constraint)')) + string="%s\n%s"%(string,fielddisplay(obj,'icefront','segments on ice front list (last column 0-> Air, 1-> Water, 2->Ice')) + + string="%s\n%s"%(string,' Rift options:') + string="%s\n%s"%(string,fielddisplay(obj,'rift_penalty_threshold','threshold for instability of mechanical constraints')) + string="%s\n%s"%(string,fielddisplay(obj,'rift_penalty_lock','number of iterations before rift penalties are locked')) + + string="%s\n%s"%(string,' Penalty options:') + string="%s\n%s"%(string,fielddisplay(obj,'penalty_factor','offset used by penalties: penalty = Kmax*10^offset')) + string="%s\n%s"%(string,fielddisplay(obj,'vertex_pairing','pairs of vertices that are penalized')) + + string="%s\n%s"%(string,' Other:') + string="%s\n%s"%(string,fielddisplay(obj,'shelf_dampening','use dampening for floating ice ? Only for Stokes model')) + string="%s\n%s"%(string,fielddisplay(obj,'stokesreconditioning','multiplier for incompressibility equation. Only for Stokes model')) + string="%s\n%s"%(string,fielddisplay(obj,'referential','local referential')) + string="%s\n%s"%(string,fielddisplay(obj,'requested_outputs','additional outputs requested')) + + return string + #}}} + + def setdefaultparameters(obj): + # {{{setdefaultparameters + #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.stokesreconditioning=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 + + return obj + #}}} Index: /issm/trunk/src/m/classes/flaim.py =================================================================== --- /issm/trunk/src/m/classes/flaim.py (revision 12329) +++ /issm/trunk/src/m/classes/flaim.py (revision 12329) @@ -0,0 +1,50 @@ +#module imports +from fielddisplay import fielddisplay + +class flaim: + #properties + def __init__(self): + # {{{ Properties + self.targets = '' + self.tracks = '' + self.flightreqs = {} + self.criterion = float('NaN') + self.gridsatequator = 200000 + self.usevalueordering = True + self.split_antimeridian = True + self.solution = '' + self.quality = 0 + self.path_optimize = False + self.opt_ndir = 1 + self.opt_dist = 25 + self.opt_niter = 30000 + #}}} + def __repr__(obj): + # {{{ Displa + string=' FLAIM - Flight Line Adaptation using Ice sheet Modeling:' + + string="%s\n\n%s"%(string,' Input:') + string="%s\n%s"%(string,fielddisplay(obj,'targets' ,'name of kml output targets file ')) + string="%s\n%s"%(string,fielddisplay(obj,'tracks' ,'name of kml input tracks file ')) + string="%s\n%s"%(string,fielddisplay(obj,'flightreqs' ,'structure of kml flight requirements (not used yet)')) + string="%s\n%s"%(string,fielddisplay(obj,'criterion' ,'element or nodal criterion for flight path evaluation (metric)')) + + string="%s\n\n%s"%(string,' Arguments:') + string="%s\n%s"%(string,fielddisplay(obj,'gridsatequator' ,'number of grids at equator (determines resolution)')) + string="%s\n%s"%(string,fielddisplay(obj,'usevalueordering' ,'flag to consider target values for flight path evaluation')) + string="%s\n%s"%(string,fielddisplay(obj,'split_antimeridian' ,'flag to split polygons on the antimeridian')) + + string="%s\n\n%s"%(string,' Optimization:') + string="%s\n%s"%(string,fielddisplay(obj,'path_optimize' ,'optimize? (default false)')) + string="%s\n%s"%(string,fielddisplay(obj,'opt_ndir' ,['number of directions to test when moving a point. If this value = 1, a random direction is tested.',\ + 'A value > 1 results in directions equally spaced from [0, 2*PI] being tested.',\ + 'For example, 4 would result in directions [0, PI/2, PI, 3PI/2].'])) + string="%s\n%s"%(string,fielddisplay(obj,'opt_dist' ,'specifies the distance in km (default 25) to move a randomly selected path point on each iteration')) + string="%s\n%s"%(string,fielddisplay(obj,'opt_niter' ,['number of iterations (default 30,000) to run for flightplan optimization',\ + 'i.e. the number of times to randomly select a point and move it.'])) + + string="%s\n\n%s"%(string,' Output:') + string="%s\n%s"%(string,fielddisplay(obj,'solution' ,'name of kml solution file')) + string="%s\n%s"%(string,fielddisplay(obj,'quality' ,'quality of kml solution')) + return string + #}}} Index: /issm/trunk/src/m/classes/flowequation.py =================================================================== --- /issm/trunk/src/m/classes/flowequation.py (revision 12329) +++ /issm/trunk/src/m/classes/flowequation.py (revision 12329) @@ -0,0 +1,41 @@ +#module imports +from fielddisplay import fielddisplay + +class flowequation: + #properties + def __init__(self): + # {{{ Properties + + self.ismacayealpattyn = 0; + self.ishutter = 0; + self.isstokes = 0; + self.vertex_equation = float('NaN') + self.element_equation = float('NaN') + self.bordermacayeal = float('NaN') + self.borderpattyn = float('NaN') + self.borderstokes = float('NaN') + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + string=' flow equation parameters:' + + string="%s\n\n%s"%(string,fielddisplay(obj,'ismacayealpattyn','is the macayeal or pattyn approximation used ?')) + string="%s\n%s"%(string,fielddisplay(obj,'ishutter','is the shallow ice approximation used ?')) + string="%s\n%s"%(string,fielddisplay(obj,'isstokes','are the Full-Stokes equations used ?')) + string="%s\n%s"%(string,fielddisplay(obj,'vertex_equation','flow equation for each vertex')) + string="%s\n%s"%(string,fielddisplay(obj,'element_equation','flow equation for each element')) + string="%s\n%s"%(string,fielddisplay(obj,'bordermacayeal','vertices on MacAyeal''s border (for tiling)')) + string="%s\n%s"%(string,fielddisplay(obj,'borderpattyn','vertices on Pattyn''s border (for tiling)')) + string="%s\n%s"%(string,fielddisplay(obj,'borderstokes','vertices on Stokes'' border (for tiling)')) + return string + #}}} + + def setdefaultparameters(obj): + # {{{setdefaultparameters + return obj + #}}} + Index: /issm/trunk/src/m/classes/friction.py =================================================================== --- /issm/trunk/src/m/classes/friction.py (revision 12329) +++ /issm/trunk/src/m/classes/friction.py (revision 12329) @@ -0,0 +1,29 @@ +#module imports +from fielddisplay import fielddisplay + +class friction: + #properties + def __init__(self): + # {{{ Properties + self.coefficient = float('NaN') + self.p = float('NaN') + self.q = float('NaN') + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + string="Sigma= drag^2 * Neff ^r * u ^s, with Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p" + string="%s\n\n%s"%(string,fielddisplay(obj,"coefficient","friction coefficient [SI]")) + string="%s\n%s"%(string,fielddisplay(obj,"p","p exponent")) + string="%s\n%s"%(string,fielddisplay(obj,"q","q exponent")) + return string + #}}} + + def setdefaultparameters(obj): + # {{{setdefaultparameters + return obj + #}}} + Index: /issm/trunk/src/m/classes/geometry.py =================================================================== --- /issm/trunk/src/m/classes/geometry.py (revision 12329) +++ /issm/trunk/src/m/classes/geometry.py (revision 12329) @@ -0,0 +1,35 @@ +#module imports +from fielddisplay import fielddisplay + +class geometry: + #properties + def __init__(self): + # {{{ Properties + self.surface = float('NaN') + self.thickness = float('NaN') + self.bed = float('NaN') + self.bathymetry = float('NaN') + self.hydrostatic_ratio = float('NaN') + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + + string=" geometry parameters:" + + string="%s\n\n%s"%(string,fielddisplay(obj,'surface','surface elevation')) + string="%s\n%s"%(string,fielddisplay(obj,'thickness','ice thickness')) + string="%s\n%s"%(string,fielddisplay(obj,'bed','bed elevation')) + string="%s\n%s"%(string,fielddisplay(obj,'bathymetry','bathymetry elevation')) + string="%s\n%s"%(string,fielddisplay(obj,'hydrostatic_ratio','coefficient for ice shelves'' thickness correction: hydrostatic_ratio H_obs+ (1-hydrostatic_ratio) H_hydro')) + return string + #}}} + + def setdefaultparameters(obj): + # {{{setdefaultparameters + return obj + #}}} + Index: /issm/trunk/src/m/classes/groundingline.m =================================================================== --- /issm/trunk/src/m/classes/groundingline.m (revision 12328) +++ /issm/trunk/src/m/classes/groundingline.m (revision 12329) @@ -37,5 +37,5 @@ end pos=find(md.mask.vertexongroundedice); - if any(md.geometry.bed(pos)-md.geometry.bathymetry(pos)), + if any(abs(md.geometry.bed(pos)-md.geometry.bathymetry(pos))>10^-10), checkmessage(['bathymetry not equal to bed on grounded ice !']); end Index: /issm/trunk/src/m/classes/groundingline.py =================================================================== --- /issm/trunk/src/m/classes/groundingline.py (revision 12329) +++ /issm/trunk/src/m/classes/groundingline.py (revision 12329) @@ -0,0 +1,34 @@ +#module imports +from fielddisplay import fielddisplay + +class groundingline: + #properties + def __init__(self): + # {{{ Properties + self.migration='' + self.melting_rate=float('NaN') + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + string=' grounding line solution parameters:' + + string="%s\n\n%s"%(string,fielddisplay(obj,'migration','type of grounding line migration: ''SoftMigration'',''AgressiveMigration'' or ''None''')) + string="%s\n%s"%(string,fielddisplay(obj,'melting_rate','melting rate applied when previously grounded parts start floating')) + return string + #}}} + def setdefaultparameters(obj): + # {{{setdefaultparameters + + #Type of migration + obj.migration='None' + + #basal melting rate correction: + obj.melting_rate=0; + + return obj + #}}} + Index: /issm/trunk/src/m/classes/hydrology.py =================================================================== --- /issm/trunk/src/m/classes/hydrology.py (revision 12329) +++ /issm/trunk/src/m/classes/hydrology.py (revision 12329) @@ -0,0 +1,50 @@ +#module imports +from fielddisplay import fielddisplay + +class hydrology: + #properties + def __init__(self): + # {{{ Properties + self.spcwatercolumn = float('NaN') + self.n = 0 + self.CR = 0 + self.p = 0 + self.q = 0 + self.kn = 0 + self.stabilization = 0 + + #set defaults + self.setdefaultparameters() + + + #}}} + def __repr__(obj): + # {{{ Display + + string=' hydrology solution parameters:' + string="%s\n\n%s"%(string,fielddisplay(obj,'spcwatercolumn','water thickness constraints (NaN means no constraint)')) + string="%s\n%s"%(string,fielddisplay(obj,'n','Manning roughness coefficient')) + string="%s\n%s"%(string,fielddisplay(obj,'CR','tortuosity parameter')) + string="%s\n%s"%(string,fielddisplay(obj,'p','dimensionless exponent in Manning velocity formula')) + string="%s\n%s"%(string,fielddisplay(obj,'q','dimensionless exponent in Manning velocity formula')) + string="%s\n%s"%(string,fielddisplay(obj,'kn','parameter in effective pressure formula')) + string="%s\n%s"%(string,fielddisplay(obj,'stabilization','artificial diffusivity (default is 1). can be more than 1 to increase diffusivity.')) + return string + #}}} + + def setdefaultparameters(obj): + # {{{setdefaultparameters + + #Parameters from Johnson's 2002 thesis, section 3.5.4 + obj.n=.02 + obj.CR=.01 + obj.p=2 + obj.q=1 + obj.kn=0 + + #Type of stabilization to use 0:nothing 1:artificial_diffusivity + obj.stabilization=1 + + return obj + #}}} + Index: /issm/trunk/src/m/classes/initialization.py =================================================================== --- /issm/trunk/src/m/classes/initialization.py (revision 12329) +++ /issm/trunk/src/m/classes/initialization.py (revision 12329) @@ -0,0 +1,42 @@ +#module imports +from fielddisplay import fielddisplay + +class initialization: + #properties + def __init__(self): + # {{{ Properties + + self.vx = float('NaN') + self.vy = float('NaN') + self.vz = float('NaN') + self.vel = float('NaN') + self.pressure = float('NaN') + self.temperature = float('NaN') + self.watercolumn = float('NaN') + self.waterfraction = float('NaN') + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + string=' initial field values:' + + string="%s\n%s"%(string,fielddisplay(obj,'vx','x component of velocity')) + string="%s\n%s"%(string,fielddisplay(obj,'vy','y component of velocity')) + string="%s\n%s"%(string,fielddisplay(obj,'vz','z component of velocity')) + string="%s\n%s"%(string,fielddisplay(obj,'vel','velocity norm')) + string="%s\n%s"%(string,fielddisplay(obj,'pressure','pressure field')) + string="%s\n%s"%(string,fielddisplay(obj,'temperature','temperature in Kelvins')) + string="%s\n%s"%(string,fielddisplay(obj,'watercolumn','thickness of subglacial water')) + string="%s\n%s"%(string,fielddisplay(obj,'waterfraction','fraction of water in the ice')) + + return string + #}}} + + def setdefaultparameters(obj): + # {{{setdefaultparameters + return obj + #}}} + Index: /issm/trunk/src/m/classes/inversion.py =================================================================== --- /issm/trunk/src/m/classes/inversion.py (revision 12329) +++ /issm/trunk/src/m/classes/inversion.py (revision 12329) @@ -0,0 +1,104 @@ +#module imports +from fielddisplay import fielddisplay + +class inversion: + #properties + def __init__(self): + # {{{ Properties + self.iscontrol = 0 + self.tao = 0 + self.incomplete_adjoint = 0 + self.control_parameters = float('NaN') + self.nsteps = 0 + self.maxiter_per_step = float('NaN') + self.cost_functions = float('NaN') + self.cost_functions_coefficients = float('NaN') + self.gradient_scaling = float('NaN') + self.cost_function_threshold = 0 + self.min_parameters = float('NaN') + self.max_parameters = float('NaN') + self.step_threshold = float('NaN') + self.gradient_only = 0 + self.vx_obs = float('NaN') + self.vy_obs = float('NaN') + self.vz_obs = float('NaN') + self.vel_obs = float('NaN') + self.thickness_obs = float('NaN') + #}}} + def __repr__(obj): + # {{{ Display + string='\n Inversion parameters:' + string="%s\n%s"%(string,fielddisplay(obj,'iscontrol','is inversion activated?')) + string="%s\n%s"%(string,fielddisplay(obj,'incomplete_adjoint','do we assume linear viscosity?')) + string="%s\n%s"%(string,fielddisplay(obj,'control_parameters','parameter where inverse control is carried out; ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''}')) + string="%s\n%s"%(string,fielddisplay(obj,'nsteps','number of optimization searches')) + string="%s\n%s"%(string,fielddisplay(obj,'cost_functions','indicate the type of response for each optimization step')) + string="%s\n%s"%(string,fielddisplay(obj,'cost_functions_coefficients','cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter')) + string="%s\n%s"%(string,fielddisplay(obj,'cost_function_threshold','misfit convergence criterion. Default is 1%, NaN if not applied')) + string="%s\n%s"%(string,fielddisplay(obj,'maxiter_per_step','maximum iterations during each optimization step')) + string="%s\n%s"%(string,fielddisplay(obj,'gradient_scaling','scaling factor on gradient direction during optimization, for each optimization step')) + string="%s\n%s"%(string,fielddisplay(obj,'step_threshold','decrease threshold for misfit, default is 30%')) + string="%s\n%s"%(string,fielddisplay(obj,'min_parameters','absolute minimum acceptable value of the inversed parameter on each vertex')) + string="%s\n%s"%(string,fielddisplay(obj,'max_parameters','absolute maximum acceptable value of the inversed parameter on each vertex')) + string="%s\n%s"%(string,fielddisplay(obj,'gradient_only','stop control method solution at gradient')) + string="%s\n%s"%(string,fielddisplay(obj,'vx_obs','observed velocity x component [m/a]')) + string="%s\n%s"%(string,fielddisplay(obj,'vy_obs','observed velocity y component [m/a]')) + string="%s\n%s"%(string,fielddisplay(obj,'vel_obs','observed velocity magnitude [m/a]')) + string="%s\n%s"%(string,fielddisplay(obj,'thickness_obs','observed thickness [m]')) + string="%s\n%s"%(string,'Available cost functions:') + string="%s\n%s"%(string,' 101: SurfaceAbsVelMisfit') + string="%s\n%s"%(string,' 102: SurfaceRelVelMisfit') + string="%s\n%s"%(string,' 103: SurfaceLogVelMisfit') + string="%s\n%s"%(string,' 104: SurfaceLogVxVyMisfit') + string="%s\n%s"%(string,' 105: SurfaceAverageVelMisfit') + string="%s\n%s"%(string,' 201: ThicknessAbsMisfit') + string="%s\n%s"%(string,' 501: DragCoefficientAbsGradient') + string="%s\n%s"%(string,' 502: RheologyBbarAbsGradient') + string="%s\n%s"%(string,' 503: ThicknessAbsGradient') + return string + #}}} + + def setdefaultparameters(obj): + # {{{setdefaultparameters + + #default is incomplete adjoint for now + obj.incomplete_adjoint=1 + + #parameter to be inferred by control methods (only + #drag and B are supported yet) + obj.control_parameters=['FrictionCoefficient'] + + #number of steps in the control methods + obj.nsteps=20 + + #maximum number of iteration in the optimization algorithm for + #each step + obj.maxiter_per_step=20*ones(obj.nsteps) + + #the inversed parameter is updated as follows: + #new_par=old_par + gradient_scaling(n)*C*gradient with C in [0 1]; + #usually the gradient_scaling must be of the order of magnitude of the + #inversed parameter (10^8 for B, 50 for drag) and can be decreased + #after the first iterations + obj.gradient_scaling=50*ones(obj.nsteps) + + #several responses can be used: + obj.cost_functions=101*ones(obj.nsteps) + + #step_threshold is used to speed up control method. When + #misfit(1)/misfit(0) < obj.step_threshold, we go directly to + #the next step + obj.step_threshold=.7*ones(obj.nsteps) #30 per cent decrement + + #stop control solution at the gradient computation and return it? + obj.gradient_only=0 + + #cost_function_threshold is a criteria to stop the control methods. + #if J[n]-J[n-1]/J[n] < criteria, the control run stops + #NaN if not applied + obj.cost_function_threshold=NaN #not activated + + return obj + #}}} + + Index: /issm/trunk/src/m/classes/mask.py =================================================================== --- /issm/trunk/src/m/classes/mask.py (revision 12329) +++ /issm/trunk/src/m/classes/mask.py (revision 12329) @@ -0,0 +1,36 @@ +#module imports +from fielddisplay import fielddisplay + +class mask: + #properties + def __init__(self): + # {{{ Properties + self.elementonfloatingice = float('NaN') + self.elementongroundedice = float('NaN') + self.elementonwater = float('NaN') + self.vertexonfloatingice = float('NaN') + self.vertexongroundedice = float('NaN') + self.vertexonwater = float('NaN') + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + + string=""; + string="%s\n%s"%(string,fielddisplay(obj,"elementonfloatingice","element on floating ice flags list")) + string="%s\n%s"%(string,fielddisplay(obj,"vertexonfloatingice","vertex on floating ice flags list")) + string="%s\n%s"%(string,fielddisplay(obj,"elementongroundedice","element on grounded ice list")) + string="%s\n%s"%(string,fielddisplay(obj,"vertexongroundedice","vertex on grounded ice flags list")) + string="%s\n%s"%(string,fielddisplay(obj,"elementonwater","element on water flags list")) + string="%s\n%s"%(string,fielddisplay(obj,"vertexonwater","vertex on water flags list")) + return string + #}}} + + def setdefaultparameters(obj): + # {{{setdefaultparameters + return obj + #}}} + Index: /issm/trunk/src/m/classes/materials.py =================================================================== --- /issm/trunk/src/m/classes/materials.py (revision 12329) +++ /issm/trunk/src/m/classes/materials.py (revision 12329) @@ -0,0 +1,80 @@ +#module imports +from fielddisplay import fielddisplay + +class materials: + #properties + def __init__(self): + # {{{ Properties + self.rho_ice = 0; + self.rho_water = 0; + self.mu_water = 0; + self.heatcapacity = 0; + self.latentheat = 0; + self.thermalconductivity = 0; + self.meltingpoint = 0; + self.beta = 0; + self.mixed_layer_capacity = 0; + self.thermal_exchange_velocity = 0; + self.rheology_B = float('NaN') + self.rheology_n = float('NaN') + self.rheology_law = ""; + + self.setdefaultparameters() + #}}} + def __repr__(obj): + # {{{ Display + string=" Materials:" + + string="%s\n\n%s"%(string,fielddisplay(obj,"rho_ice","ice density [kg/m^3]")) + string="%s\n%s"%(string,fielddisplay(obj,"rho_water","water density [kg/m^3]")) + string="%s\n%s"%(string,fielddisplay(obj,"mu_water","water viscosity [N s/m^2]")) + string="%s\n%s"%(string,fielddisplay(obj,"heatcapacity","heat capacity [J/kg/K]")) + string="%s\n%s"%(string,fielddisplay(obj,"thermalconductivity","ice thermal conductivity [W/m/K]")) + string="%s\n%s"%(string,fielddisplay(obj,"meltingpoint","melting point of ice at 1atm in K")) + string="%s\n%s"%(string,fielddisplay(obj,"latentheat","latent heat of fusion [J/m^3]")) + string="%s\n%s"%(string,fielddisplay(obj,"beta","rate of change of melting point with pressure [K/Pa]")) + string="%s\n%s"%(string,fielddisplay(obj,"mixed_layer_capacity","mixed layer capacity [W/kg/K]")) + string="%s\n%s"%(string,fielddisplay(obj,"thermal_exchange_velocity","thermal exchange velocity [m/s]")) + string="%s\n%s"%(string,fielddisplay(obj,"rheology_B","flow law parameter [Pa/s^(1/n)]")) + string="%s\n%s"%(string,fielddisplay(obj,"rheology_n","Glen""s flow law exponent")) + string="%s\n%s"%(string,fielddisplay(obj,"rheology_law","law for the temperature dependance of the rheology: ""None"", ""Paterson"" or ""Arrhenius""")) + + return string + #}}} + def setdefaultparameters(obj): + # {{{setdefaultparameters + #ice density (kg/m^3) + obj.rho_ice=917 + + #water density (kg/m^3) + obj.rho_water=1023 + + #water viscosity (N.s/m^2) + obj.mu_water=0.001787 + + #ice heat capacity cp (J/kg/K) + obj.heatcapacity=2093 + + #ice latent heat of fusion L (J/kg) + obj.latentheat=3.34*10**5 + + #ice thermal conductivity (W/m/K) + obj.thermalconductivity=2.4 + + #the melting point of ice at 1 atmosphere of pressure in K + obj.meltingpoint=273.15 + + #rate of change of melting point with pressure (K/Pa) + obj.beta=9.8*10**-8 + + #mixed layer (ice-water interface) heat capacity (J/kg/K) + obj.mixed_layer_capacity=3974 + + #thermal exchange velocity (ice-water interface) (m/s) + obj.thermal_exchange_velocity=1.00*10**-4 + + #Rheology law: what is the temperature dependence of B with T + #available: none, paterson and arrhenius + obj.rheology_law='Paterson' + return obj + #}}} Index: /issm/trunk/src/m/classes/mesh.py =================================================================== --- /issm/trunk/src/m/classes/mesh.py (revision 12329) +++ /issm/trunk/src/m/classes/mesh.py (revision 12329) @@ -0,0 +1,125 @@ +#module imports +from fielddisplay import fielddisplay + +class mesh: + #properties + def __init__(self): + # {{{ Properties + self.x = float('NaN'); + self.y = float('NaN'); + self.z = float('NaN'); + self.elements = float('NaN'); + self.dimension = 0; + self.numberoflayers = 0; + self.numberofelements = 0; + self.numberofvertices = 0; + self.numberofedges = 0; + + self.lat = float('NaN'); + self.long = float('NaN'); + self.hemisphere = float('NaN'); + + self.elementonbed = float('NaN'); + self.elementonsurface = float('NaN'); + self.vertexonbed = float('NaN'); + self.vertexonsurface = float('NaN'); + self.lowerelements = float('NaN'); + self.lowervertex = float('NaN'); + self.upperelements = float('NaN'); + self.uppervertex = float('NaN'); + self.vertexonboundary = float('NaN'); + + self.edges = float('NaN'); + self.segments = float('NaN'); + self.segmentmarkers = float('NaN'); + self.vertexconnectivity = float('NaN'); + self.elementconnectivity = float('NaN'); + self.average_vertex_connectivity = 0; + + self.x2d = float('NaN'); + self.y2d = float('NaN'); + self.elements2d = float('NaN'); + self.numberofvertices2d = 0; + self.numberofelements2d = 0; + + self.extractedvertices = float('NaN'); + self.extractedelements = float('NaN'); + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + + if obj.dimension==3: + string="\n%s"%(" Elements and vertices of the original 2d mesh:") + + string="%s\n%s"%(string,fielddisplay(obj,"numberofelements2d","number of elements")) + string="%s\n%s"%(string,fielddisplay(obj,"numberofvertices2d","number of vertices")) + string="%s\n%s"%(string,fielddisplay(obj,"elements2d","index into (x,y,z), coordinates of the vertices")) + string="%s\n%s"%(string,fielddisplay(obj,"x2d","vertices x coordinate")) + string="%s\n%s"%(string,fielddisplay(obj,"y2d","vertices y coordinate")) + + string="%s\n%s" %(string,"Elements and vertices of the extruded 3d mesh:") + else: + string="\n%s"%(" Elements and vertices:") + + string="%s\n%s"%(string,fielddisplay(obj,"numberofelements","number of elements")) + + + + + string="%s\n%s"%(string,fielddisplay(obj,"numberofvertices","number of vertices")) + string="%s\n%s"%(string,fielddisplay(obj,"elements","index into (x,y,z), coordinates of the vertices")) + string="%s\n%s"%(string,fielddisplay(obj,"x","vertices x coordinate")) + string="%s\n%s"%(string,fielddisplay(obj,"y","vertices y coordinate")) + string="%s\n%s"%(string,fielddisplay(obj,"z","vertices z coordinate")) + string="%s\n%s"%(string,fielddisplay(obj,"edges","edges of the 2d mesh (vertex1 vertex2 element1 element2)")) + string="%s\n%s"%(string,fielddisplay(obj,"numberofedges","number of edges of the 2d mesh")) + + string="%s%s"%(string,"\n Properties:") + + string="%s\n%s"%(string,fielddisplay(obj,"dimension","mesh dimension (2d or 3d)")) + string="%s\n%s"%(string,fielddisplay(obj,"numberoflayers","number of extrusion layers")) + string="%s\n%s"%(string,fielddisplay(obj,"vertexonbed","lower vertices flags list")) + string="%s\n%s"%(string,fielddisplay(obj,"elementonbed","lower elements flags list")) + string="%s\n%s"%(string,fielddisplay(obj,"vertexonsurface","upper vertices flags list")) + string="%s\n%s"%(string,fielddisplay(obj,"elementonsurface","upper elements flags list")) + string="%s\n%s"%(string,fielddisplay(obj,"uppervertex","upper vertex list (NaN for vertex on the upper surface)")) + string="%s\n%s"%(string,fielddisplay(obj,"upperelements","upper element list (NaN for element on the upper layer)")) + string="%s\n%s"%(string,fielddisplay(obj,"lowervertex","lower vertex list (NaN for vertex on the lower surface)")) + string="%s\n%s"%(string,fielddisplay(obj,"lowerelements","lower element list (NaN for element on the lower layer")) + string="%s\n%s"%(string,fielddisplay(obj,"vertexonboundary","vertices on the boundary of the domain flag list")) + + string="%s\n%s"%(string,fielddisplay(obj,"segments","edges on domain boundary (vertex1 vertex2 element)")) + string="%s\n%s"%(string,fielddisplay(obj,"segmentmarkers","number associated to each segment")) + string="%s\n%s"%(string,fielddisplay(obj,"vertexconnectivity","list of vertices connected to vertex_i")) + string="%s\n%s"%(string,fielddisplay(obj,"elementconnectivity","list of vertices connected to element_i")) + string="%s\n%s"%(string,fielddisplay(obj,"average_vertex_connectivity","average number of vertices connected to one vertex")) + + string="%s%s"%(string,"\n Extracted model:") + + string="%s\n%s"%(string,fielddisplay(obj,"extractedvertices","vertices extracted from the model")) + string="%s\n%s"%(string,fielddisplay(obj,"extractedelements","elements extracted from the model")) + + string="%s%s"%(string,"\n Projection:") + string="%s\n%s"%(string,fielddisplay(obj,"lat","vertices latitude")) + string="%s\n%s"%(string,fielddisplay(obj,"long","vertices longitude")) + string="%s\n%s"%(string,fielddisplay(obj,"hemisphere","Indicate hemisphere ""n"" or ""s"" ")) + return string + #}}} + + def setdefaultparameters(obj): + # {{{setdefaultparameters + + #the connectivity is the avergaded number of nodes linked to a + #given node through an edge. This connectivity is used to initially + #allocate memory to the stiffness matrix. A value of 16 seems to + #give a good memory/time ration. This value can be checked in + #trunk/test/Miscellaneous/runme.m + obj.average_vertex_connectivity=25 + + return obj + #}}} + Index: /issm/trunk/src/m/classes/miscellaneous.py =================================================================== --- /issm/trunk/src/m/classes/miscellaneous.py (revision 12329) +++ /issm/trunk/src/m/classes/miscellaneous.py (revision 12329) @@ -0,0 +1,30 @@ +#module imports +from fielddisplay import fielddisplay + +class miscellaneous: + #properties + def __init__(self): + # {{{ Properties + self.notes = '' + self.name = '' + self.dummy = {} + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + string=' miscellaneous parameters:' + + string="%s\n\n%s"%(string,fielddisplay(obj,'notes','notes in a cell of strings')) + string="%s\n%s"%(string,fielddisplay(obj,'name','model name')) + string="%s\n%s"%(string,fielddisplay(obj,'dummy','empty field to store some data')) + return string + #}}} + + def setdefaultparameters(obj): + # {{{setdefaultparameters + return obj + #}}} + Index: /issm/trunk/src/m/classes/model.py =================================================================== --- /issm/trunk/src/m/classes/model.py (revision 12329) +++ /issm/trunk/src/m/classes/model.py (revision 12329) @@ -0,0 +1,114 @@ +#module imports {{{ +from mesh import mesh +from mask import mask +from geometry import geometry +from constants import constants +from surfaceforcings import surfaceforcings +from basalforcings import basalforcings +from materials import materials +from friction import friction +from flowequation import flowequation +from timestepping import timestepping +from initialization import initialization +from rifts import rifts +from debug import debug +from verbose import verbose +from settings import settings +from solver import solver +from none import none +from balancethickness import balancethickness +from diagnostic import diagnostic +from groundingline import groundingline +from hydrology import hydrology +from prognostic import prognostic +from thermal import thermal +from steadystate import steadystate +from transient import transient +from autodiff import autodiff +from flaim import flaim +from inversion import inversion +from qmu import qmu +from radaroverlay import radaroverlay +from miscellaneous import miscellaneous +from private import private +#}}} +class model: + #properties + def __init__(self): + # {{{ Properties + self.mesh = mesh() + self.mask = mask() + self.geometry = geometry() + self.constants = constants() + self.surfaceforcings = surfaceforcings() + self.basalforcings = basalforcings() + self.materials = materials() + self.friction = friction() + self.flowequation = flowequation() + self.timestepping = timestepping() + self.initialization = initialization() + self.rifts = rifts() + + self.debug = debug() + self.verbose = verbose() + self.settings = settings() + self.solver = solver() + self.cluster = none() + + self.balancethickness = balancethickness() + self.diagnostic = diagnostic() + self.groundingline = groundingline() + self.hydrology = hydrology() + self.prognostic = prognostic() + self.thermal = thermal() + self.steadystate = steadystate() + self.transient = transient() + + self.autodiff = autodiff() + self.flaim = flaim() + self.inversion = inversion() + self.qmu = qmu() + + self.results = []; + self.radaroverlay = radaroverlay() + self.miscellaneous = miscellaneous() + self.private = private() + #}}} + def __repr__(obj): + # {{{ Display + + #print "Here %s the number: %d" % ("is", 37) + string="%19s: %-22s -- %s" % ("mesh","[%s,%s]" % ("1x1",obj.mesh.__class__.__name__),"mesh properties") + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("mask","[%s,%s]" % ("1x1",obj.mask.__class__.__name__),"defines grounded and floating elements")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("geometry","[%s,%s]" % ("1x1",obj.geometry.__class__.__name__),"surface elevation, bedrock topography, ice thickness,...")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("constants","[%s,%s]" % ("1x1",obj.constants.__class__.__name__),"physical constants")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("surfaceforcings","[%s,%s]" % ("1x1",obj.surfaceforcings.__class__.__name__),"surface forcings")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("basalforcings","[%s,%s]" % ("1x1",obj.basalforcings.__class__.__name__),"bed forcings")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("materials","[%s,%s]" % ("1x1",obj.materials.__class__.__name__),"material properties")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("friction","[%s,%s]" % ("1x1",obj.friction.__class__.__name__),"basal friction/drag properties")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("flowequation","[%s,%s]" % ("1x1",obj.flowequation.__class__.__name__),"flow equations")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("timestepping","[%s,%s]" % ("1x1",obj.timestepping.__class__.__name__),"time stepping for transient models")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("initialization","[%s,%s]" % ("1x1",obj.initialization.__class__.__name__),"initial guess/state")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("rifts","[%s,%s]" % ("1x1",obj.rifts.__class__.__name__),"rifts properties'")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("debug","[%s,%s]" % ("1x1",obj.debug.__class__.__name__),"debugging tools (valgrind, gprof")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("verbose","[%s,%s]" % ("1x1",obj.verbose.__class__.__name__),"verbosity level in solve")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("settings","[%s,%s]" % ("1x1",obj.settings.__class__.__name__),"settings properties")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("solver","[%s,%s]" % ("1x1",obj.solver.__class__.__name__),"PETSc options for each solution'")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("cluster","[%s,%s]" % ("1x1",obj.cluster.__class__.__name__),"cluster parameters (number of cpus...)")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("balancethickness","[%s,%s]" % ("1x1",obj.balancethickness.__class__.__name__),"parameters for balancethickness solution")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("diagnostic","[%s,%s]" % ("1x1",obj.diagnostic.__class__.__name__),"parameters for diagnostic solution")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("groundingline","[%s,%s]" % ("1x1",obj.groundingline.__class__.__name__),"parameters for groundingline solution")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("hydrology","[%s,%s]" % ("1x1",obj.hydrology.__class__.__name__),"parameters for hydrology solution")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("prognostic","[%s,%s]" % ("1x1",obj.prognostic.__class__.__name__),"parameters for prognostic solution")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("thermal","[%s,%s]" % ("1x1",obj.thermal.__class__.__name__),"parameters for thermal solution")) + 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" % ("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")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("inversion","[%s,%s]" % ("1x1",obj.inversion.__class__.__name__),"parameters for inverse methods")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("qmu","[%s,%s]" % ("1x1",obj.qmu.__class__.__name__),"dakota properties")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("results","[%s,%s]" % ("1x1",obj.results.__class__.__name__),"model results'")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("radaroverlay","[%s,%s]" % ("1x1",obj.radaroverlay.__class__.__name__),"radar image for plot overlay")) + string="%s\n%s" % (string,"%19s: %-22s -- %s" % ("miscellaneous","[%s,%s]" % ("1x1",obj.miscellaneous.__class__.__name__),"miscellaneous fields")) + return string; + #}}} Index: /issm/trunk/src/m/classes/model/model.m =================================================================== --- /issm/trunk/src/m/classes/model/model.m (revision 12328) +++ /issm/trunk/src/m/classes/model/model.m (revision 12329) @@ -379,5 +379,5 @@ md.solver = addoptions(md.solver,DiagnosticVertAnalysisEnum,iluasmoptions); end - md.cluster = none(); + md.cluster = generic(); md.balancethickness = balancethickness(); md.diagnostic = diagnostic(); Index: /issm/trunk/src/m/classes/modellist.m =================================================================== --- /issm/trunk/src/m/classes/modellist.m (revision 12328) +++ /issm/trunk/src/m/classes/modellist.m (revision 12329) @@ -7,5 +7,5 @@ properties (SetAccess=public) models = cell(0,1); - cluster = none(); + cluster = generic(); end methods Index: /issm/trunk/src/m/classes/pairoptions.py =================================================================== --- /issm/trunk/src/m/classes/pairoptions.py (revision 12329) +++ /issm/trunk/src/m/classes/pairoptions.py (revision 12329) @@ -0,0 +1,38 @@ +class pairoptions: + #properties + def __init__(self,*args): + # {{{ Properties + if len(args)%2==1: + raise RuntimeError('pairoption error message: an even number of options is required') + + #create a pairoption object + if len(args)==0: + self.list=[] + else: + self.list=[] + for i in range(int(round(len(args)/2))): + if isinstance(args[2*i],str): + self.list.append([args[2*i],args[2*i+1]]) + else: + #option is not a string, ignore it + print("%s%i%s"%('buildlist info: option number ',i,' is not a string, it will be ignored')) + continue + + #}}} + def __repr__(obj): + # {{{ Display + if not obj.list: + string=' list: empty' + else: + string=" list: (%i)"%(len(obj.list)) + for i in range(len(obj.list)): + if isinstance(obj.list[i][1],str): + string2=" field: %-10s value: '%s'"%(obj.list[i][0],obj.list[i][1]) + elif isinstance(obj.list[i][1],float): + string2=" field: %-10s value: %g"%(obj.list[i][0],obj.list[i][1]) + elif isinstance(obj.list[i][1],int): + string2=" field: %-10s value: %i"%(obj.list[i][0],obj.list[i][1]) + else: + string2=" field: %-10s value: (%i)"%(len(obj.list[i][1])) + string="%s\n%s"%(string,string2) + return string Index: /issm/trunk/src/m/classes/private.py =================================================================== --- /issm/trunk/src/m/classes/private.py (revision 12329) +++ /issm/trunk/src/m/classes/private.py (revision 12329) @@ -0,0 +1,29 @@ +#module imports +from fielddisplay import fielddisplay + +class private: + #properties + def __init__(self): + # {{{ Properties + self.runtimename = '' + self.bamg = {} + self.solution = ''; + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + string=' private parameters: do not change' + string="%s\n%s"%(string,fielddisplay(obj,'runtimename','name of the run launched')) + string="%s\n%s"%(string,fielddisplay(obj,'bamg','structure with mesh properties construced if bamg is used to mesh the domain')) + string="%s\n%s"%(string,fielddisplay(obj,'solution','type of solution launched')) + return string + #}}} + + def setdefaultparameters(obj): + # {{{setdefaultparameters + return obj + #}}} + Index: /issm/trunk/src/m/classes/prognostic.py =================================================================== --- /issm/trunk/src/m/classes/prognostic.py (revision 12329) +++ /issm/trunk/src/m/classes/prognostic.py (revision 12329) @@ -0,0 +1,47 @@ +#module imports +from fielddisplay import fielddisplay + +class prognostic: + #properties + def __init__(self): + # {{{ Properties + self.spcthickness = float('NaN') + self.min_thickness = 0 + self.hydrostatic_adjustment = 0 + self.stabilization = 0 + self.vertex_pairing = float('NaN') + self.penalty_factor = 0 + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + string=' Prognostic solution parameters:' + string="%s\n\n%s"%(string,fielddisplay(obj,'spcthickness','thickness constraints (NaN means no constraint)')) + string="%s\n%s"%(string,fielddisplay(obj,'min_thickness','minimum ice thickness allowed')) + string="%s\n%s"%(string,fielddisplay(obj,'hydrostatic_adjustment','adjustment of ice shelves surface and bed elevations: ''Incremental'' or ''Absolute'' ')) + string="%s\n%s"%(string,fielddisplay(obj,'stabilization','0->no, 1->artificial_diffusivity, 2->streamline upwinding, 3->discontinuous Galerkin')) + + return string + #}}} + + def setdefaultparameters(obj): + # {{{setdefaultparameters + + #Type of stabilization to use 0:nothing 1:artificial_diffusivity 3:Discontinuous Galerkin + obj.stabilization=1 + + #Factor applied to compute the penalties kappa=max(stiffness matrix)*10^penalty_factor + obj.penalty_factor=3 + + #Minimum ice thickness that can be used + obj.min_thickness=1 + + #Hydrostatic adjustment + obj.hydrostatic_adjustment='Absolute' + + return obj + #}}} + Index: /issm/trunk/src/m/classes/qmu.py =================================================================== --- /issm/trunk/src/m/classes/qmu.py (revision 12329) +++ /issm/trunk/src/m/classes/qmu.py (revision 12329) @@ -0,0 +1,39 @@ +#module imports +from fielddisplay import fielddisplay + +class qmu: + #properties + def __init__(self): + # {{{ Properties + self.isdakota = 0 + self.variables = {} + self.responses = {} + self.method = {} + self.params = {} + self.results = {} + self.partition = float('NaN') + self.numberofpartitions = 0 + self.numberofresponses = 0 + self.variabledescriptors = [] + self.responsedescriptors = [] + self.mass_flux_profile_directory = float('NaN') + self.mass_flux_profiles = float('NaN') + self.mass_flux_segments = [] + self.adjacency = float('NaN') + self.vertex_weight = float('NaN') + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + string=" qmu parameters: not implemented yet!" + return string + #}}} + + def setdefaultparameters(obj): + # {{{setdefaultparameters + return obj + #}}} + Index: /issm/trunk/src/m/classes/radaroverlay.py =================================================================== --- /issm/trunk/src/m/classes/radaroverlay.py (revision 12329) +++ /issm/trunk/src/m/classes/radaroverlay.py (revision 12329) @@ -0,0 +1,29 @@ +#module imports +from fielddisplay import fielddisplay + +class radaroverlay: + #properties + def __init__(self): + # {{{ Properties + self.pwr = float('NaN') + self.x = float('NaN') + self.y = float('NaN') + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + string=' radaroverlay parameters:' + string="%s\n\n%s"%(string,fielddisplay(obj,'pwr','radar power image (matrix)')) + string="%s\n%s"%(string,fielddisplay(obj,'x','corresponding x coordinates')) + string="%s\n%s"%(string,fielddisplay(obj,'y','corresponding y coordinates')) + return string + #}}} + + def setdefaultparameters(obj): + # {{{setdefaultparameters + return obj + #}}} + Index: /issm/trunk/src/m/classes/rifts.py =================================================================== --- /issm/trunk/src/m/classes/rifts.py (revision 12329) +++ /issm/trunk/src/m/classes/rifts.py (revision 12329) @@ -0,0 +1,28 @@ +#module imports +from fielddisplay import fielddisplay + +class rifts: + #properties + def __init__(self): + # {{{ Properties + self.riftstruct = float('NaN') + self.riftproperties = float('NaN') + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + string=' rifts parameters:' + + string="%s\n\n%s"%(string,fielddisplay(obj,'riftstruct','structure containing all rift information (vertices coordinates, segments, type of melange, ...)')) + string="%s\n%s"%(string,fielddisplay(obj,'riftproperties','')) + return string + #}}} + + def setdefaultparameters(obj): + # {{{setdefaultparameters + return obj + #}}} + Index: /issm/trunk/src/m/classes/settings.py =================================================================== --- /issm/trunk/src/m/classes/settings.py (revision 12329) +++ /issm/trunk/src/m/classes/settings.py (revision 12329) @@ -0,0 +1,53 @@ +#module imports +from fielddisplay import fielddisplay + +class settings: + #properties + def __init__(self): + # {{{ Properties + self.io_gather = 0 + self.lowmem = 0 + self.results_as_patches = 0 + self.output_frequency = 0 + self.waitonlock = 0 + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + string=" general settings parameters:" + + string="%s\n%s"%(string,fielddisplay(obj,"io_gather","I/O gathering strategy for result outputs (default 1)")) + string="%s\n%s"%(string,fielddisplay(obj,"lowmem","is the memory limited ? (0 or 1)")) + string="%s\n%s"%(string,fielddisplay(obj,"results_as_patches","provide results as patches for each element (0 or 1)")) + string="%s\n%s"%(string,fielddisplay(obj,"output_frequency","frequency at which results are saved in all solutions with multiple time_steps")) + string="%s\n%s"%(string,fielddisplay(obj,"waitonlock","maximum number of minutes to wait for batch results, or return 0")) + return string + #}}} + + def setdefaultparameters(obj): + # {{{setdefaultparameters + + #are we short in memory ? (0 faster but requires more memory) + obj.lowmem=0 + + #i/o: + obj.io_gather=1 + + #results frequency by default every step + obj.output_frequency=1 + + #do not use patches by default (difficult to plot) + obj.results_as_patches=0 + + #this option can be activated to load automatically the results + #onto the model after a parallel run by waiting for the lock file + #N minutes that is generated once the solution has converged + #0 to desactivate + obj.waitonlock=float('Inf') + + return obj + #}}} + Index: /issm/trunk/src/m/classes/solver.m =================================================================== --- /issm/trunk/src/m/classes/solver.m (revision 12328) +++ /issm/trunk/src/m/classes/solver.m (revision 12329) @@ -18,5 +18,8 @@ end end % }}} - function obj = addoptions(obj,analysis,solveroptions) % {{{1 + function obj = addoptions(obj,analysis,varargin) % {{{1 + % Usage example: + % md.solver=addoptions(md.solver,DiagnosticHorizAnalysisEnum,stokesoptions()); + % md.solver=addoptions(md.solver,DiagnosticHorizAnalysisEnum); %Convert analysis from enum to string @@ -29,5 +32,5 @@ %Add solver options to analysis - obj.(analysis) = solveroptions; + if nargin==3, obj.(analysis) = varargin{1}; end end %}}} Index: /issm/trunk/src/m/classes/solver.py =================================================================== --- /issm/trunk/src/m/classes/solver.py (revision 12329) +++ /issm/trunk/src/m/classes/solver.py (revision 12329) @@ -0,0 +1,49 @@ +#module imports {{{ +import fielddisplay +import ismumps +from mumpsoptions import * +from iluasmoptions import * +#}}} +class solver: + #properties + def __init__(self): + # {{{ Properties + if ismumps: + self.options=[["NoneAnalysis",mumpsoptions()]] + else: + self.options=[["NoneAnalysis",iluasmoptions()]] + #}}} + def __repr__(obj): + # {{{ Display + + string2=" solver parameters:" + for i in range(len(obj.options)): + option=obj.options[i] + analysis=option[0] + ioptions=option[1] + + string="" + for i in range(len(ioptions)): + option=ioptions[i] + if not option: + #do nothing + pass + elif len(option)==1: + #this option has only one argument + string="%s%s%s"%(string," -",option[0]) + elif len(option)==2: + #option with value. value can be string or scalar + if isinstance(option[1],float): + string="%s%s%s%s%s"%(string," -",option[0]," ","%g"%(option[1])) + elif isinstance(option[1],str): + string="%s%s%s%s%s"%(string," -",option[0]," ",option[1]) + elif isinstance(option[1],int): + string="%s%s%s%s%s"%(string," -",option[0]," ","%i"%(option[1])) + else: + raise RuntimeError("%s%s%s"%("PetscString error: option #","%i"%(i)," is not well formatted")) + else: + raise RuntimeError("%s%s%s"%("PetscString error: option #","%i"%(i)," is not well formatted")) + + string2="%s\n%s"%(string2," %s -> '%s'"%(analysis,string)) + return string2 + #}}} Index: /issm/trunk/src/m/classes/steadystate.py =================================================================== --- /issm/trunk/src/m/classes/steadystate.py (revision 12329) +++ /issm/trunk/src/m/classes/steadystate.py (revision 12329) @@ -0,0 +1,36 @@ +#module imports +from fielddisplay import fielddisplay + +class steadystate: + #properties + def __init__(self): + # {{{ Properties + self.reltol = 0 + self.maxiter = 0 + self.requested_outputs = float('NaN') + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + string=' steadystate solution parameters:' + string="%s\n%s"%(string,fielddisplay(obj,'reltol','relative tolerance criterion')) + string="%s\n%s"%(string,fielddisplay(obj,'maxiter','maximum number of iterations')) + string="%s\n%s"%(string,fielddisplay(obj,'requested_outputs','additional requested outputs')) + return string + #}}} + + def setdefaultparameters(obj): + # {{{setdefaultparameters + + #maximum of steady state iterations + obj.maxiter=100 + + #Relative tolerance for the steadystate convertgence + obj.reltol=0.01 + + return obj + #}}} + Index: /issm/trunk/src/m/classes/surfaceforcings.py =================================================================== --- /issm/trunk/src/m/classes/surfaceforcings.py (revision 12329) +++ /issm/trunk/src/m/classes/surfaceforcings.py (revision 12329) @@ -0,0 +1,29 @@ +#module imports +from fielddisplay import fielddisplay + +class surfaceforcings: + #properties + def __init__(self): + # {{{ Properties + self.precipitation = float('NaN') + self.mass_balance = float('NaN') + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + string=" surface forcings parameters:" + + string="%s\n\n%s"%(string,fielddisplay(obj,'precipitation','surface precipitation [m/yr water eq]')) + string="%s\n%s"%(string,fielddisplay(obj,'mass_balance','surface mass balance [m/yr ice eq]')) + + return string + #}}} + + def setdefaultparameters(obj): + # {{{setdefaultparameters + return obj + #}}} + Index: /issm/trunk/src/m/classes/thermal.py =================================================================== --- /issm/trunk/src/m/classes/thermal.py (revision 12329) +++ /issm/trunk/src/m/classes/thermal.py (revision 12329) @@ -0,0 +1,52 @@ +#module imports +from fielddisplay import fielddisplay + +class thermal: + #properties + def __init__(self): + # {{{ Properties + self.spctemperature = float('NaN') + self.penalty_threshold = 0 + self.stabilization = 0 + self.maxiter = 0 + self.penalty_lock = 0 + self.penalty_factor = 0 + self.isenthalpy = 0 + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + string=' Thermal solution parameters:' + string="%s\n\n%s"%(string,fielddisplay(obj,'spctemperature','temperature constraints (NaN means no constraint)')) + string="%s\n%s"%(string,fielddisplay(obj,'stabilization','0->no, 1->artificial_diffusivity, 2->SUPG')) + string="%s\n%s"%(string,fielddisplay(obj,'maxiter','maximum number of non linear iterations')) + string="%s\n%s"%(string,fielddisplay(obj,'penalty_lock','stabilize unstable thermal constraints that keep zigzagging after n iteration (default is 0, no stabilization)')) + string="%s\n%s"%(string,fielddisplay(obj,'penalty_threshold','threshold to declare convergence of thermal solution (default is 0)')) + string="%s\n%s"%(string,fielddisplay(obj,'isenthalpy','use an enthalpy formulation to include temperate ice (default is 0)')) + return string + #}}} + + def setdefaultparameters(obj): + # {{{setdefaultparameters + + #Number of unstable constraints acceptable + obj.penalty_threshold=0 + + #Type of stabilization used + obj.stabilization=1 + + #Maximum number of iterations + obj.maxiter=100 + + #factor used to compute the values of the penalties: kappa=max(stiffness matrix)*10^penalty_factor + obj.penalty_factor=3 + + #Should we use cold ice (default) or enthalpy formulation + obj.isenthalpy=0 + + return obj + #}}} + Index: /issm/trunk/src/m/classes/timestepping.py =================================================================== --- /issm/trunk/src/m/classes/timestepping.py (revision 12329) +++ /issm/trunk/src/m/classes/timestepping.py (revision 12329) @@ -0,0 +1,41 @@ +#module imports +from fielddisplay import fielddisplay + +class timestepping: + #properties + def __init__(self): + # {{{ Properties + self.time_step = 0 + self.final_time = 0 + self.time_adapt = 0 + self.cfl_coefficient = 0 + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + string=" timestepping parameters:" + string="%s\n\n%s"%(string,fielddisplay(obj,"time_step","length of time steps [yrs]")) + string="%s\n%s"%(string,fielddisplay(obj,"final_time","final time to stop the simulation [yrs]")) + string="%s\n%s"%(string,fielddisplay(obj,"time_adapt","use cfl condition to define time step ? (0 or 1) ")) + string="%s\n%s"%(string,fielddisplay(obj,"cfl_coefficient","coefficient applied to cfl condition")) + return string + #}}} + + def setdefaultparameters(obj): + # {{{setdefaultparameters + + #time between 2 time steps + obj.time_step=1/2 + + #final time + obj.final_time=10*obj.time_step + + #time adaptation? + obj.time_adapt=0 + obj.cfl_coefficient=.5 + + return obj + #}}} Index: /issm/trunk/src/m/classes/transient.py =================================================================== --- /issm/trunk/src/m/classes/transient.py (revision 12329) +++ /issm/trunk/src/m/classes/transient.py (revision 12329) @@ -0,0 +1,40 @@ +#module imports +from fielddisplay import fielddisplay + +class transient: + #properties + def __init__(self): + # {{{ Properties + self.isprognostic = 0 + self.isdiagnostic = 0 + self.isthermal = 0 + self.isgroundingline = 0 + self.requested_outputs = float('NaN') + + #set defaults + self.setdefaultparameters() + + #}}} + def __repr__(obj): + # {{{ Display + string=' transient solution parameters:' + string="%s\n%s"%(string,fielddisplay(obj,'isprognostic','indicates if a prognostic solution is used in the transient')) + string="%s\n%s"%(string,fielddisplay(obj,'isthermal','indicates if a thermal solution is used in the transient')) + string="%s\n%s"%(string,fielddisplay(obj,'isdiagnostic','indicates if a diagnostic solution is used in the transient')) + string="%s\n%s"%(string,fielddisplay(obj,'isgroundingline','indicates if a groundingline migration is used in the transient')) + string="%s\n%s"%(string,fielddisplay(obj,'requested_outputs','list of additional outputs requested')) + return string + #}}} + + def setdefaultparameters(obj): + # {{{setdefaultparameters + + #full analysis: Diagnostic, Prognostic and Thermal but no groundingline migration for now + obj.isprognostic=1 + obj.isdiagnostic=1 + obj.isthermal=1 + obj.isgroundingline=0 + + return obj + #}}} + Index: /issm/trunk/src/m/classes/verbose.py =================================================================== --- /issm/trunk/src/m/classes/verbose.py (revision 12329) +++ /issm/trunk/src/m/classes/verbose.py (revision 12329) @@ -0,0 +1,27 @@ +#module imports +from fielddisplay import fielddisplay + +class verbose: + #properties + def __init__(self): + # {{{ Properties + self.mprocessor = False + self.module = False + self.solution = False + self.solver = False + self.convergence = False + self.control = False + self.qmu = False + #}}} + def __repr__(obj): + # {{{ Display + string="%s%s%s\n\n"%("class '",obj.__class__.__name__,"'=") + string="%s%s\n"%(string," %15s : %s"%("mprocessor",obj.mprocessor)) + string="%s%s\n"%(string," %15s : %s"%("module",obj.module)) + string="%s%s\n"%(string," %15s : %s"%("solution",obj.solution)) + string="%s%s\n"%(string," %15s : %s"%("solver",obj.solver)) + string="%s%s\n"%(string," %15s : %s"%("convergence",obj.convergence)) + string="%s%s\n"%(string," %15s : %s"%("control",obj.control)) + string="%s%s\n"%(string," %15s : %s"%("qmu",obj.qmu)) + return string + #}}} Index: /issm/trunk/src/m/enum/MaxIterationConvergenceFlagEnum.m =================================================================== --- /issm/trunk/src/m/enum/MaxIterationConvergenceFlagEnum.m (revision 12329) +++ /issm/trunk/src/m/enum/MaxIterationConvergenceFlagEnum.m (revision 12329) @@ -0,0 +1,11 @@ +function macro=MaxIterationConvergenceFlagEnum() +%MAXITERATIONCONVERGENCEFLAGENUM - Enum of MaxIterationConvergenceFlag +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=MaxIterationConvergenceFlagEnum() + +macro=StringToEnum('MaxIterationConvergenceFlag'); Index: /issm/trunk/src/m/enum/MaximumNumberOfEnums.m =================================================================== --- /issm/trunk/src/m/enum/MaximumNumberOfEnums.m (revision 12329) +++ /issm/trunk/src/m/enum/MaximumNumberOfEnums.m (revision 12329) @@ -0,0 +1,11 @@ +function macro=MaximumNumberOfEnums() +%MAXIMUMNUMBEROFENUMS - Enum of MaximumNumberOfEnums +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=MaximumNumberOfEnums() + +macro=440; Index: /issm/trunk/src/m/kml/README.txt =================================================================== --- /issm/trunk/src/m/kml/README.txt (revision 12328) +++ /issm/trunk/src/m/kml/README.txt (revision 12329) @@ -61,5 +61,5 @@ There are some other utilities that are used in the construction of topological tables for the kml writing. -nodeconnectivity.m - create a node connectivity table (nnodes x mxepg+1) +kmlnodeconnectivity.m - create a node connectivity table (nnodes x mxepg+1) edgeadjacency.m - create an edge adjacency array (elems x edges) edgeperimeter.m - create an edge perimeter (edgeper x 2) and element perimeter (edgeper x 1) list Index: /issm/trunk/src/m/kml/kml_part_edges.m =================================================================== --- /issm/trunk/src/m/kml/kml_part_edges.m (revision 12328) +++ /issm/trunk/src/m/kml/kml_part_edges.m (revision 12329) @@ -150,5 +150,5 @@ elemp=md.mesh.elements(irow,:); epartp=epart(irow,:); - nodeconp=nodeconnectivity(elemp,md.mesh.numberofvertices); + nodeconp=kmlnodeconnectivity(elemp,md.mesh.numberofvertices); [edgeadjp]=edgeadjacency(elemp,nodeconp); [edgeper,elemper,iloop]=edgeperimeter(elemp,nodeconp,edgeadjp); Index: /issm/trunk/src/m/kml/kml_partitions.m =================================================================== --- /issm/trunk/src/m/kml/kml_partitions.m (revision 12328) +++ /issm/trunk/src/m/kml/kml_partitions.m (revision 12329) @@ -151,5 +151,5 @@ elemp=md.mesh.elements(irow,:); epartp=epart(irow,:); - nodeconp=nodeconnectivity(elemp,md.mesh.numberofvertices); + nodeconp=kmlnodeconnectivity(elemp,md.mesh.numberofvertices); [edgeadjp]=edgeadjacency(elemp,nodeconp); [edgeper,elemper,iloop]=edgeperimeter(elemp,nodeconp,edgeadjp); Index: /issm/trunk/src/m/kml/kml_unsh_edges.m =================================================================== --- /issm/trunk/src/m/kml/kml_unsh_edges.m (revision 12328) +++ /issm/trunk/src/m/kml/kml_unsh_edges.m (revision 12329) @@ -74,5 +74,5 @@ if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ... md.qmu.numberofpartitions - [edgeadj]=edgeadjacency(md.mesh.elements,md.nodeconnectivity); + [edgeadj]=edgeadjacency(md.mesh.elements,md.kmlnodeconnectivity); [icol,irow]=find(edgeadj'==0); edgeuns=zeros(length(irow),2); Index: /issm/trunk/src/m/kml/kmlnodeconnectivity.m =================================================================== --- /issm/trunk/src/m/kml/kmlnodeconnectivity.m (revision 12329) +++ /issm/trunk/src/m/kml/kmlnodeconnectivity.m (revision 12329) @@ -0,0 +1,58 @@ +% +% create a node connectivity table for the elements in the model. +% +% [nodecon]=edgeadjacency(elem,nnodes,mxepg) +% +% where the required input is: +% elem (numeric, element connectivity array (elems x nodes)) +% +% and the required output is: +% nodecon (numeric, node connectivity array (nnodes x mxepg+1)) +% +% the optional input is: +% nnodes (numeric, number of nodes) +% mxepg (numeric, max elements per node) +% +function [nodecon]=kmlnodeconnectivity(elem,nnodes,mxepg) + +if ~nargin + help kmlnodeconnectivity + return +end + +if ~exist('nnodes','var') || isempty(nnodes) + nnodes=max(max(elem)); +end +if ~exist('mxepg','var') || isempty(mxepg) + mxepg=25; +end + +%% create the node connectivity array + +nodecon=zeros(nnodes,mxepg+1); + +% loop over the elements + +for i=1:size(elem,1) + +% loop over the nodes for each element + + for j=1:size(elem,2) + if elem(i,j) + nodecon(elem(i,j),nodecon(elem(i,j),end)+1)=i; + nodecon(elem(i,j),end)=nodecon(elem(i,j),end)+1; + end + end +end + +%% sort the node connectivity array + +% loop over the nodes + +for i=1:size(nodecon,1) + if (nodecon(i,end) > 1) + nodecon(i,1:nodecon(i,end))=sort(nodecon(i,1:nodecon(i,end))); + end +end + +end Index: sm/trunk/src/m/kml/nodeconnectivity.m =================================================================== --- /issm/trunk/src/m/kml/nodeconnectivity.m (revision 12328) +++ (revision ) @@ -1,58 +1,0 @@ -% -% create a node connectivity table for the elements in the model. -% -% [nodecon]=edgeadjacency(elem,nnodes,mxepg) -% -% where the required input is: -% elem (numeric, element connectivity array (elems x nodes)) -% -% and the required output is: -% nodecon (numeric, node connectivity array (nnodes x mxepg+1)) -% -% the optional input is: -% nnodes (numeric, number of nodes) -% mxepg (numeric, max elements per node) -% -function [nodecon]=nodeconnectivity(elem,nnodes,mxepg) - -if ~nargin - help nodeconnectivity - return -end - -if ~exist('nnodes','var') || isempty(nnodes) - nnodes=max(max(elem)); -end -if ~exist('mxepg','var') || isempty(mxepg) - mxepg=25; -end - -%% create the node connectivity array - -nodecon=zeros(nnodes,mxepg+1); - -% loop over the elements - -for i=1:size(elem,1) - -% loop over the nodes for each element - - for j=1:size(elem,2) - if elem(i,j) - nodecon(elem(i,j),nodecon(elem(i,j),end)+1)=i; - nodecon(elem(i,j),end)=nodecon(elem(i,j),end)+1; - end - end -end - -%% sort the node connectivity array - -% loop over the nodes - -for i=1:size(nodecon,1) - if (nodecon(i,end) > 1) - nodecon(i,1:nodecon(i,end))=sort(nodecon(i,1:nodecon(i,end))); - end -end - -end Index: /issm/trunk/src/m/model/AnalysisConfiguration.m =================================================================== --- /issm/trunk/src/m/model/AnalysisConfiguration.m (revision 12329) +++ /issm/trunk/src/m/model/AnalysisConfiguration.m (revision 12329) @@ -0,0 +1,58 @@ +function [analyses,numanalyses]=AnalysisConfiguration(solutiontype), +%ANALYSISCONFIGURATION - return type of analyses, number of analyses +% +% Usage: +% [analyses, numanalyses]=AnalysisConfiguration(solutiontype); + + + +switch solutiontype, + + case DiagnosticSolutionEnum, + numanalyses=5; + analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum]; + + case SteadystateSolutionEnum, + numanalyses=7; + analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum;ThermalAnalysisEnum;MeltingAnalysisEnum]; + + case ThermalSolutionEnum, + numanalyses=2; + analyses=[ThermalAnalysisEnum;MeltingAnalysisEnum]; + + case EnthalpySolutionEnum, + numanalyses=1; + analyses=[EnthalpyAnalysisEnum]; + + case PrognosticSolutionEnum, + numanalyses=1; + analyses=[PrognosticAnalysisEnum]; + + case BalancethicknessSolutionEnum, + numanalyses=1; + analyses=[BalancethicknessAnalysisEnum]; + + case SurfaceSlopeSolutionEnum, + numanalyses=1; + analyses=[SurfaceSlopeAnalysisEnum]; + + case BedSlopeSolutionEnum, + numanalyses=1; + analyses=[BedSlopeAnalysisEnum]; + + case TransientSolutionEnum, + numanalyses=9; + analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum;ThermalAnalysisEnum;MeltingAnalysisEnum;EnthalpyAnalysisEnum;PrognosticAnalysisEnum]; + + case FlaimSolutionEnum, + numanalyses=1; + analyses=[FlaimAnalysisEnum]; + + case HydrologySolutionEnum, + numanalyses=3; + analyses=[BedSlopeAnalysisEnum;SurfaceSlopeAnalysisEnum;HydrologyAnalysisEnum]; + + otherwise + error('%s%s%s',' solution type: ',EnumToString(solutiontype),' not supported yet!'); + +end Index: /issm/trunk/src/m/model/MatlabProcessPatch.m =================================================================== --- /issm/trunk/src/m/model/MatlabProcessPatch.m (revision 12329) +++ /issm/trunk/src/m/model/MatlabProcessPatch.m (revision 12329) @@ -0,0 +1,65 @@ +function structure=MatlabProcessPatch(structure); +%PROCESSPATCH - create a structure from a patch +% +% Usage: +% Result=ProcessPatch(Result); + +%return if there is no fiel Patch +if (~isfield(structure,'Patch')), + return; +end + +%loop over steps +for i=1:length(structure), + + %Get Patch for current step + Patch=structure(i).Patch; + numvertices=structure(i).PatchVertices; + + %check that Patch is not empty + if length(Patch)==0 continue; end + + %Get number of fields; + fields=unique(Patch(:,1)); + steps=unique(Patch(:,2)); + + %parse steps + for j=1:length(steps), + + posstep=find(Patch(:,2)==steps(j)); + + %Take all the lines of the Patch for this timestep + temporarypatch=Patch(posstep,:); + time=temporarypatch(1,3); + step=temporarypatch(1,2); + + %parse fields + for i=1:length(fields), + + %get name + fieldname=EnumToString(fields(i)); + + %get line positions + pos=find(temporarypatch(:,1)==fields(i)); + + %Fill Result structure + structure(step).steps=step; + structure(step).time=time; + structure(step).(fieldname).element=temporarypatch(pos,4); + structure(step).(fieldname).interpolation=temporarypatch(pos,5); + structure(step).(fieldname).index=temporarypatch(pos,6:5+numvertices); + if structure(step).(fieldname).interpolation==P1Enum, + structure(step).(fieldname).value=temporarypatch(pos,6+numvertices:end); + end + if structure(step).(fieldname).interpolation==P0Enum, + structure(step).(fieldname).value=temporarypatch(pos,6+numvertices); + end + + end + end +end + +%remove fields +structure=rmfield(structure,'Patch'); +structure=rmfield(structure,'PatchVertices'); +structure=rmfield(structure,'PatchNodes'); Index: /issm/trunk/src/m/model/addnote.py =================================================================== --- /issm/trunk/src/m/model/addnote.py (revision 12329) +++ /issm/trunk/src/m/model/addnote.py (revision 12329) @@ -0,0 +1,29 @@ +def addnote(md, string): + + # Local Variables: md, string, i, notes, miscellaneous, newnotes + # Function calls: ischar, nargout, cell, nargin, length, addnote, error + #ADDNOTE - add a note to the existing model notes field + # + # Usage: + # md=addnote(md,string); + # + # Example: + # md=addnote(md,'Pine Island, Geometry of 2007'); + + if not isinstance(string,basetring): + print 'addnote error message: second input argument should be a string' + return [] + + notes = md.miscellaneous.notes + + if isinstance(notes,basestring): + newnotes=[notes,string] + else: + newnotes=[]; + for i in range(len(notes)): + newnotes=newnotes+notes[i] + + newnotes=newnotes+nodes; + + md.miscellaneous.notes = newnotes + return md Index: /issm/trunk/src/m/model/collapse.m =================================================================== --- /issm/trunk/src/m/model/collapse.m (revision 12328) +++ /issm/trunk/src/m/model/collapse.m (revision 12329) @@ -30,6 +30,6 @@ if ~isnan(md.inversion.vel_obs), md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers); end; if ~isnan(md.inversion.cost_functions_coefficients), md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers); end; -if ~isnan(md.inversion.min_parameters), md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers); end; -if ~isnan(md.inversion.max_parameters), md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers); end; +if numel(md.inversion.min_parameters)>1, md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers); end; +if numel(md.inversion.max_parameters)>1, md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers); end; if ~isnan(md.surfaceforcings.mass_balance), md.surfaceforcings.mass_balance=project2d(md,md.surfaceforcings.mass_balance,md.mesh.numberoflayers); @@ -101,6 +101,6 @@ %lat long -md.mesh.lat=project2d(md,md.mesh.lat,1); -md.mesh.long=project2d(md,md.mesh.long,1); +if numel(md.mesh.lat) ==md.mesh.numberofvertices, md.mesh.lat=project2d(md,md.mesh.lat,1); end +if numel(md.mesh.long)==md.mesh.numberofvertices, md.mesh.long=project2d(md,md.mesh.long,1); end %Initialize with the 2d mesh Index: /issm/trunk/src/m/model/marshall.m =================================================================== --- /issm/trunk/src/m/model/marshall.m (revision 12328) +++ /issm/trunk/src/m/model/marshall.m (revision 12329) @@ -15,4 +15,7 @@ error(['marshall error message: could not open ' [md.miscellaneous.name '.bin'],' file for binary writing']); end + +%First, write MaximumNumberOfEnum to make sure that the Enums are synchronized +WriteData(fid,'enum',MaximumNumberOfEnums(),'data',true,'format','Boolean'); %Go through all model fields: check that it is a class and call checkconsistency Index: sm/trunk/src/m/model/mesh/meshnodensity.m =================================================================== --- /issm/trunk/src/m/model/mesh/meshnodensity.m (revision 12328) +++ (revision ) @@ -1,71 +1,0 @@ -function md=meshnodensity(md,domainname,varargin) -%MESH - create model mesh -% -% This routine creates a model mesh using TriMeshNoDensity and a domain outline -% where md is a @model object, domainname is the name of an Argus domain outline file, -% Riftname is an optional argument (Argus domain outline) describing rifts. -% The difference with mesh.m is that the resolution of the mesh follows that of the domain -% outline and the riftoutline -% -% Usage: -% md=meshnodensity(md,domainname) -% or md=meshnodensity(md,domainname,riftname); -% -% Examples: -% md=meshnodensity(md,'DomainOutline.exp'); -% md=meshnodensity(md,'DomainOutline.exp','Rifts.exp'); - -if (nargin==2), - riftname=''; -end -if (nargin==3), - riftname=varargin{1}; -end - -%Mesh using TriMeshNoDensity -if strcmp(riftname,''), - [md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers]=TriMeshNoDensity(domainname); -else - [elements,x,y,segments,segmentmarkers]=TriMeshNoDensity(domainname,riftname); - - %check that all the created nodes belong to at least one element - orphan=find(~ismember([1:length(x)],sort(unique(elements(:))))); - for i=1:length(orphan), - %get rid of the orphan node i - %update x and y - x=[x(1:orphan(i)-(i-1)-1); x(orphan(i)-(i-1)+1:end)]; - y=[y(1:orphan(i)-(i-1)-1); y(orphan(i)-(i-1)+1:end)]; - %update elements - pos=find(elements>orphan(i)-(i-1)); - elements(pos)=elements(pos)-1; - %update segments - pos1=find(segments(:,1)>orphan(i)-(i-1)); - pos2=find(segments(:,2)>orphan(i)-(i-1)); - segments(pos1,1)=segments(pos1,1)-1; - segments(pos2,2)=segments(pos2,2)-1; - end - - %plug into md - md.mesh.x=x; - md.mesh.y=y; - md.mesh.elements=elements; - md.mesh.segments=segments; - md.mesh.segmentmarkers=segmentmarkers; -end - -%Fill in rest of fields: -md.mesh.numberofelements=length(md.mesh.elements); -md.mesh.numberofvertices=length(md.mesh.x); -md.mesh.z=zeros(md.mesh.numberofvertices,1); -md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; -md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1); -md.mesh.vertexonsurface=ones(md.mesh.numberofvertices,1); -md.mesh.elementonbed=ones(md.mesh.numberofelements,1); -md.mesh.elementonsurface=ones(md.mesh.numberofelements,1); - -%Now, build the connectivity tables for this mesh. -md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices); -md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity); - -%type of model -md.mesh.dimension=2; Index: sm/trunk/src/m/model/mesh/meshrefine.m =================================================================== --- /issm/trunk/src/m/model/mesh/meshrefine.m (revision 12328) +++ (revision ) @@ -1,40 +1,0 @@ -function md=meshrefine(md,areas) -%MESHREFINE: refined the mesh from a model, according to an area metric. -% -% Usage: -% md=meshrefine(md,metric) - -%some checks on list of arguments -if ((nargin~=2) | (nargout~=1)), - meshrefineusage(); - error('meshrefine error message'); -end -if ( (isempty(areas)) | (length(areas)~=md.mesh.numberofelements) | (length(find(isnan(areas))))), - meshrefineusage(); - error('meshrefine error message'); -end - -%Refine using TriMeshRefine -[md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers]=TriMeshRefine(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers,areas,'yes'); - -%Fill in rest of fields: -md.mesh.numberofelements=length(md.mesh.elements); -md.mesh.numberofvertices=length(md.mesh.x); -md.mesh.z=zeros(md.mesh.numberofvertices,1); -md.mesh.vertexonboundary=zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1; -md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1); -md.mesh.vertexonsurface=ones(md.mesh.numberofvertices,1); -md.mesh.elementonbed=ones(md.mesh.numberofelements,1); -md.mesh.elementonsurface=ones(md.mesh.numberofelements,1); - -%Now, build the connectivity tables for this mesh. -md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices); -md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity); - -%type of model -md.mesh.dimension=2; -end - -function meshrefineusage(), -disp('usage: md=meshrefine(md,areas)'); -end Index: /issm/trunk/src/m/model/mesh/triangle.m =================================================================== --- /issm/trunk/src/m/model/mesh/triangle.m (revision 12328) +++ /issm/trunk/src/m/model/mesh/triangle.m (revision 12329) @@ -38,33 +38,30 @@ %Mesh using TriMesh -if strcmp(riftname,''), - [md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers]=TriMesh(domainname,area,true); -else - [elements,x,y,segments,segmentmarkers]=TriMeshRifts(domainname,riftname,area,'yes'); +[elements,x,y,segments,segmentmarkers]=TriMesh(domainname,riftname,area); - %check that all the created nodes belong to at least one element - orphan=find(~ismember([1:length(x)],sort(unique(elements(:))))); - for i=1:length(orphan), - %get rid of the orphan node i - %update x and y - x=[x(1:orphan(i)-(i-1)-1); x(orphan(i)-(i-1)+1:end)]; - y=[y(1:orphan(i)-(i-1)-1); y(orphan(i)-(i-1)+1:end)]; - %update elements - pos=find(elements>orphan(i)-(i-1)); - elements(pos)=elements(pos)-1; - %update segments - pos1=find(segments(:,1)>orphan(i)-(i-1)); - pos2=find(segments(:,2)>orphan(i)-(i-1)); - segments(pos1,1)=segments(pos1,1)-1; - segments(pos2,2)=segments(pos2,2)-1; - end +%check that all the created nodes belong to at least one element +orphan=find(~ismember([1:length(x)],sort(unique(elements(:))))); +for i=1:length(orphan), + disp('WARNING: removing orphans'); + %get rid of the orphan node i + %update x and y + x=[x(1:orphan(i)-(i-1)-1); x(orphan(i)-(i-1)+1:end)]; + y=[y(1:orphan(i)-(i-1)-1); y(orphan(i)-(i-1)+1:end)]; + %update elements + pos=find(elements>orphan(i)-(i-1)); + elements(pos)=elements(pos)-1; + %update segments + pos1=find(segments(:,1)>orphan(i)-(i-1)); + pos2=find(segments(:,2)>orphan(i)-(i-1)); + segments(pos1,1)=segments(pos1,1)-1; + segments(pos2,2)=segments(pos2,2)-1; +end - %plug into md - md.mesh.x=x; - md.mesh.y=y; - md.mesh.elements=elements; - md.mesh.segments=segments; - md.mesh.segmentmarkers=segmentmarkers; -end +%plug into md +md.mesh.x=x; +md.mesh.y=y; +md.mesh.elements=elements; +md.mesh.segments=segments; +md.mesh.segmentmarkers=segmentmarkers; %Fill in rest of fields: Index: /issm/trunk/src/m/model/mesh/triangle.py =================================================================== --- /issm/trunk/src/m/model/mesh/triangle.py (revision 12329) +++ /issm/trunk/src/m/model/mesh/triangle.py (revision 12329) @@ -0,0 +1,56 @@ +from numpy import * +import TriMesh as tm +import NodeConnectivity as nc +import ElementConnectivity as ec + +def triangle(md, domainname, resolution,riftname=''): + #TRIANGLE - create model mesh using the triangle package + # + # This routine creates a model mesh using TriMesh and a domain outline, to within a certain resolution + # where md is a @model object, domainname is the name of an Argus domain outline file, + # and resolution is a characteristic length for the mesh (same unit as the domain outline + # unit). Riftname is an optional argument (Argus domain outline) describing rifts. + # + # Usage: + # md=triangle(md,domainname,resolution) + # or md=triangle(md,domainname, resolution, riftname) + # + # Examples: + # md=triangle(md,'DomainOutline.exp',1000); + # md=triangle(md,'DomainOutline.exp',1000,'Rifts.exp'); + + + #Figure out a characteristic area. Resolution is a node oriented concept (ex a 1000m resolution node would + #be made of 1000*1000 area squares). + + #Check that mesh was not already run, and warn user: + if md.mesh.numberofelements != 0.: + choice = input('This model already has a mesh. Are you sure you want to go ahead? (y/n)') + if choice != 'y': + print 'no meshing done ... exiting' + return [] + + area = resolution**2. + + #Mesh using TriMesh + [md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers]=tm.TriMesh(domainname,riftname,area) + + + #Fill in rest of fields: + md.mesh.numberofelements = len(md.mesh.elements) + md.mesh.numberofvertices = len(md.mesh.x) + md.mesh.z = zeros(md.mesh.numberofvertices) + md.mesh.vertexonboundary = zeros(md.mesh.numberofvertices) + md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1] = 1. + md.mesh.vertexonbed = ones(md.mesh.numberofvertices) + md.mesh.vertexonsurface = ones(md.mesh.numberofvertices) + md.mesh.elementonbed = ones(md.mesh.numberofelements) + md.mesh.elementonsurface = ones(md.mesh.numberofelements) + + #Now, build the connectivity tables for this mesh. + [md.mesh.vertexconnectivity]= nc.NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices) + [md.mesh.elementconnectivity] = ec.ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity) + + #type of model + md.mesh.dimension = 2. + return md Index: /issm/trunk/src/m/model/modis.m =================================================================== --- /issm/trunk/src/m/model/modis.m (revision 12328) +++ /issm/trunk/src/m/model/modis.m (revision 12329) @@ -14,9 +14,9 @@ %Get path to gdal binaries -path_gdal=[issmtier() '/externalpackages/gdal/install/bin/']; +path_gdal=[issmdir() '/externalpackages/gdal/install/bin/']; %Was gdal compiled? if ~exist([path_gdal 'gdal_translate']), - error(['modis error message: GDAL library needs to be compiled to use this routine. Compile GDAL in ' issmtier() '/externalpackages/gdal to use this routine.']); + error(['modis error message: GDAL library needs to be compiled to use this routine. Compile GDAL in ' issmdir() '/externalpackages/gdal to use this routine.']); end Index: /issm/trunk/src/m/model/parameterization/parameterize.py =================================================================== --- /issm/trunk/src/m/model/parameterization/parameterize.py (revision 12329) +++ /issm/trunk/src/m/model/parameterization/parameterize.py (revision 12329) @@ -0,0 +1,27 @@ +import os +def parameterize(md,parametername): + #PARAMETERIZE - parameterize a model + # + # from a parameter matlab file, start filling in all the model fields that were not + # filled in by the mesh.py and setmask.py model methods. + # Warning: the parameter file must be able to be run in Python + # + # Usage: + # md=parameterize(md,parametername) + # + # Example: + # md=parameterize(md,'Square.par'); + + #some checks + if not os.path.isfile(parametername): + print 'parameterize error message: file '+parametername+' not found.' + return [] + + #Try and run parameter file. + execfile(parametername) + + #ame and notes + if len(md.miscellaneous.name)==0: + md.miscellaneous.name=os.path.basename(parametername) + + md=addnote(md,'Model created by using parameter file: '+parametername+' on: '+str(datetime.datetime.now())) Index: /issm/trunk/src/m/model/partition/AreaAverageOntoPartition.m =================================================================== --- /issm/trunk/src/m/model/partition/AreaAverageOntoPartition.m (revision 12328) +++ /issm/trunk/src/m/model/partition/AreaAverageOntoPartition.m (revision 12329) @@ -29,5 +29,5 @@ %finally, project vector: vector=project2d(md3d,vector,layer); - md.qmu.partition=project2d(md3d,md3d.part,layer); + md.qmu.partition=project2d(md3d,md3d.qmu.partition,layer); end Index: /issm/trunk/src/m/model/petscversion.m =================================================================== --- /issm/trunk/src/m/model/petscversion.m (revision 12328) +++ /issm/trunk/src/m/model/petscversion.m (revision 12329) @@ -8,5 +8,5 @@ PETSC_VERSION=3; -configfile=[issmtier() '/config.h']; +configfile=[issmdir() '/config.h']; if ~exist(configfile,'file'), error(['File ' configfile ' not found. ISSM has not been configured yet!']); Index: /issm/trunk/src/m/model/plot/applyoptions.m =================================================================== --- /issm/trunk/src/m/model/plot/applyoptions.m (revision 12328) +++ /issm/trunk/src/m/model/plot/applyoptions.m (revision 12329) @@ -219,4 +219,7 @@ warning on MATLAB:log:logOfZero; set(c,'YTickLabel',labels); + end + if exist(options,'cbYLim'); + set(c,'YLim',getfieldvalue(options,'cbYLim')); end if exist(options,'colorbartitle'), @@ -289,5 +292,5 @@ textpositioni=textposition{i}; textrotationi=textrotation{i}; - h=text(textpositioni(1),textpositioni(2),textstringi,'FontSize',textsizei,'FontWeight',textweighti,'Color',textcolori,'Rotation',textrotationi); + h=text(textpositioni(1),textpositioni(2),10,textstringi,'FontSize',textsizei,'FontWeight',textweighti,'Color',textcolori,'Rotation',textrotationi); set(h,'Clipping','on'); %prevent text from appearing outside of the box end @@ -346,4 +349,5 @@ if exist(options,'axispos'), Axis=getfieldvalue(options,'axispos'); + hold on set(gca,'pos',Axis); end @@ -411,18 +415,20 @@ %box off if strcmpi(md.mesh.hemisphere,'n') | strcmpi(md.mesh.hemisphere,'north'), - A=expread([ issmdir() '/projects/Exp/GreenlandBoxFront.exp']); + A=expread([ jplsvn() '/projects/Exp/GreenlandBoxFront.exp']); [A.x A.y]=ll2xy(A.x,A.y,+1,45,70); + A.x = A.x(1:30:end); + A.y = A.y(1:30:end); elseif strcmpi(md.mesh.hemisphere,'s') | strcmpi(md.mesh.hemisphere,'south'), - A=expread([ issmdir() '/projects/Exp/Antarctica.exp']); + A=expread([ jplsvn() '/projects/Exp/Antarctica.exp']); else error('applyoptions error message: hemisphere not defined'); end - Ax=[min(A.x) max(A.x)]; - Ay=[min(A.y) max(A.y)]; + offset=3*10^4; + Ax=[min(A.x)-offset max(A.x)+offset]; + Ay=[min(A.y)-offset max(A.y)+offset]; %if we are zooming on a basin, don't take the mesh for the boundaries! if exist(options,'basin'), [mdx mdy]=basinzoom(options); else - offset=3*10^4; mdx=[min(md.mesh.x)-offset max(md.mesh.x)+offset]; mdy=[min(md.mesh.y)-offset max(md.mesh.y)+offset]; @@ -430,5 +436,5 @@ line(A.x,A.y,ones(size(A.x)),'color','b'); patch([Ax(1) Ax(2) Ax(2) Ax(1) Ax(1)],[Ay(1) Ay(1) Ay(2) Ay(2) Ay(1)],[1 1 1],'EdgeColor',[0 0 0],'LineWidth',1,'FaceLighting','none') - patch( [mdx(1) mdx(2) mdx(2) mdx(1)],[mdy(1) mdy(1) mdy(2) mdy(2)],ones(4,1),'EdgeColor',[0 0 0],'FaceColor','r','FaceAlpha',0.5) + patch([mdx(1) mdx(2) mdx(2) mdx(1)],[mdy(1) mdy(1) mdy(2) mdy(2)],ones(4,1),'EdgeColor',[0 0 0],'FaceColor','r','FaceAlpha',0.5) colorbar('off'); %back to main gca @@ -438,13 +444,20 @@ %flag edges of a partition if exist(options,'partitionedges') -[xsegments ysegments]=flagedges(md.mesh.elements,md.mesh.x,md.mesh.y,md.qmu.partition); -xsegments=xsegments*getfieldvalue(options,'unit',1); -ysegments=ysegments*getfieldvalue(options,'unit',1); -color=getfieldvalue(options,'partitionedgescolor','r-'); -linewidth=getfieldvalue(options,'linewidth',1); -hold on; -for i=1:length(xsegments), - plot(xsegments(i,:),ysegments(i,:),color,'LineWidth',linewidth); -end + [xsegments ysegments]=flagedges(md.mesh.elements,md.mesh.x,md.mesh.y,md.qmu.partition); + xsegments=xsegments*getfieldvalue(options,'unit',1); + ysegments=ysegments*getfieldvalue(options,'unit',1); + color=getfieldvalue(options,'partitionedgescolor','r-'); + linewidth=getfieldvalue(options,'linewidth',1); + hold on; + for i=1:length(xsegments), + plot(xsegments(i,:),ysegments(i,:),color,'LineWidth',linewidth); + end +end + +%Scatter +if exist(options,'scatter') + data=getfieldvalue(options,'scatter'); + hold on + plot_scatter(data(:,1),data(:,2),getfieldvalue(options,'scattersize',3),data(:,3),options); end Index: /issm/trunk/src/m/model/plot/northarrow.m =================================================================== --- /issm/trunk/src/m/model/plot/northarrow.m (revision 12328) +++ /issm/trunk/src/m/model/plot/northarrow.m (revision 12329) @@ -16,5 +16,5 @@ structure(6)=16; %default fontsize elseif length(structure)>6 - error('plotmodel error message: to many input arguments for northarrow: [x0 y0 length [ratio [width]]]'); + error('plotmodel error message: to many input arguments for northarrow: [x0 y0 length [ratio width fontsize]]'); end @@ -72,5 +72,5 @@ %Text North -xN=max([A(1) D(1) E(1) F(1) G(1)])+ratio/3*lengtharrow; +xN=max([A(1) D(1) E(1) F(1) G(1)])+ratio/3*abs(lengtharrow); yN=mean([A(2) F(2) H(2)]); text(xN,yN,'North','FontSize',fontsize,'FontWeight','b'); Index: /issm/trunk/src/m/model/plot/plot_gridded.m =================================================================== --- /issm/trunk/src/m/model/plot/plot_gridded.m (revision 12328) +++ /issm/trunk/src/m/model/plot/plot_gridded.m (revision 12329) @@ -40,4 +40,5 @@ data_max=max(data_grid(:)); end +options=changefieldvalue(options,'cbYLim',[data_min data_max]); if whitepos==1, white =data_max + (data_max-data_min)/55; Index: /issm/trunk/src/m/model/plot/plot_overlay.m =================================================================== --- /issm/trunk/src/m/model/plot/plot_overlay.m (revision 12328) +++ /issm/trunk/src/m/model/plot/plot_overlay.m (revision 12329) @@ -96,5 +96,5 @@ %h_data(find(data_grid>data_mean))=1; h_data=1*ones(size(data_grid)); - h_data(find(data_grid>data_mean))=0.7; + h_data(find(data_gridxlim(1) & md.mesh.yylim(1)); + flag=prod(flag_nodes(md.mesh.elements),2); + else: + #ok, flag elements + flag=ContourToMesh(md.mesh.elements[:,0:3],md.mesh.x,md.mesh.y,region,'element',1); + + if invert: + flag=~flag; + + elif isinstance(region,nparray): + if len(region)!=md.mesh.numberofelements: + print FlagElements.__doc__ + print 'Flaglist for region must be of same size as number of elements in model' + return [] + flag=region; + else: + print 'Invalid region option' + return [] + + return flag; Index: sm/trunk/src/m/utils/ImageProcessing/gradient_perso.m =================================================================== --- /issm/trunk/src/m/utils/ImageProcessing/gradient_perso.m (revision 12328) +++ (revision ) @@ -1,91 +1,0 @@ -function [ax,ay ]=gradient_perso(a,dx,dy,n); -%GRADIENT_PERSO - gradient computation -% -% INPUT a,dx,dy,n where a is the scalar value, dx and dy the spacing of one pixel -% in a and n is the width wanted in the gradient computation, in pixels. -% -% Usage: -% [ax,ay ]=gradient_perso(a,dx,dy,n) - -s=size(a); -ax=zeros(s(1),s(2)); -ay=zeros(s(1),s(2)); - -for k=n+1:s(1)-n, - if mod(k,10)==0, - disp(k/s(1)*100); - end - - for j=n+1:s(2)-n, - if isnan(a(k,j)), - ax(k,j)=NaN; - ay(k,j)=NaN; - else - temp=a(k,j); - temp2=a(k,j); - - count=1; - while ~isnan(a(k,j+count)), - temp=[temp a(k,j+count)]; - count=count+1; - if count>n, - count=count-1; - break; - end - - end - count=1; - while ~isnan(a(k,j-count)), - temp=[a(k,j-count) temp]; - count=count+1; - if count>n, - count=count-1; - break; - end - - end - - count=1; - while ~isnan(a(k+count,j)), - temp2=[temp2 a(k+count,j)]; - count=count+1; - if count>n, - count=count-1; - break; - end - - end - count=1; - while ~isnan(a(k-count,j)), - temp2=[a(k-count,j) temp2]; - count=count+1; - if count>n, - count=count-1; - break; - end - - end - - if length(temp)==1, - ax(k,j)=NaN; - else - ax(k,j)=(temp(length(temp))-temp(1))/(length(temp)-1)/dx; - end - - if length(temp2)==1, - ay(k,j)=NaN; - else - ay(k,j)=(temp2(length(temp2))-temp2(1))/(length(temp2)-1)/dy; - end - - end - end -end - - - - - - - - Index: sm/trunk/src/m/utils/ImageProcessing/transp.m =================================================================== --- /issm/trunk/src/m/utils/ImageProcessing/transp.m (revision 12328) +++ (revision ) @@ -1,10 +1,0 @@ -bg = uint8(255.*rand(100, 100, 3)); -imview(bg); - -fg(:, :, 1) = uint8(255 .*(cumsum(ones(100, 100)) ./ 100));% -fg(:, :, 2) = uint8(zeros(100)); -fg(:, :, 3) = flipud(fg(:, :, 1)); -imview(fg); - -res=immerge(bg, fg, .3); -imview(res); Index: /issm/trunk/src/m/utils/Meca/paterson.py =================================================================== --- /issm/trunk/src/m/utils/Meca/paterson.py (revision 12329) +++ /issm/trunk/src/m/utils/Meca/paterson.py (revision 12329) @@ -0,0 +1,50 @@ +from numpy import * + +def paterson(temperature): + + # Local Variables: pos11, pos5, pos10, temperature, pos, T, pos8, pos9, pos6, pos7, pos4, rigidity, pos2, pos3, pos1 + # Function calls: length, zeros, argwhere, paterson, error + #PATERSON - figure out the rigidity of ice for a given temperature + # + # rigidigty (in s^(1/3)Pa) is the flow law paramter in the flow law sigma=B*e(1/3) (Paterson, p97). + # temperature is in Kelvin degrees + # + # Usage: + # rigidity=paterson(temperature) + + pos=argwhere(temperature<0.) + if len(pos): + print 'input temperature should be in Kelvin (positive)' + return [] + + T = temperature-273.15 + #The routine below is equivalent to: + # n=3; T=temperature-273; + # %From paterson, + # Temp=[0;-2;-5;-10;-15;-20;-25;-30;-35;-40;-45;-50]; + # A=[6.8*10^-15;2.4*10^-15;1.6*10^-15;4.9*10^-16;2.9*10^-16;1.7*10^-16;9.4* + # 10^-17;5.1*10^-17;2.7*10^-17;1.4*10^-17;7.3*10^-18;3.6*10^-18];;%s-1(kPa-3) + # %Convert into rigidity B + # B=A.^(-1/n)*10^3; %s^(1/3)Pa + # %Now, do a cubic fit between Temp and B: + # fittedmodel=fit(Temp,B,'cubicspline'); + # rigidity=fittedmodel(temperature); + + rigidity=zeros(len(T)) + pos1=argwhere(T<=-45); rigidity[pos1]=10**8*(-0.000292866376675*(T[pos1]+50)**3+ 0.011672640664130*(T[pos1]+50)**2 -0.325004442485481*(T[pos1]+50)+ 6.524779401948101) + pos2=argwhere(logical_and(-45<=T,T<-40)); rigidity[pos2]=10**8*(-0.000292866376675*(T[pos2]+45)**3+ 0.007279645014004*(T[pos2]+45)**2 -0.230243014094813*(T[pos2]+45)+ 5.154964909039554) + pos3=argwhere(logical_and(-40<=T,T<-35)); rigidity[pos3]=10**8*(0.000072737147457*(T[pos3]+40)**3+ 0.002886649363879*(T[pos3]+40)**2 -0.179411542205399*(T[pos3]+40)+ 4.149132666831214) + pos4=argwhere(logical_and(-35<=T,T<-30)); rigidity[pos4]=10**8*(-0.000086144770023*(T[pos4]+35)**3+ 0.003977706575736*(T[pos4]+35)**2 -0.145089762507325*(T[pos4]+35)+ 3.333333333333331) + pos5=argwhere(logical_and(-30<=T,T<-25)); rigidity[pos5]=10**8*(-0.000043984685769*(T[pos5]+30)**3+ 0.002685535025386*(T[pos5]+30)**2 -0.111773554501713*(T[pos5]+30)+ 2.696559088937191) + pos6=argwhere(logical_and(-25<=T,T<-20)); rigidity[pos6]=10**8*(-0.000029799523463*(T[pos6]+25)**3+ 0.002025764738854*(T[pos6]+25)**2 -0.088217055680511*(T[pos6]+25)+ 2.199331606342181) + pos7=argwhere(logical_and(-20<=T,T<-15)); rigidity[pos7]=10**8*(0.000136920904777*(T[pos7]+20)**3+ 0.001578771886910*(T[pos7]+20)**2 -0.070194372551690*(T[pos7]+20)+ 1.805165505978111) + pos8=argwhere(logical_and(-15<=T,T<-10)); rigidity[pos8]=10**8*(-0.000899763781026*(T[pos8]+15)**3+ 0.003632585458564*(T[pos8]+15)**2 -0.044137585824322*(T[pos8]+15)+ 1.510778053489523) + pos9=argwhere(logical_and(-10<=T,T<-5)); rigidity[pos9]=10**8*(0.001676964325070*(T[pos9]+10)**3- 0.009863871256831*(T[pos9]+10)**2 -0.075294014815659*(T[pos9]+10)+ 1.268434288203714) + pos10=argwhere(logical_and(-5<=T,T<-2)); rigidity[pos10]=10**8*(-0.003748937622487*(T[pos10]+5)**3+0.015290593619213*(T[pos10]+5)**2 -0.048160403003748*(T[pos10]+5)+ 0.854987973338348) + pos11=argwhere(-2<=T); rigidity[pos11]=10**8*(-0.003748937622488*(T[pos11]+2)**3-0.018449844983174*(T[pos11]+2)**2 -0.057638157095631*(T[pos11]+2)+ 0.746900791092860) + + #Now make sure that rigidity is positive + pos=argwhere(rigidity<0); rigidity[pos]=1**6 + + return rigidity + Index: /issm/trunk/src/m/utils/Miscellaneous/issmdoc.m =================================================================== --- /issm/trunk/src/m/utils/Miscellaneous/issmdoc.m (revision 12328) +++ /issm/trunk/src/m/utils/Miscellaneous/issmdoc.m (revision 12329) @@ -2,9 +2,9 @@ %First get ISSM tier: -ISSM_TIER=issmtier; +ISSM_DIR=issmdir; disp(sprintf('\n%s',' A comprehensive documentation is available on http://issm.jpl.nasa.gov')); disp(sprintf('\n%s',' Example: how to create a square ice shelf')); -disp(sprintf('%s',' go to ',ISSM_TIER,'/examples/SquareIceshelf')); +disp(sprintf('%s',' go to ',ISSM_DIR,'/examples/SquareIceshelf')); disp(sprintf('%s',' md=model; %creates a new empty model structure')); disp(sprintf('%s',' md=triangle(md,''DomainOutline.exp'',50000); %creates a mesh of the domain outline with a resolution of 50000m')); Index: /issm/trunk/src/m/utils/OS/ismumps.m =================================================================== --- /issm/trunk/src/m/utils/OS/ismumps.m (revision 12328) +++ /issm/trunk/src/m/utils/OS/ismumps.m (revision 12329) @@ -6,5 +6,5 @@ -configfile=[issmtier() '/bin/config.h']; %should find it in the install target +configfile=[issmdir() '/bin/config.h']; %should find it in the install target if ~exist(configfile,'file'), error(['File ' configfile ' not found. ISSM has not been configured yet!']); Index: /issm/trunk/src/m/utils/OS/ispetsc.m =================================================================== --- /issm/trunk/src/m/utils/OS/ispetsc.m (revision 12328) +++ /issm/trunk/src/m/utils/OS/ispetsc.m (revision 12329) @@ -6,5 +6,5 @@ -configfile=[issmtier() '/bin/config.h']; %should find it in the install target +configfile=[issmdir() '/bin/config.h']; %should find it in the install target if ~exist(configfile,'file'), error(['File ' configfile ' not found. ISSM has not been configured yet!']); Index: /issm/trunk/src/m/utils/OS/issmscpin.m =================================================================== --- /issm/trunk/src/m/utils/OS/issmscpin.m (revision 12328) +++ /issm/trunk/src/m/utils/OS/issmscpin.m (revision 12329) @@ -32,10 +32,10 @@ %use the putty project pscp.exe: it should be in the path. - %get ISSM_TIER variable - [status,ISSM_TIER]=system('echo [%ISSM_TIER_WIN%]'); + %get ISSM_DIR variable + [status,ISSM_DIR]=system('echo [%ISSM_DIR_WIN%]'); if status, - error('scpin error message: could not find ISSM_TIER_WIN envirnoment variable'); + error('scpin error message: could not find ISSM_DIR_WIN envirnoment variable'); end - ISSM_TIER=ISSM_TIER(2:end-2); + ISSM_DIR=ISSM_DIR(2:end-2); username=input('Username: (quoted string) '); @@ -43,5 +43,5 @@ for i=1:numel(packages), - [status,result]=system([ISSM_TIER '/externalpackages/ssh/pscp.exe -l "' username '" -pw "' key '" ' host ':' path '/' packages{i} ' ./']); + [status,result]=system([ISSM_DIR '/externalpackages/ssh/pscp.exe -l "' username '" -pw "' key '" ' host ':' path '/' packages{i} ' ./']); if status, error('scpin error message: could not call putty pscp'); Index: /issm/trunk/src/m/utils/OS/issmscpout.m =================================================================== --- /issm/trunk/src/m/utils/OS/issmscpout.m (revision 12328) +++ /issm/trunk/src/m/utils/OS/issmscpout.m (revision 12329) @@ -23,10 +23,10 @@ %use the putty project pscp.exe: it should be in the path. - %get ISSM_TIER variable - [status,ISSM_TIER]=system('echo [%ISSM_TIER_WIN%]'); + %get ISSM_DIR variable + [status,ISSM_DIR]=system('echo [%ISSM_DIR_WIN%]'); if status, - error('scpout error message: could not find ISSM_TIER_WIN envirnoment variable'); + error('scpout error message: could not find ISSM_DIR_WIN envirnoment variable'); end - ISSM_TIER=ISSM_TIER(2:end-2); + ISSM_DIR=ISSM_DIR(2:end-2); username=input('Username: (quoted string) '); @@ -34,5 +34,5 @@ for i=1:numel(packages), - [status,result]=system([ISSM_TIER '/externalpackages/ssh/pscp.exe -l "' username '" -pw "' key '" ' packages{i} ' ' host ':' path]); + [status,result]=system([ISSM_DIR '/externalpackages/ssh/pscp.exe -l "' username '" -pw "' key '" ' packages{i} ' ' host ':' path]); if status, error('scpout error message: could not call putty pscp'); Index: /issm/trunk/src/m/utils/OS/issmssh.m =================================================================== --- /issm/trunk/src/m/utils/OS/issmssh.m (revision 12328) +++ /issm/trunk/src/m/utils/OS/issmssh.m (revision 12329) @@ -15,15 +15,15 @@ %use the putty project plink.exe: it should be in the path. - %get ISSM_TIER variable - [status,ISSM_TIER]=system('echo [%ISSM_TIER_WIN%]'); + %get ISSM_DIR variable + [status,ISSM_DIR]=system('echo [%ISSM_DIR_WIN%]'); if status, - error('issmssh error message: could not find ISSM_TIER_WIN envirnoment variable'); + error('issmssh error message: could not find ISSM_DIR_WIN envirnoment variable'); end - ISSM_TIER=ISSM_TIER(2:end-2); + ISSM_DIR=ISSM_DIR(2:end-2); username=input('Username: (quoted string) '); key=input('Key: (quoted string) '); - system([ISSM_TIER '/externalpackages/ssh/plink.exe -ssh -l "' username '" -pw "' key '" ' host ' "' command '"']); + system([ISSM_DIR '/externalpackages/ssh/plink.exe -ssh -l "' username '" -pw "' key '" ' host ' "' command '"']); else Index: /issm/trunk/src/m/utils/Shell/flaimdir.m =================================================================== --- /issm/trunk/src/m/utils/Shell/flaimdir.m (revision 12328) +++ /issm/trunk/src/m/utils/Shell/flaimdir.m (revision 12329) @@ -5,3 +5,3 @@ % FLAIM_DIR=flaimdir() -FLAIM_DIR=[issmtier '/externalpackages/flaim/install']; +FLAIM_DIR=[issmdir '/externalpackages/flaim/install']; Index: /issm/trunk/src/m/utils/Shell/issmdir.m =================================================================== --- /issm/trunk/src/m/utils/Shell/issmdir.m (revision 12328) +++ /issm/trunk/src/m/utils/Shell/issmdir.m (revision 12329) @@ -1,4 +1,4 @@ function ISSM_DIR=issmdir() -%ISSMDIR - Get ISSM_DIR environment variable contents. +%ISSMDIR - Get ISSM_DIR environment variable % % Usage: @@ -9,4 +9,7 @@ else ISSM_DIR =getenv('ISSM_DIR_WIN'); + if strcmpi(ISSM_DIR(end),'/') | strcmpi(ISSM_DIR(end),'\'), + ISSM_DIR = ISSM_DIR(1:end-1); %shave off the last '/' + end end Index: sm/trunk/src/m/utils/Shell/issmtier.m =================================================================== --- /issm/trunk/src/m/utils/Shell/issmtier.m (revision 12328) +++ (revision ) @@ -1,18 +1,0 @@ -function ISSM_TIER=issmtier() -%ISSMTIER - Get ISSM_TIER environment variable contents. -% -% Usage: -% ISSM_TIER=issmtier() - -if ~ispc, - ISSM_TIER =getenv('ISSM_TIER'); -else - ISSM_TIER =getenv('ISSM_TIER_WIN'); - if strcmpi(ISSM_TIER(end),'/') | strcmpi(ISSM_TIER(end),'\'), - ISSM_TIER = ISSM_TIER(1:end-1); %shave off the last '/' - end -end - -if (isempty(ISSM_TIER)), - error('issmdir error message: ''ISSM_TIER'' environment variable is empty! You should define ISSM_TIER in your .cshrc or .bashrc!'); -end Index: /issm/trunk/src/m/utils/Shell/jplsvn.m =================================================================== --- /issm/trunk/src/m/utils/Shell/jplsvn.m (revision 12329) +++ /issm/trunk/src/m/utils/Shell/jplsvn.m (revision 12329) @@ -0,0 +1,15 @@ +function JPL_SVN=jplsvn() +%ISSMDIR - Get JPL_SVN environment variable +% +% Usage: +% JPL_SVN=jplsvn() + +if ~ispc, + JPL_SVN =getenv('JPL_SVN'); +else + JPL_SVN =getenv('JPL_SVN_WIN'); +end + +if (isempty(JPL_SVN)), + error('jplsvn error message: ''JPL_SVN'' environment variable is empty! You should define JPL_SVN in your .cshrc or .bashrc'); +end Index: /issm/trunk/src/m/utils/Shell/ucisvn.m =================================================================== --- /issm/trunk/src/m/utils/Shell/ucisvn.m (revision 12329) +++ /issm/trunk/src/m/utils/Shell/ucisvn.m (revision 12329) @@ -0,0 +1,15 @@ +function UCI_SVN=ucisvn() +%ISSMDIR - Get UCI_SVN environment variable +% +% Usage: +% UCI_SVN=ucisvn() + +if ~ispc, + UCI_SVN =getenv('UCI_SVN'); +else + UCI_SVN =getenv('UCI_SVN_WIN'); +end + +if (isempty(UCI_SVN)), + error('ucisvn error message: ''UCI_SVN'' environment variable is empty! You should define UCI_SVN in your .cshrc or .bashrc'); +end