source: issm/trunk/src/m/classes/fourierlove.m

Last change on this file was 27035, checked in by Mathieu Morlighem, 3 years ago

merged trunk-jpl and trunk for revision 27033
File size: 10.5 KB
Line 
1%FOURIERLOVE class definition
2%
3% Usage:
4% md.love=fourierlove();
5
6classdef fourierlove
7 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;
29
30 end
31 methods (Static)
32 function self = loadobj(self) % {{{
33 % This function is directly called by matlab when a model object is
34 % loaded. Update old properties here
35 end% }}}
36 end
37 methods
38 function self = fourierlove(varargin) % {{{
39 switch nargin
40 case 0
41 self=setdefaultparameters(self);
42 otherwise
43 error('constructor not supported');
44 end
45 end % }}}
46 function self = setdefaultparameters(self) % {{{
47 %we setup an elastic love number computation by default.
48 self.nfreq=1;
49 self.frequencies=[0]; %Hz
50 self.sh_nmax=256; % .35 degree, 40 km at the equator.
51 self.sh_nmin=1;
52 % work on matlab script for computing g0 for given Earth's structure.
53 self.g0=9.81; % m/s^2;
54 self.r0=6371*1e3; %m;
55 self.mu0=1e11; % Pa
56 self.Gravitational_Constant=6.67259e-11; % m^3 kg^-1 s^-2
57 self.chandler_wobble=0;
58 self.allow_layer_deletion=1;
59 self.underflow_tol=1e-16; %threshold of deep to surface love number ratio to trigger the deletion of layer
60 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
61 self.integration_steps_per_layer=100;
62 self.istemporal=0;
63 self.n_temporal_iterations=8;
64 self.time=[0]; %s
65 self.love_kernels=0;
66 self.forcing_type = 11; % surface loading
67 self.inner_core_boundary=1;
68 self.core_mantle_boundary=2;
69 self.complex_computation=0;
70 end % }}}
71 function disp(self) % {{{
72 fielddisplay(self,'nfreq','number of frequencies sampled (default: 1, elastic) [Hz]');
73 fielddisplay(self,'frequencies','frequencies sampled (convention defaults to 0 for the elastic case) [Hz]');
74 fielddisplay(self,'sh_nmax','maximum spherical harmonic degree (default: 256, .35 deg, or 40 km at equator)');
75 fielddisplay(self,'sh_nmin','minimum spherical harmonic degree (default: 1)');
76 fielddisplay(self,'g0','adimensioning constant for gravity (default: 10) [m/s^2]');
77 fielddisplay(self,'r0','adimensioning constant for radius (default: 6371*10^3) [m]');
78 fielddisplay(self,'mu0','adimensioning constant for stress (default: 10^11) [Pa]');
79 fielddisplay(self,'Gravitational_Constant','Newtonian constant of gravitation (default: 6.67259e-11 [m^3 kg^-1 s^-2])');
80 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 yet)');
81 fielddisplay(self,'allow_layer_deletion','allow for migration of the integration boundary with increasing spherical harmonics degree (default: 1)');
82 fielddisplay(self,'underflow_tol','threshold of deep to surface love number ratio to trigger the deletion of layers (default: 1e-16)');
83 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)');
84 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)');
85 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'});
86 fielddisplay(self,'n_temporal_iterations','max number of iterations in the inverse Laplace transform. Also the number of spectral samples per time step requested (default: 8)');
87 fielddisplay(self,'time','time vector for deformation if istemporal (default: 0) [s]');
88 fielddisplay(self,'love_kernels','compute love numbers at depth? (default: 0)');
89 fielddisplay(self,'forcing_type',{'integer indicating the nature and depth of the forcing for the Love number calculation (default: 11):','1: Inner core boundary -- Volumic Potential','2: Inner core boundary -- Pressure','3: Inner core boundary -- Loading','4: Inner core boundary -- Tangential traction','5: Core mantle boundary -- Volumic Potential','6: Core mantle boundary -- Pressure','7: Core mantle boundary -- Loading','8: Core mantle boundary -- Tangential traction','9: Surface -- Volumic Potential','10: Surface -- Pressure','11: Surface -- Loading','12: Surface -- Tangential traction '});
90 fielddisplay(self,'inner_core_boundary','interface index in materials.radius locating forcing. Only used for forcing_type 1--4 (default: 1)');
91 fielddisplay(self,'core_mantle_boundary','interface index in materials.radius locating forcing. Only used for forcing_type 5--8 (default: 2)');
92
93 end % }}}
94 function md = checkconsistency(self,md,solution,analyses) % {{{
95
96 if ~ismember('LoveAnalysis',analyses), return; end
97
98 md = checkfield(md,'fieldname','love.nfreq','NaN',1,'Inf',1,'numel',1,'>',0);
99 md = checkfield(md,'fieldname','love.frequencies','NaN',1,'Inf',1,'numel',md.love.nfreq);
100 md = checkfield(md,'fieldname','love.sh_nmax','NaN',1,'Inf',1,'numel',1,'>',0);
101 md = checkfield(md,'fieldname','love.sh_nmin','NaN',1,'Inf',1,'numel',1,'>',0);
102 md = checkfield(md,'fieldname','love.g0','NaN',1,'Inf',1,'numel',1,'>',0);
103 md = checkfield(md,'fieldname','love.r0','NaN',1,'Inf',1,'numel',1,'>',0);
104 md = checkfield(md,'fieldname','love.mu0','NaN',1,'Inf',1,'numel',1,'>',0);
105 md = checkfield(md,'fieldname','love.Gravitational_Constant','NaN',1,'Inf',1,'numel',1,'>',0);
106 md = checkfield(md,'fieldname','love.chandler_wobble','values',[0 1]);
107 md = checkfield(md,'fieldname','love.allow_layer_deletion','values',[0 1]);
108 md = checkfield(md,'fieldname','love.underflow_tol','NaN',1,'Inf',1,'numel',1,'>',0);
109 md = checkfield(md,'fieldname','love.pw_threshold','NaN',1,'Inf',1,'numel',1,'>',0);
110 md = checkfield(md,'fieldname','love.integration_steps_per_layer','NaN',1,'Inf',1,'numel',1,'>',0);
111 md = checkfield(md,'fieldname','love.love_kernels','values',[0 1]);
112 md = checkfield(md,'fieldname','love.forcing_type','NaN',1,'Inf',1,'numel',1,'>',0, '<=', 12);
113 md = checkfield(md,'fieldname','love.complex_computation','NaN',1,'Inf',1,'numel',1,'values',[0 1]);
114
115 md = checkfield(md,'fieldname','love.istemporal','values',[0 1]);
116 if md.love.istemporal==1
117 md = checkfield(md,'fieldname','love.n_temporal_iterations','NaN',1,'Inf',1,'numel',1,'>',0);
118 md = checkfield(md,'fieldname','love.time','NaN',1,'Inf',1,'numel',md.love.nfreq/2/md.love.n_temporal_iterations);
119 end
120 if md.love.sh_nmin<=1 & (md.love.forcing_type==1 || md.love.forcing_type==5 || md.love.forcing_type==9)
121 error(['Degree 1 not supported for forcing type ' num2str(md.love.forcing_type) '. Use sh_min>=2 for this kind of calculation.'])
122 end
123
124 if md.love.chandler_wobble==1
125 disp('Warning, Chandler Wobble in Love number calculator has not been validated yet');
126 end
127
128 %need 'litho' material:
129 if ~isa(md.materials,'materials') | ~sum(strcmpi(md.materials.nature,'litho'))
130 error('Need a ''litho'' material to run a Fourier Love number analysis');
131 end
132
133 mat=find(strcmpi(md.materials.nature,'litho'));
134 if md.love.forcing_type<=4
135 md = checkfield(md,'fieldname','love.inner_core_boundary','NaN',1,'Inf',1,'numel',1,'>',0, '<=', md.materials(mat).numlayers);
136 elseif md.love.forcing_type<=8
137 md = checkfield(md,'fieldname','love.core_mantle_boundary','NaN',1,'Inf',1,'numel',1,'>',0, '<=', md.materials(mat).numlayers);
138 end
139
140 end % }}}
141 function marshall(self,prefix,md,fid) % {{{
142
143 WriteData(fid,prefix,'object',self,'fieldname','nfreq','format','Integer');
144 WriteData(fid,prefix,'object',self,'fieldname','frequencies','format','DoubleMat','mattype',3);
145 WriteData(fid,prefix,'object',self,'fieldname','sh_nmax','format','Integer');
146 WriteData(fid,prefix,'object',self,'fieldname','sh_nmin','format','Integer');
147 WriteData(fid,prefix,'object',self,'fieldname','g0','format','Double');
148 WriteData(fid,prefix,'object',self,'fieldname','r0','format','Double');
149 WriteData(fid,prefix,'object',self,'fieldname','mu0','format','Double');
150 WriteData(fid,prefix,'object',self,'fieldname','Gravitational_Constant','format','Double');
151 WriteData(fid,prefix,'object',self,'fieldname','chandler_wobble','format','Boolean');
152 WriteData(fid,prefix,'object',self,'fieldname','allow_layer_deletion','format','Boolean');
153 WriteData(fid,prefix,'object',self,'fieldname','underflow_tol','format','Double');
154 WriteData(fid,prefix,'object',self,'fieldname','pw_threshold','format','Double');
155 WriteData(fid,prefix,'object',self,'fieldname','integration_steps_per_layer','format','Integer');
156 WriteData(fid,prefix,'object',self,'fieldname','istemporal','format','Boolean');
157 WriteData(fid,prefix,'object',self,'fieldname','n_temporal_iterations','format','Integer');
158 WriteData(fid,prefix,'object',self,'fieldname','complex_computation','format','Boolean');
159 %note: no need to marshall the time vector, we have frequencies
160 WriteData(fid,prefix,'object',self,'fieldname','love_kernels','format','Boolean');
161 WriteData(fid,prefix,'object',self,'fieldname','forcing_type','format','Integer');
162 WriteData(fid,prefix,'object',self,'fieldname','inner_core_boundary','format','Integer');
163 WriteData(fid,prefix,'object',self,'fieldname','core_mantle_boundary','format','Integer');
164
165 end % }}}
166 function self = extrude(self,md) % {{{
167 end % }}}
168 function savemodeljs(self,fid,modelname) % {{{
169 error('not implemented yet!');
170 end % }}}
171 function self=build_frequencies_from_time(self) % {{{
172 if ~self.istemporal
173 error('cannot build frequencies for temporal love numbers if love.istemporal==0')
174 end
175 disp('Temporal love numbers: Overriding md.love.nfreq and md.love.frequencies');
176 self.nfreq=length(self.time)*2*self.n_temporal_iterations;
177 self.frequencies=zeros(self.nfreq,1);
178 for i=1:length(self.time)
179 for j=1:2*self.n_temporal_iterations
180 if self.time(i)==0
181 self.frequencies((i-1)*2*self.n_temporal_iterations +j) =0; % convention to avoid marshalling infinite numbers
182 else
183 self.frequencies((i-1)*2*self.n_temporal_iterations +j) = j*log(2)/self.time(i)/2/pi;
184 end
185 end
186 end
187 end % }}}
188 end
189end
Note: See TracBrowser for help on using the repository browser.