Changeset 25161


Ignore:
Timestamp:
06/25/20 21:09:13 (5 years ago)
Author:
jdquinn
Message:

BUG: Marshalling of solid earth requested outputs.

Location:
issm/trunk-jpl
Files:
12 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/m/boundaryconditions/getlovenumbers.py

    r25158 r25161  
    1004710047    #cut off series at maxdeg
    1004810048    love_numbers = np.delete(love_numbers, range(maxdeg + 1, len(love_numbers)), axis=0)
    10049     # print(len(love_numbers))
    10050     # print(love_numbers[-1,:])
    10051     # exit()
    1005210049
    1005310050    #retrieve right type

  • issm/trunk-jpl/src/m/classes/friction.m

    r24666 r25161  
    11%FRICTION class definition
    22%
    3 %   Usage:
    4 %      friction=friction();
     3%       Usage:
     4%               friction=friction();
    55
    66classdef friction

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

    r24668 r25161  
     1from checkfield import checkfield
    12from fielddisplay import fielddisplay
    23from project3d import project3d
    3 from checkfield import checkfield
    44from WriteData import WriteData
    55
    66
    77class friction(object):
    8     """
     8    '''
    99    FRICTION class definition
    1010
    11        Usage:
    12           friction = friction()
    13     """
     11        Usage:
     12            friction = friction()
     13    '''
    1414
    1515    def __init__(self):  # {{{
    … …  
    2222        #set defaults
    2323        self.setdefaultparameters()
    24         self.requested_outputs = []
     24        #self.requested_outputs = []
    2525    #}}}
    2626
    … …  
    3434        string = "%s\n%s" % (string, fielddisplay(self, 'effective_pressure', 'Effective Pressure for the forcing if not coupled [Pa]'))
    3535        string = "%s\n%s" % (string, fielddisplay(self, 'effective_pressure_limit', 'Neff do not allow to fall below a certain limit: effective_pressure_limit * rho_ice * g * thickness (default 0)'))
    36         string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
     36        #string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
    3737        return string
    3838    #}}}
    … …  
    5151
    5252    def setdefaultparameters(self):  # {{{
    53         self.requested_outputs = ['default']
     53        #self.requested_outputs = ['default']
    5454        self.effective_pressure_limit = 0
    5555        return self
    … …  
    6262
    6363    def checkconsistency(self, md, solution, analyses):  # {{{
    64 
    6564        #Early return
    6665        if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses:
    … …  
    7675        elif self.coupling > 4:
    7776            raise ValueError('md.friction.coupling larger than 4, not supported yet')
    78         md = checkfield(md, 'fieldname', 'friction.requested_outputs', 'stringrow', 1)
     77        #md = checkfield(md, 'fieldname', 'friction.requested_outputs', 'stringrow', 1)
    7978        return md
    8079    # }}}
    … …  
    9493            raise ValueError('md.friction.coupling larger than 4, not supported yet')
    9594
    96     #process requested outputs
    97         outputs = self.requested_outputs
    98         indices = [i for i, x in enumerate(outputs) if x == 'default']
    99         if len(indices) > 0:
    100             outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:]
    101             outputs = outputscopy
    102         WriteData(fid, prefix, 'data', outputs, 'name', 'md.friction.requested_outputs', 'format', 'StringArray')
     95    # #process requested outputs
     96    #     outputs = self.requested_outputs
     97    #     indices = [i for i, x in enumerate(outputs) if x == 'default']
     98    #     if len(indices) > 0:
     99    #         outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:]
     100    #         outputs = outputscopy
     101    #     WriteData(fid, prefix, 'data', outputs, 'name', 'md.friction.requested_outputs', 'format', 'StringArray')
    103102    # }}}

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

    r25158 r25161  
    148148                'initialization',
    149149                'rifts',
     150                'dsl',
    150151                'solidearth',
    151                 'dsl',
    152152                'debug',
    153153                'verbose',

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

    r25160 r25161  
    8888        outputs = self.requested_outputs
    8989        pos = np.where(np.asarray(outputs) == 'default')[0]
    90         for i in pos:
    91             outputs[i] = ''  #remove 'default' from outputs
    92             outputs.extend(self.defaultoutputs(md))  #add defaults
    93 
     90        if len(pos):
     91            outputs = np.delete(outputs, pos) #remove 'default' from outputs
     92            outputs = np.append(outputs, self.defaultoutputs(md)) #add defaults
    9493        WriteData(fid, prefix, 'data', outputs, 'name', 'md.solidearth.requested_outputs', 'format', 'StringArray')
    9594    #}}}

  • issm/trunk-jpl/src/m/coordsystems/gmtmask.py

    r25011 r25161  
    1 from MatlabFuncs import *
    2 from model import *
    3 import numpy as np
    41from os import getenv, putenv
    52import subprocess
    63
     4import numpy as np
    75
    8 def gmtmask(lat, long, *varargin):
    9     '''GMTMASK - figure out which lat, long points are on the ocean
     6from MatlabFuncs import *
     7from model import *
    108
    11     Usage:
    12       mask.ocean = gmtmask(md.mesh.lat, md.mesh.long)
     9
     10def gmtmask(lat, long, *args):
     11    '''
     12    GMTMASK - figure out which lat, long points are on the ocean
     13
     14        Usage:
     15            mask.ocean = gmtmask(md.mesh.lat, md.mesh.long)
    1316    '''
    1417    lenlat = len(lat)
    … …  
    1619
    1720    #are we doing a recursive call?
    18     if len(varargin) == 3:
     21    if len(args) == 3:
    1922        recursive = 1
    2023    else:
    … …  
    2427        print(('             recursing: num vertices  #' + str(lenlat)))
    2528    else:
    26         print(('gmtmask: num vertices  #' + str(lenlat)))
     29        print(('gmtmask: num vertices ' + str(lenlat)))
    2730
    28     #Check lat and long size is not more than 50, 000 If so, recursively call gmtmask:
    29 
     31    #Check lat and long size is not more than 50,000. If so, recursively call gmtmask:
    3032    if lenlat > 50000:
    3133        for i in range(int(ceil(lenlat / 50000))):

  • issm/trunk-jpl/src/m/solve/marshall.m

    r21097 r25161  
    3434
    3535        %Marshall current object
    36         %disp(['marshalling ' field '...']);
     36        %disp(['marshalling ' field '...']); %Uncomment for debugging
    3737        marshall(md.(field),['md.' field],md,fid);
    3838end

  • issm/trunk-jpl/src/m/solve/marshall.py

    r25158 r25161  
    1313    '''
    1414    if md.verbose.solution:
    15         print(("marshalling file '%s.bin'." % md.miscellaneous.name))
     15        print("marshalling file {}.bin".format(md.miscellaneous.name))
    1616
    1717    #open file for binary writing

  • issm/trunk-jpl/src/m/solve/solve.m

    r24541 r25161  
    151151        return;
    152152end
     153
    153154%wait on lock
    154155if isnan(md.settings.waitonlock),

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

    r25147 r25161  
    1515late=sum(md.mesh.lat(md.mesh.elements),2)/3;
    1616longe=sum(md.mesh.long(md.mesh.elements),2)/3;
    17 pos=find(late <-80);
     17pos=find(late < -80);
    1818md.solidearth.surfaceload.icethicknesschange(pos)=-100;
    1919%greenland
    20 pos=find(late > 70 &  late < 80 & longe>-60 & longe<-30);
     20pos=find(late>70 & late<80 & longe>-60 & longe<-30);
    2121md.solidearth.surfaceload.icethicknesschange(pos)=-100;
    2222
    … …  
    2828mask=gmtmask(md.mesh.lat,md.mesh.long);
    2929icemask=ones(md.mesh.numberofvertices,1);
    30 pos=find(mask==0);  icemask(pos)=-1;
    31 pos=find(sum(mask(md.mesh.elements),2)<3);   icemask(md.mesh.elements(pos,:))=-1;
     30pos=find(mask==0);
     31icemask(pos)=-1;
     32pos=find(sum(mask(md.mesh.elements),2)<3);
     33icemask(md.mesh.elements(pos,:))=-1;
    3234md.mask.ice_levelset=icemask;
    3335md.mask.ocean_levelset=-icemask;

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

    r25147 r25161  
    33%mesh earth:
    44md=model;
    5 
    6 md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',1000.); %500 km resolution mesh
     5md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',1000.); %1000 km resolution mesh
    76
    87%parameterize slr solution:
    9 %slr loading:  {{{
     8%solidearth loading:  {{{
    109md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1);
    1110md.solidearth.sealevel=zeros(md.mesh.numberofvertices,1);
    … …  
    2726mask=gmtmask(md.mesh.lat,md.mesh.long);
    2827icemask=ones(md.mesh.numberofvertices,1);
    29 pos=find(mask==0);  icemask(pos)=-1;
    30 pos=find(sum(mask(md.mesh.elements),2)<3);   icemask(md.mesh.elements(pos,:))=-1;
     28pos=find(mask==0);
     29icemask(pos)=-1;
     30pos=find(sum(mask(md.mesh.elements),2)<3);
     31icemask(md.mesh.elements(pos,:))=-1;
    3132md.mask.ice_levelset=icemask;
    3233md.mask.ocean_levelset=-icemask;
    … …  
    4142% }}}
    4243
    43 % use model representation of ocea area (not the ture area)
     44% use model representation of ocen area (not the true area)
    4445md.solidearth.settings.ocean_area_scaling = 0;
    4546
    … …  
    6566
    6667%eustatic + rigid + elastic run:
    67 md.solidearth.settings.rigid=1; md.solidearth.settings.elastic=1; md.solidearth.settings.rotation=0;
     68md.solidearth.settings.rigid=1;
     69md.solidearth.settings.elastic=1;
     70md.solidearth.settings.rotation=0;
    6871md.cluster=generic('name',oshostname(),'np',3);
    6972%md.verbose=verbose('111111111');
    … …  
    7275
    7376%eustatic + rigid + elastic + rotation run:
    74 md.solidearth.settings.rigid=1; md.solidearth.settings.elastic=1; md.solidearth.settings.rotation=1;
     77md.solidearth.settings.rigid=1;
     78md.solidearth.settings.elastic=1; md.solidearth.settings.rotation=1;
    7579md.cluster=generic('name',oshostname(),'np',3);
    7680%md.verbose=verbose('111111111');

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

    r25158 r25161  
    1414#mesh earth:
    1515md = model()
    16 md.mesh = gmshplanet('radius', 6.371012 * 1e3, 'resolution', 1000.)  #500 km resolution mesh
     16md.mesh = gmshplanet('radius', 6.371012 * 1e3, 'resolution', 1000.)  #1000 km resolution mesh
    1717
    1818#parameterize solidearth solution:
    1919#solidearth loading:  {{{
    20 md.solidearth.deltathickness = np.zeros((md.mesh.numberofelements, ))
     20md.solidearth.surfaceload.icethicknesschange = np.zeros((md.mesh.numberofelements, ))
    2121md.solidearth.sealevel = np.zeros((md.mesh.numberofvertices, ))
    2222md.dsl.global_average_thermosteric_sea_level_change=np.zeros((2, ))
    … …  
    3030longe = sum(md.mesh.long[md.mesh.elements - 1], 1) / 3
    3131pos = np.intersect1d(np.array(np.where(late < -75)), np.array(np.where(longe < 0)))
    32 md.solidearth.deltathickness[pos] = -1
     32md.solidearth.surfaceload.icethicknesschange[pos] = -1
    3333
    3434#elastic loading from love numbers:
    … …  
    3838#mask:  {{{
    3939mask = gmtmask(md.mesh.lat, md.mesh.long)
    40 icemask = np.ones((md.mesh.numberofvertices, ))
    41 pos = np.where(mask == 0)
    42 #pos[0] because np.where(mask = 0) returns a 2d array, the latter parts of which are all array / s of 0s
    43 icemask[pos[0]] = -1
    44 pos = np.where(sum(mask[md.mesh.elements - 1], 1) < 3)
     40icemask = np.ones(md.mesh.numberofvertices)
     41pos = np.where(mask == 0)[0]
     42icemask[pos] = -1
     43pos = np.where(np.sum(mask[md.mesh.elements - 1], axis=1) < 3)
    4544icemask[md.mesh.elements[pos, :] - 1] = -1
    4645md.mask.ice_levelset = icemask
    4746md.mask.ocean_levelset = -icemask
    4847
    49 #make sure that the ice level set is all inclusive:
    50 md.mask.land_levelset = np.zeros((md.mesh.numberofvertices, ))
    51 md.mask.ocean_levelset = -np.ones((md.mesh.numberofvertices, ))
    52 
    53 #make sure that the elements that have loads are fully grounded:
    54 pos = np.nonzero(md.solidearth.deltathickness)[0]
     48#make sure that the elements that have loads are fully grounded
     49pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0]
    5550md.mask.ocean_levelset[md.mesh.elements[pos, :] - 1] = 1
    5651
    5752#make sure wherever there is an ice load, that the mask is set to ice:
    58 icemask[md.mesh.elements[pos, :] - 1] = -1
    59 md.mask.ice_levelset = icemask
     53#pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0] # Do we need to do this twice?
     54md.mask.ice_levelset[md.mesh.elements[pos, :] - 1] = 1
    6055# }}}
    6156
    62 # use model representation of ocea area (not the ture area)
    63 md.solidearth.ocean_area_scaling = 0
     57# use model representation of ocen area (not the true area)
     58md.solidearth.settings.ocean_area_scaling = 0
    6459
    6560#geometry
    6661di = md.materials.rho_ice / md.materials.rho_water
    67 md.geometry.thickness = np.ones((md.mesh.numberofvertices, ))
    68 md.geometry.surface = (1 - di) * np.zeros((md.mesh.numberofvertices, ))
     62md.geometry.thickness = np.ones(md.mesh.numberofvertices)
     63md.geometry.surface = (1 - di) * np.zeros(md.mesh.numberofvertices)
    6964md.geometry.base = md.geometry.surface - md.geometry.thickness
    7065md.geometry.bed = md.geometry.base
    … …  
    7873md.miscellaneous.name = 'test2003'
    7974
    80 #New stuff
    81 md.solidearth.spcthickness = np.nan * np.ones((md.mesh.numberofvertices, ))
    82 md.solidearth.Ngia = np.zeros((md.mesh.numberofvertices, ))
    83 md.solidearth.Ugia = np.zeros((md.mesh.numberofvertices, ))
    84 md.solidearth.hydro_rate = np.zeros((md.mesh.numberofvertices, ))
    85 
    8675#Solution parameters
    87 md.solidearth.reltol = float('NaN')
    88 md.solidearth.abstol = 1e-3
    89 md.solidearth.geodetic = 1
     76md.solidearth.settings.reltol = np.nan
     77md.solidearth.settings.abstol = 1e-3
     78md.solidearth.settings.computesealevelchange = 1
    9079
    9180#eustatic + rigid + elastic run:
    92 md.solidearth.rigid = 1
    93 md.solidearth.elastic = 1
    94 md.solidearth.rotation = 0
     81md.solidearth.settings.rigid = 1
     82md.solidearth.settings.elastic = 1
     83md.solidearth.settings.rotation = 0
    9584md.cluster = generic('name', gethostname(), 'np', 3)
    9685#md.verbose = verbose('111111111')
    97 #print md.calving
    98 #print md.gia
    99 #print md.love
    100 #print md.esa
    101 #print md.autodiff
    10286md = solve(md, 'Sealevelrise')
    10387SnoRotation = md.results.SealevelriseSolution.Sealevel
    10488
    10589#eustatic + rigid + elastic + rotation run:
    106 md.solidearth.rigid = 1
    107 md.solidearth.elastic = 1
    108 md.solidearth.rotation = 1
     90md.solidearth.settings.rigid = 1
     91md.solidearth.settings.elastic = 1
     92md.solidearth.settings.rotation = 1
    10993md.cluster = generic('name', gethostname(), 'np', 3)
    11094#md.verbose = verbose('111111111')
Note: See TracChangeset for help on using the changeset viewer.