Changeset 22305


Ignore:
Timestamp:
12/22/17 08:45:11 (7 years ago)
Author:
Mathieu Morlighem
Message:

BUG: revert to previous version

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r22304 r22305  
    1616                                self.nature=varargin;
    1717                        end
    18 
     18                       
     19                        %check this is acceptable:
     20                        for i=1:length(self.nature),
     21                                if ~(strcmpi(self.nature{i},'litho') | strcmpi(self.nature{i},'ice')),
     22                                        error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')');
     23                                end
     24                        end
     25                       
    1926                        %start filling in the dynamic fields:
    2027                        for i=1:length(self.nature),
    21             list = listoffields(self,self.nature{i});
    22             for j=1:size(list,1),
    23                                         %Add property to materials
    24                self.addprop(list{j,1});
    25                                         %Set default value according to table
    26                                         self.(list{j,1}) = list{j,2};
    27                                 end
    28                         end
    29 
     28                                nat=self.nature{i};
     29                                switch nat
     30                                case 'ice'
     31                                        self.addprop('rho_ice');
     32                                        self.addprop('rho_water');
     33                                        self.addprop('rho_freshwater');
     34                                        self.addprop('mu_water');
     35                                        self.addprop('heatcapacity');
     36                                        self.addprop('latentheat');
     37                                        self.addprop('thermalconductivity');
     38                                        self.addprop('temperateiceconductivity');
     39                                    self.addprop('meltingpoint');
     40                                        self.addprop('beta');
     41                                        self.addprop('mixed_layer_capacity');
     42                                        self.addprop('thermal_exchange_velocity');
     43                                        self.addprop('rheology_B');
     44                                        self.addprop('rheology_n');
     45                                        self.addprop('rheology_law');
     46                                case 'litho'
     47                                        self.addprop('numlayers');
     48                                    self.addprop('radius');
     49                                        self.addprop('viscosity');
     50                                        self.addprop('lame_lambda');
     51                                        self.addprop('lame_mu');
     52                                        self.addprop('burgers_viscosity');
     53                                        self.addprop('burgers_mu');
     54                                        self.addprop('isburgers');
     55                                        self.addprop('density');
     56                                        self.addprop('issolid');
     57                                otherwise
     58                                        error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')');
     59                                end
     60                        end
     61                        %set default parameters:
     62                        self.setdefaultparameters();
     63
     64                end % }}}
     65                function self = setdefaultparameters(self) % {{{
     66
     67                        for i=1:length(self.nature),
     68                                nat=self.nature{i};
     69                                switch nat
     70                                case 'ice'
     71                                        %ice density (kg/m^3)
     72                                        self.rho_ice=917.;
     73
     74                                        %ocean water density (kg/m^3)
     75                                        self.rho_water=1023.;
     76
     77                                        %fresh water density (kg/m^3)
     78                                        self.rho_freshwater=1000.;
     79
     80                                        %water viscosity (N.s/m^2)
     81                                        self.mu_water=0.001787; 
     82
     83                                        %ice heat capacity cp (J/kg/K)
     84                                        self.heatcapacity=2093.;
     85
     86                                        %ice latent heat of fusion L (J/kg)
     87                                        self.latentheat=3.34*10^5;
     88
     89                                        %ice thermal conductivity (W/m/K)
     90                                        self.thermalconductivity=2.4;
     91                                       
     92                                        %wet ice thermal conductivity (W/m/K)
     93                                        self.temperateiceconductivity=.24;
     94
     95                                        %the melting point of ice at 1 atmosphere of pressure in K
     96                                        self.meltingpoint=273.15;
     97
     98                                        %rate of change of melting point with pressure (K/Pa)
     99                                        self.beta=9.8*10^-8;
     100
     101                                        %mixed layer (ice-water interface) heat capacity (J/kg/K)
     102                                        self.mixed_layer_capacity=3974.;
     103
     104                                        %thermal exchange velocity (ice-water interface) (m/s)
     105                                        self.thermal_exchange_velocity=1.00*10^-4;
     106
     107                                        %Rheology law: what is the temperature dependence of B with T
     108                                        %available: none, paterson and arrhenius
     109                                        self.rheology_law='Paterson';
     110
     111                                case 'litho'
     112                                        %we default to a configuration that enables running GIA solutions using giacaron and/or giaivins.
     113                                        self.numlayers=2;
     114
     115                                        %center of the earth (approximation, must not be 0), then the lab (lithosphere/asthenosphere boundary) then the surface
     116                                        %(with 1d3 to avoid numerical singularities)
     117                                        self.radius=[1e3;6278*1e3;6378*1e3];
     118
     119                                        self.viscosity=[1e21;1e40]; %mantle and lithosphere viscosity (respectively) [Pa.s]
     120                                        self.lame_mu=[1.45*1e11;6.7*10^10];  % (Pa) %lithosphere and mantle shear modulus (respectively) [Pa]
     121                                        self.lame_lambda=self.lame_mu;  % (Pa) %mantle and lithosphere lamba parameter (respectively) [Pa]
     122                                        self.burgers_viscosity=[NaN;NaN];
     123                                        self.burgers_mu=[NaN;NaN];
     124                                        self.isburgers=[0;0];
     125                                        self.density=[5.51*1e3;5.50*1e3];  % (Pa) %mantle and lithosphere density [kg/m^3]
     126                                        self.issolid=[1;1]; % is layer solid or liquid.
     127
     128                                otherwise
     129                                        error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')');
     130                                end
     131                        end
    30132                end % }}}
    31133                function disp(self) % {{{
    … …  
    33135
    34136                        for i=1:length(self.nature),
    35                                 list = listoffields(self,self.nature{i});
    36                                 disp(sprintf('\n      %s:',upper(self.nature{i})));
    37                                 for j=1:size(list,1),
    38                                         fielddisplay(self,list{j,1},list{j,3});
     137                                nat=self.nature{i};
     138                                switch nat
     139                                case 'ice'
     140                                        disp(sprintf('   \nIce:'));
     141                                        fielddisplay(self,'rho_ice','ice density [kg/m^3]');
     142                                        fielddisplay(self,'rho_water','ocean water density [kg/m^3]');
     143                                        fielddisplay(self,'rho_freshwater','fresh water density [kg/m^3]');
     144                                        fielddisplay(self,'mu_water','water viscosity [N s/m^2]');
     145                                        fielddisplay(self,'heatcapacity','heat capacity [J/kg/K]');
     146                                        fielddisplay(self,'thermalconductivity',['ice thermal conductivity [W/m/K]']);
     147                                        fielddisplay(self,'temperateiceconductivity','temperate ice thermal conductivity [W/m/K]');
     148                                        fielddisplay(self,'meltingpoint','melting point of ice at 1atm in K');
     149                                        fielddisplay(self,'latentheat','latent heat of fusion [J/kg]');
     150                                        fielddisplay(self,'beta','rate of change of melting point with pressure [K/Pa]');
     151                                        fielddisplay(self,'mixed_layer_capacity','mixed layer capacity [W/kg/K]');
     152                                        fielddisplay(self,'thermal_exchange_velocity','thermal exchange velocity [m/s]');
     153                                        fielddisplay(self,'rheology_B','flow law parameter [Pa/s^(1/n)]');
     154                                        fielddisplay(self,'rheology_n','Glen''s flow law exponent');
     155                                        fielddisplay(self,'rheology_law',['law for the temperature dependance of the rheology: ''None'', ''BuddJacka'', Cuffey'', ''CuffeyTemperate'', ''Paterson'', ''Arrhenius'' or ''LliboutryDuval''']);
     156                                case 'litho'
     157                                        disp(sprintf('   \nLitho:'));
     158                                        fielddisplay(self,'numlayers','number of layers (default 2)');
     159                                        fielddisplay(self,'radius','array describing the radius for each interface (numlayers+1) [m]');
     160                                        fielddisplay(self,'viscosity','array describing each layer''s viscosity (numlayers) [Pa.s]');
     161                                        fielddisplay(self,'lame_lambda','array describing the lame lambda parameter (numlayers) [Pa]');
     162                                        fielddisplay(self,'lame_mu','array describing the shear modulus for each layers (numlayers) [Pa]');
     163                                        fielddisplay(self,'burgers_viscosity','array describing each layer''s transient viscosity, only for Burgers rheologies  (numlayers) [Pa.s]');
     164                                        fielddisplay(self,'burgers_mu','array describing each layer''s transient shear modulus, only for Burgers rheologies  (numlayers) [Pa]');
     165                                        fielddisplay(self,'isburgers','array describing whether we adopt a MaxWell (0) or Burgers (1) rheology (default 0)');
     166                                        fielddisplay(self,'density','array describing each layer''s density (numlayers) [kg/m^3]');
     167                                        fielddisplay(self,'issolid','array describing whether the layer is solid or liquid (default 1) (numlayers)');
     168                                otherwise
     169                                        error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')');
    39170                                end
    40171                        end
    … …  
    42173                function md = checkconsistency(self,md,solution,analyses) % {{{
    43174
    44                         if isprop(self,'rho_ice'), md = checkfield(md,'fieldname','materials.rho_ice','>',0); end
    45                         if isprop(self,'rho_water'), md = checkfield(md,'fieldname','materials.rho_water','>',0); end
    46                         if isprop(self,'rho_freshwater'), md = checkfield(md,'fieldname','materials.rho_freshwater','>',0); end
    47                         if isprop(self,'mu_water'), md = checkfield(md,'fieldname','materials.mu_water','>',0); end
    48                         if isprop(self,'rheology_B'), md = checkfield(md,'fieldname','materials.rheology_B','>',0,'timeseries',1,'NaN',1,'Inf',1); end
    49                         if isprop(self,'rheology_n'), md = checkfield(md,'fieldname','materials.rheology_n','>',0,'size',[md.mesh.numberofelements 1]); end
    50                         if isprop(self,'rheology_law'), md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'BuddJacka' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval'}); end
    51                         if ~ismember('LoveAnalysis',analyses), return; end
    52                         if isprop(self,'numlayers'), md = checkfield(md,'fieldname','materials.numlayers','NaN',1,'Inf',1,'>',0,'numel',1); end
    53                         if isprop(self,'radius'), md = checkfield(md,'fieldname','materials.radius','NaN',1,'Inf',1,'size',[md.materials.numlayers+1 1],'>',0); end
    54                         if isprop(self,'lame_mu'), md = checkfield(md,'fieldname','materials.lame_mu','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0); end
    55                         if isprop(self,'lame_lambda'), md = checkfield(md,'fieldname','materials.lame_lambda','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0); end
    56                         if isprop(self,'issolid'), md = checkfield(md,'fieldname','materials.issolid','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0,'<',2); end
    57                         if isprop(self,'density'), md = checkfield(md,'fieldname','materials.density','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>',0); end
    58                         if isprop(self,'viscosity'), md = checkfield(md,'fieldname','materials.viscosity','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0); end
    59                         if isprop(self,'isburgers'), md = checkfield(md,'fieldname','materials.isburgers','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0,'<=',1); end
    60                         if isprop(self,'burgers_viscosity'), md = checkfield(md,'fieldname','materials.burgers_viscosity','Inf',1,'size',[md.materials.numlayers 1],'>=',0); end
    61                         if isprop(self,'burgers_mu'), md = checkfield(md,'fieldname','materials.burgers_mu','Inf',1,'size',[md.materials.numlayers 1],'>=',0); end
    62 
    63                         for i=1:md.materials.numlayers,
    64                                 if md.materials.isburgers(i) & (isnan(md.materials.burgers_viscosity(i) | isnan(md.materials.burgers_mu(i)))),
    65                                         error('materials checkconsistency error message: Litho burgers_viscosity or burgers_mu has NaN values, inconsistent with isburgers choice');
    66                                 end
    67                         end
    68                         if md.materials.issolid(1)==0 | md.materials.lame_mu(1)==0
    69                                 error('First layer must be solid (issolid(1) > 0 AND lame_mu(1) > 0). Add a weak inner core if necessary.');
    70                         end
    71                         ind=find(md.materials.issolid==0);
    72                         if sum(ismember(diff(ind),1)>=1) %if there are at least two consecutive indices that contain issolid = 0
    73                                 error(['Fluid layers detected at layers #', num2str(ind'), ', but having 2 or more adjacent fluid layers is not supported yet. Consider merging them.'])
     175                        for i=1:length(self.nature),
     176                                nat=self.nature{i};
     177                                switch nat
     178                                case 'ice'
     179                                        md = checkfield(md,'fieldname','materials.rho_ice','>',0);
     180                                        md = checkfield(md,'fieldname','materials.rho_water','>',0);
     181                                        md = checkfield(md,'fieldname','materials.rho_freshwater','>',0);
     182                                        md = checkfield(md,'fieldname','materials.mu_water','>',0);
     183                                        md = checkfield(md,'fieldname','materials.rheology_B','>',0,'timeseries',1,'NaN',1,'Inf',1);
     184                                        md = checkfield(md,'fieldname','materials.rheology_n','>',0,'size',[md.mesh.numberofelements 1]);
     185                                        md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'BuddJacka' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval'});
     186                                case 'litho'
     187                                        if ~ismember('LoveAnalysis',analyses), return; end
     188                                        md = checkfield(md,'fieldname','materials.numlayers','NaN',1,'Inf',1,'>',0,'numel',1);
     189                                        md = checkfield(md,'fieldname','materials.radius','NaN',1,'Inf',1,'size',[md.materials.numlayers+1 1],'>',0);
     190                                        md = checkfield(md,'fieldname','materials.lame_mu','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0);
     191                                        md = checkfield(md,'fieldname','materials.lame_lambda','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0);
     192                                        md = checkfield(md,'fieldname','materials.issolid','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0,'<',2);
     193                                        md = checkfield(md,'fieldname','materials.density','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>',0);
     194                                        md = checkfield(md,'fieldname','materials.viscosity','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0);
     195                                        md = checkfield(md,'fieldname','materials.isburgers','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>=',0,'<=',1);
     196                                        md = checkfield(md,'fieldname','materials.burgers_viscosity','Inf',1,'size',[md.materials.numlayers 1],'>=',0);
     197                                        md = checkfield(md,'fieldname','materials.burgers_mu','Inf',1,'size',[md.materials.numlayers 1],'>=',0);
     198
     199                                        for i=1:md.materials.numlayers,
     200                                                if md.materials.isburgers(i) & (isnan(md.materials.burgers_viscosity(i) | isnan(md.materials.burgers_mu(i)))),
     201                                                        error('materials checkconsistency error message: Litho burgers_viscosity or burgers_mu has NaN values, inconsistent with isburgers choice');
     202                                                end
     203                                        end
     204                                        if md.materials.issolid(1)==0 | md.materials.lame_mu(1)==0
     205                                                        error('First layer must be solid (issolid(1) > 0 AND lame_mu(1) > 0). Add a weak inner core if necessary.');
     206                                        end
     207                                        ind=find(md.materials.issolid==0);
     208                                        if sum(ismember(diff(ind),1)>=1) %if there are at least two consecutive indices that contain issolid = 0
     209                                                        error(['Fluid layers detected at layers #', num2str(ind'), ', but having 2 or more adjacent fluid layers is not supported yet. Consider merging them.'])
     210                                        end
     211                                otherwise
     212                                        error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')');
     213                                end
    74214                        end
    75215
    … …  
    81221                        WriteData(fid,prefix,'name','md.materials.nature','data',naturetointeger(self.nature),'format','IntMat','mattype',3);
    82222
    83                         if isprop(self,'rho_ice'),                   WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_ice','format','Double'); end
    84                         if isprop(self,'rho_water'),                 WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_water','format','Double'); end
    85                         if isprop(self,'rho_freshwater'),            WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_freshwater','format','Double'); end
    86                         if isprop(self,'mu_water'),                  WriteData(fid,prefix,'object',self,'class','materials','fieldname','mu_water','format','Double'); end
    87                         if isprop(self,'heatcapacity'),              WriteData(fid,prefix,'object',self,'class','materials','fieldname','heatcapacity','format','Double'); end
    88                         if isprop(self,'latentheat'),                WriteData(fid,prefix,'object',self,'class','materials','fieldname','latentheat','format','Double'); end
    89                         if isprop(self,'thermalconductivity'),       WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermalconductivity','format','Double'); end
    90                         if isprop(self,'temperateiceconductivity'),  WriteData(fid,prefix,'object',self,'class','materials','fieldname','temperateiceconductivity','format','Double'); end
    91                         if isprop(self,'meltingpoint'),              WriteData(fid,prefix,'object',self,'class','materials','fieldname','meltingpoint','format','Double'); end
    92                         if isprop(self,'beta'),                      WriteData(fid,prefix,'object',self,'class','materials','fieldname','beta','format','Double'); end
    93                         if isprop(self,'mixed_layer_capacity'),      WriteData(fid,prefix,'object',self,'class','materials','fieldname','mixed_layer_capacity','format','Double'); end
    94                         if isprop(self,'thermal_exchange_velocity'), WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermal_exchange_velocity','format','Double'); end
    95                         if isprop(self,'rheology_B'),                WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_B','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); end
    96                         if isprop(self,'rheology_n'),                WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_n','format','DoubleMat','mattype',2); end
    97                         if isprop(self,'rheology_law'),              Write Data(fid,prefix,'data',self.rheology_law,'nae','md.materials.rheology_law','format','String'); end
    98                         if isprop(self,'numlayers'),                 WriteData(fid,prefix,'object',self,'class','materials','fieldname','numlayers','format','Integer'); end
    99                         if isprop(self,'radius'),                    WriteData(fid,prefix,'object',self,'class','materials','fieldname','radius','format','DoubleMat','mattype',3); end
    100                         if isprop(self,'lame_mu'),                   WriteData(fid,prefix,'object',self,'class','materials','fieldname','lame_mu','format','DoubleMat','mattype',3); end
    101                         if isprop(self,'lame_lambda'),               WriteData(fid,prefix,'object',self,'class','materials','fieldname','lame_lambda','format','DoubleMat','mattype',3); end
    102                         if isprop(self,'issolid'),                   WriteData(fid,prefix,'object',self,'class','materials','fieldname','issolid','format','DoubleMat','mattype',3); end
    103                         if isprop(self,'density'),                   WriteData(fid,prefix,'object',self,'class','materials','fieldname','density','format','DoubleMat','mattype',3); end
    104                         if isprop(self,'viscosity'),                 WriteData(fid,prefix,'object',self,'class','materials','fieldname','viscosity','format','DoubleMat','mattype',3); end
    105                         if isprop(self,'isburgers'),                 WriteData(fid,prefix,'object',self,'class','materials','fieldname','isburgers','format','DoubleMat','mattype',3); end
    106                         if isprop(self,'burgers_viscosity'),         WriteData(fid,prefix,'object',self,'class','materials','fieldname','burgers_viscosity','format','DoubleMat','mattype',3); end
    107                         if isprop(self,'burgers_mu'),                WriteData(fid,prefix,'object',self,'class','materials','fieldname','burgers_mu','format','DoubleMat','mattype',3); end
    108         end % }}}
     223                        for i=1:length(self.nature),
     224                                nat=self.nature{i};
     225                                switch nat
     226                                case 'ice'
     227                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_ice','format','Double');
     228                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_water','format','Double');
     229                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_freshwater','format','Double');
     230                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','mu_water','format','Double');
     231                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','heatcapacity','format','Double');
     232                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','latentheat','format','Double');
     233                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermalconductivity','format','Double');
     234                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','temperateiceconductivity','format','Double');
     235                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','meltingpoint','format','Double');
     236                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','beta','format','Double');
     237                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','mixed_layer_capacity','format','Double');
     238                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermal_exchange_velocity','format','Double');
     239                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_B','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
     240                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_n','format','DoubleMat','mattype',2);
     241                                        WriteData(fid,prefix,'data',self.rheology_law,'name','md.materials.rheology_law','format','String');
     242                                case 'litho'
     243                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','numlayers','format','Integer');
     244                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','radius','format','DoubleMat','mattype',3);
     245                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','lame_mu','format','DoubleMat','mattype',3);
     246                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','lame_lambda','format','DoubleMat','mattype',3);
     247                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','issolid','format','DoubleMat','mattype',3);
     248                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','density','format','DoubleMat','mattype',3);
     249                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','viscosity','format','DoubleMat','mattype',3);
     250                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','isburgers','format','DoubleMat','mattype',3);
     251                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','burgers_viscosity','format','DoubleMat','mattype',3);
     252                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','burgers_mu','format','DoubleMat','mattype',3);
     253                                otherwise
     254                                        error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')');
     255                                end
     256                        end
     257                end % }}}
    109258                function self = extrude(self,md) % {{{
    110                         if isprop(self,'rheology_B'), self.rheology_B=project3d(md,'vector',self.rheology_B,'type','node'); end
    111                         if isprop(self,'rheology_n'), self.rheology_n=project3d(md,'vector',self.rheology_n,'type','element'); end
     259                        for i=1:length(self.nature),
     260                                nat=self.nature{i};
     261                                switch nat
     262                                case 'ice'
     263                                        self.rheology_B=project3d(md,'vector',self.rheology_B,'type','node');
     264                                        self.rheology_n=project3d(md,'vector',self.rheology_n,'type','element');
     265                                end
     266                        end
    112267                end % }}}
    113268                function savemodeljs(self,fid,modelname) % {{{
    114 
    115                         proplist = properties(self);
    116 
    117                         for i=1:numel(proplist),
    118                                 if numel(self.(proplist{i}))==1,
    119                                         writejsdouble(fid,[modelname '.materials.' proplist{i}],self.(proplist{i}));
    120                                 elseif size(self.(proplist{i}),2)==1,
    121                                         writejs1Darray(fid,[modelname '.materials.' proplist{i}],self.(proplist{i}));
    122                                 else
    123                                         error('not supported yet');
    124                                 end
    125                         end
    126                 end % }}}
    127                 function list = listoffields(self,nat) % {{{
    128 
    129                         switch nat
    130                                 case 'ice'
    131                                         list={...
    132                                                 'rho_ice'                   , 917., 'ice density [kg/m^3]'; ...
    133                                                 'rho_water'                 ,1023., 'ocean water density [kg/m^3]'; ...
    134                                                 'rho_freshwater'            ,1000., 'fresh water density [kg/m^3]';...
    135                                                 'mu_water'                  ,0.001787, 'water viscosity [N s/m^2]';...
    136                                                 'heatcapacity'              ,2093., 'heat capacity [J/kg/K]';...
    137                                                 'latentheat'                ,3.34e5, 'latent heat of fusion [J/kg]';...
    138                                                 'thermalconductivity'       ,2.4, 'ice thermal conductivity [W/m/K]';...
    139                                                 'temperateiceconductivity'  ,.24, 'temperate ice thermal conductivity [W/m/K]';...
    140                                                 'meltingpoint'              ,273.15, 'melting point of ice at 1 atm in K';...
    141                                                 'beta'                      ,9.8e-8 'rate of change of melting point with pressure [K/Pa]';...
    142                                                 'mixed_layer_capacity'      ,3974., 'mixed layer capacity [W/kg/K]'; ...
    143                                                 'thermal_exchange_velocity' ,1e-4 'thermal exchange velocity [m/s]'; ...
    144                                                 'rheology_B'                ,NaN, 'flow law parameter [Pa/s^(1/n)]'; ...
    145                                                 'rheology_n'                ,NaN, 'Glen''s flow law exponent'; ...
    146                                                 'rheology_law'              ,'Paterson','law for the temperature dependance of the rheology: ''None'' , ''BuddJacka'' , Cuffey'' , ''CuffeyTemperate'' , ''Paterson'' , ''Arrhenius'' or ''LliboutryDuval'''};
    147                                 case 'damageice'
    148                                         list = listoffields(self,'ice'),
    149                                 case 'enhancedice'
    150                                         list = listoffields(self,'ice');
    151                                         list(end+1,:) = {'rheology_E',NaN,'enhancement factor'};
    152                                 case 'estarice'
    153                                         list = listoffields(self,'ice');
    154                                         list(end+1,:) = {'rheology_Ec',NaN,'compressive enhancement factor'};
    155                                         list(end+1,:) = {'rheology_Es',NaN,'shear enhancement factor'};
    156                                         list(find(strcmp(list(:,1),'rheology_n')),:) = [];
    157                                 case 'litho'
    158                                         %Note for the radius:
    159                                         %   center of the earth (approximation, must not be 0), then the lab (lithosphere/asthenosphere boundary) then the surface
    160                                         %   (with 1d3 to avoid numerical singularities)
    161                                         %Note: here we consider by default two layers: mantle and lithosphere
    162                                         list={...
    163                                                 'numlayers'         ,2, 'number of layers (default 2)'; ...
    164                                                 'radius'            ,[1e3;6278*1e3;6378*1e3] 'array describing the radius for each interface (numlayers+1) [m]'; ...
    165                                                 'viscosity'         ,[1e21;1e40], 'array describing each layer''s viscosity (numlayers) [Pa.s]'; ...
    166                                                 'lame_lambda'       ,[1.45e11;6.7e10], 'array describing the lame lambda parameter (numlayers) [Pa]'; ...
    167                                                 'lame_mu'           ,[1.45e11;6.7e10], 'array describing the shear modulus for each layers (numlayers) [Pa]'; ...
    168                                                 'burgers_viscosity' ,[NaN;NaN], 'array describing each layer''s transient viscosity, only for Burgers rheologies  (numlayers) [Pa.s]',...
    169                                                 'burgers_mu'        ,[0;0], 'array describing each layer''s transient shear modulus, only for Burgers rheologies  (numlayers) [Pa]'; ...
    170                                                 'isburgers'         ,0, 'array describing whether we adopt a MaxWell (0) or Burgers (1) rheology (default 0)'; ...
    171                                                 'density'           ,[5.51*1e3;5.50*1e3], 'array describing each layer''s density (numlayers) [kg/m^3]'; ...
    172                                                 'issolid'           ,[1;1], 'array describing whether the layer is solid or liquid (default 1) (numlayers)'
    173                                         }
    174                                 otherwise
    175                                         error('nature of the material not supported yet!');
    176                         end
    177                 end % }}}
     269       
     270                        for i=1:length(self.nature),
     271                                nat=self.nature{i};
     272                                switch nat
     273                                case 'ice'
     274                                        writejsdouble(fid,[modelname '.materials.rho_ice'],self.rho_ice);
     275                                        writejsdouble(fid,[modelname '.materials.rho_water'],self.rho_water);
     276                                        writejsdouble(fid,[modelname '.materials.rho_freshwater'],self.rho_freshwater);
     277                                        writejsdouble(fid,[modelname '.materials.mu_water'],self.mu_water);
     278                                        writejsdouble(fid,[modelname '.materials.heatcapacity'],self.heatcapacity);
     279                                        writejsdouble(fid,[modelname '.materials.latentheat'],self.latentheat);
     280                                        writejsdouble(fid,[modelname '.materials.thermalconductivity'],self.thermalconductivity);
     281                                        writejsdouble(fid,[modelname '.materials.temperateiceconductivity'],self.temperateiceconductivity);
     282                                        writejsdouble(fid,[modelname '.materials.meltingpoint'],self.meltingpoint);
     283                                        writejsdouble(fid,[modelname '.materials.beta'],self.beta);
     284                                        writejsdouble(fid,[modelname '.materials.mixed_layer_capacity'],self.mixed_layer_capacity);
     285                                        writejsdouble(fid,[modelname '.materials.thermal_exchange_velocity'],self.thermal_exchange_velocity);
     286                                        writejsdouble(fid,[modelname '.materials.mixed_layer_capacity'],self.mixed_layer_capacity);
     287                                        writejs1Darray(fid,[modelname '.materials.rheology_B'],self.rheology_B);
     288                                        writejs1Darray(fid,[modelname '.materials.rheology_n'],self.rheology_n);
     289                                        writejsstring(fid,[modelname '.materials.rheology_law'],self.rheology_law);
     290                                case 'litho'
     291                                        writejsdouble(fid,[modelname '.materials.numlayers'],self.numlayers);
     292                                        writejsdouble(fid,[modelname '.materials.radius'],self.radius);
     293                                        writejsdouble(fid,[modelname '.materials.lame_mu'],self.lame_mu);
     294                                        writejsdouble(fid,[modelname '.materials.lame_lambda'],self.lame_lambda);
     295                                        writejsdouble(fid,[modelname '.materials.issolid'],self.issolid);
     296                                        writejsdouble(fid,[modelname '.materials.density'],self.density);
     297                                        writejsdouble(fid,[modelname '.materials.viscosity'],self.viscosity);
     298                                        writejsdouble(fid,[modelname '.materials.isburgers'],self.isburgers);
     299                                        writejsdouble(fid,[modelname '.materials.burgers_viscosity'],self.burgers_viscosity);
     300                                        writejsdouble(fid,[modelname '.materials.burgers_mu'],self.burgers_mu);
     301                                otherwise
     302                                        error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')');
     303                                end
     304                        end
     305
     306
     307
     308                end % }}}
     309        end
     310end
    178311                function intnat = naturetointeger(strnat) % {{{
    179312                        intnat=zeros(length(strnat),1);
    180313                        for i=1:length(strnat),
    181314                                switch strnat{i},
    182                                         case 'damageice'
    183                                                 intnat(i)=1;
    184                                         case 'estar'
    185                                                 intnat(i)=2;
    186                                         case 'ice'
    187                                                 intnat(i)=3;
    188                                         case 'enhancedice'
    189                                                 intnat(i)=4;
    190                                         case 'litho'
    191                                                 intnat(i)=5;
    192                                         otherwise
    193                                                 error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')');
    194                                 end
    195                         end
    196                 end % }}}
    197         end
    198 end
     315                                case 'damageice'
     316                                        intnat(i)=1;
     317                                case 'estar'
     318                                        intnat(i)=2;
     319                                case 'ice'
     320                                        intnat(i)=3;
     321                                case 'enhancedice'
     322                                        intnat(i)=4;
     323                                case 'litho'
     324                                        intnat(i)=5;
     325                                otherwise
     326                                        error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')');
     327                                end
     328                        end
     329                end % }}}
Note: See TracChangeset for help on using the changeset viewer.