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

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

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