Changeset 22304


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

CHG: added P1xP4 for HO

Location:
issm/trunk-jpl/src/m/classes
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r21049 r22304  
    1111
    1212                omega             = NaN;
     13                slopex            = NaN;
     14                slopey            = NaN;
    1315        end
    1416        methods
    … …  
    5355                        WriteData(fid,prefix,'object',self,'fieldname','stabilization','format','Integer');
    5456
     57                        WriteData(fid,prefix,'object',self,'fieldname','slopex','format','DoubleMat','mattype',1);
     58                        WriteData(fid,prefix,'object',self,'fieldname','slopey','format','DoubleMat','mattype',1);
    5559                        WriteData(fid,prefix,'object',self,'fieldname','omega','format','DoubleMat','mattype',1);
    5660                end % }}}

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

    r21049 r22304  
    9999                        md = checkfield(md,'fieldname','flowequation.isFS','numel',[1],'values',[0 1]);
    100100                        md = checkfield(md,'fieldname','flowequation.fe_SSA','values',{'P1','P1bubble','P1bubblecondensed','P2','P2bubble'});
    101                         md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',{'P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P2bubble','P1xP3','P2xP4'});
     101                        md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',{'P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P2bubble','P1xP3','P1xP4','P2xP4'});
    102102                        md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',{'P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','LATaylorHood','XTaylorHood','OneLayerP4z','CrouzeixRaviart','LACrouzeixRaviart'});
    103103                        md = checkfield(md,'fieldname','flowequation.augmented_lagrangian_r','numel',[1],'>=',0.);

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

    r22004 r22304  
    1616                                self.nature=varargin;
    1717                        end
    18                        
    19                         %check this is acceptable:
     18
     19                        %start filling in the dynamic fields:
    2020                        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'')');
     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};
    2327                                end
    2428                        end
    25                        
    26                         %start filling in the dynamic fields:
    27                         for i=1:length(self.nature),
    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();
    6329
    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
    13230                end % }}}
    13331                function disp(self) % {{{
    … …  
    13533
    13634                        for i=1:length(self.nature),
    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'')');
     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});
    17039                                end
    17140                        end
    … …  
    17342                function md = checkconsistency(self,md,solution,analyses) % {{{
    17443
    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);
     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
    19862
    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'')');
     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');
    21366                                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.'])
    21474                        end
    21575
    … …  
    22181                        WriteData(fid,prefix,'name','md.materials.nature','data',naturetointeger(self.nature),'format','IntMat','mattype',3);
    22282
    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'')');
     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 % }}}
     109                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
     112                end % }}}
     113                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');
    255124                                end
    256125                        end
    257126                end % }}}
    258                 function self = extrude(self,md) % {{{
    259                         for i=1:length(self.nature),
    260                                 nat=self.nature{i};
    261                                 switch nat
     127                function list = listoffields(self,nat) % {{{
     128
     129                        switch nat
    262130                                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
     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!');
    266176                        end
    267177                end % }}}
    268                 function savemodeljs(self,fid,modelname) % {{{
    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
    310 end
    311178                function intnat = naturetointeger(strnat) % {{{
    312179                        intnat=zeros(length(strnat),1);
    313180                        for i=1:length(strnat),
    314181                                switch strnat{i},
    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'')');
     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'')');
    327194                                end
    328195                        end
    329196                end % }}}
     197        end
     198end
Note: See TracChangeset for help on using the changeset viewer.