Changeset 26840


Ignore:
Timestamp:
01/30/22 15:46:47 (3 years ago)
Author:
jdquinn
Message:

CHG: MATLAB -> Python for SE modules and tests; clean up

Location:
issm/trunk-jpl
Files:
20 edited

Legend:

Unmodified
Added
Removed

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

    r26800 r26840  
    66classdef fourierlove
    77        properties (SetAccess=public)
    8                 nfreq                       = 0;
    9                 frequencies                 = 0;
    10                 sh_nmax                     = 0;
    11                 sh_nmin                     = 0;
    12                 g0                          = 0;
    13                 r0                          = 0;
    14                 mu0                         = 0;
    15                 Gravitational_Constant      = 0;
    16                 chandler_wobble             = 0;
    17                 allow_layer_deletion        = 0;
    18                 underflow_tol               = 0;
    19                 pw_threshold                = 0;
    20                 integration_steps_per_layer = 0;
    21                 istemporal                  = 0;
    22                 n_temporal_iterations       = 0;
    23                 time                        = 0;
    24                 love_kernels                = 0;
    25                 forcing_type                = 0;
    26                 inner_core_boundary         = 0;
    27                 core_mantle_boundary        = 0;
    28                 complex_computation         = 0;
     8                nfreq                                           = 0;
     9                frequencies                                     = 0;
     10                sh_nmax                                         = 0;
     11                sh_nmin                                         = 0;
     12                g0                                                      = 0;
     13                r0                                                      = 0;
     14                mu0                                                     = 0;
     15                Gravitational_Constant          = 0;
     16                chandler_wobble                         = 0;
     17                allow_layer_deletion            = 0;
     18                underflow_tol                           = 0;
     19                pw_threshold                            = 0;
     20                integration_steps_per_layer     = 0;
     21                istemporal                                      = 0;
     22                n_temporal_iterations           = 0;
     23                time                                            = 0;
     24                love_kernels                            = 0;
     25                forcing_type                            = 0;
     26                inner_core_boundary                     = 0;
     27                core_mantle_boundary            = 0;
     28                complex_computation                     = 0;
    2929
    3030        end
    … …  
    180180                                        if self.time(i)==0
    181181                                                self.frequencies((i-1)*2*self.n_temporal_iterations +j) =0; % convention to avoid marshalling infinite numbers
    182                                         else                                           
     182                                        else
    183183                                                self.frequencies((i-1)*2*self.n_temporal_iterations +j) = j*log(2)/self.time(i)/2/pi;
    184184                                        end

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

    r26358 r26840  
    2222        self.mu0 = 0
    2323        self.Gravitational_Constant = 0
     24        self.chandler_wobble = 0
    2425        self.allow_layer_deletion = 0
    2526        self.underflow_tol = 0
     27        self.pw_threshold = 0
    2628        self.integration_steps_per_layer = 0
    2729        self.istemporal = 0
    … …  
    5153        s += '{}\n'.format(fielddisplay(self, 'allow_layer_deletion', 'allow for migration of the integration boundary with increasing spherical harmonics degree (default: 1)'))
    5254        s += '{}\n'.format(fielddisplay(self, 'Gravitational_Constant', 'Newtonian constant of gravitation (default: 6.67259e-11 [m^3 kg^-1 s^-2])'))
     55        s += '{}\n'.format(fielddisplay(self, 'chandler_wobble', 'includes the inertial terms for the chandler wobble in the rotational feedback love numbers, only for forcing_type=11 (default: 0) (/!\\ 1 has not been validated yet)'))
    5356        s += '{}\n'.format(fielddisplay(self, 'allow_layer_deletion', 'allow for migration of the integration boundary with increasing spherical harmonics degree (default: 1)'))
    54         s += '{}\n'.format(fielddisplay(self, 'underflow_tol', 'threshold of deep to surface love number ratio to trigger the deletion of layers (default: 1e-16)'))
     57        s += '{}\n'.format(fielddisplay(self, 'pw_threshold', 'if relative variation across frequencies is smaller than this ratio, the post-widder transform for time-dependent love numbers is bypassed (default (1e-3)'))
     58        s += '{}\n'.format(fielddisplay(self, 'chandler_wobble', 'includes the inertial terms for the chandler wobble in the rotational feedback love numbers, only for forcing_type=11 (default: 0) (/!\\ 1 is untested)'))
    5559        s += '{}\n'.format(fielddisplay(self, 'integration_steps_per_layer', 'number of radial steps to propagate the yi system from the bottom to the top of each layer (default: 100)'))
    5660        s += '{}\n'.format(fielddisplay(self, 'istemporal', ['1 for time-dependent love numbers, 0 for frequency-dependent or elastic love numbers (default: 0)', 'If 1: use fourierlove function build_frequencies_from_time to meet consistency']))
    … …  
    8892        self.mu0 = 1e11 # Pa
    8993        self.Gravitational_Constant = 6.67259e-11 # m^3 kg^-1 s^-2
     94        self.chandler_wobble = 0
    9095        self.allow_layer_deletion = 1
    9196        self.underflow_tol = 1e-16 # Threshold of deep to surface love number ratio to trigger the deletion of layer
     97        self.pw_threshold = 1e-3 # If relative variation across frequencies is smaller than this ratio, the post-widder transform for time-dependent love numbers is bypassed
    9298        self.integration_steps_per_layer = 100
    9399        self.istemporal = 0
    … …  
    113119        md = checkfield(md, 'fieldname', 'love.mu0', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
    114120        md = checkfield(md, 'fieldname', 'love.Gravitational_Constant', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     121        md = checkfield(md, 'fieldname', 'love.chandler_wobble', 'values', [0, 1])
    115122        md = checkfield(md, 'fieldname', 'love.allow_layer_deletion', 'values', [0, 1])
    116123        md = checkfield(md, 'fieldname', 'love.underflow_tol', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
     124        md = checkfield(md, 'fieldname', 'love.pw_threshold', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
    117125        md = checkfield(md, 'fieldname', 'love.integration_steps_per_layer', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
    118126        md = checkfield(md, 'fieldname', 'love.love_kernels', 'values', [0, 1])
    … …  
    127135            raise RuntimeError('Degree 1 not supported for forcing type {}. Use sh_min >= 2 for this kind of calculation.'.format(md.love.forcing_type))
    128136
     137        if md.love.chandler_wobble  == 1:
     138            print('Warning: Chandler wobble in Love number calculator has not been validated yet')
     139
    129140        # Need 'litho' material
    130         if not hasattr(md.materials, 'materials') or 'litho' not in md.materials.nature:
     141        if md.materials.__class__.__name__ != 'materials' or 'litho' not in md.materials.nature:
    131142            raise RuntimeError('Need a \'litho\' material to run a Fourier Love number analysis')
    132143
    133         mat = np.where(np.array(nature) == 'litho')[0]
     144        mat = np.where(np.array(md.materials.nature) == 'litho')[0]
    134145        if md.love.forcing_type <= 4:
    135146            md = checkfield(md, 'fieldname', 'love.inner_core_boundary', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0, '<=', md.materials[mat].numlayers)
    … …  
    149160        WriteData(fid, prefix, 'object', self, 'fieldname', 'mu0', 'format', 'Double')
    150161        WriteData(fid, prefix, 'object', self, 'fieldname', 'Gravitational_Constant', 'format', 'Double')
     162        WriteData(fid, prefix, 'object', self, 'fieldname', 'chandler_wobble', 'format', 'Boolean')
    151163        WriteData(fid, prefix, 'object', self, 'fieldname', 'allow_layer_deletion', 'format', 'Boolean')
    152164        WriteData(fid, prefix, 'object', self, 'fieldname', 'underflow_tol', 'format', 'Double')
     165        WriteData(fid, prefix, 'object', self, 'fieldname', 'pw_threshold', 'format', 'Double')
    153166        WriteData(fid, prefix, 'object', self, 'fieldname', 'integration_steps_per_layer', 'format', 'Integer')
    154167        WriteData(fid, prefix, 'object', self, 'fieldname', 'istemporal', 'format', 'Boolean')
    … …  
    171184        print('Temporal love numbers: Overriding md.love.nfreq and md.love.frequencies')
    172185        self.nfreq = len(self.time) * 2 * self.n_temporal_iterations
    173         self.frequencies = np.zeros((self.nfreq, 1))
     186        self.frequencies = np.zeros((self.nfreq,))
    174187        for i in range(len(self.time)):
    175188            for j in range(2 * self.n_temporal_iterations):
    176                 self.frequencies[(i - 1) * 2 * self.n_temporal_iterations + j] = j * np.log(2) / self.time[i] / 2 / np.pi
     189                if self.time[i] == 0:
     190                    self.frequencies[(i - 1) * 2 * self.n_temporal_iterations + j] = 0 # Convention to avoid marshalling infinite numbers
     191                else:
     192                    self.frequencies[(i - 1) * 2 * self.n_temporal_iterations + j] = j * np.log(2) / self.time[i] / 2 / np.pi
     193        return self
    177194    #}}}

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

    r26358 r26840  
    6161
    6262    def marshall(self, prefix, md, fid):  # {{{
    63         length_thickness = len(self.thickness)
     63        length_thickness = 1 if np.isnan(self.thickness) else len(self.thickness)
    6464        if (length_thickness == md.mesh.numberofvertices) or (length_thickness == md.mesh.numberofvertices + 1):
    6565            WriteData(fid, prefix, 'object', self, 'fieldname', 'thickness', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)

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

    r26755 r26840  
    5151
    5252    def marshall(self, prefix, md, fid):  # {{{
    53         WriteData(fid, prefix, 'name', 'md.hydrology.model', 'data', 2, 'format', 'Integer')
     53        WriteData(fid, prefix, 'name', 'md.hydrology.model', 'data', 6, 'format', 'Integer')
    5454        WriteData(fid, prefix, 'object', self, 'fieldname', 'spcwatercolumn', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
    5555        outputs = self.requested_outputs

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

    r26828 r26840  
    66from WriteData import WriteData
    77
    8 class lovenumbers(object):
     8
     9class lovenumbers(object):  #{{{
    910    """LOVENUMBERS class definition
    1011
    … …  
    1819
    1920    def __init__(self, *args):  #{{{
    20         # Regular love numbers
    21         self.h = []   # Provided by PREM model
    22         self.k = []   # idem
    23         self.l = []   # idem
     21        # Loading love numbers
     22        self.h = [] # Provided by PREM model
     23        self.k = [] # idem
     24        self.l = [] # idem
    2425
    2526        # Tidal love numbers for computing rotational feedback
    … …  
    2728        self.tk = []
    2829        self.tl = []
    29         self.tk2secular = 0  # deg 2 secular number
     30        self.tk2secular = 0 # deg 2 secular number
     31        self.pmtf_colinear = []
     32        self.pmtf_ortho = []
    3033        pmtf_colinear   = []
    3134        pmtf_ortho      = []
    … …  
    4043        self.setdefaultparameters(maxdeg, referenceframe)
    4144    #}}}
     45
    4246    def __repr__(self):  #{{{
    4347        s = '   lovenumbers parameters:\n'
    … …  
    5357        s += '{}\n'.format(fielddisplay(self, 'istime', 'time (default: 1) or frequency love numbers (0)'))
    5458        s += '{}\n'.format(fielddisplay(self, 'timefreq', 'time/frequency vector (yr or 1/yr)'))
     59        s += '{}\n'.format(fielddisplay(self, 'pmtf_colinear', 'Colinear component of the Polar Motion Transfer Function (e.g. x-motion due to x-component perturbation of the inertia tensor)'))
     60        s += '{}\n'.format(fielddisplay(self, 'pmtf_ortho', 'Orthogonal component of the Polar Motion Transfer Function (couples x and y components, only used for Chandler Wobble)'))
    5561        return s
    5662    #}}}
     63
    5764    def setdefaultparameters(self, maxdeg, referenceframe):  #{{{
    5865        # Initialize love numbers
    … …  
    7279            self.pmtf_ortho= 0.0
    7380
     81        self.pmtf_colinear = np.array([0.0]).reshape(-1, 1)
     82        self.pmtf_ortho = np.array([0.0]).reshape(-1, 1)
     83        if maxdeg >= 2:
     84            self.pmtf_colinear = ((1.0 + self.k[2, :]) / (1.0 - self.tk[2, :] / self.tk2secular)).reshape(-1, 1) # Valid only for elastic regime, not viscous. Also neglects chandler wobble.
     85            self.pmtf_ortho = np.array([0.0]).reshape(-1, 1)
    7486        # Time
    7587        self.istime = 1 # Temporal love numbers by default
    … …  
    7789        return self
    7890    #}}}
     91
    7992    def checkconsistency(self, md, solution, analyses):  #{{{
    8093        if ('SealevelchangeAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc):
    … …  
    98111
    99112        ntf = len(self.timefreq)
    100         if (np.shape(self.h)[1] != ntf or np.shape(self.k)[1] != ntf or np.shape(self.l)[1] != ntf or np.shape(self.th)[1] != ntf or np.shape(self.tk)[1] != ntf or np.shape(self.tl)[1] != ntf):
     113        if (np.shape(self.h)[1] != ntf or np.shape(self.k)[1] != ntf or np.shape(self.l)[1] != ntf or np.shape(self.th)[1] != ntf or np.shape(self.tk)[1] != ntf or np.shape(self.tl)[1] != ntf or np.shape(self.pmtf_colinear)[1] != ntf or np.shape(self.pmtf_ortho)[1] != ntf):
    101114            raise ValueError('lovenumbers error message: love numbers should have as many time/frequency steps as the time/frequency vector')
    102115
     116        if self.istime and self.timefreq[0] != 0:
     117            raise ValueError('temporal love numbers must start with elastic response, i.e. timefreq[0] = 0')
    103118        return md
    104119    #}}}
     120
    105121    def defaultoutputs(self, md):  #{{{
    106122        return[]
    107123    #}}}
     124
    108125    def marshall(self, prefix, md, fid):  #{{{
    109126        WriteData(fid, prefix, 'object', self, 'fieldname', 'h', 'name', 'md.solidearth.lovenumbers.h', 'format', 'DoubleMat', 'mattype', 1)
    … …  
    114131        WriteData(fid, prefix, 'object', self, 'fieldname', 'tk', 'name', 'md.solidearth.lovenumbers.tk', 'format', 'DoubleMat', 'mattype', 1)
    115132        WriteData(fid, prefix, 'object', self, 'fieldname', 'tl', 'name', 'md.solidearth.lovenumbers.tl', 'format', 'DoubleMat', 'mattype', 1)
     133        WriteData(fid, prefix, 'object', self, 'fieldname', 'pmtf_colinear', 'name', 'md.solidearth.lovenumbers.pmtf_colinear', 'format', 'DoubleMat', 'mattype', 1)
     134        WriteData(fid, prefix, 'object', self, 'fieldname', 'pmtf_ortho', 'name', 'md.solidearth.lovenumbers.pmtf_ortho', 'format', 'DoubleMat', 'mattype', 1)
    116135        WriteData(fid, prefix, 'object', self, 'data', self.tk2secular, 'fieldname', 'lovenumbers.tk2secular', 'format', 'Double')
    117136        WriteData(fid, prefix, 'object', self, 'fieldname', 'pmtf_colinear','name','md.solidearth.lovenumbers.pmtf_colinear','format','DoubleMat','mattype',1);
    … …  
    125144        WriteData(fid, prefix, 'object', self, 'fieldname', 'timefreq', 'name', 'md.solidearth.lovenumbers.timefreq', 'format', 'DoubleMat', 'mattype', 1, 'scale', scale);
    126145    #}}}
     146 
    127147    def extrude(self, md):  #{{{
    128148        return

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

    r26358 r26840  
    250250                    raise RuntimeError('First layer must be solid (issolid[0] > 0 AND lame_mu[0] > 0). Add a weak inner core if necessary.')
    251251                ind = np.where(md.materials.issolid == 0)[0]
    252                 if np.sum(np.in1d(np.diff(ind),1) >= 1): # If there are at least two consecutive indices that contain issolid = 0
    253                     raise RuntimeError('Fluid layers detected at layers #' + indices + ', but having 2 or more adjacent fluid layers is not supported yet. Consider merging them.')
     252                #if np.sum(np.in1d(np.diff(ind),1) >= 1): # If there are at least two consecutive indices that contain issolid = 0
     253                #    raise RuntimeError('Fluid layers detected at layers #' + indices + ', but having 2 or more adjacent fluid layers is not supported yet. Consider merging them.')
    254254
    255255            elif nat == 'hydro':

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

    r25817 r26840  
    143143
    144144    def __getattr__(self, key):  #{{{
    145         if len(self.steps) == 1:
    146             return getattr(self.steps[0], key)
    147         else:
    148             raise Exception('<results>.<solution> error: Currently, can only get a field if we are not working with a transient solution.')
     145        # NOTE: Currently only returning value from first frame of transient solution (see retrieval of md.results.TransientSolution.BedGRD in test2051.py for justification)
     146        return getattr(self.steps[0], key)
     147
     148        # Original code
     149        # if len(self.steps) == 1:
     150        #     return getattr(self.steps[0], key)
     151        # else:
     152        #     raise Exception('<results>.<solution> error: Currently, can only get a field if we are not working with a transient solution.')
    149153    #}}}
    150154

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

    r26358 r26840  
    120120            if md.mesh.__class__.__name__ == 'mesh3dsurface':
    121121                if self.grdmodel == 2:
    122                     raise RuntimeException('model requires a 2D mesh to run gia Ivins computations (change mesh from mesh3dsurface to mesh2d)')
     122                    raise Exception('model requires a 2D mesh to run gia Ivins computations (change mesh from mesh3dsurface to mesh2d)')
    123123            else:
    124124                if self.grdmodel == 1:
    125                     raise RuntimeException('model requires a 3D surface mesh to run GRD computations (change mesh from mesh2d to mesh3dsurface)')
     125                    raise Exception('model requires a 3D surface mesh to run GRD computations (change mesh from mesh2d to mesh3dsurface)')
    126126            if self.sealevelloading and not self.grdocean:
    127                 raise RuntimeException('solidearthsettings checkconsistency error message: need grdocean on if sealevelloading flag is set')
     127                raise Exception('solidearthsettings checkconsistency error message: need grdocean on if sealevelloading flag is set')
    128128
    129129        if self.compute_bp_grd and not md.solidearth.settings.isgrd:
    130             raise RuntimeException('solidearthsettings checkconsistency error message; if bottom pressure grd patterns are requested, solidearth settings class should have isgrd flag on')
     130            raise Exception('solidearthsettings checkconsistency error message; if bottom pressure grd patterns are requested, solidearth settings class should have isgrd flag on')
    131131
    132132        return md

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

    r26358 r26840  
    9797md.transient.ismasstransport = 1
    9898md.transient.isslc = 1
    99 md.solidearth.requested_outputs = ['Sealevel', 'Bed']
     99md.solidearth.requested_outputs = ['Sealevel', 'Bed', 'SealevelBarystaticIceLoad', 'SealevelBarystaticIceArea', 'SealevelBarystaticIceWeights']
    100100
    101101# Max number of iterations reverted back to 10 (i.e., the original default value)

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

    r26358 r26840  
    465465    'DeltaIceThickness',
    466466    'Sealevel',
    467     'SealevelUGrd',
    468     'SealevelchangeBarystaticMask',
    469     'SealevelchangeBarystaticOceanMask',
     467    'Bed',
     468    'SealevelBarystaticIceMask',
     469    'SealevelBarystaticOceanMask',
    470470]
    471471md = solve(md, 'Transient')

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

    r26809 r26840  
    113113moi_yz = sum(-loadice.*areaice.*rad_e^2.*sin(lat).*cos(lat).*sin(lon));
    114114moi_zz = sum(-loadice.*areaice.*rad_e^2.*(-1.0/3.0+sin(lat).^2));
    115 theoretical_value_check=[moixz/moi_xz moiyz/moi_yz moizz/moi_zz] % should yield [1.0 1.0 1.0]
     115theoretical_value_check=[moixz/moi_xz moiyz/moi_yz moizz/moi_zz]; % should yield [1.0 1.0 1.0]
    116116% }}}
    117117

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

    r26638 r26840  
    115115eus = md.results.TransientSolution.Bslc
    116116slc = md.results.TransientSolution.Sealevel
    117 moixz = md.results.TransientSolution.SealevelInertiaTensorXZ / (1 / (1 - tide_love_k2 / tide_love_k2secular) * (1 + load_love_k2) / (moi_p - moi_e))
    118 moiyz = md.results.TransientSolution.SealevelInertiaTensorYZ / (1 / (1 - tide_love_k2 / tide_love_k2secular) * (1 + load_love_k2) / (moi_p - moi_e))
    119 moizz = md.results.TransientSolution.SealevelInertiaTensorZZ / ( -(1 + load_love_k2) / moi_p)
     117moixz = md.results.TransientSolution.SealevelchangePolarMotionX / (1 / (1 - tide_love_k2 / tide_love_k2secular) * (1 + load_love_k2) / (moi_p - moi_e))
     118moiyz = md.results.TransientSolution.SealevelchangePolarMotionY / (1 / (1 - tide_love_k2 / tide_love_k2secular) * (1 + load_love_k2) / (moi_p - moi_e))
     119moizz = md.results.TransientSolution.SealevelchangePolarMotionZ / ( -(1 + load_love_k2) / moi_p)
    120120
    121121areaice = md.results.TransientSolution.SealevelBarystaticIceArea
     122areaice[np.isnan(areaice)] = 0
     123print(np.isnan(areaice))
     124print(np.sum(areaice))
    122125loadice = md.results.TransientSolution.SealevelBarystaticIceLoad
    123 
    124 # analytical moi = > just checking FOR ICE only!!! {{{
    125 # ...have to mute ** slc induced MOI in Tria.cpp**prior to the comparison
    126126rad_e = md.solidearth.planetradius
    127127
    … …  
    130130moi_xz = sum(-loadice * areaice * pow(rad_e, 2) * np.sin(lat) * np.cos(lat) * np.cos(lon))
    131131moi_yz = sum(-loadice * areaice * pow(rad_e, 2) * np.sin(lat) * np.cos(lat) * np.sin(lon))
    132 moi_zz = sum(-loadice * areaice * pow(rad_e, 2) * (1 - np.sin(lat) ** 2))
    133 theoretical_value_check = [moixz / moi_xz, moiyz / moi_yz, moizz / moi_zz]
    134 print('\ntheoretical_value_check =\n')
    135 print('\t{}\n'.format(theoretical_value_check))
     132moi_zz = sum(-loadice * areaice * pow(rad_e, 2) * (-1.0 / 3.0 + np.sin(lat) ** 2))
     133theoretical_value_check = [moixz / moi_xz, moiyz / moi_yz, moizz / moi_zz] # Should yield [1.0, 1.0, 1.0]
    136134# }}}
    137135

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

    r26800 r26840  
    8080
    8181%Fields and tolerances to track changes
    82 field_names     ={'U_AB2dA1','URate_AB2dA1','U_AB2dA2','URate_AB2dA2','U_AB2dA3','URate_AB2dA3'};
    83 field_tolerances={1e-13,1e-13,1e-13,1e-13,1e-13,1e-13};
    84 field_values={U_AB2dA1,URate_AB2dA1,U_AB2dA2,URate_AB2dA2,U_AB2dA3,URate_AB2dA3};
     82field_names     ={'U_AB2dA1','U_AB2dA2','U_AB2dA3'};
     83field_tolerances={1e-13,1e-13,1e-13};
     84field_values={U_AB2dA1,U_AB2dA2,U_AB2dA3};
     85% field_names     ={'U_AB2dA1','URate_AB2dA1','U_AB2dA2','URate_AB2dA2','U_AB2dA3','URate_AB2dA3'};
     86% field_tolerances={1e-13,1e-13,1e-13,1e-13,1e-13,1e-13};
     87% field_values={U_AB2dA1,URate_AB2dA1,U_AB2dA2,URate_AB2dA2,U_AB2dA3,URate_AB2dA3};
    8588

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

    r25300 r26840  
    1616md = parameterize(md, '../Par/GiaIvinsBenchmarksAB.py')
    1717
    18 # indicate what you want to compute
    19 md.gia.cross_section_shape = 1  # for square-edged x-section
     18# GIA Ivins, 2 layer model
     19md.solidearth.settings.grdmodel = 2
     20md.solidearth.settings.isgrd = 1
     21md.initialization.sealevel = np.zeros((md.mesh.numberofvertices,))
    2022
    21 # evaluation time (termed start_time)
    22 md.timestepping.start_time = 2002100  # after 2 kyr of deglaciation
    23 md.timestepping.final_time = 2500000  # 2,500 kyr
     23# Indicate what you want to compute
     24md.solidearth.settings.cross_section_shape = 1 # For square-edged x-section
     25md.solidearth.settings.grdocean = 0 # Do not compute sea level, only deformation
     26md.solidearth.settings.sealevelloading = 0 # Do not compute sea level, only deformation
    2427
    25 # define loading history
    26 md.geometry.thickness = np.array([
     28# Evaluation time (termed start_time)
     29md.timestepping.time_step = 2002100 # after 2 kyr of deglaciation
     30md.timestepping.start_time = -md.timestepping.time_step # Need one time step before t = 0 to generate a thickness change at t = 0
     31md.timestepping.final_time = 2500000 # 2,500 kyr
     32
     33# Define loading history
     34md.masstransport.spcthickness = np.array([
    2735    np.append(md.geometry.thickness * 0.0, 0.0),
    2836    np.append(md.geometry.thickness, 1000),
    2937    np.append(md.geometry.thickness, 2000000),
    3038    np.append(md.geometry.thickness * 0.0, 2000100),
    31     np.append(md.geometry.thickness * 0.0, md.timestepping.start_time)
     39    np.append(md.geometry.thickness * 0.0, md.timestepping.start_time + 2 * md.timestepping.time_step)
    3240    ]).T
    3341
    34 # find out the elements that have zero loads throughout the loading history
     42md.geometry.bed = np.zeros((md.mesh.numberofvertices,))
     43
     44# Find out the elements that have zero loads throughout the loading history
    3545pos = np.where(np.abs(md.geometry.thickness[0:-2, :].sum(axis=1)) == 0)[0]
    36 md.mask.ice_levelset[pos] = 1 # no ice
     46md.mask.ice_levelset[pos] = 1 # No ice
     47
     48# Physics
     49md.transient.issmb = 0
     50md.transient.isstressbalance = 0
     51md.transient.isthermal = 0
     52md.transient.ismasstransport = 1
     53md.transient.isslc = 1
    3754
    3855md.cluster = generic('name', gethostname(), 'np', 3)
    3956md.verbose = verbose('1111111')
     57md.solidearth.requested_outputs = ['Sealevel', 'BedGRD']
    4058
    41 # solve for GIA deflection
    42 md = solve(md, 'Gia')
     59# Solve for GIA deflection
     60md = solve(md, 'Transient')
    4361
    4462# Test Name: GiaIvinsBenchmarksAB2dA1
    45 U_AB2dA1 = md.results.GiaSolution.UGia
    46 URate_AB2dA1 = md.results.GiaSolution.UGiaRate
     63U_AB2dA1 = md.results.TransientSolution.BedGRD
     64#URate_AB2dA1 = md.results.TransientSolution.UGiaRate
    4765
    4866# Test Name: GiaIvinsBenchmarksAB2dA2
    4967# different evaluation time # {{{
    50 md.timestepping.start_time = 2005100 # after 5 kyr of deglaciation
    51 md.geometry.thickness[-1, -1] = md.timestepping.start_time
     68md.timestepping.time_step = 2005100 # After 5 kyr of deglaciation
     69md.timestepping.start_time = -md.timestepping.time_step # Need one time step before t = 0 to generate a thickness change at t = 0
     70md.masstransport.spcthickness[-1, -1] = md.timestepping.start_time + 2 * md.timestepping.time_step
    5271
    53 md = solve(md, 'Gia')
     72md = solve(md, 'Transient')
    5473
    55 U_AB2dA2 = md.results.GiaSolution.UGia
    56 URate_AB2dA2 = md.results.GiaSolution.UGiaRate
     74U_AB2dA2 = md.results.TransientSolution.BedGRD
     75#URate_AB2dA2 = md.results.TransientSolution.BedGRDRate
    5776# }}}
    5877
    5978# Test Name: GiaIvinsBenchmarksAB2dA3
    6079# different evaluation time # {{{
    61 md.timestepping.start_time = 2010100 # after 10 kyr of deglaciation
    62 md.geometry.thickness[-1, -1] = md.timestepping.start_time
     80md.timestepping.time_step = 2010100 # After 10 kyr of deglaciation
     81md.timestepping.start_time = -md.timestepping.time_step # Need one time step before t = 0 to generate a thickness change at t = 0
     82md.masstransport.spcthickness[-1, -1] = md.timestepping.start_time + 2 * md.timestepping.time_step
    6383
    64 md = solve(md, 'Gia')
     84md = solve(md, 'Transient')
    6585
    66 U_AB2dA3 = md.results.GiaSolution.UGia
    67 URate_AB2dA3 = md.results.GiaSolution.UGiaRate
     86U_AB2dA3 = md.results.TransientSolution.BedGRD
     87#URate_AB2dA3 = md.results.TransientSolution.UGiaRate
    6888# }}}
    6989
    7090# Fields and tolerances to track changes
    71 field_names = ['U_AB2dA1','URate_AB2dA1','U_AB2dA2','URate_AB2dA2','U_AB2dA3','URate_AB2dA3']
    72 field_tolerances = [1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13]
    73 field_values = [U_AB2dA1, URate_AB2dA1, U_AB2dA2, URate_AB2dA2, U_AB2dA3, URate_AB2dA3]
     91field_names = ['U_AB2dA1','U_AB2dA2','U_AB2dA3']
     92field_tolerances = [1e-13, 1e-13, 1e-13]
     93field_values = [U_AB2dA1, U_AB2dA2, U_AB2dA3]
     94# field_names = ['U_AB2dA1','URate_AB2dA1','U_AB2dA2','URate_AB2dA2','U_AB2dA3','URate_AB2dA3']
     95# field_tolerances = [1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13]
     96# field_values = [U_AB2dA1, URate_AB2dA1, U_AB2dA2, URate_AB2dA2, U_AB2dA3, URate_AB2dA3]

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

    r26638 r26840  
    1515
    1616md = model()
    17 md.cluster = generic('name', gethostname(), 'np', 1)
     17md.cluster = generic('name', gethostname(), 'np', 8)
    1818
    19 # Set validation=1 for comparing against the Spada benchark
     19# Set validation=1 for comparing against the Spada benchmark
    2020validation = 0
    2121
    … …  
    2525
    2626md.verbose = verbose('all')
     27md.verbose = verbose('1111111111111111')
    2728cst = 365.25 * 24 * 3600 * 1000
    2829
    … …  
    4849
    4950md.love.allow_layer_deletion = 1
    50 md.love.frequencies = (np.array([0]) * 2 * np.pi).reshape(-1, 1) / cst
     51md.love.frequencies = np.vstack(([0], np.logspace(-6, 3, 1000).reshape(-1, 1) / cst))
    5152md.love.nfreq = len(md.love.frequencies)
    5253md.love.sh_nmin = 1
    53 md.love.sh_nmax = 256
     54md.love.sh_nmax = 1000
    5455md.love.underflow_tol = 1e-20
     56md.love.pw_threshold = 1e-3
    5557md.love.Gravitational_Constant = 6.6732e-11
    56 md.love.integration_steps_per_layer = 200
     58md.love.integration_steps_per_layer = 100
     59md.love.allow_layer_deletion = 1
     60md.love.forcing_type = 11
     61md.love.chandler_wobble = 0
     62md.love.complex_computation = 0
    5763
    5864md.love.istemporal = 1
    5965md.love.n_temporal_iterations = 8
    60 #md.love.time = np.logspace(-4, 5, 2).reshape(-1, 1) * cst
    61 md.love.time = np.logspace(-1, 2, 50).reshape(-1, 1) * cst
     66#md.love.time = np.logspace(-6, 5, 2).reshape(-1, 1) * cst
     67md.love.time = np.vstack(([0], np.logspace(-3, 5, 99).reshape(-1, 1) * cst))
     68
     69#md.love.time = np.linspace(1/12, 10, 10 * 12).reshape(-1, 1) * cst / 1e3
    6270md.love.love_kernels = 1
    6371if md.love.istemporal:
    64     md.love = md.love.build_frequencies_from_time
     72    md.love = md.love.build_frequencies_from_time()
    6573
    6674md = solve(md, 'lv')
    6775
    68 ht2 = md.results.LoveSolution.LoveHr
    69 lt2 = md.results.LoveSolution.LoveLr
    70 kt2 = md.results.LoveSolution.LoveKr
     76ht = md.results.LoveSolution.LoveHt.reshape(-1, 1)
     77lt = md.results.LoveSolution.LoveLt.reshape(-1, 1)
     78kt = md.results.LoveSolution.LoveKt.reshape(-1, 1)
    7179t = md.love.time / cst * 1e3
    7280
    73 #Fields and tolerances to track changes
    74 #loading love numbers
     81# Fields and tolerances to track changes
     82# loading love numbers
    7583field_names = ['LoveH_loading_elastic', 'LoveK_loading_elastic', 'LoveL_loading_elastic']
    7684field_tolerances = [2.0e-8, 2.0e-8, 2.0e-8]
    7785field_values = [
    78     np.array(md.results.LoveSolution.LoveHr)[:, 0],
    79     np.array(md.results.LoveSolution.LoveKr)[:, 0],
    80     np.array(md.results.LoveSolution.LoveLr)[:, 0]
     86    np.array(md.results.LoveSolution.LoveHt)[:, 0],
     87    np.array(md.results.LoveSolution.LoveKt)[:, 0],
     88    np.array(md.results.LoveSolution.LoveLt)[:, 0]
    8189]
    82 
    83 # Validate elastic loading solutions against the Spada benchmark {{{
    8490
    8591# TODO:
    8692# - Implement read from file and comparison
    8793# - Implement plot
    88 #
    89 
    90 #}}}
    91 
    92 md.love.frequencies = (np.array([1e-3, 1e-2, 1e-1, 1, -1e-3, -1e-2, -1e-1, -1]) * 2 * np.np.pi).reshape(-1, 1) / cst
    93 md.love.nfreq = len(md.love.frequencies)
    94 md.love.sh_nmax = 256
    95 md.materials.burgers_mu = md.materials.lame_mu
    96 md.materials.burgers_viscosity = md.materials.viscosity
    97 
    98 md = solve(md, 'lv')
    99 
    100 #Fields and tolerances to track changes
    101 field_names += ['LoveH_loading_realpart', 'LoveK_loading_realpart', 'LoveL_loading_realpart', 'LoveH_loading_imagpart', 'LoveK_loading_imagpart', 'LoveL_loading_imagpart']
    102 field_tolerances += [5e-7, 5e-7, 5e-7, 5e-7, 5e-7, 5e-7]
    103 field_values += [
    104     np.array(md.results.LoveSolution.LoveHr),
    105     np.array(md.results.LoveSolution.LoveKr),
    106     np.array(md.results.LoveSolution.LoveLr),
    107     np.array(md.results.LoveSolution.LoveHi),
    108     np.array(md.results.LoveSolution.LoveKi),
    109     np.array(md.results.LoveSolution.LoveLi)
    110 ]
    111 
    112 md.love.forcing_type = 9
    113 md.love.sh_nmin = 2
    114 md.love.frequencies = ((np.array([0, 1e-3, 1e-2, 1e-1, 1, -1e-3, -1e-2, -1e-1, -1]) * 2 * np.pi).reshape(-1, 1) / cst)
    115 md.love.nfreq = len(md.love.frequencies)
    116 md = solve(md, 'lv')
    117 
    118 # Validate elastic tidal solutions against the Spada benchmark #{{{
    119 
    120 # TODO:
    121 # - Implement read from file and comparison
    122 # - Implement plot
    123 #
    124 
    125 #}}}
    126 
    127 #tidal love numbers, check
    128 field_names += ['LoveH_tidal_elastic', 'LoveK_tidal_elastic', 'LoveL_tidal_elastic', 'LoveH_tidal_realpart', 'LoveK_tidal_realpart', 'LoveL_tidal_realpart', 'LoveH_tidal_imagpart', 'LoveK_tidal_imagpart', 'LoveL_tidal_imagpart']
    129 field_tolerances += [8e-6, 8e-6, 8e-6, 8e-6, 8e-6, 8e-6, 8e-6, 8e-6, 8e-6]
    130 field_values += [
    131     np.array(md.results.LoveSolution.LoveHr)[:, 0],
    132     np.array(md.results.LoveSolution.LoveKr)[:, 0],
    133     np.array(md.results.LoveSolution.LoveLr)[:, 0],
    134     np.array(md.results.LoveSolution.LoveHr)[:, 1:],
    135     np.array(md.results.LoveSolution.LoveKr)[:, 1:],
    136     np.array(md.results.LoveSolution.LoveLr)[:, 1:],
    137     np.array(md.results.LoveSolution.LoveHi)[:, 1:],
    138     np.array(md.results.LoveSolution.LoveKi)[:, 1:],
    139     np.array(md.results.LoveSolution.LoveLi)[:, 1:]
    140 ]
    141 
    142 # Many layers PREM-based model
    143 #data = load('../Data/PREM_500layers')
    144 #md.love.sh_nmin = 1
    145 #md.materials.radius = data(2:end - 2, 1)
    146 #md.materials.density = data(3:end - 2, 2)
    147 #md.materials.lame_lambda = data(3:end - 2, 3)
    148 #md.materials.lame_mu = data(3:end - 2, 4)
    149 #md.materials.issolid = data(3:end - 2, 4) > 0
    150 #ind = find(md.materials.issolid = 0)
    151 #md.materials.density(ind(1))=sum((md.materials.radius(ind + 1).^3 - md.materials.radius(ind).^3). * md.materials.density(ind)) / (md.materials.radius(ind(end) + 1).^3 - md.materials.radius(ind(1) + 1).^3)
    152 #md.materials.lame_lambda(ind(1))=sum((md.materials.radius(ind + 1).^3 - md.materials.radius(ind).^3). * md.materials.lame_lambda(ind)) / (md.materials.radius(ind(end) + 1).^3 - md.materials.radius(ind(1) + 1).^3)
    153 #md.materials.lame_mu(ind(1))=sum((md.materials.radius(ind + 1).^3 - md.materials.radius(ind).^3). * md.materials.lame_mu(ind)) / (md.materials.radius(ind(end) + 1).^3 - md.materials.radius(ind(1)).^3)
    154 #md.materials.radius(ind(2:end) + 1) = []
    155 #md.materials.density(ind(2:end)) = []
    156 #md.materials.lame_lambda(ind(2:end)) = []
    157 #md.materials.lame_mu(ind(2:end)) = []
    158 #md.materials.issolid(ind(2:end)) = []
    159 #md.materials.viscosity = 10.^interp1([0 3479e3 3480e3 3680e3 5720e3 5800e3 6270e3 6371e3], log10([1e8 1e8 5e21 1e23 1e22 1e20 1e21 1e40]), md.materials.radius(2:end), 'PCHIP')
    160 #md.materials.viscosity = md.materials.viscosity. * md.materials.issolid
    161 #md.materials.burgers_mu = md.materials.lame_mu
    162 #md.materials.burgers_viscosity = md.materials.viscosity
    163 #md.materials.rheologymodel = md.materials.issolid * 0
    164 #md.love.forcing_type = 11
    165 #md.materials.numlayers = len(md.materials.viscosity)
    166 #md = solve(md, 'lv')
    167 #
    168 #field_names = [field_names, 'LoveH_loadingVSS96_elastic', 'LoveK_loadingVSS96_elastic', 'LoveL_loadingVSS96_elastic', 'LoveH_loadingVSS96_realpart', 'LoveK_loadingVSS96_realpart', 'LoveL_loadingVSS96_realpart', 'LoveH_loadingVSS96_imagpart', 'LoveK_loadingVSS96_imagpart', 'LoveL_loadingVSS96_imagpart']
    169 #field_tolerances = [field_tolerances, 4.3e-9, 4.3e-9, 4.3e-9, 4.3e-9, 4.3e-9, 4.3e-9, 4.3e-9, 4.3e-9, 4.3e-9]
    170 #field_values = [field_values,
    171 #       (md.results.LoveSolution.LoveHr[:][0]),
    172 #       (md.results.LoveSolution.LoveKr[:][0]),
    173 #       (md.results.LoveSolution.LoveLr[:][0]),
    174 #       (md.results.LoveSolution.LoveHr[:][1:]),
    175 #       (md.results.LoveSolution.LoveKr[:][1:]),
    176 #       (md.results.LoveSolution.LoveLr[:][1:]),
    177 #       (md.results.LoveSolution.LoveHi[:][1:]),
    178 #       (md.results.LoveSolution.LoveKi[:][1:]),
    179 #       (md.results.LoveSolution.LoveLi[:][1:]),
    180 #       ]
    181 
    182 # Model VSS96 from Vermeersen, L.L.A., Sabadini, R. & Spada, G., 1996a. Analytical visco-elastic relaxation models, Geophys. Res. Lett., 23, 697-700.
    183 md.materials.radius = np.array([10, 1222.5, 3480., 3600., 3630.5, 3700., 3900., 4000.,
    184                                 4200., 4300., 4500., 4600., 4800., 4900., 5100., 5200.,
    185                                 5400., 5500., 5600.5, 5650., 5701., 5736., 5771.5,
    186                                 5821., 5951., 5970.5, 6016., 6061., 6150.5, 6151.5,
    187                                 6251., 6371.]).reshape(-1, 1) * 1e3
    188 md.materials.lame_mu = np.array([1e-5, 0., 2.933, 2.8990002, 2.8550003, 2.7340002, 2.675,
    189                                 2.559, 2.502, 2.388, 2.331, 2.215, 2.157, 2.039, 1.979,
    190                                 1.8560001, 1.794, 1.73, 1.639, 1.2390001, 1.224, 1.21,
    191                                 1.128, 0.97700006, 0.906, 0.79, 0.773, 0.741, 0.656, 0.665,
    192                                 0.602]).reshape(-1, 1) * 1e11
    193 md.materials.density = np.array([10925., 10925., 5506.42, 5491.45, 5456.57, 5357.06,
    194                                 5307.24, 5207.13, 5156.69, 5054.69, 5002.99, 4897.83,
    195                                 4844.22, 4734.6, 4678.44, 4563.07, 4503.72, 4443.16,
    196                                 4412.41, 3992.14, 3983.99, 3975.84, 3912.82, 3786.78,
    197                                 3723.78, 3516.39, 3489.51, 3435.78, 3359.5, 3367.1,
    198                                 3184.3]).reshape(-1, 1)
    199 md.materials.viscosity = np.array([0., 0., 7.999999999999999e+21, 8.5e+21,
    200                                    8.999999999999999e+21, 3.e+22, 4.e+22,
    201                                    5.0000000000000004e+22, 6.e+22,
    202                                    5.0000000000000004e+22, 4.5e+22, 3.e+22,
    203                                    2.5000000000000002e+22, 1.7999999999999998e+22,
    204                                    1.3e+22, 7.999999999999999e+21, 6.999999999999999e+21,
    205                                    6.5e+21, 6.e+21, 5.5e+21, 5.e+21, 4.4999999999999995e+21,
    206                                    3.9999999999999995e+21, 2.5e+21,
    207                                    1.9999999999999997e+21, 1.5e+21, 9.999999999999999e+20,
    208                                    6.e+20, 5.5000000000000007e+20, 2.e+20,
    209                                    1.E40]).reshape(-1, 1)
    210 md.materials.lame_lambda = np.array(md.materials.lame_mu) * 0 + 5e14
    211 md.materials.issolid = np.ones(len(md.materials.lame_mu)).reshape(-1, 1)
    212 md.materials.issolid[1] = 0
    213 md.materials.numlayers = len(md.materials.lame_mu)
    214 md.materials.burgers_mu = md.materials.lame_mu
    215 md.materials.burgers_viscosity = md.materials.viscosity
    216 md.materials.rheologymodel = md.materials.issolid * 0
    217 md.love.forcing_type = 11
    218 md.love.sh_nmin = 1
    219 md.love.sh_nmax = 100
    220 md = solve(md, 'lv')
    221 md.love.frequencies = (np.array([0, 1e-3, 1e-2, 1, -1e-3, -1e-2, -1]) * 2 * np.pi).reshape(-1, 1) / cst
    222 md.love.nfreq = len(md.love.frequencies)
    223 
    224 field_names += [
    225     'LoveH_loadingVSS96_elastic',
    226     'LoveK_loadingVSS96_elastic',
    227     'LoveL_loadingVSS96_elastic',
    228     'LoveH_loadingVSS96_realpart',
    229     'LoveK_loadingVSS96_realpart',
    230     'LoveL_loadingVSS96_realpart',
    231     'LoveH_loadingVSS96_imagpart',
    232     'LoveK_loadingVSS96_imagpart',
    233     'LoveL_loadingVSS96_imagpart'
    234 ]
    235 field_tolerances += [2e-6, 2e-6, 2e-6, 2e-6, 2e-6, 2e-6, 2e-6, 2e-6, 2e-6]
    236 field_values += [
    237     np.array(md.results.LoveSolution.LoveHr)[:, 0],
    238     np.array(md.results.LoveSolution.LoveKr)[:, 0],
    239     np.array(md.results.LoveSolution.LoveLr)[:, 0],
    240     np.array(md.results.LoveSolution.LoveHr)[:, 1:],
    241     np.array(md.results.LoveSolution.LoveKr)[:, 1:],
    242     np.array(md.results.LoveSolution.LoveLr)[:, 1:],
    243     np.array(md.results.LoveSolution.LoveHi)[:, 1:],
    244     np.array(md.results.LoveSolution.LoveKi)[:, 1:],
    245     np.array(md.results.LoveSolution.LoveLi)[:, 1:]
    246 ]
     94# - Implement validation of elastic loading solutions against the Spada benchmark

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

    r26800 r26840  
    88
    99% for volumetric potential
    10         md=model();
    11         md.groundingline.migration='None';
    12        
    13         md.materials=materials('litho');
    14         cst=365.25*24*3600*1000;
    15 
    16         md.materials.numlayers = 40;
    17         md.love.forcing_type = 9;
    18 
    19         md.materials.density=zeros(md.materials.numlayers,1)+5511;
    20         md.materials.lame_mu=zeros(md.materials.numlayers,1)+0.75e11;
    21         md.materials.lame_lambda=zeros(md.materials.numlayers,1)+5e17;
    22         md.materials.issolid=ones(md.materials.numlayers,1);
    23         md.materials.rheologymodel=zeros(md.materials.numlayers,1);
    24 
    25         %the following isn't used here but needs to have arrays of consistent size with the rest of the materials
    26         md.materials.viscosity=zeros(md.materials.numlayers,1)+1e21;
    27         md.materials.burgers_mu=md.materials.lame_mu/3;
    28         md.materials.burgers_viscosity=md.materials.viscosity/10;
    29         md.materials.ebm_alpha= ones(md.materials.numlayers,1)*.9;
    30         md.materials.ebm_delta= ones(md.materials.numlayers,1)*0.2;
    31         md.materials.ebm_taul= ones(md.materials.numlayers,1)*54*60; %54min
    32         md.materials.ebm_tauh= ones(md.materials.numlayers,1)*18.6*cst/1e3; %18.6yr
    33 
    34 
    35         md.materials.radius =  linspace(10e3,6371e3,md.materials.numlayers+1)';
    36         md.love.g0 = 9.8134357285509388; % directly grabbed from fourierlovesolver for this particular case.
    37 
    38         md.love.allow_layer_deletion=1;
    39         md.love.frequencies=0;
    40         md.love.nfreq=1;
    41         md.love.istemporal=0;
    42 
    43         md.love.sh_nmin = 2;
    44         md.love.sh_nmax = 200;
    45         md.love.love_kernels=1;
    46 
    47         md.miscellaneous.name='kernels';
    48         md.cluster=generic('name',oshostname(),'np',3);
    49         md.verbose=verbose('111111101');
    50        
    51         md=solve(md,'lv');
    52 
    53         % save yi's for all layers except for the inner-most one, at select degrees.
    54         degrees = [2 20 200];  % we archive solutions for degrees 2, 20, 200
    55         kernels=reshape(md.results.LoveSolution.LoveKernels, [md.love.sh_nmax+1 md.love.nfreq md.materials.numlayers+1 6]);
    56         % extract love kernels {{{
    57         % degree 2.
    58         y1_tidal_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,1));
    59         y2_tidal_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,2));
    60         y3_tidal_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,3));
    61         y4_tidal_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,4));
    62         y5_tidal_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,5));
    63         y6_tidal_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,6));
    64 
    65         % degree 20.
    66         y1_tidal_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,1));
    67         y2_tidal_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,2));
    68         y3_tidal_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,3));
    69         y4_tidal_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,4));
    70         y5_tidal_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,5));
    71         y6_tidal_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,6));
    72 
    73         % degree 200.
    74         y1_tidal_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,1));
    75         y2_tidal_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,2));
    76         y3_tidal_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,3));
    77         y4_tidal_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,4));
    78         y5_tidal_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,5));
    79         y6_tidal_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,6));
    80         % }}}
     10md=model();
     11md.groundingline.migration='None';
     12
     13md.materials=materials('litho');
     14cst=365.25*24*3600*1000;
     15
     16md.materials.numlayers = 40;
     17md.love.forcing_type = 9;
     18
     19md.materials.density=zeros(md.materials.numlayers,1)+5511;
     20md.materials.lame_mu=zeros(md.materials.numlayers,1)+0.75e11;
     21md.materials.lame_lambda=zeros(md.materials.numlayers,1)+5e17;
     22md.materials.issolid=ones(md.materials.numlayers,1);
     23md.materials.rheologymodel=zeros(md.materials.numlayers,1);
     24
     25%the following isn't used here but needs to have arrays of consistent size with the rest of the materials
     26md.materials.viscosity=zeros(md.materials.numlayers,1)+1e21;
     27md.materials.burgers_mu=md.materials.lame_mu/3;
     28md.materials.burgers_viscosity=md.materials.viscosity/10;
     29md.materials.ebm_alpha= ones(md.materials.numlayers,1)*0.9;
     30md.materials.ebm_delta= ones(md.materials.numlayers,1)*0.2;
     31md.materials.ebm_taul= ones(md.materials.numlayers,1)*54*60; %54min
     32md.materials.ebm_tauh= ones(md.materials.numlayers,1)*18.6*cst/1e3; %18.6yr
     33
     34md.materials.radius =  linspace(10e3,6371e3,md.materials.numlayers+1)';
     35md.love.g0 = 9.8134357285509388; % directly grabbed from fourierlovesolver for this particular case.
     36
     37md.love.allow_layer_deletion=1;
     38md.love.frequencies=0;
     39md.love.nfreq=1;
     40md.love.istemporal=0;
     41
     42md.love.sh_nmin = 2;
     43md.love.sh_nmax = 200;
     44md.love.love_kernels=1;
     45
     46md.miscellaneous.name='kernels';
     47md.cluster=generic('name',oshostname(),'np',3);
     48md.verbose=verbose('111111101');
     49
     50md=solve(md,'lv');
     51
     52% save yi's for all layers except for the inner-most one, at select degrees.
     53degrees = [2 20 200];  % we archive solutions for degrees 2, 20, 200
     54kernels=reshape(md.results.LoveSolution.LoveKernels, [md.love.sh_nmax+1 md.love.nfreq md.materials.numlayers+1 6]);
     55
     56% extract love kernels {{{
     57% degree 2.
     58y1_tidal_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,1));
     59y2_tidal_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,2));
     60y3_tidal_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,3));
     61y4_tidal_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,4));
     62y5_tidal_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,5));
     63y6_tidal_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,6));
     64
     65% degree 20.
     66y1_tidal_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,1));
     67y2_tidal_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,2));
     68y3_tidal_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,3));
     69y4_tidal_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,4));
     70y5_tidal_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,5));
     71y6_tidal_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,6));
     72
     73% degree 200.
     74y1_tidal_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,1));
     75y2_tidal_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,2));
     76y3_tidal_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,3));
     77y4_tidal_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,4));
     78y5_tidal_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,5));
     79y6_tidal_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,6));
     80% }}}
    8181
    8282% validate tidal potential solutions against the analytic solutions. {{{
    … …  
    176176        legend(axes6,'n=2','n=4','n=8','n=16','n=32');
    177177        %export_fig('/Users/adhikari/issm/trunk-jpl/src/m/contrib/adhikari/issm_vs_analytic_loading_homogeneous.pdf');
    178        
    179178else
    180179        %
    … …  
    183182
    184183% for surface load.
    185 
    186         md.love.forcing_type = 11;
    187         md=solve(md,'lv');
    188         kernels=reshape(md.results.LoveSolution.LoveKernels, [md.love.sh_nmax+1 md.love.nfreq md.materials.numlayers+1 6]);
    189 
    190         % extract love kernels {{{
    191         % degree 2.
    192         y1_loading_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,1));
    193         y2_loading_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,2));
    194         y3_loading_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,3));
    195         y4_loading_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,4));
    196         y5_loading_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,5));
    197         y6_loading_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,6));
    198 
    199         % degree 20.
    200         y1_loading_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,1));
    201         y2_loading_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,2));
    202         y3_loading_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,3));
    203         y4_loading_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,4));
    204         y5_loading_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,5));
    205         y6_loading_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,6));
    206 
    207         % degree 200.
    208         y1_loading_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,1));
    209         y2_loading_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,2));
    210         y3_loading_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,3));
    211         y4_loading_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,4));
    212         y5_loading_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,5));
    213         y6_loading_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,6));
    214         % }}}
     184md.love.forcing_type = 11;
     185md=solve(md,'lv');
     186kernels=reshape(md.results.LoveSolution.LoveKernels, [md.love.sh_nmax+1 md.love.nfreq md.materials.numlayers+1 6]);
     187
     188% extract love kernels {{{
     189% degree 2.
     190y1_loading_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,1));
     191y2_loading_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,2));
     192y3_loading_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,3));
     193y4_loading_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,4));
     194y5_loading_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,5));
     195y6_loading_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,6));
     196
     197% degree 20.
     198y1_loading_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,1));
     199y2_loading_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,2));
     200y3_loading_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,3));
     201y4_loading_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,4));
     202y5_loading_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,5));
     203y6_loading_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,6));
     204
     205% degree 200.
     206y1_loading_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,1));
     207y2_loading_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,2));
     208y3_loading_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,3));
     209y4_loading_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,4));
     210y5_loading_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,5));
     211y6_loading_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,6));
     212% }}}
    215213
    216214% validate loading solutions against the analytic solutions. {{{
    … …  
    303301        legend(axes6,'n=2','n=4','n=8','n=16','n=32');
    304302        %export_fig('/Users/adhikari/issm/trunk-jpl/src/m/contrib/adhikari/issm_vs_analytic_tidal_homogeneous.pdf');
    305        
    306303else
    307304        %

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

    r25342 r26840  
    1616# For volumetric potential
    1717md = model()
     18md.groundingline.migration = 'None'
    1819
    1920md.materials = materials('litho')
     21cst = 365.25 * 24 * 3600 * 1000
    2022
    2123md.materials.numlayers = 40
    … …  
    2426md.materials.density = np.zeros((md.materials.numlayers, 1)) + 5511
    2527md.materials.lame_mu = np.zeros((md.materials.numlayers, 1)) + 0.75e11
     28md.materials.lame_lambda = np.zeros((md.materials.numlayers, 1)) + 5e17
     29md.materials.issolid = np.ones((md.materials.numlayers, 1))
     30md.materials.rheologymodel = np.zeros((md.materials.numlayers, 1))
     31
     32# The following isn't used here but needs to hhave arrays of consistent size with the rest of the materials
    2633md.materials.viscosity = np.zeros((md.materials.numlayers, 1)) + 1e21
    27 md.materials.lame_lambda = np.zeros((md.materials.numlayers, 1)) + 0.5e17
    28 md.materials.issolid = np.ones((md.materials.numlayers, 1))
    29 md.materials.isburgers = np.zeros((md.materials.numlayers, 1))
    3034md.materials.burgers_mu = md.materials.lame_mu / 3
    3135md.materials.burgers_viscosity = md.materials.viscosity / 10
     36md.materials.ebm_alpha = np.ones((md.materials.numlayers, 1)) * 0.9
     37md.materials.ebm_delta = np.ones((md.materials.numlayers, 1)) * 0.2
     38md.materials.ebm_taul = np.ones((md.materials.numlayers, 1)) * 54 * 60 # 54 min
     39md.materials.ebm_tauh = np.ones((md.materials.numlayers, 1)) * 18.6 * cst / 1e3 # 18.6 yr
    3240
    3341md.materials.radius = np.linspace(10e3, 6371e3, md.materials.numlayers + 1).reshape(-1, 1)
    … …  
    3745md.love.frequencies = 0
    3846md.love.nfreq = 1
     47md.love.istemporal = 0
    3948
    4049md.love.sh_nmin = 2
    … …  
    4352
    4453md.miscellaneous.name = 'kernels'
    45 md.cluster = generic('name', gethostname(), 'np', 1)
     54md.cluster = generic('name', gethostname(), 'np', 3)
    4655md.verbose = verbose('111111101')
    4756
    … …  
    4958
    5059# Save yi's for all layers except for the inner-most one, at select degrees.
    51 degrees = [2, 20, 200] # we archive solutions for degrees 2, 20, 200
     60degrees = [1, 19, 199] # we archive solutions for degrees 2, 20, 200
     61kernels = np.reshape(md.results.LoveSolution.LoveKernels, (md.love.sh_nmax + 1, md.love.nfreq, md.materials.numlayers + 1, 6), order='F')
    5262
    5363# Extract love kernels #{{{
    5464# degree 2
    55 y1_tidal_degree002 = md.results.LoveSolution.LoveKernelsReal[degrees[0],0,1:,0].squeeze()
    56 y2_tidal_degree002 = md.results.LoveSolution.LoveKernelsReal[degrees[0],0,1:,1].squeeze()
    57 y3_tidal_degree002 = md.results.LoveSolution.LoveKernelsReal[degrees[0],0,1:,2].squeeze()
    58 y4_tidal_degree002 = md.results.LoveSolution.LoveKernelsReal[degrees[0],0,1:,3].squeeze()
    59 y5_tidal_degree002 = md.results.LoveSolution.LoveKernelsReal[degrees[0],0,1:,4].squeeze()
    60 y6_tidal_degree002 = md.results.LoveSolution.LoveKernelsReal[degrees[0],0,1:,5].squeeze()
     65y1_tidal_degree002 = kernels[degrees[0]+1,0,1:,0].squeeze()
     66y2_tidal_degree002 = kernels[degrees[0]+1,0,1:,1].squeeze()
     67y3_tidal_degree002 = kernels[degrees[0]+1,0,1:,2].squeeze()
     68y4_tidal_degree002 = kernels[degrees[0]+1,0,1:,3].squeeze()
     69y5_tidal_degree002 = kernels[degrees[0]+1,0,1:,4].squeeze()
     70y6_tidal_degree002 = kernels[degrees[0]+1,0,1:,5].squeeze()
    6171
    6272# degree 20
    63 y1_tidal_degree020 = md.results.LoveSolution.LoveKernelsReal[degrees[1],0,1:,0].squeeze()
    64 y2_tidal_degree020 = md.results.LoveSolution.LoveKernelsReal[degrees[1],0,1:,1].squeeze()
    65 y3_tidal_degree020 = md.results.LoveSolution.LoveKernelsReal[degrees[1],0,1:,2].squeeze()
    66 y4_tidal_degree020 = md.results.LoveSolution.LoveKernelsReal[degrees[1],0,1:,3].squeeze()
    67 y5_tidal_degree020 = md.results.LoveSolution.LoveKernelsReal[degrees[1],0,1:,4].squeeze()
    68 y6_tidal_degree020 = md.results.LoveSolution.LoveKernelsReal[degrees[1],0,1:,5].squeeze()
     73y1_tidal_degree020 = kernels[degrees[1]+1,0,1:,0].squeeze()
     74y2_tidal_degree020 = kernels[degrees[1]+1,0,1:,1].squeeze()
     75y3_tidal_degree020 = kernels[degrees[1]+1,0,1:,2].squeeze()
     76y4_tidal_degree020 = kernels[degrees[1]+1,0,1:,3].squeeze()
     77y5_tidal_degree020 = kernels[degrees[1]+1,0,1:,4].squeeze()
     78y6_tidal_degree020 = kernels[degrees[1]+1,0,1:,5].squeeze()
    6979
    7080# degree 200
    71 y1_tidal_degree200 = md.results.LoveSolution.LoveKernelsReal[degrees[2],0,1:,0].squeeze()
    72 y2_tidal_degree200 = md.results.LoveSolution.LoveKernelsReal[degrees[2],0,1:,1].squeeze()
    73 y3_tidal_degree200 = md.results.LoveSolution.LoveKernelsReal[degrees[2],0,1:,2].squeeze()
    74 y4_tidal_degree200 = md.results.LoveSolution.LoveKernelsReal[degrees[2],0,1:,3].squeeze()
    75 y5_tidal_degree200 = md.results.LoveSolution.LoveKernelsReal[degrees[2],0,1:,4].squeeze()
    76 y6_tidal_degree200 = md.results.LoveSolution.LoveKernelsReal[degrees[2],0,1:,5].squeeze()
     81y1_tidal_degree200 = kernels[degrees[2]+1,0,1:,0].squeeze()
     82y2_tidal_degree200 = kernels[degrees[2]+1,0,1:,1].squeeze()
     83y3_tidal_degree200 = kernels[degrees[2]+1,0,1:,2].squeeze()
     84y4_tidal_degree200 = kernels[degrees[2]+1,0,1:,3].squeeze()
     85y5_tidal_degree200 = kernels[degrees[2]+1,0,1:,4].squeeze()
     86y6_tidal_degree200 = kernels[degrees[2]+1,0,1:,5].squeeze()
    7787#}}}
    7888
    7989# For surface load
    8090md.love.forcing_type = 11
    81 
    8291md = solve(md,'lv')
     92kernels = np.reshape(md.results.LoveSolution.LoveKernels, (md.love.sh_nmax + 1, md.love.nfreq, md.materials.numlayers + 1, 6), order='F')
    8393
    8494# Extract love kernels #{{{
    8595# degree 2
    86 y1_loading_degree002 = md.results.LoveSolution.LoveKernelsReal[degrees[0],0,1:,0].squeeze()
    87 y2_loading_degree002 = md.results.LoveSolution.LoveKernelsReal[degrees[0],0,1:,1].squeeze()
    88 y3_loading_degree002 = md.results.LoveSolution.LoveKernelsReal[degrees[0],0,1:,2].squeeze()
    89 y4_loading_degree002 = md.results.LoveSolution.LoveKernelsReal[degrees[0],0,1:,3].squeeze()
    90 y5_loading_degree002 = md.results.LoveSolution.LoveKernelsReal[degrees[0],0,1:,4].squeeze()
    91 y6_loading_degree002 = md.results.LoveSolution.LoveKernelsReal[degrees[0],0,1:,5].squeeze()
     96y1_loading_degree002 = kernels[degrees[0]+1,0,1:,0].squeeze()
     97y2_loading_degree002 = kernels[degrees[0]+1,0,1:,1].squeeze()
     98y3_loading_degree002 = kernels[degrees[0]+1,0,1:,2].squeeze()
     99y4_loading_degree002 = kernels[degrees[0]+1,0,1:,3].squeeze()
     100y5_loading_degree002 = kernels[degrees[0]+1,0,1:,4].squeeze()
     101y6_loading_degree002 = kernels[degrees[0]+1,0,1:,5].squeeze()
    92102
    93103# degree 20
    94 y1_loading_degree020 = md.results.LoveSolution.LoveKernelsReal[degrees[1],0,1:,0].squeeze()
    95 y2_loading_degree020 = md.results.LoveSolution.LoveKernelsReal[degrees[1],0,1:,1].squeeze()
    96 y3_loading_degree020 = md.results.LoveSolution.LoveKernelsReal[degrees[1],0,1:,2].squeeze()
    97 y4_loading_degree020 = md.results.LoveSolution.LoveKernelsReal[degrees[1],0,1:,3].squeeze()
    98 y5_loading_degree020 = md.results.LoveSolution.LoveKernelsReal[degrees[1],0,1:,4].squeeze()
    99 y6_loading_degree020 = md.results.LoveSolution.LoveKernelsReal[degrees[1],0,1:,5].squeeze()
     104y1_loading_degree020 = kernels[degrees[1]+1,0,1:,0].squeeze()
     105y2_loading_degree020 = kernels[degrees[1]+1,0,1:,1].squeeze()
     106y3_loading_degree020 = kernels[degrees[1]+1,0,1:,2].squeeze()
     107y4_loading_degree020 = kernels[degrees[1]+1,0,1:,3].squeeze()
     108y5_loading_degree020 = kernels[degrees[1]+1,0,1:,4].squeeze()
     109y6_loading_degree020 = kernels[degrees[1]+1,0,1:,5].squeeze()
    100110
    101111# degree 200
    102 y1_loading_degree200 = md.results.LoveSolution.LoveKernelsReal[degrees[2],0,1:,0].squeeze()
    103 y2_loading_degree200 = md.results.LoveSolution.LoveKernelsReal[degrees[2],0,1:,1].squeeze()
    104 y3_loading_degree200 = md.results.LoveSolution.LoveKernelsReal[degrees[2],0,1:,2].squeeze()
    105 y4_loading_degree200 = md.results.LoveSolution.LoveKernelsReal[degrees[2],0,1:,3].squeeze()
    106 y5_loading_degree200 = md.results.LoveSolution.LoveKernelsReal[degrees[2],0,1:,4].squeeze()
    107 y6_loading_degree200 = md.results.LoveSolution.LoveKernelsReal[degrees[2],0,1:,5].squeeze()
     112y1_loading_degree200 = kernels[degrees[2]+1,0,1:,0].squeeze()
     113y2_loading_degree200 = kernels[degrees[2]+1,0,1:,1].squeeze()
     114y3_loading_degree200 = kernels[degrees[2]+1,0,1:,2].squeeze()
     115y4_loading_degree200 = kernels[degrees[2]+1,0,1:,3].squeeze()
     116y5_loading_degree200 = kernels[degrees[2]+1,0,1:,4].squeeze()
     117y6_loading_degree200 = kernels[degrees[2]+1,0,1:,5].squeeze()
    108118#}}}
    109119

  • issm/trunk-jpl/test/Par/GiaIvinsBenchmarksAB.py

    r25300 r26840  
    66
    77from arch import *
    8 from giaivins import giaivins
    98from InterpFromMeshToMesh2d import *
     9from materials import *
    1010from paterson import *
     11from SetIceSheetBC import *
    1112from verbose import *
    12 from SetIceSheetBC import *
    1313
    1414
    … …  
    2828md.geometry.surface = md.geometry.thickness + md.geometry.base.reshape(-1, 1)  #would otherwise create a 91x91 matrix
    2929
    30 #Ice density used for benchmarking, not 917 kg/m^3
    31 md.materials.rho_ice = 1000  #kg/m^3
     30# GIA parameters specific to Experiments A and B
     31md.materials = materials('litho', 'ice')
     32md.materials.radius = [10, 6271e3, 6371e3] # 100km lithosphere, the rest is mantle material
     33md.materials.density = [3.38e3, 3.36e3]
     34md.materials.lame_u = [1.45e11, 6.7e10]
     35md.materials.viscosity = [1e21, 1e40]
    3236
    33 #GIA parameters specific to Experiments A and B
    34 md.gia=giaivins();
    35 md.gia.mantle_viscosity = 1e21 * np.ones((md.mesh.numberofvertices, 1))  #in Pa.s
    36 md.gia.lithosphere_thickness = 100 * np.ones((md.mesh.numberofvertices, 1))  #in km
    37 md.materials.lithosphere_shear_modulus = 6.7 * 1e10  #in Pa
    38 md.materials.lithosphere_density = 3.36  #in g/cm^3
    39 md.materials.mantle_shear_modulus = 1.45 * 1e11  #in Pa
    40 md.materials.mantle_density = 3.38  #in g/cm^3
     37# Ice density used for benchmarking, not 917 kg/m^3
     38md.materials.rho_ice = 1000 # kg/m^3
    4139
    42 #Initial velocity
     40# Initial velocity
    4341x = archread('../Data/SquareSheetConstrained.arch', 'x')
    4442y = archread('../Data/SquareSheetConstrained.arch', 'y')
    … …  
    5755md.initialization.pressure = np.zeros((md.mesh.numberofvertices, 1))
    5856
    59 #Materials
     57# Materials
    6058md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices, 1))
    6159md.materials.rheology_B = paterson(md.initialization.temperature)
    6260md.materials.rheology_n = 3. * np.ones((md.mesh.numberofelements, 1))
    6361
    64 #Friction
     62# Friction
    6563md.friction.coefficient = 20. * np.ones((md.mesh.numberofvertices, 1))
    6664md.friction.coefficient[np.where(md.mask.ocean_levelset < 0.)] = 0.
    … …  
    6866md.friction.q = np.ones((md.mesh.numberofelements, 1))
    6967
    70 #Numerical parameters
     68# Numerical parameters
    7169md.groundingline.migration = 'None'
    7270md.masstransport.stabilization = 1
    … …  
    8179md.timestepping.final_time = 3.
    8280
    83 #Boundary conditions:
     81# Boundary conditions:
    8482md = SetIceSheetBC(md)
    8583
    86 #Change name so that no test have the same name
     84# Change name so that no test have the same name
    8785if len(inspect.stack()) > 2:
    8886    md.miscellaneous.name = os.path.basename(inspect.stack()[2][1]).split('.')[0]
Note: See TracChangeset for help on using the changeset viewer.