Changeset 25183


Ignore:
Timestamp:
06/29/20 16:10:31 (5 years ago)
Author:
jdquinn
Message:

BUG: Corrected remaining tests (Love number retrieval)

Location:
issm/trunk-jpl
Files:
8 edited

Legend:

Unmodified
Added
Removed

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

    r25166 r25183  
    1006910069    elif frame == 'CF': #from Blewitt, 2003, JGR
    1007010070        if type == 'loadingverticaldisplacement':
    10071             series[0] = -0.269
     10071            series[1] = -0.269
    1007210072        elif type == 'loadinggravitationalpotential':
    10073             series[0] = 0.021
     10073            series[1] = 0.021
    1007410074        elif type == 'loadinghorizontaldisplacement':
    10075             series[0] = 0.134
     10075            series[1] = 0.134
    1007610076    else:
    1007710077        raise Exception("love_numbers error message: unknown reference frame: {}".format(frame))

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

    r25166 r25183  
    7878        s += '{}\n'.format(fielddisplay(self, 'tk', 'tidal load Love number (deg 2)'))
    7979        s += '{}\n'.format(fielddisplay(self, 'tk2secular', 'secular fluid Love number'))
     80
     81        return s
    8082    #}}}
    8183

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

    r25147 r25183  
    33
    44%mesh ais: {{{
    5         md=model();
    6         md=triangle(md,'../Exp/Ais.exp',200000); % max element size
     5md=model();
     6md=triangle(md,'../Exp/Ais.exp',200000); % max element size
    77% }}}
    88%define load: {{{
    9         md.esa.deltathickness=zeros(md.mesh.numberofelements,1);
    10         disc_radius=500; % km
    11         index=md.mesh.elements;
    12         x_element=mean(md.mesh.x(index),2)-1.0e6;
    13         y_element=mean(md.mesh.y(index),2)-1.0e6;
    14         rad_dist=sqrt(x_element.^2+y_element.^2)/1000;  % radial distance in km
    15         md.esa.deltathickness(rad_dist<=disc_radius)=-1;   % 1 m water withdrawl
     9md.esa.deltathickness=zeros(md.mesh.numberofelements,1);
     10disc_radius=500; % km
     11index=md.mesh.elements;
     12x_element=mean(md.mesh.x(index),2)-1.0e6;
     13y_element=mean(md.mesh.y(index),2)-1.0e6;
     14rad_dist=sqrt(x_element.^2+y_element.^2)/1000;  % radial distance in km
     15pos=find(rad_dist<=disc_radius);
     16md.esa.deltathickness(pos)=-1;   % 1 m water withdrawl
    1617% }}}
    1718%read in love numbers:{{{
    … …  
    1920% }}}
    2021%mask:  {{{
    21         %make sure wherever there is an ice load, that the mask is set to ice:
    22         md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
    23         pos=find(md.esa.deltathickness); md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
     22%make sure wherever there is an ice load, that the mask is set to ice:
     23md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
     24pos=find(md.esa.deltathickness);
     25md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
    2426
    25         %is ice grounded?
    26         md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1);
    27         pos=find(md.mask.ice_levelset<=0); md.mask.ocean_levelset(pos)=1;
     27%is ice grounded?
     28md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1);
     29pos=find(md.mask.ice_levelset<=0);
     30md.mask.ocean_levelset(pos)=1;
    2831% }}}
    2932%geometry:  {{{
    30         di=md.materials.rho_ice/md.materials.rho_water;
    31         md.geometry.thickness=ones(md.mesh.numberofvertices,1);
    32         md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1);
    33         md.geometry.base=md.geometry.surface-md.geometry.thickness;
    34         md.geometry.bed=md.geometry.base;
     33di=md.materials.rho_ice/md.materials.rho_water;
     34md.geometry.thickness=ones(md.mesh.numberofvertices,1);
     35md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1);
     36md.geometry.base=md.geometry.surface-md.geometry.thickness;
     37md.geometry.bed=md.geometry.base;
    3538% }}}
    3639%materials:  {{{
    37         md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
    38         md.materials.rheology_B=paterson(md.initialization.temperature);
    39         md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     40md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
     41md.materials.rheology_B=paterson(md.initialization.temperature);
     42md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
    4043% }}}
    4144%additional parameters, miscellaneous: {{{
    42         md.miscellaneous.name='test2111';
    43         md.esa.degacc=0.01;
    44         md.esa.hemisphere = -1;
     45md.miscellaneous.name='test2111';
     46md.esa.degacc=0.01;
     47md.esa.hemisphere = -1;
    4548% }}}
    4649

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

    r25158 r25183  
    1717# }}}
    1818#define load: {{{
    19 md.esa.deltathickness = np.zeros((md.mesh.numberofelements, ))
     19md.esa.deltathickness = np.zeros((md.mesh.numberofelements, 1))
    2020disc_radius = 500  # km
    2121index = md.mesh.elements
    22 x_element = np.mean(md.mesh.x[index - 1], 1) - 1.0e6
    23 y_element = np.mean(md.mesh.y[index - 1], 1) - 1.0e6
     22x_element = md.mesh.x[index - 1].mean(axis=1) - 1.0e6
     23y_element = md.mesh.y[index - 1].mean(axis=1) - 1.0e6
    2424rad_dist = np.sqrt(x_element**2 + y_element**2) / 1000  # radial distance in km
    25 md.esa.deltathickness[np.where(rad_dist <= disc_radius)] = -1  # 1 m water withdrawl
     25pos = np.where(rad_dist <= disc_radius)[0]
     26md.esa.deltathickness[pos] = -1  # 1 m water withdrawl
    2627# }}}
    2728#love numbers: {{{
    … …  
    3031#mask:  {{{
    3132#make sure wherever there is an ice load, that the mask is set to ice:
    32 md.mask.ice_levelset = np.ones((md.mesh.numberofvertices, ))
    33 pos = np.where(md.esa.deltathickness)
     33md.mask.ice_levelset = np.ones((md.mesh.numberofvertices, 1))
     34pos = np.where(md.esa.deltathickness)[0]
    3435md.mask.ice_levelset[md.mesh.elements[pos, :]] = -1
    3536
    3637#is ice grounded?
    37 md.mask.ocean_levelset = -np.ones((md.mesh.numberofvertices, ))
    38 pos = np.where(md.mask.ice_levelset <= 0)
     38md.mask.ocean_levelset = -np.ones((md.mesh.numberofvertices, 1))
     39pos = np.where(md.mask.ice_levelset <= 0)[0]
    3940md.mask.ocean_levelset[pos] = 1
    4041# }}}
    4142#geometry:  {{{
    4243di = md.materials.rho_ice / md.materials.rho_water
    43 md.geometry.thickness = np.ones((md.mesh.numberofvertices, ))
    44 md.geometry.surface = (1 - di) * np.zeros((md.mesh.numberofvertices, ))
     44md.geometry.thickness = np.ones((md.mesh.numberofvertices, 1))
     45md.geometry.surface = (1 - di) * np.zeros((md.mesh.numberofvertices, 1))
    4546md.geometry.base = md.geometry.surface - md.geometry.thickness
    4647md.geometry.bed = md.geometry.base
    4748# }}}
    4849#materials:  {{{
    49 md.initialization.temperature = 273.25 * np.ones((md.mesh.numberofvertices, ))
     50md.initialization.temperature = 273.25 * np.ones((md.mesh.numberofvertices, 1))
    5051md.materials.rheology_B = paterson(md.initialization.temperature)
    51 md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements, ))
     52md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements, 1))
    5253# }}}
    5354#additional parameters, miscellaneous: {{{

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

    r25147 r25183  
    33
    44%mesh ais: {{{
    5         md=model();
    6         md=triangle(md,'../Exp/Ais.exp',100000); % max element size
     5md=model();
     6md=triangle(md,'../Exp/Ais.exp',100000); % max element size
    77% }}}
    88%define load: {{{
    9         md.esa.deltathickness=zeros(md.mesh.numberofelements,1);
    10         disc_radius=500; % km
    11         index=md.mesh.elements;
    12         x_element=mean(md.mesh.x(index),2);
    13         y_element=mean(md.mesh.y(index),2);
    14         rad_dist=sqrt(x_element.^2+y_element.^2)/1000;  % radial distance in km
    15         md.esa.deltathickness(rad_dist<=disc_radius)=-1;   % 1 m water withdrawl
     9md.esa.deltathickness=zeros(md.mesh.numberofelements,1);
     10disc_radius=500; % km
     11index=md.mesh.elements;
     12x_element=mean(md.mesh.x(index),2);
     13y_element=mean(md.mesh.y(index),2);
     14rad_dist=sqrt(x_element.^2+y_element.^2)/1000;  % radial distance in km
     15pos=find(rad_dist<=disc_radius);
     16md.esa.deltathickness(pos)=-1;   % 1 m water withdrawl
    1617% }}}
    1718%read in love numbers:{{{
    … …  
    1920% }}}
    2021%mask:  {{{
    21         %make sure wherever there is an ice load, that the mask is set to ice:
    22         md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
    23         pos=find(md.esa.deltathickness); md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
     22%make sure wherever there is an ice load, that the mask is set to ice:
     23md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
     24pos=find(md.esa.deltathickness); md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
    2425
    25         %is ice grounded?
    26         md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1);
    27         pos=find(md.mask.ice_levelset<=0); md.mask.ocean_levelset(pos)=1;
     26%is ice grounded?
     27md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1);
     28pos=find(md.mask.ice_levelset<=0); md.mask.ocean_levelset(pos)=1;
    2829% }}}
    2930%geometry:  {{{
    30         di=md.materials.rho_ice/md.materials.rho_water;
    31         md.geometry.thickness=ones(md.mesh.numberofvertices,1);
    32         md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1);
    33         md.geometry.base=md.geometry.surface-md.geometry.thickness;
    34         md.geometry.bed=md.geometry.base;
     31di=md.materials.rho_ice/md.materials.rho_water;
     32md.geometry.thickness=ones(md.mesh.numberofvertices,1);
     33md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1);
     34md.geometry.base=md.geometry.surface-md.geometry.thickness;
     35md.geometry.bed=md.geometry.base;
    3536% }}}
    3637%materials:  {{{
    37         md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
    38         md.materials.rheology_B=paterson(md.initialization.temperature);
    39         md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     38md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
     39md.materials.rheology_B=paterson(md.initialization.temperature);
     40md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
    4041% }}}
    4142%additional parameters, miscellaneous: {{{
    42         md.miscellaneous.name='test2112';
    43         md.esa.degacc=0.01;
    44         md.esa.hemisphere = -1;
     43md.miscellaneous.name='test2112';
     44md.esa.degacc=0.01;
     45md.esa.hemisphere = -1;
    4546% }}}
    4647

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

    r25158 r25183  
    2020disc_radius = 500  # km
    2121index = md.mesh.elements
    22 x_element = np.mean(md.mesh.x[index - 1], 1)
    23 y_element = np.mean(md.mesh.y[index - 1], 1)
     22x_element = md.mesh.x[index - 1].mean(axis=1)
     23y_element = md.mesh.y[index - 1].mean(axis=1)
    2424rad_dist = np.sqrt(x_element**2 + y_element**2) / 1000  # radial distance in km
    25 md.esa.deltathickness[np.where(rad_dist <= disc_radius)] = -1  # 1 m water withdrawl
     25pos = np.where(rad_dist <= disc_radius)[0]
     26md.esa.deltathickness[pos] = -1  # 1 m water withdrawl
    2627# }}}
    2728#love numbers: {{{
    … …  
    3031#mask:  {{{
    3132#make sure wherever there is an ice load, that the mask is set to ice:
    32 md.mask.ice_levelset = np.ones((md.mesh.numberofvertices, ))
    33 pos = np.where(md.esa.deltathickness)
     33md.mask.ice_levelset = np.ones((md.mesh.numberofvertices, 1))
     34pos = np.nonzero(md.esa.deltathickness)[0]
    3435md.mask.ice_levelset[md.mesh.elements[pos, :]] = -1
    3536
    3637#is ice grounded?
    37 md.mask.ocean_levelset = -np.ones((md.mesh.numberofvertices, ))
    38 pos = np.where(md.mask.ice_levelset <= 0)
     38md.mask.ocean_levelset = -np.ones((md.mesh.numberofvertices, 1))
     39pos = np.where(md.mask.ice_levelset <= 0)[0]
    3940md.mask.ocean_levelset[pos] = 1
    4041# }}}
    4142#geometry:  {{{
    4243di = md.materials.rho_ice / md.materials.rho_water
    43 md.geometry.thickness = np.ones((md.mesh.numberofvertices, ))
    44 md.geometry.surface = (1 - di) * np.zeros((md.mesh.numberofvertices, ))
     44md.geometry.thickness = np.ones((md.mesh.numberofvertices, 1))
     45md.geometry.surface = (1 - di) * np.zeros((md.mesh.numberofvertices, 1))
    4546md.geometry.base = md.geometry.surface - md.geometry.thickness
    4647md.geometry.bed = md.geometry.base
    4748# }}}
    4849#materials:  {{{
    49 md.initialization.temperature = 273.25 * np.ones((md.mesh.numberofvertices, ))
     50md.initialization.temperature = 273.25 * np.ones((md.mesh.numberofvertices, 1))
    5051md.materials.rheology_B = paterson(md.initialization.temperature)
    51 md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements, ))
     52md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements, 1))
    5253# }}}
    5354#additional parameters, miscellaneous: {{{

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

    r25147 r25183  
    44
    55%mesh ais: {{{
    6         md=model();
    7         md=triangle(md,'../Exp/Ais.exp',200000); % max element size
     6md=model();
     7md=triangle(md,'../Exp/Ais.exp',200000); % max element size
    88% }}}
    99%define load: {{{
    10         md.esa.deltathickness=zeros(md.mesh.numberofelements,1);
    11         disc_radius=500; % km
    12         index=md.mesh.elements;
    13         x_element=mean(md.mesh.x(index),2)-1.0e6;
    14         y_element=mean(md.mesh.y(index),2)+1.0e6;
    15         rad_dist=sqrt(x_element.^2+y_element.^2)/1000;  % radial distance in km
    16         md.esa.deltathickness(rad_dist<=disc_radius)=-1;   % 1 m water withdrawl
     10md.esa.deltathickness=zeros(md.mesh.numberofelements,1);
     11disc_radius=500; % km
     12index=md.mesh.elements;
     13x_element=mean(md.mesh.x(index),2)-1.0e6;
     14y_element=mean(md.mesh.y(index),2)+1.0e6;
     15rad_dist=sqrt(x_element.^2+y_element.^2)/1000;  % radial distance in km
     16pos=find(rad_dist<=disc_radius);
     17md.esa.deltathickness(pos)=-1;   % 1 m water withdrawl
    1718% }}}
    1819%read in love numbers:{{{
    … …  
    2021% }}}
    2122%mask:  {{{
    22         %make sure wherever there is an ice load, that the mask is set to ice:
    23         md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
    24         pos=find(md.esa.deltathickness); md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
     23%make sure wherever there is an ice load, that the mask is set to ice:
     24md.mask.ice_levelset=ones(md.mesh.numberofvertices,1);
     25pos=find(md.esa.deltathickness); md.mask.ice_levelset(md.mesh.elements(pos,:))=-1;
    2526
    26         %is ice grounded?
    27         md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1);
    28         pos=find(md.mask.ice_levelset<=0); md.mask.ocean_levelset(pos)=1;
     27%is ice grounded?
     28md.mask.ocean_levelset=-ones(md.mesh.numberofvertices,1);
     29pos=find(md.mask.ice_levelset<=0); md.mask.ocean_levelset(pos)=1;
    2930% }}}
    3031%geometry:  {{{
    31         di=md.materials.rho_ice/md.materials.rho_water;
    32         md.geometry.thickness=ones(md.mesh.numberofvertices,1);
    33         md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1);
    34         md.geometry.base=md.geometry.surface-md.geometry.thickness;
    35         md.geometry.bed=md.geometry.base;
     32di=md.materials.rho_ice/md.materials.rho_water;
     33md.geometry.thickness=ones(md.mesh.numberofvertices,1);
     34md.geometry.surface=(1-di)*zeros(md.mesh.numberofvertices,1);
     35md.geometry.base=md.geometry.surface-md.geometry.thickness;
     36md.geometry.bed=md.geometry.base;
    3637% }}}
    3738%materials:  {{{
    38         md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
    39         md.materials.rheology_B=paterson(md.initialization.temperature);
    40         md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     39md.initialization.temperature=273.25*ones(md.mesh.numberofvertices,1);
     40md.materials.rheology_B=paterson(md.initialization.temperature);
     41md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
    4142% }}}
    4243%additional parameters, miscellaneous: {{{
    43         md.miscellaneous.name='test2113';
    44         md.esa.degacc=0.01;
    45         md.esa.hemisphere = 1; % AIS is placed in Northern Hemisphere
     44md.miscellaneous.name='test2113';
     45md.esa.degacc=0.01;
     46md.esa.hemisphere = 1; % AIS is placed in Northern Hemisphere
    4647% }}}
    4748

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

    r25158 r25183  
    1818# }}}
    1919#define load: {{{
    20 md.esa.deltathickness = np.zeros((md.mesh.numberofelements, ))
     20md.esa.deltathickness = np.zeros((md.mesh.numberofelements, 1))
    2121disc_radius = 500  # km
    2222index = md.mesh.elements
    … …  
    2424y_element = np.mean(md.mesh.y[index - 1], 1) + 1.0e6
    2525rad_dist = np.sqrt(x_element**2 + y_element**2) / 1000  # radial distance in km
    26 md.esa.deltathickness[np.where(rad_dist <= disc_radius)] = -1  # 1 m water withdrawl
     26pos = np.where(rad_dist <= disc_radius)[0]
     27md.esa.deltathickness[pos] = -1  # 1 m water withdrawl
    2728# }}}
    2829#love numbers: {{{
    … …  
    3132#mask:  {{{
    3233#make sure wherever there is an ice load, that the mask is set to ice:
    33 md.mask.ice_levelset = np.ones((md.mesh.numberofvertices, ))
    34 pos = np.where(md.esa.deltathickness)
     34md.mask.ice_levelset = np.ones((md.mesh.numberofvertices, 1))
     35pos = np.nonzero(md.esa.deltathickness)[0]
    3536md.mask.ice_levelset[md.mesh.elements[pos, :]] = -1
    3637
    3738#is ice grounded?
    38 md.mask.ocean_levelset = -np.ones((md.mesh.numberofvertices, ))
    39 pos = np.where(md.mask.ice_levelset <= 0)
     39md.mask.ocean_levelset = -np.ones((md.mesh.numberofvertices, 1))
     40pos = np.where(md.mask.ice_levelset <= 0)[0]
    4041md.mask.ocean_levelset[pos] = 1
    4142# }}}
    4243#geometry:  {{{
    4344di = md.materials.rho_ice / md.materials.rho_water
    44 md.geometry.thickness = np.ones((md.mesh.numberofvertices, ))
    45 md.geometry.surface = (1 - di) * np.zeros((md.mesh.numberofvertices, ))
     45md.geometry.thickness = np.ones((md.mesh.numberofvertices, 1))
     46md.geometry.surface = (1 - di) * np.zeros((md.mesh.numberofvertices, 1))
    4647md.geometry.base = md.geometry.surface - md.geometry.thickness
    4748md.geometry.bed = md.geometry.base
    4849# }}}
    4950#materials:  {{{
    50 md.initialization.temperature = 273.25 * np.ones((md.mesh.numberofvertices, ))
     51md.initialization.temperature = 273.25 * np.ones((md.mesh.numberofvertices, 1))
    5152md.materials.rheology_B = paterson(md.initialization.temperature)
    52 md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements, ))
     53md.materials.rheology_n = 3 * np.ones((md.mesh.numberofelements, 1))
    5354# }}}
    5455#additional parameters, miscellaneous: {{{
Note: See TracChangeset for help on using the changeset viewer.