Index: /issm/trunk-jpl/src/m/classes/giamme.py =================================================================== --- /issm/trunk-jpl/src/m/classes/giamme.py (revision 25164) +++ /issm/trunk-jpl/src/m/classes/giamme.py (revision 25165) @@ -42,5 +42,5 @@ def checkconsistency(self, md, solution, analyses): # {{{ - if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and md.transient.istidallove == 0): + if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and md.transient.isslr == 0): return md Index: /issm/trunk-jpl/src/m/classes/model.py =================================================================== --- /issm/trunk-jpl/src/m/classes/model.py (revision 25164) +++ /issm/trunk-jpl/src/m/classes/model.py (revision 25165) @@ -688,6 +688,5 @@ md.materials.extrude(md) - if md.damage.isdamage == 1: - md.damage.extrude(md) + md.damage.extrude(md) md.mask.extrude(md) md.qmu.extrude(md) Index: /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py =================================================================== --- /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py (revision 25164) +++ /issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py (revision 25165) @@ -1,53 +1,15 @@ -def AnalysisConfiguration(solutiontype): #{{{ - """ - ANALYSISCONFIGURATION - return type of analyses, number of analyses - - Usage: - [analyses] = AnalysisConfiguration(solutiontype) - """ - - if solutiontype == 'StressbalanceSolution': - analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis'] - elif solutiontype == 'SteadystateSolution': - analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis'] - elif solutiontype == 'ThermalSolution': - analyses = ['EnthalpyAnalysis', 'ThermalAnalysis', 'MeltingAnalysis'] - elif solutiontype == 'MasstransportSolution': - analyses = ['MasstransportAnalysis'] - elif solutiontype == 'BalancethicknessSolution': - analyses = ['BalancethicknessAnalysis'] - elif solutiontype == 'SurfaceSlopeSolution': - analyses = ['L2ProjectionBaseAnalysis'] - elif solutiontype == 'BalancevelocitySolution': - analyses = ['BalancevelocityAnalysis'] - elif solutiontype == 'BedSlopeSolution': - analyses = ['L2ProjectionBaseAnalysis'] - elif solutiontype == 'GiaSolution': - analyses = ['GiaIvinsAnalysis'] - elif solutiontype == 'LoveSolution': - analyses = ['LoveAnalysis'] - elif solutiontype == 'TransientSolution': - analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis', 'MasstransportAnalysis', 'HydrologyShaktiAnalysis', 'HydrologyGladsAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis'] - elif solutiontype == 'HydrologySolution': - analyses = ['L2ProjectionBaseAnalysis', 'HydrologyShreveAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis'] - elif 'DamageEvolutionSolution': - analyses = ['DamageEvolutionAnalysis'] - - else: - raise TypeError("solution type: '%s' not supported yet!" % solutiontype) - return analyses - #}}} - - -def ismodelselfconsistent(md): - """ +def ismodelselfconsistent(md): #{{{ + ''' ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem. Usage: ismodelselfconsistent(md), - """ + ''' #initialize consistency as true md.private.isconsistent = True + + # print(md.mesh.z) + # exit() #Get solution and associated analyses @@ -64,5 +26,5 @@ #Check that current field is an object if not hasattr(getattr(md, field), 'checkconsistency'): - md.checkmessage("field '%s' is not an object." % field) + md.checkmessage('field {} is not an object.'.format(field)) #Check consistency of the object @@ -72,2 +34,51 @@ if not md.private.isconsistent: raise RuntimeError('Model not consistent, see messages above.') +#}}} + +def AnalysisConfiguration(solutiontype): #{{{ + ''' + ANALYSISCONFIGURATION - return type of analyses, number of analyses + + Usage: + [analyses] = AnalysisConfiguration(solutiontype) + ''' + + if solutiontype == 'StressbalanceSolution': + analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis'] + elif solutiontype == 'SteadystateSolution': + analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis'] + elif solutiontype == 'ThermalSolution': + analyses = ['EnthalpyAnalysis', 'ThermalAnalysis', 'MeltingAnalysis'] + elif solutiontype == 'MasstransportSolution': + analyses = ['MasstransportAnalysis'] + elif solutiontype == 'BalancethicknessSolution': + analyses = ['BalancethicknessAnalysis'] + elif solutiontype == 'Balancethickness2Solution': + analyses = ['Balancethickness2Analysis'] + elif solutiontype == 'BalancethicknessSoftSolution': + analyses = ['BalancethicknessAnalysis'] + elif solutiontype == 'BalancevelocitySolution': + analyses = ['BalancevelocityAnalysis'] + elif solutiontype == 'SurfaceSlopeSolution': + analyses = ['L2ProjectionBaseAnalysis'] + elif solutiontype == 'BedSlopeSolution': + analyses = ['L2ProjectionBaseAnalysis'] + elif solutiontype == 'GiaSolution': + analyses = ['GiaIvinsAnalysis'] + elif solutiontype == 'LoveSolution': + analyses = ['LoveAnalysis'] + elif solutiontype == 'EsaSolution': + analyses = ['EsaAnalysis'] + elif solutiontype == 'TransientSolution': + analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis', 'MasstransportAnalysis', 'HydrologyShaktiAnalysis', 'HydrologyGladsAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis', 'SealevelriseAnalysis'] + elif solutiontype == 'SealevelriseSolution': + analyses = ['SealevelriseAnalysis'] + elif solutiontype == 'HydrologySolution': + analyses = ['L2ProjectionBaseAnalysis', 'HydrologyShreveAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis'] + elif 'DamageEvolutionSolution': + analyses = ['DamageEvolutionAnalysis'] + else: + raise TypeError('solution type: {} not supported yet!'.format(solutiontype)) + + return analyses +#}}} Index: /issm/trunk-jpl/src/m/extrusion/project3d.m =================================================================== --- /issm/trunk-jpl/src/m/extrusion/project3d.m (revision 25164) +++ /issm/trunk-jpl/src/m/extrusion/project3d.m (revision 25165) @@ -38,8 +38,7 @@ elseif strcmpi(type,'node'), - %Initialize 3d vector if size(vector2d,1)==md.mesh.numberofvertices2d - projected_vector=paddingvalue*ones(md.mesh.numberofvertices, size(vector2d,2)); + projected_vector=paddingvalue*ones(md.mesh.numberofvertices, size(vector2d,2)); elseif size(vector2d,1)==md.mesh.numberofvertices2d+1 projected_vector=paddingvalue*ones(md.mesh.numberofvertices+1,size(vector2d,2)); Index: /issm/trunk-jpl/src/m/extrusion/project3d.py =================================================================== --- /issm/trunk-jpl/src/m/extrusion/project3d.py (revision 25164) +++ /issm/trunk-jpl/src/m/extrusion/project3d.py (revision 25165) @@ -1,26 +1,32 @@ import numpy as np + from pairoptions import pairoptions def project3d(md, *args): - """ + ''' PROJECT3D - vertically project a vector from 2d mesh - vertically project a vector from 2d mesh (split in noncoll and coll areas) into a 3d mesh. - This vector can be a node vector of size (md.mesh.numberofvertices2d, N / A) or an - element vector of size (md.mesh.numberofelements2d, N / A). - arguments: - 'vector': 2d vector - 'type': 'element' or 'node'. - options: - 'layer' a layer number where vector should keep its values. If not specified, all layers adopt the - value of the 2d vector. - 'padding': default to 0 (value adopted by other 3d layers not being projected + vertically project a vector from 2d mesh (split in noncoll and coll + areas) into a 3d mesh. + This vector can be a node vector of size (md.mesh.numberofvertices2d, + N/A) or an element vector of size (md.mesh.numberofelements2d, N/A). - Examples: - extruded_vector = project3d(md, 'vector', vector2d, 'type', 'node', 'layer', 1, 'padding', NaN) - extruded_vector = project3d(md, 'vector', vector2d, 'type', 'element', 'padding', 0) - extruded_vector = project3d(md, 'vector', vector2d, 'type', 'node') - """ + arguments: + 'vector': 2d vector + 'type': 'element' or 'node' + + options: + 'layer' a layer number where vector should keep its values. If + not specified, all layers adopt the value of the 2d + vector. + 'padding': default to 0 (value adopted by other 3d layers not + being projected. + + Examples: + extruded_vector = project3d(md, 'vector', vector2d, 'type', 'node', 'layer', 1, 'padding', NaN) + extruded_vector = project3d(md, 'vector', vector2d, 'type', 'element', 'padding', 0) + extruded_vector = project3d(md, 'vector', vector2d, 'type', 'node') + ''' #some regular checks @@ -37,10 +43,9 @@ paddingvalue = options.getfieldvalue('padding', 0) #0 by default - vector1d = False - if isinstance(vector2d, np.ndarray) and np.ndim(vector2d) == 1: - vector1d = True - vector2d = vector2d.reshape(-1, ) + #Handle special case where vector2d is single element (differs from representation in MATLAB) + if isinstance(vector2d, (bool, int, float)): + projected_vector = vector2d - if isinstance(vector2d, (bool, int, float)) or np.size(vector2d) == 1: + if np.size(vector2d) == 1: projected_vector = vector2d @@ -56,5 +61,5 @@ else: raise TypeError("vector length not supported") - #Fill in + #Fill in if layer == 0: for i in range(md.mesh.numberoflayers): @@ -71,5 +76,5 @@ else: raise TypeError("vector length not supported") - #Fill in + #Fill in if layer == 0: for i in range(md.mesh.numberoflayers): @@ -89,5 +94,5 @@ else: raise TypeError("vector length not supported") - #Fill in + #Fill in if layer == 0: for i in range(md.mesh.numberoflayers - 1): @@ -104,5 +109,5 @@ else: raise TypeError("vector length not supported") - #Fill in + #Fill in if layer == 0: for i in range(md.mesh.numberoflayers - 1): @@ -114,6 +119,3 @@ raise TypeError("project3d error message: unknown projection type") - if vector1d: - projected_vector = projected_vector.reshape(-1, ) - return projected_vector Index: /issm/trunk-jpl/src/m/partition/partitioner.m =================================================================== --- /issm/trunk-jpl/src/m/partition/partitioner.m (revision 25164) +++ /issm/trunk-jpl/src/m/partition/partitioner.m (revision 25165) @@ -14,4 +14,5 @@ % Usage: % partitionvector=partitioner(md,'package','chaco','npart',100,'weighting','on'); +% [partitionvector,md]=partitioner(md,'package','chaco','npart',100,'weighting','on'); % @@ -36,5 +37,5 @@ %partitioning essentially happens in 2D. So partition in 2D, then %extrude the partition vector vertically. - md3d=md; %save for later + md3d=md; %save for later md.mesh.elements=md.mesh.elements2d; md.mesh.x=md.mesh.x2d; @@ -78,5 +79,5 @@ % partition into nparts if isa(md.mesh,'mesh2d'), - part=Chaco(md.qmu.adjacency,weights,[],md.mesh.x, md.mesh.y,zeros(md.mesh.numberofvertices,1),method,npart,[])'+1; %index partitions from 1 up. like metis. + part=Chaco(md.qmu.adjacency,weights,[],md.mesh.x,md.mesh.y,zeros(md.mesh.numberofvertices,1),method,npart,[])'+1; %index partitions from 1 up. like metis. else part=Chaco(md.qmu.adjacency,weights,[],md.mesh.x, md.mesh.y,md.mesh.z,method,npart,[])'+1; %index partitions from 1 up. like metis. Index: /issm/trunk-jpl/src/m/partition/partitioner.py =================================================================== --- /issm/trunk-jpl/src/m/partition/partitioner.py (revision 25164) +++ /issm/trunk-jpl/src/m/partition/partitioner.py (revision 25165) @@ -1,7 +1,8 @@ +import copy + import numpy as np from adjacency import * from Chaco import * -import copy from mesh2d import * from MeshPartition import * @@ -10,41 +11,40 @@ -def partitioner(md, *varargin): +def partitioner(md, *args): ''' PARTITIONER - partition mesh - List of options to partitioner: + List of options to partitioner: + package: 'chaco', 'metis', or 'scotch' + npart: number of partitions + weighting: 'on' or 'off': default off + section: 1 by defaults(1=bisection, 2=quadrisection, 3=octasection) + recomputeadjacency: 'on' by default (set to 'off' to compute existing one) + type: 'node' or 'element' partition vector (default to 'node') + Output: partitionvector: the partition vector - package: 'chaco', 'metis', or 'scotch' - npart: number of partitions - weighting: 'on' or 'off': default off - section: 1 by defaults(1=bisection, 2=quadrisection, 3=octasection) - recomputeadjacency: 'on' by default (set to 'off' to compute existing one) - type: 'node' or 'element' partition vector (default to 'node') - Output: partitionvector: the partition vector - - Usage: - partitionvector = partitioner(md, 'package', 'chaco', 'npart', 100, 'weighting', 'on') + Usage: + partitionvector, md = partitioner(md, 'package', 'chaco', 'npart', 100, 'weighting', 'on') ''' #get options: - options = pairoptions(*varargin) + options = pairoptions(*args) - #get options: + #set defaults + package = options.getfieldvalue('package', 'chaco') + npart = options.getfieldvalue('npart', 10) + weighting = options.getfieldvalue('weighting', 'on') section = options.getfieldvalue('section', 1) - weighting = options.getfieldvalue('weighting', 'on') - package = options.getfieldvalue('package', 'chaco') #default to chaco - npart = options.getfieldvalue('npart', 10) # default to 10 - recomputeadjacency = options.getfieldvalue('recomputeadjacency', 'on') # default to on - vectortype = options.getfieldvalue('type', 'node') #default to node + recomputeadjacency = options.getfieldvalue('recomputeadjacency', 'on') + vectortype = options.getfieldvalue('type', 'node') # Python only: short-circuit if vectortype == 'element' and not package == 'linear': - raise RuntimeError('partitioner error message: package %s does not allow element partitions.' % package) + raise RuntimeError('partitioner error message: package {} does not allow element partitions.'.format(package)) - if(md.mesh.dimension() == 3): + if md.mesh.dimension() == 3: #partitioning essentially happens in 2D. So partition in 2D, then #extrude the partition vector vertically. - md3d = copy.deepcopy(md) + md3d = copy.deepcopy(md) # save for later md.mesh.elements = md.mesh.elements2d md.mesh.x = md.mesh.x2d @@ -64,8 +64,8 @@ if package == 'chaco': if vectortype == 'element': - raise RuntimeError('partitioner error message: package %s does not allow element partitions.' % package) + raise RuntimeError('partitioner error message: package {} does not allow element partitions.'.format(package)) else: - # default method (from chaco.m) - method = np.array([1, 1, 0, 0, 1, 1, 50, 0, 0.001, 7654321]) + # default method (from chaco.m) + method = np.array([1, 1, 0, 0, 1, 1, 50, 0, 0.001, 7654321]).conj().transpose() method[0] = 3 # global method (3 = inertial (geometric)) method[2] = 0 # vertex weights (0 = off, 1 = on) @@ -81,12 +81,11 @@ weights = [] - method = method.reshape(-1, 1) # transpose to 1x10 instead of 10 + #method = method.reshape(-1, 1) # transpose to 1x10 instead of 10 # partition into nparts if isinstance(md.mesh, mesh2d): - part = np.array(Chaco(md.qmu.adjacency, weights, np.array([]), md.mesh.x, md.mesh.y, np.zeros((md.mesh.numberofvertices, )), method, npart, np.array([]))).T + 1 #index partitions from 1 up. like metis. + part = Chaco(md.qmu.adjacency, weights, np.array([]), md.mesh.x, md.mesh.y, np.zeros((md.mesh.numberofvertices, 1)), method, npart, np.array([]))[0].conj().transpose() + 1 #index partitions from 1 up. like metis. else: - part = np.array(Chaco(md.qmu.adjacency, weights, np.array([]), md.mesh.x, md.mesh.y, md.mesh.z, method, npart, np.array([]))).T + 1 #index partitions from 1 up. like metis. - + part = Chaco(md.qmu.adjacency, weights, np.array([]), md.mesh.x, md.mesh.y, md.mesh.z, method, npart, np.array([]))[0].conj().transpose() + 1 #index partitions from 1 up. like metis. elif package == 'scotch': if vectortype == 'element': @@ -116,5 +115,5 @@ else: - raise RuntimeError('partitioner error message: could not find %s partitioner' % package) + raise RuntimeError('partitioner error message: could not find {} partitioner'.format(package)) #extrude if we are in 3D: @@ -122,11 +121,12 @@ md3d.qmu.vertex_weight = md.qmu.vertex_weight md3d.qmu.adjacency = md.qmu.adjacency - md = md3d + md = copy.deepcopy(md3d) if vectortype == 'element': - part = project3d(md, 'vector', np.squeeze(part), 'type', 'element') + part = project3d(md, 'vector', part.conj().transpose(), 'type', 'element') else: - part = project3d(md, 'vector', np.squeeze(part), 'type', 'node') + part = project3d(md, 'vector', part.conj().transpose(), 'type', 'node') - part = part.reshape(-1, 1) + if part.shape[0] == 1: + part = part.conj().transpose() # Output Index: /issm/trunk-jpl/test/NightlyRun/test445.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test445.m (revision 25164) +++ /issm/trunk-jpl/test/NightlyRun/test445.m (revision 25165) @@ -16,5 +16,7 @@ %dakota version -version=IssmConfig('_DAKOTA_VERSION_'); version=version(1:3); version=str2num(version); +version=IssmConfig('_DAKOTA_VERSION_'); +version=version(1:3); +version=str2num(version); %partitioning Index: /issm/trunk-jpl/test/NightlyRun/test445.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test445.py (revision 25164) +++ /issm/trunk-jpl/test/NightlyRun/test445.py (revision 25165) @@ -1,15 +1,18 @@ #Test Name: SquareSheetShelfSteaEnthalpyHO3dDakotaSampNeff +from os import getcwd +from socket import gethostname + import numpy as np -from os import getcwd + +from ContourToMesh import * +from dmeth_params_set import * from model import * -from socket import gethostname +from parameterize import * +from partitioner import * +from setflowequation import * +from setmask import * +from solve import * from triangle import * -from setmask import * -from parameterize import * -from setflowequation import * -from solve import * -from partitioner import * -from dmeth_params_set import * -from ContourToMesh import * + #model not consistent: equality thickness = surface-base violated @@ -24,6 +27,6 @@ md.thermal.isenthalpy = 1 md.thermal.isdynamicbasalspc = 1 -md.initialization.waterfraction = np.zeros((md.mesh.numberofvertices)) -md.initialization.watercolumn = np.zeros((md.mesh.numberofvertices)) +md.initialization.waterfraction = np.zeros((md.mesh.numberofvertices, 1)) +md.initialization.watercolumn = np.zeros((md.mesh.numberofvertices, 1)) md.friction.coupling = 3 @@ -36,5 +39,6 @@ #partitioning npart = 10 -partition = partitioner(md, 'package', 'chaco', 'npart', npart, 'weighting', 'on')[0] - 1 +partition, md = partitioner(md, 'package', 'chaco', 'npart', npart, 'weighting', 'on') +partition -= 1 md.qmu.isdakota = 1 @@ -69,5 +73,10 @@ #method md.qmu.method = dakota_method.dakota_method('nond_samp') -md.qmu.method = dmeth_params_set(md.qmu.method, 'seed', 1234, 'samples', 20, 'sample_type', 'random') +md.qmu.method = dmeth_params_set( + md.qmu.method, + 'seed', 1234, + 'samples', 20, + 'sample_type', 'random' + ) #parameters