Changeset 25165


Ignore:
Timestamp:
06/26/20 15:31:01 (5 years ago)
Author:
jdquinn
Message:

BUG: When partitioning 3d object, need to assign second return value back into md!

Location:
issm/trunk-jpl
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/m/classes/giamme.py

    r25158 r25165  
    4242
    4343    def checkconsistency(self, md, solution, analyses): # {{{
    44         if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and md.transient.istidallove == 0):
     44        if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and md.transient.isslr == 0):
    4545            return md
    4646

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

    r25164 r25165  
    688688
    689689        md.materials.extrude(md)
    690         if md.damage.isdamage == 1:
    691             md.damage.extrude(md)
     690        md.damage.extrude(md)
    692691        md.mask.extrude(md)
    693692        md.qmu.extrude(md)

  • issm/trunk-jpl/src/m/consistency/ismodelselfconsistent.py

    r24240 r25165  
    1 def AnalysisConfiguration(solutiontype):  #{{{
    2     """
    3     ANALYSISCONFIGURATION - return type of analyses, number of analyses
    4 
    5             Usage:
    6                     [analyses] = AnalysisConfiguration(solutiontype)
    7     """
    8 
    9     if solutiontype == 'StressbalanceSolution':
    10         analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis']
    11     elif solutiontype == 'SteadystateSolution':
    12         analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis']
    13     elif solutiontype == 'ThermalSolution':
    14         analyses = ['EnthalpyAnalysis', 'ThermalAnalysis', 'MeltingAnalysis']
    15     elif solutiontype == 'MasstransportSolution':
    16         analyses = ['MasstransportAnalysis']
    17     elif solutiontype == 'BalancethicknessSolution':
    18         analyses = ['BalancethicknessAnalysis']
    19     elif solutiontype == 'SurfaceSlopeSolution':
    20         analyses = ['L2ProjectionBaseAnalysis']
    21     elif solutiontype == 'BalancevelocitySolution':
    22         analyses = ['BalancevelocityAnalysis']
    23     elif solutiontype == 'BedSlopeSolution':
    24         analyses = ['L2ProjectionBaseAnalysis']
    25     elif solutiontype == 'GiaSolution':
    26         analyses = ['GiaIvinsAnalysis']
    27     elif solutiontype == 'LoveSolution':
    28         analyses = ['LoveAnalysis']
    29     elif solutiontype == 'TransientSolution':
    30         analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis', 'MasstransportAnalysis', 'HydrologyShaktiAnalysis', 'HydrologyGladsAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis']
    31     elif solutiontype == 'HydrologySolution':
    32         analyses = ['L2ProjectionBaseAnalysis', 'HydrologyShreveAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis']
    33     elif 'DamageEvolutionSolution':
    34         analyses = ['DamageEvolutionAnalysis']
    35 
    36     else:
    37         raise TypeError("solution type: '%s' not supported yet!" % solutiontype)
    38     return analyses
    39     #}}}
    40 
    41 
    42 def ismodelselfconsistent(md):
    43     """
     1def ismodelselfconsistent(md): #{{{
     2    '''
    443    ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem.
    454
    465       Usage:
    476          ismodelselfconsistent(md),
    48     """
     7    '''
    498
    509    #initialize consistency as true
    5110    md.private.isconsistent = True
     11
     12    # print(md.mesh.z)
     13    # exit()
    5214
    5315    #Get solution and associated analyses
    … …  
    6426        #Check that current field is an object
    6527        if not hasattr(getattr(md, field), 'checkconsistency'):
    66             md.checkmessage("field '%s' is not an object." % field)
     28            md.checkmessage('field {} is not an object.'.format(field))
    6729
    6830        #Check consistency of the object
    … …  
    7234    if not md.private.isconsistent:
    7335        raise RuntimeError('Model not consistent, see messages above.')
     36#}}}
     37
     38def AnalysisConfiguration(solutiontype): #{{{
     39    '''
     40    ANALYSISCONFIGURATION - return type of analyses, number of analyses
     41
     42        Usage:
     43                [analyses] = AnalysisConfiguration(solutiontype)
     44    '''
     45
     46    if solutiontype == 'StressbalanceSolution':
     47        analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis']
     48    elif solutiontype == 'SteadystateSolution':
     49        analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis']
     50    elif solutiontype == 'ThermalSolution':
     51        analyses = ['EnthalpyAnalysis', 'ThermalAnalysis', 'MeltingAnalysis']
     52    elif solutiontype == 'MasstransportSolution':
     53        analyses = ['MasstransportAnalysis']
     54    elif solutiontype == 'BalancethicknessSolution':
     55        analyses = ['BalancethicknessAnalysis']
     56    elif solutiontype == 'Balancethickness2Solution':
     57        analyses = ['Balancethickness2Analysis']
     58    elif solutiontype == 'BalancethicknessSoftSolution':
     59        analyses = ['BalancethicknessAnalysis']
     60    elif solutiontype == 'BalancevelocitySolution':
     61        analyses = ['BalancevelocityAnalysis']
     62    elif solutiontype == 'SurfaceSlopeSolution':
     63        analyses = ['L2ProjectionBaseAnalysis']
     64    elif solutiontype == 'BedSlopeSolution':
     65        analyses = ['L2ProjectionBaseAnalysis']
     66    elif solutiontype == 'GiaSolution':
     67        analyses = ['GiaIvinsAnalysis']
     68    elif solutiontype == 'LoveSolution':
     69        analyses = ['LoveAnalysis']
     70    elif solutiontype == 'EsaSolution':
     71        analyses = ['EsaAnalysis']
     72    elif solutiontype == 'TransientSolution':
     73        analyses = ['StressbalanceAnalysis', 'StressbalanceVerticalAnalysis', 'StressbalanceSIAAnalysis', 'L2ProjectionBaseAnalysis', 'ThermalAnalysis', 'MeltingAnalysis', 'EnthalpyAnalysis', 'MasstransportAnalysis', 'HydrologyShaktiAnalysis', 'HydrologyGladsAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis', 'SealevelriseAnalysis']
     74    elif solutiontype == 'SealevelriseSolution':
     75        analyses = ['SealevelriseAnalysis']
     76    elif solutiontype == 'HydrologySolution':
     77        analyses = ['L2ProjectionBaseAnalysis', 'HydrologyShreveAnalysis', 'HydrologyDCInefficientAnalysis', 'HydrologyDCEfficientAnalysis']
     78    elif 'DamageEvolutionSolution':
     79        analyses = ['DamageEvolutionAnalysis']
     80    else:
     81        raise TypeError('solution type: {} not supported yet!'.format(solutiontype))
     82
     83    return analyses
     84#}}}

  • issm/trunk-jpl/src/m/extrusion/project3d.m

    r17686 r25165  
    3838
    3939elseif strcmpi(type,'node'),
    40 
    4140        %Initialize 3d vector
    4241        if size(vector2d,1)==md.mesh.numberofvertices2d
    43                 projected_vector=paddingvalue*ones(md.mesh.numberofvertices,  size(vector2d,2));
     42                projected_vector=paddingvalue*ones(md.mesh.numberofvertices, size(vector2d,2));
    4443        elseif size(vector2d,1)==md.mesh.numberofvertices2d+1
    4544                projected_vector=paddingvalue*ones(md.mesh.numberofvertices+1,size(vector2d,2));

  • issm/trunk-jpl/src/m/extrusion/project3d.py

    r24256 r25165  
    11import numpy as np
     2
    23from pairoptions import pairoptions
    34
    45
    56def project3d(md, *args):
    6     """
     7    '''
    78    PROJECT3D - vertically project a vector from 2d mesh
    89
    9        vertically project a vector from 2d mesh (split in noncoll and coll areas) into a 3d mesh.
    10        This vector can be a node vector of size (md.mesh.numberofvertices2d, N / A) or an
    11        element vector of size (md.mesh.numberofelements2d, N / A).
    12        arguments:
    13           'vector': 2d vector
    14           'type': 'element' or 'node'.
    15        options:
    16           'layer' a layer number where vector should keep its values. If not specified, all layers adopt the
    17                  value of the 2d vector.
    18           'padding': default to 0 (value adopted by other 3d layers not being projected
     10        vertically project a vector from 2d mesh (split in noncoll and coll
     11        areas) into a 3d mesh.
     12        This vector can be a node vector of size (md.mesh.numberofvertices2d,
     13        N/A) or an element vector of size (md.mesh.numberofelements2d, N/A).
    1914
    20        Examples:
    21           extruded_vector = project3d(md, 'vector', vector2d, 'type', 'node', 'layer', 1, 'padding', NaN)
    22           extruded_vector = project3d(md, 'vector', vector2d, 'type', 'element', 'padding', 0)
    23           extruded_vector = project3d(md, 'vector', vector2d, 'type', 'node')
    24     """
     15        arguments:
     16            'vector': 2d vector
     17            'type': 'element' or 'node'
     18
     19        options:
     20            'layer'     a layer number where vector should keep its values. If
     21                        not specified, all layers adopt the value of the 2d
     22                        vector.
     23            'padding':  default to 0 (value adopted by other 3d layers not
     24                        being projected.
     25
     26        Examples:
     27            extruded_vector = project3d(md, 'vector', vector2d, 'type', 'node', 'layer', 1, 'padding', NaN)
     28            extruded_vector = project3d(md, 'vector', vector2d, 'type', 'element', 'padding', 0)
     29            extruded_vector = project3d(md, 'vector', vector2d, 'type', 'node')
     30    '''
    2531
    2632    #some regular checks
    … …  
    3743    paddingvalue = options.getfieldvalue('padding', 0)  #0 by default
    3844
    39     vector1d = False
    40     if isinstance(vector2d, np.ndarray) and np.ndim(vector2d) == 1:
    41         vector1d = True
    42         vector2d = vector2d.reshape(-1, )
     45    #Handle special case where vector2d is single element (differs from representation in MATLAB)
     46    if isinstance(vector2d, (bool, int, float)):
     47        projected_vector = vector2d
    4348
    44     if isinstance(vector2d, (bool, int, float)) or np.size(vector2d) == 1:
     49    if np.size(vector2d) == 1:
    4550        projected_vector = vector2d
    4651
    … …  
    5661            else:
    5762                raise TypeError("vector length not supported")
    58     #Fill in
     63            #Fill in
    5964            if layer == 0:
    6065                for i in range(md.mesh.numberoflayers):
    … …  
    7176            else:
    7277                raise TypeError("vector length not supported")
    73     #Fill in
     78            #Fill in
    7479            if layer == 0:
    7580                for i in range(md.mesh.numberoflayers):
    … …  
    8994            else:
    9095                raise TypeError("vector length not supported")
    91     #Fill in
     96            #Fill in
    9297            if layer == 0:
    9398                for i in range(md.mesh.numberoflayers - 1):
    … …  
    104109            else:
    105110                raise TypeError("vector length not supported")
    106     #Fill in
     111            #Fill in
    107112            if layer == 0:
    108113                for i in range(md.mesh.numberoflayers - 1):
    … …  
    114119        raise TypeError("project3d error message: unknown projection type")
    115120
    116     if vector1d:
    117         projected_vector = projected_vector.reshape(-1, )
    118 
    119121    return projected_vector

  • issm/trunk-jpl/src/m/partition/partitioner.m

    r25109 r25165  
    1414%   Usage:
    1515%      partitionvector=partitioner(md,'package','chaco','npart',100,'weighting','on');
     16%      [partitionvector,md]=partitioner(md,'package','chaco','npart',100,'weighting','on');
    1617%
    1718
    … …  
    3637        %partitioning essentially happens in 2D. So partition in 2D, then
    3738        %extrude the partition vector vertically.
    38         md3d=md; %save  for later
     39        md3d=md; %save for later
    3940        md.mesh.elements=md.mesh.elements2d;
    4041        md.mesh.x=md.mesh.x2d;
    … …  
    7879                %  partition into nparts
    7980                if isa(md.mesh,'mesh2d'),
    80                         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.
     81                        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.
    8182                else
    8283                        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.

  • issm/trunk-jpl/src/m/partition/partitioner.py

    r25132 r25165  
     1import copy
     2
    13import numpy as np
    24
    35from adjacency import *
    46from Chaco import *
    5 import copy
    67from mesh2d import *
    78from MeshPartition import *
    … …  
    1011
    1112
    12 def partitioner(md, *varargin):
     13def partitioner(md, *args):
    1314    '''
    1415    PARTITIONER - partition mesh
    1516
    16     List of options to partitioner:
     17        List of options to partitioner:
     18            package: 'chaco', 'metis', or 'scotch'
     19            npart: number of partitions
     20            weighting: 'on' or 'off': default off
     21            section: 1 by defaults(1=bisection, 2=quadrisection, 3=octasection)
     22            recomputeadjacency: 'on' by default (set to 'off' to compute existing one)
     23            type: 'node' or 'element' partition vector (default to 'node')
     24            Output: partitionvector: the partition vector
    1725
    18     package: 'chaco', 'metis', or 'scotch'
    19     npart: number of partitions
    20     weighting: 'on' or 'off': default off
    21     section: 1 by defaults(1=bisection, 2=quadrisection, 3=octasection)
    22     recomputeadjacency: 'on' by default (set to 'off' to compute existing one)
    23     type: 'node' or 'element' partition vector (default to 'node')
    24     Output: partitionvector: the partition vector
    25 
    26     Usage:
    27         partitionvector = partitioner(md, 'package', 'chaco', 'npart', 100, 'weighting', 'on')
     26        Usage:
     27            partitionvector, md = partitioner(md, 'package', 'chaco', 'npart', 100, 'weighting', 'on')
    2828    '''
    2929
    3030    #get options:
    31     options = pairoptions(*varargin)
     31    options = pairoptions(*args)
    3232
    33     #get options:
     33    #set defaults
     34    package = options.getfieldvalue('package', 'chaco')
     35    npart = options.getfieldvalue('npart', 10)
     36    weighting = options.getfieldvalue('weighting', 'on')
    3437    section = options.getfieldvalue('section', 1)
    35     weighting = options.getfieldvalue('weighting', 'on')
    36     package = options.getfieldvalue('package', 'chaco')  #default to chaco
    37     npart = options.getfieldvalue('npart', 10)  # default to 10
    38     recomputeadjacency = options.getfieldvalue('recomputeadjacency', 'on')  # default to on
    39     vectortype = options.getfieldvalue('type', 'node')  #default to node
     38    recomputeadjacency = options.getfieldvalue('recomputeadjacency', 'on')
     39    vectortype = options.getfieldvalue('type', 'node')
    4040
    4141    # Python only: short-circuit
    4242    if vectortype == 'element' and not package == 'linear':
    43         raise RuntimeError('partitioner error message: package %s does not allow element partitions.' % package)
     43        raise RuntimeError('partitioner error message: package {} does not allow element partitions.'.format(package))
    4444
    45     if(md.mesh.dimension() == 3):
     45    if md.mesh.dimension() == 3:
    4646        #partitioning essentially happens in 2D. So partition in 2D, then
    4747        #extrude the partition vector vertically.
    48         md3d = copy.deepcopy(md)
     48        md3d = copy.deepcopy(md) # save for later
    4949        md.mesh.elements = md.mesh.elements2d
    5050        md.mesh.x = md.mesh.x2d
    … …  
    6464    if package == 'chaco':
    6565        if vectortype == 'element':
    66             raise RuntimeError('partitioner error message: package %s does not allow element partitions.' % package)
     66            raise RuntimeError('partitioner error message: package {} does not allow element partitions.'.format(package))
    6767        else:
    68             #  default method (from chaco.m)
    69             method = np.array([1, 1, 0, 0, 1, 1, 50, 0, 0.001, 7654321])
     68            # default method (from chaco.m)
     69            method = np.array([1, 1, 0, 0, 1, 1, 50, 0, 0.001, 7654321]).conj().transpose()
    7070            method[0] = 3  #  global method (3 = inertial (geometric))
    7171            method[2] = 0  #  vertex weights (0 = off, 1 = on)
    … …  
    8181                weights = []
    8282
    83             method = method.reshape(-1, 1)  # transpose to 1x10 instead of 10
     83            #method = method.reshape(-1, 1)  # transpose to 1x10 instead of 10
    8484
    8585            #  partition into nparts
    8686            if isinstance(md.mesh, mesh2d):
    87                 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.
     87                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.
    8888            else:
    89                 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.
    90 
     89                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.
    9190    elif package == 'scotch':
    9291        if vectortype == 'element':
    … …  
    116115
    117116    else:
    118         raise RuntimeError('partitioner error message: could not find %s partitioner' % package)
     117        raise RuntimeError('partitioner error message: could not find {} partitioner'.format(package))
    119118
    120119    #extrude if we are in 3D:
    … …  
    122121        md3d.qmu.vertex_weight = md.qmu.vertex_weight
    123122        md3d.qmu.adjacency = md.qmu.adjacency
    124         md = md3d
     123        md = copy.deepcopy(md3d)
    125124        if vectortype == 'element':
    126             part = project3d(md, 'vector', np.squeeze(part), 'type', 'element')
     125            part = project3d(md, 'vector', part.conj().transpose(), 'type', 'element')
    127126        else:
    128             part = project3d(md, 'vector', np.squeeze(part), 'type', 'node')
     127            part = project3d(md, 'vector', part.conj().transpose(), 'type', 'node')
    129128
    130         part = part.reshape(-1, 1)
     129    if part.shape[0] == 1:
     130        part = part.conj().transpose()
    131131
    132132    # Output

  • issm/trunk-jpl/test/NightlyRun/test445.m

    r25005 r25165  
    1616
    1717%dakota version
    18 version=IssmConfig('_DAKOTA_VERSION_'); version=version(1:3); version=str2num(version);
     18version=IssmConfig('_DAKOTA_VERSION_');
     19version=version(1:3);
     20version=str2num(version);
    1921
    2022%partitioning

  • issm/trunk-jpl/test/NightlyRun/test445.py

    r25133 r25165  
    11#Test Name: SquareSheetShelfSteaEnthalpyHO3dDakotaSampNeff
     2from os import getcwd
     3from socket import gethostname
     4
    25import numpy as np
    3 from os import getcwd
     6
     7from ContourToMesh import *
     8from dmeth_params_set import *
    49from model import *
    5 from socket import gethostname
     10from parameterize import *
     11from partitioner import *
     12from setflowequation import *
     13from setmask import *
     14from solve import *
    615from triangle import *
    7 from setmask import *
    8 from parameterize import *
    9 from setflowequation import *
    10 from solve import *
    11 from partitioner import *
    12 from dmeth_params_set import *
    13 from ContourToMesh import *
     16
    1417
    1518#model not consistent:  equality thickness = surface-base violated
    … …  
    2427md.thermal.isenthalpy = 1
    2528md.thermal.isdynamicbasalspc = 1
    26 md.initialization.waterfraction = np.zeros((md.mesh.numberofvertices))
    27 md.initialization.watercolumn = np.zeros((md.mesh.numberofvertices))
     29md.initialization.waterfraction = np.zeros((md.mesh.numberofvertices, 1))
     30md.initialization.watercolumn = np.zeros((md.mesh.numberofvertices, 1))
    2831
    2932md.friction.coupling = 3
    … …  
    3639#partitioning
    3740npart = 10
    38 partition = partitioner(md, 'package', 'chaco', 'npart', npart, 'weighting', 'on')[0] - 1
     41partition, md = partitioner(md, 'package', 'chaco', 'npart', npart, 'weighting', 'on')
     42partition -= 1
    3943md.qmu.isdakota = 1
    4044
    … …  
    6973#method
    7074md.qmu.method = dakota_method.dakota_method('nond_samp')
    71 md.qmu.method = dmeth_params_set(md.qmu.method, 'seed', 1234, 'samples', 20, 'sample_type', 'random')
     75md.qmu.method = dmeth_params_set(
     76    md.qmu.method,
     77    'seed', 1234,
     78    'samples', 20,
     79    'sample_type', 'random'
     80    )
    7281
    7382#parameters
Note: See TracChangeset for help on using the changeset viewer.