Changeset 26187


Ignore:
Timestamp:
04/08/21 14:07:35 (4 years ago)
Author:
Mathieu Morlighem
Message:

CHG: added possibility to couple with friciton schoof

Location:
issm/trunk-jpl/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

    r26155 r26187  
    921921                        break;
    922922                case 11:
     923                        iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling");
    923924                        iomodel->FetchDataToInput(inputs,elements,"md.friction.m",FrictionMEnum);
    924925                        iomodel->FetchDataToInput(inputs,elements,"md.friction.C",FrictionCEnum);
    925926                        iomodel->FetchDataToInput(inputs,elements,"md.friction.Cmax",FrictionCmaxEnum);
     927                        if(FrictionCoupling==3){
     928                                iomodel->FetchDataToInput(inputs,elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);}
     929                        else if(FrictionCoupling==4){
     930                                iomodel->FetchDataToInput(inputs,elements,"md.friction.effective_pressure",EffectivePressureEnum);
     931                        }
    926932                        break;
    927933                case 12:
    … …  
    938944                        else if(FrictionCoupling==4){
    939945                                iomodel->FetchDataToInput(inputs,elements,"md.friction.effective_pressure",EffectivePressureEnum);
    940 
    941946                        }
    942947                        break;
    … …  
    10471052                        break;
    10481053                case 11:
    1049                         parameters->AddObject(new IntParam(FrictionCouplingEnum,2));
     1054                        parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
    10501055                        parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
    10511056                        break;

  • issm/trunk-jpl/src/c/classes/IoModel.cpp

    r26073 r26187  
    616616                                if(fread(record_name,record_name_size*sizeof(char),1,fid)!=0){};
    617617                               
    618                                 _error_("error while looking in binary file. String " << record_name << " a string of size "<<record_name_size);
     618                                _error_("error while looking in binary file. String \"" << record_name << "\" is a string of size "<<record_name_size);
    619619                        }
    620620

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

    r25688 r26187  
    66classdef friction
    77        properties (SetAccess=public)
    8                 coefficient = NaN;
    9                 p           = NaN;
    10                 q           = NaN;
    11                 coupling    = 0;
    12                 effective_pressure = NaN;
     8                coefficient              = NaN;
     9                p                        = NaN;
     10                q                        = NaN;
     11                coupling                 = 0;
     12                effective_pressure       = NaN;
    1313                effective_pressure_limit = 0;
    1414        end
    … …  
    1818                        self.p=project3d(md,'vector',self.p,'type','element');
    1919                        self.q=project3d(md,'vector',self.q,'type','element');
    20                         switch self.coupling
    21                                 case 0
    22                                 case 1
    23                                 case 2
    24                                 case 3
    25                                         self.effective_pressure=project3d(md,'vector',self.effective_pressure,'type','node','layer',1);
    26                                 case 4
    27                                         self.effective_pressure=project3d(md,'vector',self.effective_pressure,'type','node','layer',1);
    28                                 otherwise
    29                                         error('not supported yet');             
     20                        if self.coupling==3 || self.coupling==4
     21                                self.effective_pressure=project3d(md,'vector',self.effective_pressure,'type','node','layer',1);
    3022                        end
    3123                end % }}}
    … …  
    5749                        md = checkfield(md,'fieldname','friction.coupling','numel',[1],'values',[0:4]);
    5850                        md = checkfield(md,'fieldname','friction.effective_pressure_limit','numel',[1],'>=',0);
    59                         switch self.coupling
    60                                 case 0
    61                                 case 1
    62                                 case 2
    63                                 case 3
    64                                         md = checkfield(md,'fieldname','friction.effective_pressure','NaN',1,'Inf',1,'timeseries',1);
    65                                 case 4
    66                                
    67                                 otherwise
    68                                         error('not supported yet');
     51         if self.coupling==3
     52            md = checkfield(md,'fieldname','friction.effective_pressure','NaN',1,'Inf',1,'timeseries',1);
    6953                        end
    7054                end % }}}
    … …  
    9276                        WriteData(fid,prefix,'class','friction','object',self,'fieldname','coupling','format','Integer');
    9377                        WriteData(fid,prefix,'object',self,'class','friction','fieldname','effective_pressure_limit','format','Double');
    94                         switch self.coupling
    95                                 case 0
    96                                 case 1
    97                                 case 2
    98                                 case 3
    99                                         WriteData(fid,prefix,'class','friction','object',self,'fieldname','effective_pressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    100                                 case 4
    101                                         WriteData(fid,prefix,'class','friction','object',self,'fieldname','effective_pressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    102                                 otherwise
    103                                         error('not supported yet');             
     78                        if self.coupling==3 || self.coupling==4
     79                                WriteData(fid,prefix,'class','friction','object',self,'fieldname','effective_pressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    10480                        end
    10581                end % }}}

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

    r25688 r26187  
    2323        self.setdefaultparameters()
    2424    #}}}
    25 
    2625    def __repr__(self):  # {{{
    2726        s = 'Basal shear stress parameters: Sigma_b = coefficient^2 * Neff ^r * |u_b|^(s - 1) * u_b,\n'
    … …  
    3534        return s
    3635    #}}}
    37 
    3836    def setdefaultparameters(self):  # {{{
    3937        self.effective_pressure_limit = 0
    4038        return self
    4139    #}}}
    42 
    4340    def extrude(self, md):  # {{{
    4441        self.coefficient = project3d(md, 'vector', self.coefficient, 'type', 'node', 'layer', 1)
    4542        self.p = project3d(md, 'vector', self.p, 'type', 'element')
    4643        self.q = project3d(md, 'vector', self.q, 'type', 'element')
    47         # If self.coupling == 0:  #doesnt work with empty loop, so just skip it?
    4844        if self.coupling in[3, 4]:
    4945            self.effective_pressure = project3d(md, 'vector', self.effective_pressure, 'type', 'node', 'layer', 1)
    50         elif self.coupling > 4:
    51             raise ValueError('md.friction.coupling larger than 4, not supported yet')
    5246        return self
    5347    #}}}
    54 
    5548    def defaultoutputs(self, md):  # {{{
    5649        list = []
    5750        return list
    5851    #}}}
    59 
    6052    def checkconsistency(self, md, solution, analyses):  # {{{
    6153        # Early return
    … …  
    7163        if self.coupling == 3:
    7264            md = checkfield(md, 'fieldname', 'friction.effective_pressure', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    73         elif self.coupling > 4:
    74             raise ValueError('not supported yet')
    7565        return md
    7666    # }}}
    77 
    7867    def marshall(self, prefix, md, fid):  # {{{
    7968        WriteData(fid, prefix, 'name', 'md.friction.law', 'data', 1, 'format', 'Integer')
    … …  
    8978        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'coupling', 'format', 'Integer')
    9079        WriteData(fid, prefix, 'object', self, 'class', 'friction', 'fieldname', 'effective_pressure_limit', 'format', 'Double')
    91         if self.coupling == 3:
     80        if self.coupling == 3 or self.coupling == 4:
    9281            WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'effective_pressure', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
    93         elif self.coupling == 4:
    94             WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'effective_pressure', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
    95         elif self.coupling > 4:
    96             raise ValueError('not supported yet')
    9782    # }}}

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

    r25514 r26187  
    66classdef frictionschoof
    77        properties (SetAccess=public)
    8                 C    = NaN;
    9                 Cmax = NaN;
    10                 m    = NaN;
     8                C                        = NaN;
     9                Cmax                     = NaN;
     10                m                        = NaN;
     11                coupling                 = 0;
     12                effective_pressure       = NaN;
    1113                effective_pressure_limit = 0;
    1214        end
    … …  
    2628                        self.Cmax = project3d(md,'vector',self.Cmax,'type','node');
    2729                        self.m    = project3d(md,'vector',self.m,'type','element');
     30                        if self.coupling==3 || self.coupling==4
     31                                self.effective_pressure=project3d(md,'vector',self.effective_pressure,'type','node','layer',1);
     32                        end
    2833                end % }}}
    2934                function self = setdefaultparameters(self) % {{{
    3035
     36         self.coupling = 0;
    3137                        self.effective_pressure_limit = 0;
    3238
    … …  
    3642                        %Early return
    3743                        if ~ismember('StressbalanceAnalysis',analyses) & ~ismember('ThermalAnalysis',analyses), return; end
    38                         md = checkfield(md,'fieldname','friction.C','timeseries',1,'NaN',1,'Inf',1,'>',0.);
     44                        md = checkfield(md,'fieldname','friction.C','timeseries',1,'NaN',1,'Inf',1,'>=',0.);
    3945                        md = checkfield(md,'fieldname','friction.Cmax','timeseries',1,'NaN',1,'Inf',1,'>',0.);
    4046                        md = checkfield(md,'fieldname','friction.m','NaN',1,'Inf',1,'>',0.,'size',[md.mesh.numberofelements,1]);
    4147                        md = checkfield(md,'fieldname','friction.effective_pressure_limit','numel',[1],'>=',0);
     48         md = checkfield(md,'fieldname','friction.coupling','numel',[1],'values',[0:4]);
     49         if self.coupling==3
     50            md = checkfield(md,'fieldname','friction.effective_pressure','NaN',1,'Inf',1,'timeseries',1);
     51         end
    4252                end % }}}
    4353                function disp(self) % {{{
    … …  
    5262                        fielddisplay(self,'Cmax','Iken''s bound (typically between 0.17 and 0.84) [SI]');
    5363                        fielddisplay(self,'m','m exponent (generally taken as m = 1/n = 1/3)');
     64                        fielddisplay(self,'effective_pressure','Effective Pressure for the forcing if not coupled [Pa]');
     65         fielddisplay(self,'coupling','Coupling flag 0: uniform sheet (negative pressure ok, default), 1: ice pressure only, 2: water pressure assuming uniform sheet (no negative pressure), 3: use provided effective_pressure, 4: use coupled model (not implemented yet)');
    5466                        fielddisplay(self,'effective_pressure_limit','Neff do not allow to fall below a certain limit: effective_pressure_limit*rho_ice*g*thickness (default 0)');
    5567                end % }}}
    … …  
    6274                        WriteData(fid,prefix,'class','friction','object',self,'fieldname','m','format','DoubleMat','mattype',2);
    6375                        WriteData(fid,prefix,'object',self,'class','friction','fieldname','effective_pressure_limit','format','Double');
     76         WriteData(fid,prefix,'class','friction','object',self,'fieldname','coupling','format','Integer');
     77         if self.coupling==3 || self.coupling==4
     78            WriteData(fid,prefix,'class','friction','object',self,'fieldname','effective_pressure','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
     79         end
    6480                end % }}}
    6581        end

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

    r25806 r26187  
    1616
    1717    def __init__(self, *args):  # {{{
    18         self.C                          = np.nan
    19         self.Cmax                       = np.nan
    20         self.m                          = np.nan
    21         self.effective_pressure_limit   = 0
     18        self.C                        = np.nan
     19        self.Cmax                     = np.nan
     20        self.m                        = np.nan
     21        self.coupling                 = 0
     22        self.effective_pressure       = np.nan
     23        self.effective_pressure_limit = 0
    2224       
    2325        nargs = len(args)
    … …  
    2931            raise Exception('constructor not supported')
    3032    #}}}
    31 
    3233    def __repr__(self):  # {{{
    3334        # See Brondex et al. 2019 https://www.the-cryosphere.net/13/177/2019/
    … …  
    4142        s += "{}\n".format(fielddisplay(self, 'Cmax', 'Iken\'s bound (typically between 0.17 and 0.84) [SI]'))
    4243        s += "{}\n".format(fielddisplay(self, 'm', 'm exponent (generally taken as m = 1/n = 1/3)'))
     44        s += '{}\n'.format(fielddisplay(self, 'coupling', 'Coupling flag 0: uniform sheet (negative pressure ok, default), 1: ice pressure only, 2: water pressure assuming uniform sheet (no negative pressure), 3: use provided effective_pressure, 4: used coupled model (not implemented yet)'))
     45        s += '{}\n'.format(fielddisplay(self, 'effective_pressure', 'Effective Pressure for the forcing if not coupled [Pa]'))
    4346        s += "{}\n".format(fielddisplay(self, 'effective_pressure_limit', 'fNeff do not allow to fall below a certain limit: effective_pressure_limit*rho_ice*g*thickness (default 0)'))
    4447        return s
    4548    #}}}
    46 
    4749    def setdefaultparameters(self):  # {{{
    4850        self.effective_pressure_limit = 0
    4951        return self
    5052    #}}}
    51 
    5253    def extrude(self, md):  # {{{
    5354        self.C = project3d(md, 'vector', self.C, 'type', 'node')
    5455        self.Cmax = project3d(md, 'vector', self.Cmax, 'type', 'node')
    5556        self.m = project3d(md, 'vector', self.m, 'type', 'element')
     57        if self.coupling in[3, 4]:
     58            self.effective_pressure = project3d(md, 'vector', self.effective_pressure, 'type', 'node', 'layer', 1)
    5659        return self
    5760    #}}}
    58 
    5961    def checkconsistency(self, md, solution, analyses):  # {{{
    6062        # Early return
    … …  
    6567        md = checkfield(md, 'fieldname', 'friction.m', 'NaN', 1, 'Inf', 1, '>', 0., 'size', [md.mesh.numberofelements, 1])
    6668        md = checkfield(md, 'fieldname', 'friction.effective_pressure_limit', 'numel', [1], '>=', 0)
     69        md = checkfield(md, 'fieldname', 'friction.coupling', 'numel', [1], 'values', [0, 1, 2, 3, 4])
     70        if self.coupling == 3:
     71            md = checkfield(md, 'fieldname', 'friction.effective_pressure', 'NaN', 1, 'Inf', 1, 'timeseries', 1)
    6772        return md
    6873    # }}}
    69 
    7074    def marshall(self, prefix, md, fid):  # {{{
    7175        yts = md.constants.yts
    … …  
    7579        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'Cmax', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
    7680        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'm', 'format', 'DoubleMat', 'mattype', 2)
     81        WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'coupling', 'format', 'Integer')
    7782        WriteData(fid, prefix, 'object', self, 'class', 'friction', 'fieldname', 'effective_pressure_limit', 'format', 'Double')
     83        if self.coupling == 3 or self.coupling == 4:
     84            WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'effective_pressure', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
    7885    # }}}
    79 
Note: See TracChangeset for help on using the changeset viewer.