Changeset 21739


Ignore:
Timestamp:
05/23/17 00:21:35 (8 years ago)
Author:
sjohnsen
Message:

NEW: adding a friction effective pressure flag

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

Legend:

Unmodified
Added
Removed

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

    r21049 r21739  
    99                p           = NaN;
    1010                q           = NaN;
     11                coupling    = 0;        %Silje
     12                effective_pressure = NaN;       %Silje
    1113        end
    1214        methods
    … …  
    1517                        self.p=project3d(md,'vector',self.p,'type','element');
    1618                        self.q=project3d(md,'vector',self.q,'type','element');
     19                        switch self.coupling
     20                                case 0
     21                                case 1
     22                                        self.effective_pressure=project3d(md,'vector',self.effective_pressure,'type','node','layer',1);
     23                                case 2
     24                                        error('not implemented yet');
     25                                otherwise
     26                                        error('not supported yet');             
     27                        end
    1728                end % }}}
    1829                function self = friction(varargin) % {{{
    … …  
    3647                        md = checkfield(md,'fieldname','friction.q','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]);
    3748                        md = checkfield(md,'fieldname','friction.p','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]);
     49                        md = checkfield(md,'fieldname','friction.coupling','numel',[1],'values',[0 1 2]);%Silje
     50                        switch self.coupling
     51                                case 0
     52                                case 1
     53                                        md = checkfield(md,'fieldname','friction.effective_pressure','NaN',1,'Inf',1,'timeseries',1);
     54                                case 2
     55                                        error('not implemented yet');
     56                                otherwise
     57                                        error('not supported yet');             
     58                        end
    3859                end % }}}
    3960                function disp(self) % {{{
    … …  
    4263                        fielddisplay(self,'p','p exponent');
    4364                        fielddisplay(self,'q','q exponent');
     65                        fielddisplay(self,'effective_pressure','Effective Pressure for the forcing if not coupled [Pa]');%Silje
     66                        fielddisplay(self,'coupling','Coupling flag: 0 for default, 1 for forcing(provide md.friction.effective_pressure)  and 2 for coupled(not implemented yet)');%Silje
    4467                end % }}}
    4568                function marshall(self,prefix,md,fid) % {{{
    … …  
    5073                        WriteData(fid,prefix,'object',self,'fieldname','p','format','DoubleMat','mattype',2);
    5174                        WriteData(fid,prefix,'object',self,'fieldname','q','format','DoubleMat','mattype',2);
    52                        
    53 
     75                        WriteData(fid,prefix,'class','friction','object',self,'fieldname','coupling','format','Integer'); %Silje
     76                        switch self.coupling
     77                                case 0
     78                                case 1
     79                                        WriteData(fid,prefix,'class','friction','object',self,'fieldname','effective_pressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
     80                                case 2
     81                                        error('not implemented yet');
     82                                otherwise
     83                                        error('not supported yet');             
     84                        end
    5485                end % }}}
    5586                function savemodeljs(self,fid,modelname) % {{{
    … …  
    5889                        writejs1Darray(fid,[modelname '.friction.p'],self.p);
    5990                        writejs1Darray(fid,[modelname '.friction.q'],self.q);
     91                        writejs1Darray(fid,[modelname '.friction.coupling'],self.coupling);%Silje
     92                        writejs1Darray(fid,[modelname '.friction.effective_pressure'],self.effective_pressure);%Silje
    6093
    6194                end % }}}

  • issm/trunk-jpl/src/m/classes/friction.py

    r21049 r21739  
     1import numpy as np
    12from fielddisplay import fielddisplay
    23from project3d import project3d
    … …  
    1617                self.p           = float('NaN')
    1718                self.q           = float('NaN')
    18 
     19                self.coupling    = 0
     20                self.effective_pressure = float('NaN')
    1921                #set defaults
    2022                self.setdefaultparameters()
    … …  
    2224                #}}}
    2325        def __repr__(self): # {{{
    24                 string="Basal shear stress parameters: Sigma_b = coefficient^2 * Neff ^r * |u_b|^(s-1) * u_b,\n(effective stress Neff=rho_ice*g*thickness+rho_water*g*bed, r=q/p and s=1/p)"
     26                string="Basal shear stress parameters: Sigma_b = coefficient^2 * Neff ^r * |u_b|^(s-1) * u_b,\n(effective stress Neff=rho_ice*g*thickness+rho_water*g*base, r=q/p and s=1/p)"
    2527
    2628                string="%s\n%s"%(string,fielddisplay(self,"coefficient","friction coefficient [SI]"))
    2729                string="%s\n%s"%(string,fielddisplay(self,"p","p exponent"))
    2830                string="%s\n%s"%(string,fielddisplay(self,"q","q exponent"))
     31                string="%s\n%s"%(string,fielddisplay(self,'coupling','Coupling flag: 0 for default, 1 for forcing(provide md.friction.effective_pressure)  and 2 for coupled(not implemented yet)'))
     32                string="%s\n%s"%(string,fielddisplay(self,'effective_pressure','Effective Pressure for the forcing if not coupled [Pa]'))
    2933                return string
    3034                #}}}
    … …  
    3337                self.p=project3d(md,'vector',self.p,'type','element')
    3438                self.q=project3d(md,'vector',self.q,'type','element')
     39                #if self.coupling==0: #doesnt work with empty loop, so just skip it?
     40                if self.coupling==1:
     41                        self.effective_pressure=project3d(md,'vector',self.effective_pressure,'type','node','layer',1)
     42                elif self.coupling==2:
     43                        raise ValueError('coupling not supported yet')
     44                elif self.coupling > 2:
     45                        raise ValueError('md.friction.coupling larger than 2, not supported yet')       
    3546                return self
    3647        #}}}
    … …  
    4758                md = checkfield(md,'fieldname','friction.q','NaN',1,'Inf',1,'size',[md.mesh.numberofelements])
    4859                md = checkfield(md,'fieldname','friction.p','NaN',1,'Inf',1,'size',[md.mesh.numberofelements])
    49 
     60                md = checkfield(md,'fieldname','friction.coupling','numel',[1],'values',[0,1,2])
     61                if self.coupling==1:
     62                        md = checkfield(md,'fieldname','friction.effective_pressure','NaN',1,'Inf',1,'timeseries',1)
     63                elif self.coupling==2:
     64                        raise ValueError('coupling not supported yet')
     65                elif self.coupling > 2:
     66                        raise ValueError('md.friction.coupling larger than 2, not supported yet')
    5067                return md
    5168        # }}}
    … …  
    5572                WriteData(fid,prefix,'object',self,'fieldname','p','format','DoubleMat','mattype',2)
    5673                WriteData(fid,prefix,'object',self,'fieldname','q','format','DoubleMat','mattype',2)
     74                WriteData(fid,prefix,'class','friction','object',self,'fieldname','coupling','format','Integer')
     75                if self.coupling==1:
     76                        WriteData(fid,prefix,'class','friction','object',self,'fieldname','effective_pressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
     77                elif self.coupling==2:
     78                        raise ValueError('coupling not supported yet')
     79                elif self.coupling > 2:
     80                        raise ValueError('md.friction.coupling larger than 2, not supported yet')                       
    5781        # }}}

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

    r21049 r21739  
    66classdef frictionhydro
    77        properties (SetAccess=public)
    8                 Coupling           = 0;
     8                coupling           = 0;
    99                q                  = NaN;
    1010                C                  = NaN;
    … …  
    2828                        %Early return
    2929                        if ~ismember('StressbalanceAnalysis',analyses) & ~ismember('ThermalAnalysis',analyses), return; end
    30                         md = checkfield(md,'fieldname','friction.Coupling','numel',[1],'values',[0 1]);
     30                        md = checkfield(md,'fieldname','friction.coupling','numel',[1],'values',[0 1 2]);
    3131                        md = checkfield(md,'fieldname','friction.q','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]);
    3232                        md = checkfield(md,'fieldname','friction.C','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]);
    3333                        md = checkfield(md,'fieldname','friction.As','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]);
    34                         if self.Coupling==0,
    35                                 md = checkfield(md,'fieldname','friction.effective_pressure','NaN',1,'Inf',1,'timeseries',1);
     34                        switch self.coupling
     35                                case 0
     36                                case 1
     37                                        md = checkfield(md,'fieldname','friction.effective_pressure','NaN',1,'Inf',1,'timeseries',1);
     38                                case 2
     39                                        error('not implemented yet');
     40                                otherwise
     41                                        error('not supported yet');             
     42                        end
    3643            end
    3744                end % }}}
    … …  
    4047                        self.C=project3d(md,'vector',self.C,'type','element');
    4148                        self.As=project3d(md,'vector',self.As,'type','element');
    42                         if self.Coupling==0,
    43                                 self.effective_pressure=project3d(md,'vector',self.effective_pressure,'type','node','layer',1);
     49                        switch self.coupling
     50                                case 0
     51                                case 1
     52                                        self.effective_pressure=project3d(md,'vector',self.effective_pressure,'type','node','layer',1);
     53                                case 2
     54                                        error('not implemented yet');
     55                                otherwise
     56                                        error('not supported yet');             
    4457                        end
    4558          end % }}}
    4659                function disp(self) % {{{
    4760                        disp(sprintf('Effective Pressure based friction law described in Gagliardini 2007'));
    48                         fielddisplay(self,'Coupling','Coupling flag, 1 for coupling and 0 for forcing');
     61                        fielddisplay(self,'coupling','Coupling flag: 0 for default, 1 for forcing(provide md.friction.effective_pressure)  and 2 for coupled(not implemented yet)');
    4962                        fielddisplay(self,'q','friction law exponent q>=1');
    5063                        fielddisplay(self,'C','friction law max value [SI]');
    … …  
    5467                function marshall(self,prefix,md,fid) % {{{
    5568                        WriteData(fid,prefix,'name','md.friction.law','data',3,'format','Integer');
    56                         WriteData(fid,prefix,'class','friction','object',self,'fieldname','Coupling','format','Integer');
     69                        WriteData(fid,prefix,'class','friction','object',self,'fieldname','coupling','format','Integer');
    5770                        WriteData(fid,prefix,'class','friction','object',self,'fieldname','q','format','DoubleMat','mattype',2);
    5871                        WriteData(fid,prefix,'class','friction','object',self,'fieldname','C','format','DoubleMat','mattype',2);
    5972                        WriteData(fid,prefix,'class','friction','object',self,'fieldname','As','format','DoubleMat','mattype',2);
    60                         if self.Coupling==0,
    61                                 WriteData(fid,prefix,'class','friction','object',self,'fieldname','effective_pressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
     73                        switch self.coupling
     74                                case 0
     75                                case 1
     76                                        WriteData(fid,prefix,'class','friction','object',self,'fieldname','effective_pressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
     77                                case 2
     78                                        error('not implemented yet');
     79                                otherwise
     80                                        error('not supported yet');             
    6281                        end
    6382          end % }}}
Note: See TracChangeset for help on using the changeset viewer.