Changeset 21697


Ignore:
Timestamp:
05/02/17 16:29:43 (8 years ago)
Author:
Eric.Larour
Message:

CHG: finished first prototype materials.m implementation

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r21696 r21697  
    6565                function self = setdefaultparameters(self) % {{{
    6666
    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                                 %surface, then the lab (lithosphere/asthenosphere boundary) then the center of the earth
    116                                 %(with 1d3 to avoid numerical singularities)
    117                                 self.radius=[6378*1e3;6278*1e3;1e3];
    118 
    119                                 self.viscosity=[1e40;1e21]; %lithosphere and mantle viscosity (respectively) [Pa.s]
    120                                 self.lame_mu=[6.7*10^10;1.45*1e11];  % (Pa) %lithosphere and mantle shear modulus (respectively) [Pa]
    121                                 self.lame_lambda=[6.7*10^10;1.45*1e11];  % (Pa) %lithosphere and mantle lamba parameter (respectively) [Pa]
    122                                 self.burgers_viscosity=[NaN;NaN];
    123                                 self.burgers_mu=[NaN;NaN];
    124                                 self.isburgers=[false;false];
    125                                 self.density=[3.32*1e3,3.34*1e3];  % (Pa) %lithosphere and mantle density [kg/m^3]
    126                                 self.issolid=[true,true]; % 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
    132 
     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                                        %surface, then the lab (lithosphere/asthenosphere boundary) then the center of the earth
     116                                        %(with 1d3 to avoid numerical singularities)
     117                                        self.radius=[6378*1e3;6278*1e3;1e3];
     118
     119                                        self.viscosity=[1e40;1e21]; %lithosphere and mantle viscosity (respectively) [Pa.s]
     120                                        self.lame_mu=[6.7*10^10;1.45*1e11];  % (Pa) %lithosphere and mantle shear modulus (respectively) [Pa]
     121                                        self.lame_lambda=[6.7*10^10;1.45*1e11];  % (Pa) %lithosphere and mantle lamba parameter (respectively) [Pa]
     122                                        self.burgers_viscosity=[NaN;NaN];
     123                                        self.burgers_mu=[NaN;NaN];
     124                                        self.isburgers=[false;false];
     125                                        self.density=[3.32*1e3;3.34*1e3];  % (Pa) %lithosphere and mantle density [kg/m^3]
     126                                        self.issolid=[true;true]; % 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
    133132                end % }}}
    134133                function disp(self) % {{{
    … …  
    173172                end % }}}
    174173                function md = checkconsistency(self,md,solution,analyses) % {{{
    175                         md = checkfield(md,'fieldname','materials.rho_ice','>',0);
    176                         md = checkfield(md,'fieldname','materials.rho_water','>',0);
    177                         md = checkfield(md,'fieldname','materials.rho_freshwater','>',0);
    178                         md = checkfield(md,'fieldname','materials.mu_water','>',0);
    179                         md = checkfield(md,'fieldname','materials.rheology_B','>',0,'timeseries',1,'NaN',1,'Inf',1);
    180                         md = checkfield(md,'fieldname','materials.rheology_n','>',0,'size',[md.mesh.numberofelements 1]);
    181                         md = checkfield(md,'fieldname','materials.rheology_law','values',{'None' 'BuddJacka' 'Cuffey' 'CuffeyTemperate' 'Paterson' 'Arrhenius' 'LliboutryDuval'});
    182 
    183                         if ismember('GiaAnalysis',analyses),
    184                                 md = checkfield(md,'fieldname','materials.lithosphere_shear_modulus','>',0,'numel',1);
    185                                 md = checkfield(md,'fieldname','materials.lithosphere_density','>',0,'numel',1);
    186                                 md = checkfield(md,'fieldname','materials.mantle_shear_modulus','>',0,'numel',1);
    187                                 md = checkfield(md,'fieldname','materials.mantle_density','>',0,'numel',1);
    188                         end
    189                         if ismember('SealevelriseAnalysis',analyses),
    190                                 md = checkfield(md,'fieldname','materials.earth_density','>',0,'numel',1);
    191                         end
    192 
     174
     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,'<',2);
     196                                        md = checkfield(md,'fieldname','materials.burgers_viscosity','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>',0);
     197                                        md = checkfield(md,'fieldname','materials.burgers_mu','NaN',1,'Inf',1,'size',[md.materials.numlayers 1],'>',0);
     198
     199                                otherwise
     200                                        error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')');
     201                                end
     202                        end
     203
     204                end % }}}
     205                function intnat = naturetointeger(strnat) % {{{
     206                        intnat=zeros(length(strnat),1);
     207                        for i=1:length(strnat),
     208                                switch strnat{i},
     209                                case 'damageice'
     210                                        intnat(i)=1;
     211                                case 'estar'
     212                                        intnat(i)=2;
     213                                case 'ice'
     214                                        intnat(i)=3;
     215                                case 'enhancedice'
     216                                        intnat(i)=4;
     217                                case 'litho'
     218                                        intnat(i)=5;
     219                                otherwise
     220                                        error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')');
     221                                end
     222                        end
    193223                end % }}}
    194224                function marshall(self,prefix,md,fid) % {{{
    195                         WriteData(fid,prefix,'name','md.materials.type','data',3,'format','Integer');
    196                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_ice','format','Double');
    197                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_water','format','Double');
    198                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_freshwater','format','Double');
    199                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','mu_water','format','Double');
    200                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','heatcapacity','format','Double');
    201                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','latentheat','format','Double');
    202                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermalconductivity','format','Double');
    203                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','temperateiceconductivity','format','Double');
    204                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','meltingpoint','format','Double');
    205                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','beta','format','Double');
    206                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','mixed_layer_capacity','format','Double');
    207                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermal_exchange_velocity','format','Double');
    208                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_B','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    209                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_n','format','DoubleMat','mattype',2);
    210                         WriteData(fid,prefix,'data',self.rheology_law,'name','md.materials.rheology_law','format','String');
    211 
    212                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','lithosphere_shear_modulus','format','Double');
    213                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','lithosphere_density','format','Double','scale',10^3);
    214                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','mantle_shear_modulus','format','Double');
    215                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','mantle_density','format','Double','scale',10^3);
    216                         WriteData(fid,prefix,'object',self,'class','materials','fieldname','earth_density','format','Double');
     225
     226                        %1: MatdamageiceEnum 2: MatestarEnum 3: MaticeEnum 4: MatenhancediceEnum 5: MaterialsEnum
     227                        WriteData(fid,prefix,'name','md.materials.type','data',6,'format','Integer');
     228                        WriteData(fid,prefix,'name','md.materials.nature','data',naturetointeger(self.nature),'format','IntMat','mattype',1);
     229
     230                        for i=1:length(self.nature),
     231                                nat=self.nature{i};
     232                                switch nat
     233                                case 'ice'
     234                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_ice','format','Double');
     235                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_water','format','Double');
     236                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_freshwater','format','Double');
     237                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','mu_water','format','Double');
     238                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','heatcapacity','format','Double');
     239                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','latentheat','format','Double');
     240                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermalconductivity','format','Double');
     241                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','temperateiceconductivity','format','Double');
     242                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','meltingpoint','format','Double');
     243                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','beta','format','Double');
     244                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','mixed_layer_capacity','format','Double');
     245                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','thermal_exchange_velocity','format','Double');
     246                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_B','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
     247                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','rheology_n','format','DoubleMat','mattype',2);
     248                                        WriteData(fid,prefix,'data',self.rheology_law,'name','md.materials.rheology_law','format','String');
     249                                case 'litho'
     250                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','numlayers','format','Double');
     251                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','radius','format','DoubleMat','mattype',1);
     252                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','lame_mu','format','DoubleMat','mattype',1);
     253                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','lame_lambda ','format','DoubleMat','mattype',1);
     254                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','issolid','format','BooleanMat','mattype',1);
     255                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','density','format','DoubleMat','mattype',1);
     256                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','viscosity','format','DoubleMat','mattype',1);
     257                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','isburgers','format','BooleanMat','mattype',1);
     258                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','burgers_viscosity','format','DoubleMat','mattype',1);
     259                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','burgers_mu','format','DoubleMat','mattype',1);
     260                                otherwise
     261                                        error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')');
     262                                end
     263                        end
    217264                end % }}}
    218265                function self = extrude(self,md) % {{{
    219                         self.rheology_B=project3d(md,'vector',self.rheology_B,'type','node');
    220                         self.rheology_n=project3d(md,'vector',self.rheology_n,'type','element');
     266                        for i=1:length(self.nature),
     267                                nat=self.nature{i};
     268                                switch nat
     269                                case 'ice'
     270                                        self.rheology_B=project3d(md,'vector',self.rheology_B,'type','node');
     271                                        self.rheology_n=project3d(md,'vector',self.rheology_n,'type','element');
     272                                end
     273                        end
    221274                end % }}}
    222275                function savemodeljs(self,fid,modelname) % {{{
    223                
    224                         writejsdouble(fid,[modelname '.materials.rho_ice'],self.rho_ice);
    225                         writejsdouble(fid,[modelname '.materials.rho_water'],self.rho_water);
    226                         writejsdouble(fid,[modelname '.materials.rho_freshwater'],self.rho_freshwater);
    227                         writejsdouble(fid,[modelname '.materials.mu_water'],self.mu_water);
    228                         writejsdouble(fid,[modelname '.materials.heatcapacity'],self.heatcapacity);
    229                         writejsdouble(fid,[modelname '.materials.latentheat'],self.latentheat);
    230                         writejsdouble(fid,[modelname '.materials.thermalconductivity'],self.thermalconductivity);
    231                         writejsdouble(fid,[modelname '.materials.temperateiceconductivity'],self.temperateiceconductivity);
    232                         writejsdouble(fid,[modelname '.materials.meltingpoint'],self.meltingpoint);
    233                         writejsdouble(fid,[modelname '.materials.beta'],self.beta);
    234                         writejsdouble(fid,[modelname '.materials.mixed_layer_capacity'],self.mixed_layer_capacity);
    235                         writejsdouble(fid,[modelname '.materials.thermal_exchange_velocity'],self.thermal_exchange_velocity);
    236                         writejsdouble(fid,[modelname '.materials.mixed_layer_capacity'],self.mixed_layer_capacity);
    237                         writejs1Darray(fid,[modelname '.materials.rheology_B'],self.rheology_B);
    238                         writejs1Darray(fid,[modelname '.materials.rheology_n'],self.rheology_n);
    239                         writejsstring(fid,[modelname '.materials.rheology_law'],self.rheology_law);
    240                         writejsdouble(fid,[modelname '.materials.lithosphere_shear_modulus'],self.lithosphere_shear_modulus);
    241                         writejsdouble(fid,[modelname '.materials.lithosphere_density'],self.lithosphere_density);
    242                         writejsdouble(fid,[modelname '.materials.mantle_shear_modulus'],self.mantle_shear_modulus);
    243                         writejsdouble(fid,[modelname '.materials.mantle_density'],self.mantle_density);
    244                         writejsdouble(fid,[modelname '.materials.earth_density'],self.earth_density);
     276       
     277                        for i=1:length(self.nature),
     278                                nat=self.nature{i};
     279                                switch nat
     280                                case 'ice'
     281                                        writejsdouble(fid,[modelname '.materials.rho_ice'],self.rho_ice);
     282                                        writejsdouble(fid,[modelname '.materials.rho_water'],self.rho_water);
     283                                        writejsdouble(fid,[modelname '.materials.rho_freshwater'],self.rho_freshwater);
     284                                        writejsdouble(fid,[modelname '.materials.mu_water'],self.mu_water);
     285                                        writejsdouble(fid,[modelname '.materials.heatcapacity'],self.heatcapacity);
     286                                        writejsdouble(fid,[modelname '.materials.latentheat'],self.latentheat);
     287                                        writejsdouble(fid,[modelname '.materials.thermalconductivity'],self.thermalconductivity);
     288                                        writejsdouble(fid,[modelname '.materials.temperateiceconductivity'],self.temperateiceconductivity);
     289                                        writejsdouble(fid,[modelname '.materials.meltingpoint'],self.meltingpoint);
     290                                        writejsdouble(fid,[modelname '.materials.beta'],self.beta);
     291                                        writejsdouble(fid,[modelname '.materials.mixed_layer_capacity'],self.mixed_layer_capacity);
     292                                        writejsdouble(fid,[modelname '.materials.thermal_exchange_velocity'],self.thermal_exchange_velocity);
     293                                        writejsdouble(fid,[modelname '.materials.mixed_layer_capacity'],self.mixed_layer_capacity);
     294                                        writejs1Darray(fid,[modelname '.materials.rheology_B'],self.rheology_B);
     295                                        writejs1Darray(fid,[modelname '.materials.rheology_n'],self.rheology_n);
     296                                        writejsstring(fid,[modelname '.materials.rheology_law'],self.rheology_law);
     297                                case 'litho'
     298                                        writejsdouble(fid,[modelname '.materials.numlayers'],self.numlayers);
     299                                        writejsdouble(fid,[modelname '.materials.radius'],self.radius);
     300                                        writejsdouble(fid,[modelname '.materials.lame_mu'],self.lame_mu);
     301                                        writejsdouble(fid,[modelname '.materials.lame_lambda'],self.lame_lambda);
     302                                        writejsdouble(fid,[modelname '.materials.issolid'],self.issolid);
     303                                        writejsdouble(fid,[modelname '.materials.density'],self.density);
     304                                        writejsdouble(fid,[modelname '.materials.viscosity'],self.viscosity);
     305                                        writejsdouble(fid,[modelname '.materials.isburgers'],self.isburgers);
     306                                        writejsdouble(fid,[modelname '.materials.burgers_viscosity'],self.burgers_viscosity);
     307                                        writejsdouble(fid,[modelname '.materials.burgers_mu'],self.burgers_mu);
     308                                otherwise
     309                                        error('materials constructor error message: nature of the material not supported yet! (''ice'' or ''litho'')');
     310                                end
     311                        end
     312
     313
    245314
    246315                end % }}}
Note: See TracChangeset for help on using the changeset viewer.