Changeset 18816


Ignore:
Timestamp:
11/20/14 11:49:30 (10 years ago)
Author:
srebuffi
Message:

NEW : added melting rate

Location:
issm/trunk-jpl
Files:
16 edited

Legend:

Unmodified
Added
Removed

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

    r18808 r18816  
    4848                        case DefaultCalvingEnum:
    4949                                iomodel->FetchDataToInput(elements,CalvingCalvingrateEnum);
     50                                iomodel->FetchDataToInput(elements,CalvingMeltingrateEnum);
    5051                                break;
    5152                        case CalvingLevermannEnum:
    5253                                iomodel->FetchDataToInput(elements,CalvinglevermannCoeffEnum);
     54                                iomodel->FetchDataToInput(elements,CalvinglevermannMeltingrateEnum);
    5355                                break;
    5456                        default:
    … …  
    125127        IssmDouble h,hx,hy,hz;
    126128        IssmDouble vel;
    127         IssmDouble norm_dlsf, calvingrate;
     129        IssmDouble norm_dlsf, norm_calving, calvingrate, meltingrate;
    128130        IssmDouble* xyz_list = NULL;
    129131
    … …  
    151153        IssmDouble*    w        = xNew<IssmDouble>(dim);
    152154        IssmDouble*    c        = xNewZeroInit<IssmDouble>(dim);
     155        IssmDouble*    m        = xNewZeroInit<IssmDouble>(dim);
    153156        IssmDouble*    dlsf     = xNew<IssmDouble>(dim);
    154157
    … …  
    163166        Input* lsf_slopey_input   = NULL;
    164167        Input* calvingrate_input  = NULL;
     168        Input* meltingrate_input  = NULL;
    165169
    166170        /*Load velocities*/
    … …  
    186190                                if(dim==2) lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
    187191                                calvingrate_input = basalelement->GetInput(CalvingCalvingrateEnum);     _assert_(calvingrate_input);
     192                                meltingrate_input = basalelement->GetInput(CalvingMeltingrateEnum);     _assert_(meltingrate_input);
    188193                                break;
    189194                        case CalvingLevermannEnum:
    … …  
    201206                                        }
    202207                                }
     208                                meltingrate_input = basalelement->GetInput(CalvinglevermannMeltingrateEnum);     _assert_(meltingrate_input);
    203209                                break;
    204210                        default:
    … …  
    238244                                        if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
    239245                                        calvingrate_input->GetInputValue(&calvingrate,gauss);
     246                                        meltingrate_input->GetInputValue(&meltingrate,gauss);
    240247
    241248                                        norm_dlsf=0.;
    … …  
    244251
    245252                                        if(norm_dlsf>1.e-10)
    246                                          for(i=0;i<dim;i++) c[i]=calvingrate*dlsf[i]/norm_dlsf;
     253                                         for(i=0;i<dim;i++){
     254                                          c[i]=calvingrate*dlsf[i]/norm_dlsf;
     255                                          m[i]=meltingrate*dlsf[i]/norm_dlsf;
     256                                         }
    247257                                        else
    248                                          for(i=0;i<dim;i++) c[i]=0.;
     258                                         for(i=0;i<dim;i++){
     259                                                 c[i]=0.;
     260                                                 m[i]=0.;
     261                                         }
     262
    249263                                        break;
    250264                                case CalvingLevermannEnum:
    251265                                        calvingratex_input->GetInputValue(&c[0],gauss); // in 3D case, add mesh velocity
    252266                                        if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss);
     267                                        meltingrate_input->GetInputValue(&meltingrate,gauss);
     268
     269                                        norm_calving=0.;
     270                                        for(i=0;i<dim;i++) norm_calving+=pow(c[i],2);
     271                                        norm_calving=sqrt(norm_calving)+1.e-14;
     272
     273                                        for(i=0;i<dim;i++) m[i]=meltingrate*c[i]/norm_calving;
     274                                       
    253275                                        break;
    254276                                default:
    … …  
    258280
    259281                /*Levelset speed is ice velocity - calving rate*/
    260                 for(i=0;i<dim;i++) w[i]=v[i]-c[i];
     282                for(i=0;i<dim;i++) w[i]=v[i]-c[i]-m[i];
    261283
    262284                /*Compute D*/

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r18812 r18816  
    213213        CalvingLawEnum,
    214214        CalvingCalvingrateEnum,
     215        CalvingMeltingrateEnum,
    215216        CalvingLevermannEnum,
    216217        DefaultCalvingEnum,
    217218        CalvingRequestedOutputsEnum,
    218219        CalvinglevermannCoeffEnum,
     220        CalvinglevermannMeltingrateEnum,
    219221        CalvingratexEnum,
    220222        CalvingrateyEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r18812 r18816  
    221221                case CalvingLawEnum : return "CalvingLaw";
    222222                case CalvingCalvingrateEnum : return "CalvingCalvingrate";
     223                case CalvingMeltingrateEnum : return "CalvingMeltingrate";
    223224                case CalvingLevermannEnum : return "CalvingLevermann";
    224225                case DefaultCalvingEnum : return "DefaultCalving";
    225226                case CalvingRequestedOutputsEnum : return "CalvingRequestedOutputs";
    226227                case CalvinglevermannCoeffEnum : return "CalvinglevermannCoeff";
     228                case CalvinglevermannMeltingrateEnum : return "CalvinglevermannMeltingrate";
    227229                case CalvingratexEnum : return "Calvingratex";
    228230                case CalvingrateyEnum : return "Calvingratey";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r18812 r18816  
    224224              else if (strcmp(name,"CalvingLaw")==0) return CalvingLawEnum;
    225225              else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
     226              else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
    226227              else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum;
    227228              else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum;
    228229              else if (strcmp(name,"CalvingRequestedOutputs")==0) return CalvingRequestedOutputsEnum;
    229230              else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
     231              else if (strcmp(name,"CalvinglevermannMeltingrate")==0) return CalvinglevermannMeltingrateEnum;
    230232              else if (strcmp(name,"Calvingratex")==0) return CalvingratexEnum;
    231233              else if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum;
    … …  
    258260              else if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum;
    259261              else if (strcmp(name,"MeshVertexonboundary")==0) return MeshVertexonboundaryEnum;
    260               else if (strcmp(name,"MeshX")==0) return MeshXEnum;
    261               else if (strcmp(name,"MeshY")==0) return MeshYEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"MeshZ")==0) return MeshZEnum;
     265              if (strcmp(name,"MeshX")==0) return MeshXEnum;
     266              else if (strcmp(name,"MeshY")==0) return MeshYEnum;
     267              else if (strcmp(name,"MeshZ")==0) return MeshZEnum;
    266268              else if (strcmp(name,"MeshElementtype")==0) return MeshElementtypeEnum;
    267269              else if (strcmp(name,"DomainType")==0) return DomainTypeEnum;
    … …  
    381383              else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum;
    382384              else if (strcmp(name,"AnalysisType")==0) return AnalysisTypeEnum;
    383               else if (strcmp(name,"ConfigurationType")==0) return ConfigurationTypeEnum;
    384               else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum;
     388              if (strcmp(name,"ConfigurationType")==0) return ConfigurationTypeEnum;
     389              else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
     390              else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum;
    389391              else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
    390392              else if (strcmp(name,"AnalysisCounter")==0) return AnalysisCounterEnum;
    … …  
    504506              else if (strcmp(name,"MassconaxpbyNamex")==0) return MassconaxpbyNamexEnum;
    505507              else if (strcmp(name,"MassconaxpbyNamey")==0) return MassconaxpbyNameyEnum;
    506               else if (strcmp(name,"MassconaxpbyAlpha")==0) return MassconaxpbyAlphaEnum;
    507               else if (strcmp(name,"MassconaxpbyBeta")==0) return MassconaxpbyBetaEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
     511              if (strcmp(name,"MassconaxpbyAlpha")==0) return MassconaxpbyAlphaEnum;
     512              else if (strcmp(name,"MassconaxpbyBeta")==0) return MassconaxpbyBetaEnum;
     513              else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
    512514              else if (strcmp(name,"VectorParam")==0) return VectorParamEnum;
    513515              else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
    … …  
    627629              else if (strcmp(name,"StrainRateyz")==0) return StrainRateyzEnum;
    628630              else if (strcmp(name,"StrainRatezz")==0) return StrainRatezzEnum;
    629               else if (strcmp(name,"Divergence")==0) return DivergenceEnum;
    630               else if (strcmp(name,"MaxDivergence")==0) return MaxDivergenceEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"GiaCrossSectionShape")==0) return GiaCrossSectionShapeEnum;
     634              if (strcmp(name,"Divergence")==0) return DivergenceEnum;
     635              else if (strcmp(name,"MaxDivergence")==0) return MaxDivergenceEnum;
     636              else if (strcmp(name,"GiaCrossSectionShape")==0) return GiaCrossSectionShapeEnum;
    635637              else if (strcmp(name,"GiadWdt")==0) return GiadWdtEnum;
    636638              else if (strcmp(name,"GiaW")==0) return GiaWEnum;
    … …  
    750752              else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum;
    751753              else if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum;
    752               else if (strcmp(name,"Seq")==0) return SeqEnum;
    753               else if (strcmp(name,"Mpi")==0) return MpiEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"Mumps")==0) return MumpsEnum;
     757              if (strcmp(name,"Seq")==0) return SeqEnum;
     758              else if (strcmp(name,"Mpi")==0) return MpiEnum;
     759              else if (strcmp(name,"Mumps")==0) return MumpsEnum;
    758760              else if (strcmp(name,"Gsl")==0) return GslEnum;
    759761              else if (strcmp(name,"Option")==0) return OptionEnum;

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

    r18757 r18816  
    77        properties (SetAccess=public)
    88                 calvingrate            = NaN;
     9                 meltingrate            = NaN;
    910        end
    1011        methods (Static)
    … …  
    6061                        if (solution~=TransientSolutionEnum() | md.transient.iscalving==0), return; end
    6162                        md = checkfield(md,'fieldname','calving.calvingrate','NaN',1,'size',[md.mesh.numberofvertices 1],'>=',0);
     63                        md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'size',[md.mesh.numberofvertices 1],'>=',0);
    6264                end % }}}
    6365                function disp(obj) % {{{
    6466                        disp(sprintf('   Calving parameters:'));
    6567                        fielddisplay(obj,'calvingrate','calving rate at given location [m/a]');
     68                        fielddisplay(obj,'meltingrate','melting rate at given location [m/a]');
    6669                end % }}}
    6770                function marshall(obj,md,fid) % {{{
    … …  
    6972                        WriteData(fid,'enum',CalvingLawEnum(),'data',DefaultCalvingEnum(),'format','Integer');
    7073                        WriteData(fid,'object',obj,'fieldname','calvingrate','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1,'scale',1./yts);
     74                        WriteData(fid,'object',obj,'fieldname','meltingrate','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1,'scale',1./yts);
    7175                end % }}}
    7276        end

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

    r18766 r18816  
    1515        def __init__(self): # {{{
    1616                self.calvingrate            = float('NaN')
     17                self.meltingrate            = float('NaN')
    1718
    1819                #set defaults
    … …  
    2324                string='   Calving parameters:'
    2425                string="%s\n%s"%(string,fielddisplay(self,'calvingrate','calving rate at given location [m/a]'))
     26                string="%s\n%s"%(string,fielddisplay(self,'meltingrate','melting rate at given location [m/a]'))
    2527
    2628                return string
    … …  
    3739
    3840                md = checkfield(md,'fieldname','calving.calvingrate','NaN',1,'size',[md.mesh.numberofvertices],'>=',0)
     41                md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'size',[md.mesh.numberofvertices],'>=',0)
    3942                return md
    4043        # }}}
    … …  
    4548                WriteData(fid,'enum',CalvingLawEnum(),'data',DefaultCalvingEnum(),'format','Integer');
    4649                WriteData(fid,'object',self,'fieldname','calvingrate','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1,'scale',1./yts)
     50                WriteData(fid,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1,'scale',1./yts)
    4751        # }}}

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

    r18766 r18816  
    66classdef calvinglevermann
    77        properties (SetAccess=public)
    8                  coeff = NaN;
     8                 coeff       = NaN;
     9                 meltingrate = NaN;
    910        end
    1011        methods
    … …  
    3536                        if (solution~=TransientSolutionEnum() | md.transient.iscalving==0), return; end
    3637                        md = checkfield(md,'fieldname','calving.coeff','>',0,'size',[md.mesh.numberofvertices 1]);
     38                        md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'size',[md.mesh.numberofvertices 1],'>=',0);
    3739                end % }}}
    3840                function disp(obj) % {{{
    3941                        disp(sprintf('   Calving Levermann parameters:'));
    4042                        fielddisplay(obj,'coeff','proportionality coefficient in Levermann model');
     43                        fielddisplay(obj,'meltingrate','melting rate at given location [m/a]');
    4144
    4245                end % }}}
    4346                function marshall(obj,md,fid) % {{{
     47                        yts=365.0*24.0*3600.0;
    4448                        WriteData(fid,'enum',CalvingLawEnum(),'data',CalvingLevermannEnum(),'format','Integer');
    4549                        WriteData(fid,'enum',CalvinglevermannCoeffEnum(),'data',obj.coeff,'format','DoubleMat','mattype',1);
     50                        WriteData(fid,'object',obj,'fieldname','meltingrate','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1,'scale',1./yts);
    4651                end % }}}
    4752        end

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

    r18792 r18816  
    1414
    1515        def __init__(self): # {{{
    16                 self.coeff = float('NaN')
     16                self.coeff       = float('NaN')
     17                self.meltingrate = float('NaN')
    1718
    1819                #set defaults
    … …  
    2324                string='   Calving Levermann parameters:'
    2425                string="%s\n%s"%(string,fielddisplay(self,'coeff','proportionality coefficient in Levermann model'))
     26                string="%s\n%s"%(string,fielddisplay(self,'meltingrate','melting rate at given location [m/a]'))
    2527
    2628                return string
    … …  
    3840
    3941                md = checkfield(md,'fieldname','calving.coeff','size',[md.mesh.numberofvertices],'>',0)
     42                md = checkfield(md,'fieldname','calving.meltingrate','NaN',1,'size',[md.mesh.numberofvertices],'>=',0)
    4043                return md
    4144        # }}}
    4245        def marshall(self,md,fid):    # {{{
     46                yts=365.*24.*3600.
    4347                WriteData(fid,'enum',CalvingLawEnum(),'data',CalvingLevermannEnum(),'format','Integer');
    4448                WriteData(fid,'enum',CalvinglevermannCoeffEnum(),'data',self.coeff,'format','DoubleMat','mattype',1)
     49                WriteData(fid,'object',self,'fieldname','meltingrate','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1,'scale',1./yts)
    4550        # }}}

  • issm/trunk-jpl/test/NightlyRun/IdToName.m

    r18782 r18816  
    195195        case 804, name='ValleyGlacierLevelsetCalvingSIA2d';
    196196        case 805, name='ValleyGlacierLevelsetEnthCalvingHO3d';
    197         case 806, name='SquareShelfLevelsetCalvingSIA2dLevermann';
     197        case 806, name='SquareShelfLevelsetCalvingSSA2dLevermann';
     198        case 807, name='SquareShelfLevelsetMeltingSSA2d';
    198199        case 1101, name='ISMIPAHO';
    199200        case 1102, name='ISMIPAFS';

  • issm/trunk-jpl/test/NightlyRun/IdToName.py

    r18782 r18816  
    180180        804  : 'ValleyGlacierLevelsetCalvingSIA2d',     
    181181        805  : 'ValleyGlacierLevelsetEnthCalvingHO3d', 
    182         806  : 'SquareShelfLevelsetCalvingSIA2dLevermann',     
     182        806  : 'SquareShelfLevelsetCalvingSSA2dLevermann',     
     183        807  : 'SquareShelfLevelsetMeltingSSA2d',
    183184        1101 : 'ISMIPAHO',
    184185        1102 : 'ISMIPAFS',

  • issm/trunk-jpl/test/NightlyRun/test804.m

    r18760 r18816  
    1515
    1616md.calving.calvingrate=1000.*ones(md.mesh.numberofvertices,1);
     17md.calving.meltingrate=zeros(md.mesh.numberofvertices,1);
    1718
    1819md=solve(md,TransientSolutionEnum());

  • issm/trunk-jpl/test/NightlyRun/test804.py

    r18762 r18816  
    2525
    2626md.calving.calvingrate=1000.*numpy.ones((md.mesh.numberofvertices,1))
     27md.calving.meltingrate=numpy.zeros((md.mesh.numberofvertices,1))
    2728
    2829md=solve(md,TransientSolutionEnum())

  • issm/trunk-jpl/test/NightlyRun/test805.m

    r18760 r18816  
    2222
    2323md.calving.calvingrate=1000.*ones(md.mesh.numberofvertices,1);
     24md.calving.meltingrate=zeros(md.mesh.numberofvertices,1);
    2425
    2526md=solve(md,TransientSolutionEnum());

  • issm/trunk-jpl/test/NightlyRun/test805.py

    r18762 r18816  
    3232
    3333md.calving.calvingrate=1000.*numpy.ones((md.mesh.numberofvertices,1))
     34md.calving.meltingrate=numpy.zeros((md.mesh.numberofvertices,1))
    3435
    3536md=solve(md,TransientSolutionEnum())

  • issm/trunk-jpl/test/NightlyRun/test806.m

    r18782 r18816  
    2626md.calving=calvinglevermann();
    2727md.calving.coeff=4.89e13*ones(md.mesh.numberofvertices,1);
     28md.calving.meltingrate=zeros(md.mesh.numberofvertices,1);
    2829
    2930md.transient.requested_outputs={'default','StrainRateparallel','StrainRateperpendicular','Calvingratex','Calvingratey','CalvingCalvingrate'};

  • issm/trunk-jpl/test/NightlyRun/test806.py

    r18793 r18816  
    3737md.transient.iscalving=True;
    3838
    39 md.calving=calvinglevermann();
    40 md.calving.coeff=4.89e13*numpy.ones((md.mesh.numberofvertices,1));
     39md.calving=calvinglevermann()
     40md.calving.coeff=4.89e13*numpy.ones((md.mesh.numberofvertices,1))
     41md.calving.meltingrate=numpy.zeros((md.mesh.numberofvertices,1))
    4142
    4243md.transient.requested_outputs=['default','StrainRateparallel','StrainRateperpendicular','Calvingratex','Calvingratey','CalvingCalvingrate']
Note: See TracChangeset for help on using the changeset viewer.