| [18296] | 1 | Index: ../trunk-jpl/src/m/mech/backstressfrominversion.m
|
|---|
| 2 | ===================================================================
|
|---|
| 3 | --- ../trunk-jpl/src/m/mech/backstressfrominversion.m (revision 18010)
|
|---|
| 4 | +++ ../trunk-jpl/src/m/mech/backstressfrominversion.m (revision 18011)
|
|---|
| 5 | @@ -22,9 +22,8 @@
|
|---|
| 6 | % 'xy': x and y axes same as in polar stereographic projection
|
|---|
| 7 | %
|
|---|
| 8 | % Return values:
|
|---|
| 9 | -% 'backstress' is the inferred backstress necessary to balance the
|
|---|
| 10 | -% analytical solution (keeping damage within its appropriate limits, e.g. D
|
|---|
| 11 | -% in [0,1]).
|
|---|
| 12 | +% 'backstress' is the inferred backstress based on the analytical
|
|---|
| 13 | +% solution for ice shelf creep
|
|---|
| 14 | %
|
|---|
| 15 | % Usage:
|
|---|
| 16 | % backstress=backstressfrominversion(md,options)
|
|---|
| 17 | Index: ../trunk-jpl/src/m/mech/calcbackstress.m
|
|---|
| 18 | ===================================================================
|
|---|
| 19 | --- ../trunk-jpl/src/m/mech/calcbackstress.m (revision 0)
|
|---|
| 20 | +++ ../trunk-jpl/src/m/mech/calcbackstress.m (revision 18011)
|
|---|
| 21 | @@ -0,0 +1,64 @@
|
|---|
| 22 | +function backstress=calcbackstress(md,varargin)
|
|---|
| 23 | +%BACKSTRESSFROMINVERSION - compute ice shelf backstress
|
|---|
| 24 | +%
|
|---|
| 25 | +% This routine computes backstress based on the analytical formalism of
|
|---|
| 26 | +% Thomas (1973) and Borstad et al. (2013) based on the ice rigidity,
|
|---|
| 27 | +% thickness, the densities of ice and seawater, and (optionally)
|
|---|
| 28 | +% damage. Strain rates must also be included, either from observed or
|
|---|
| 29 | +% modeled velocities.
|
|---|
| 30 | +
|
|---|
| 31 | +% Available options:
|
|---|
| 32 | +% - 'smoothing' : the amount of smoothing to be applied to the strain rate data.
|
|---|
| 33 | +% Type 'help averaging' for more information on its
|
|---|
| 34 | +% usage. Defaults to 0.
|
|---|
| 35 | +% - 'coordsys' : coordinate system for calculating the strain rate
|
|---|
| 36 | +% components. Must be one of:
|
|---|
| 37 | +% 'longitudinal': x axis aligned along a flowline at every point (default)
|
|---|
| 38 | +% 'principal': x axis aligned along maximum principal strain rate
|
|---|
| 39 | +% at every point
|
|---|
| 40 | +% 'xy': x and y axes same as in polar stereographic projection
|
|---|
| 41 | +%
|
|---|
| 42 | +% Return values:
|
|---|
| 43 | +% 'backstress' is the inferred backstress based on the analytical
|
|---|
| 44 | +% solution for ice shelf creep
|
|---|
| 45 | +%
|
|---|
| 46 | +% Usage:
|
|---|
| 47 | +% backstress=calcbackstress(md,options)
|
|---|
| 48 | +%
|
|---|
| 49 | +% Example:
|
|---|
| 50 | +% backstress=backstressfrominversion(md,'smoothing',2,'coordsys','longitudinal');
|
|---|
| 51 | +
|
|---|
| 52 | +% check inputs
|
|---|
| 53 | +if (nargin<1),
|
|---|
| 54 | + help backstressfrominversion
|
|---|
| 55 | + error('bad usage');
|
|---|
| 56 | +end
|
|---|
| 57 | +if isempty(fieldnames(md.results)),
|
|---|
| 58 | + error(['md.results.strainrate is not present. Calculate using md=mechanicalproperties(md,vx,vy)']);
|
|---|
| 59 | +end
|
|---|
| 60 | +if dimension(md.mesh)~=2,
|
|---|
| 61 | + error('only 2d model supported currently');
|
|---|
| 62 | +end
|
|---|
| 63 | +if any(md.flowequation.element_equation~=2),
|
|---|
| 64 | + disp('Warning: the model has some non SSA elements. These will be treated like SSA elements');
|
|---|
| 65 | +end
|
|---|
| 66 | +
|
|---|
| 67 | +% process options
|
|---|
| 68 | +options = pairoptions(varargin{:});
|
|---|
| 69 | +smoothing = getfieldvalue(options,'smoothing',0);
|
|---|
| 70 | +coordsys = getfieldvalue(options,'coordsys','longitudinal');
|
|---|
| 71 | +tempmask = getfieldvalue(options,'tempmask',false);
|
|---|
| 72 | +
|
|---|
| 73 | +T=0.5*md.materials.rho_ice*md.constants.g*(1-md.materials.rho_ice/md.materials.rho_water)*md.geometry.thickness;
|
|---|
| 74 | +n=averaging(md,md.materials.rheology_n,0);
|
|---|
| 75 | +B=md.materials.rheology_B;
|
|---|
| 76 | +if md.damage.isdamage,
|
|---|
| 77 | + D=md.damage.D
|
|---|
| 78 | +else
|
|---|
| 79 | + D=0.
|
|---|
| 80 | +
|
|---|
| 81 | +[a0,b0,theta0,ex0]=thomasparams(md,'eq','Thomas','smoothing',smoothing,'coordsys',coordsys);
|
|---|
| 82 | +
|
|---|
| 83 | +% analytical backstress solution
|
|---|
| 84 | +backstress=T-(1.-D).*B.*sign(ex0).*(2+a0).*abs(ex0).^(1./n)./((1+a0+a0.^2+b0.^2).^((n-1)/2./n));
|
|---|
| 85 | +backstress(find(backstress<0))=0;
|
|---|
| 86 | Index: ../trunk-jpl/src/m/mech/backstressfrominversion.py
|
|---|
| 87 | ===================================================================
|
|---|
| 88 | --- ../trunk-jpl/src/m/mech/backstressfrominversion.py (revision 18010)
|
|---|
| 89 | +++ ../trunk-jpl/src/m/mech/backstressfrominversion.py (revision 18011)
|
|---|
| 90 | @@ -7,7 +7,7 @@
|
|---|
| 91 | Compute ice shelf backstress from inversion results.
|
|---|
| 92 |
|
|---|
| 93 | This routine computes backstress based on the analytical formalism of
|
|---|
| 94 | - Thomas (1973) and Borstad et al. (2013, The Cryoshpere). The model
|
|---|
| 95 | + Thomas (1973) and Borstad et al. (2013, The Cryosphere). The model
|
|---|
| 96 | must contain inversion results for ice rigidity. Strain rates must
|
|---|
| 97 | also be included, either from observed or modeled velocities. Ice
|
|---|
| 98 | rigidity B is assumed to be parameterized by the ice temperature in
|
|---|
| 99 | @@ -28,9 +28,8 @@
|
|---|
| 100 | 'xy': x and y axes same as in polar stereographic projection
|
|---|
| 101 |
|
|---|
| 102 | Return values:
|
|---|
| 103 | - 'backstress' is the inferred backstress necessary to balance the
|
|---|
| 104 | - analytical solution (keeping damage within its appropriate limits, e.g. D
|
|---|
| 105 | - in [0,1]).
|
|---|
| 106 | + 'backstress' is the inferred backstress based on the analytical
|
|---|
| 107 | + solution for ice shelf creep
|
|---|
| 108 |
|
|---|
| 109 | Usage:
|
|---|
| 110 | backstress=backstressfrominversion(md,options)
|
|---|
| 111 | Index: ../trunk-jpl/src/m/mech/calcbackstress.py
|
|---|
| 112 | ===================================================================
|
|---|
| 113 | --- ../trunk-jpl/src/m/mech/calcbackstress.py (revision 0)
|
|---|
| 114 | +++ ../trunk-jpl/src/m/mech/calcbackstress.py (revision 18011)
|
|---|
| 115 | @@ -0,0 +1,66 @@
|
|---|
| 116 | +import numpy as npy
|
|---|
| 117 | +from averaging import averaging
|
|---|
| 118 | +from thomasparams import thomasparams
|
|---|
| 119 | +
|
|---|
| 120 | +def calcbackstress(md,**kwargs):
|
|---|
| 121 | + '''
|
|---|
| 122 | + Compute ice shelf backstress.
|
|---|
| 123 | +
|
|---|
| 124 | + This routine computes backstress based on the analytical formalism of
|
|---|
| 125 | + Thomas (1973) and Borstad et al. (2013, The Cryosphere) based on the
|
|---|
| 126 | + ice rigidity, thickness, the densities of ice and seawater, and
|
|---|
| 127 | + (optionally) damage. Strain rates must also be included, either from
|
|---|
| 128 | + observed or modeled velocities.
|
|---|
| 129 | +
|
|---|
| 130 | + Available options:
|
|---|
| 131 | + - 'smoothing' : the amount of smoothing to be applied to the strain rate data.
|
|---|
| 132 | + Type 'help averaging' for more information on its
|
|---|
| 133 | + usage. Defaults to 0.
|
|---|
| 134 | + - 'coordsys' : coordinate system for calculating the strain rate
|
|---|
| 135 | + components. Must be one of:
|
|---|
| 136 | + 'longitudinal': x axis aligned along a flowline at every point (default)
|
|---|
| 137 | + 'principal': x axis aligned along maximum principal strain rate
|
|---|
| 138 | + at every point
|
|---|
| 139 | + 'xy': x and y axes same as in polar stereographic projection
|
|---|
| 140 | +
|
|---|
| 141 | + Return values:
|
|---|
| 142 | + 'backstress' is the inferred backstress based on the analytical
|
|---|
| 143 | + solution for ice shelf creep
|
|---|
| 144 | +
|
|---|
| 145 | + Usage:
|
|---|
| 146 | + backstress=calcbackstress(md,options)
|
|---|
| 147 | +
|
|---|
| 148 | + Example:
|
|---|
| 149 | + backstress=calcbackstress(md,'smoothing',2,'coordsys','longitudinal')
|
|---|
| 150 | + '''
|
|---|
| 151 | +
|
|---|
| 152 | + # unpack kwargs
|
|---|
| 153 | + smoothing=kwargs.pop('smoothing',0)
|
|---|
| 154 | + if 'smoothing' in kwargs: del kwargs['smoothing']
|
|---|
| 155 | + coordsys=kwargs.pop('coordsys','longitudinal')
|
|---|
| 156 | + if 'coordsys' in kwargs: del kwargs['coordsys']
|
|---|
| 157 | + assert len(kwargs)==0, 'error, unexpected or misspelled kwargs'
|
|---|
| 158 | +
|
|---|
| 159 | + # some checks
|
|---|
| 160 | + if not hasattr(md.results,'strainrate'):
|
|---|
| 161 | + raise StandardError('md.results.strainrate not present. Calculate using md=mechanicalproperties(md,vx,vy)')
|
|---|
| 162 | + if not '2d' in md.mesh.__doc__:
|
|---|
| 163 | + raise StandardError('only 2d (planview) model supported currently')
|
|---|
| 164 | + if any(md.flowequation.element_equation!=2):
|
|---|
| 165 | + raise StandardError('Warning: the model has some non-SSA elements. These will be treated like SSA elements')
|
|---|
| 166 | +
|
|---|
| 167 | + T=0.5*md.materials.rho_ice*md.constants.g*(1-md.materials.rho_ice/md.materials.rho_water)*md.geometry.thickness
|
|---|
| 168 | + n=averaging(md,md.materials.rheology_n,0)
|
|---|
| 169 | + B=md.materials.rheology_B
|
|---|
| 170 | + if md.damage.isdamage:
|
|---|
| 171 | + D=md.damage.D
|
|---|
| 172 | + else:
|
|---|
| 173 | + D=0.
|
|---|
| 174 | +
|
|---|
| 175 | + a0,b0,theta0,ex0=thomasparams(md,eq='Thomas',smoothing=smoothing,coordsys=coordsys)
|
|---|
| 176 | +
|
|---|
| 177 | + # analytical backstress solution
|
|---|
| 178 | + backstress=T-(1.-D)*B*npy.sign(ex0)*(2+a0)*npy.abs(ex0)**(1./n)/((1+a0+a0**2+b0**2)**((n-1.)/2./n))
|
|---|
| 179 | + backstress[npy.nonzero(backstress<0)]=0
|
|---|
| 180 | +
|
|---|
| 181 | + return backstress
|
|---|
![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
