Changeset 17364


Ignore:
Timestamp:
02/27/14 16:38:27 (11 years ago)
Author:
jbondzio
Message:

ADD: Calvingrate as model field; LevelsetAnalysis::SetDistanceOnIntersectedElements; MINOR: code cleanup

Location:
issm/trunk-jpl
Files:
1 added
12 edited

Legend:

Unmodified
Added
Removed

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

    r17282 r17364  
    1515/*}}}*/
    1616void LevelsetAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
    17         int    stabilization,finiteelement;
     17        int  finiteelement;
    1818
    1919        /*Finite element type*/
    … …  
    3333        iomodel->FetchDataToInput(elements,VxEnum);
    3434        iomodel->FetchDataToInput(elements,VyEnum);
     35        iomodel->FetchDataToInput(elements,MasstransportCalvingrateEnum);
    3536       
    3637}
    … …  
    9394        /*Intermediaries */
    9495        const int  dim = 2; // solve for LSF in horizontal plane only
     96        int i, row, col;
    9597        IssmDouble kappa;
    9698        IssmDouble Jdet, dt, D_scalar;
    9799        IssmDouble h,hx,hy,hz;
    98         IssmDouble vel,vx,vy,bx,by;
     100        IssmDouble vel,v[dim];
     101        IssmDouble calvingrate, c[dim];
     102        IssmDouble dlsf[dim], norm_dlsf, normal[dim];
    99103        IssmDouble* xyz_list = NULL;
    100104
    … …  
    107111        IssmDouble*    B        = xNew<IssmDouble>(dim*numnodes);
    108112        IssmDouble*    Bprime   = xNew<IssmDouble>(dim*numnodes);
    109         IssmDouble     D[dim][dim], K[dim][dim];
     113        IssmDouble     D[dim][dim];
    110114
    111115        /*Retrieve all inputs and parameters*/
    … …  
    114118        Input* vx_input  = element->GetInput(VxEnum);     _assert_(vx_input);
    115119        Input* vy_input  = element->GetInput(VyEnum);     _assert_(vy_input);
     120        Input* lsf_slopex_input  = element->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
     121        Input* lsf_slopey_input  = element->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
     122        Input* calvingrate_input  = element->GetInput(MasstransportCalvingrateEnum);     _assert_(calvingrate_input);
    116123       
    117124        /* Start  looping on the number of gaussian points: */
    … …  
    136143                GetB(B,element,xyz_list,gauss);
    137144                GetBprime(Bprime,element,xyz_list,gauss);
    138                 vx_input->GetInputValue(&vx,gauss); // in 3D case, add mesh velocity
    139                 vy_input->GetInputValue(&vy,gauss);
    140                 bx=0.; // horizontal mass change velocities (melt/refreeze/calving) FIXME: insert values from model here
    141                 by=0.;
    142                 D[0][0]=D_scalar*(vx+bx);
    143                 D[0][1]=0.;
    144                 D[1][0]=0.;
    145                 D[1][1]=D_scalar*(vy+by);
     145                vx_input->GetInputValue(&v[0],gauss); // in 3D case, add mesh velocity
     146                vy_input->GetInputValue(&v[1],gauss);
     147                lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
     148                lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
     149                calvingrate_input->GetInputValue(&calvingrate,gauss);
     150
     151                norm_dlsf=0.;
     152                for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
     153                norm_dlsf=sqrt(norm_dlsf);
     154
     155                if(norm_dlsf>1.e-10)
     156                        for(i=0;i<dim;i++) c[i]=calvingrate*dlsf[i]/norm_dlsf;
     157                else
     158                        for(i=0;i<dim;i++) c[i]=0.;
     159
     160                for(row=0;row<dim;row++)
     161                        for(col=0;col<dim;col++)
     162                                D[row][col]=((row==col)?D_scalar*(v[row]-c[row]):0.);
     163
    146164                TripleMultiply(B,dim,numnodes,1,
    147165                                        &D[0][0],dim,dim,0,
    … …  
    151169                /* Stabilization */
    152170                int stabilization=2;
     171                vel=0.;
     172                for(i=0;i<dim;i++) vel+=pow(v[i],2);
     173                vel=sqrt(vel)+1.e-14;
    153174                switch(stabilization){
    154175                        case 0:
    … …  
    158179                                /* Artificial Diffusion */
    159180                                element->ElementSizes(&hx,&hy,&hz);
    160                                 vel=sqrt(vx*vx + vy*vy) + 1e-14;
    161                                 h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) ); //FIXME: is this correct?
     181                                h=sqrt( pow(hx*v[0]/vel,2) + pow(hy*v[1]/vel,2) ); //FIXME: is this correct?
    162182
    163183                                kappa=h*vel/2.; //FIXME: insert suitable value for kappa
    164                                 //GetBprime(Bprime,element,xyz_list,gauss); // recalculation of Bprime needed?
    165                                 D[0][0]=D_scalar*kappa;
    166                                 D[0][1]=0.;
    167                                 D[1][0]=0.;
    168                                 D[1][1]=D_scalar*kappa;
     184                                for(row=0;row<dim;row++)
     185                                        for(col=0;col<dim;col++)
     186                                                D[row][col]=((row==col)?D_scalar*kappa:0.);
     187
    169188                                TripleMultiply(Bprime,dim,numnodes,1,
    170189                                                        &D[0][0],dim,dim,0,
    … …  
    175194                                /* Streamline Upwinding */
    176195                                element->ElementSizes(&hx,&hy,&hz);
    177                                 vel=sqrt(vx*vx + vy*vy )+1.e-14;
    178                                 h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) );
    179                                 K[0][0]=h/(2.*vel)*vx*vx;  K[0][1]=h/(2.*vel)*vx*vy;
    180                                 K[1][0]=h/(2.*vel)*vy*vx;  K[1][1]=h/(2.*vel)*vy*vy;
    181                                 for(int i=0;i<dim;i++) for(int j=0;j<dim;j++) K[i][j] = D_scalar*K[i][j];
    182 
    183                                 //GetBprime(Bprime,element,xyz_list,gauss);
     196                                h=sqrt( pow(hx*v[0]/vel,2) + pow(hy*v[1]/vel,2) );
     197                                for(row=0;row<dim;row++)
     198                                        for(col=0;col<dim;col++)
     199                                                D[row][col] = D_scalar*h/(2.*vel)*v[row]*v[col];
     200
    184201                                TripleMultiply(Bprime,dim,numnodes,1,
    185                                                         &K[0][0],dim,dim,0,
     202                                                        &D[0][0],dim,dim,0,
    186203                                                        Bprime,dim,numnodes,0,
    187204                                                        &Ke->values[0],1);
    … …  
    203220
    204221        /*Intermediaries */
    205         int i, ig;
     222        const int dim = 2;
     223        int i, ig, k;
    206224        IssmDouble  Jdet,dt;
    207         IssmDouble  phi;
     225        IssmDouble  lsf;
    208226        IssmDouble* xyz_list = NULL;
    209227       
    … …  
    230248                        element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    231249                        element->NodalFunctions(basis,gauss);
    232                         levelset_input->GetInputValue(&phi,gauss);
    233                         for(i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*phi*basis[i];
     250
     251                        /* old function value */
     252                        levelset_input->GetInputValue(&lsf,gauss);
     253                        for(i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*lsf*basis[i];
    234254                }
    235255
    … …  
    320340
    321341}/*}}}*/
    322 
     342void LevelsetAnalysis::SetDistanceOnIntersectedElements(FemModel* femmodel){/*{{{*/
     343
     344        /* Intermediaries */
     345        int i,k;
     346        IssmDouble dmaxp=0.,dmaxm=0,val=0.;
     347
     348        /*Initialize vector with number of vertices*/
     349        int numvertices=femmodel->vertices->NumberOfVertices();
     350        Element* element;
     351
     352        Vector<IssmDouble>* vec_dist_zerolevelset = NULL;
     353        GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexEnum);
     354       
     355        /* set NaN on elements intersected by zero levelset */
     356        for(i=0;i<femmodel->elements->Size();i++){
     357                element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
     358                if(element->IsZeroLevelset(MaskIceLevelsetEnum))
     359                        for(k=0;k<element->GetNumberOfVertices();k++)
     360                                vec_dist_zerolevelset->SetValue(element->vertices[k]->Sid(),NAN,INS_VAL);
     361        }
     362
     363        /* set distance on elements intersected by zero levelset */
     364        for(i=0;i<femmodel->elements->Size();i++){
     365                element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
     366                if(element->IsZeroLevelset(MaskIceLevelsetEnum)){
     367                        SetDistanceToZeroLevelsetElement(vec_dist_zerolevelset, element);
     368
     369                        /* Get maximum distance to interface along vertices */
     370                        for(k=0;k<element->GetNumberOfVertices();k++){
     371                                        vec_dist_zerolevelset->GetValue(&val,element->vertices[k]->Sid());
     372                                        if((val>0.) && (val>dmaxp))
     373                                                 dmaxp=val;
     374                                        else if((val<0.) && (val<dmaxm))
     375                                                 dmaxm=val;
     376                        }
     377                }
     378        }
     379
     380        /* set all none intersected vertices to max/min distance */
     381        for(i=0;i<numvertices;i++){
     382                vec_dist_zerolevelset->GetValue(&val,i);
     383                if(val==1.) //FIXME: improve check
     384                        vec_dist_zerolevelset->SetValue(i,3.*dmaxp,INS_VAL);
     385                else if(val==-1.)
     386                        vec_dist_zerolevelset->SetValue(i,3.*dmaxm,INS_VAL);
     387        }
     388
     389        /*Assemble vector and serialize */
     390        vec_dist_zerolevelset->Assemble();
     391        IssmDouble* dist_zerolevelset=vec_dist_zerolevelset->ToMPISerial();
     392        InputUpdateFromVectorx(femmodel,dist_zerolevelset,MaskIceLevelsetEnum,VertexSIdEnum);
     393
     394        /*Clean up and return*/
     395        delete vec_dist_zerolevelset;
     396        delete dist_zerolevelset;
     397}/*}}}*/
     398void LevelsetAnalysis::SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_signed_dist, Element* element){/*{{{*/
     399
     400        if(!element->IsZeroLevelset(MaskIceLevelsetEnum))
     401                return;
     402
     403        /* Intermediaries */
     404        const int dim=3;
     405        int i,d;
     406        int numvertices=element->GetNumberOfVertices();
     407        IssmDouble s0[dim], s1[dim], v[dim];
     408        IssmDouble dist,lsf_old;
     409
     410        IssmDouble* lsf = xNew<IssmDouble>(numvertices);
     411        IssmDouble* sign_lsf = xNew<IssmDouble>(numvertices);
     412        IssmDouble* signed_dist = xNew<IssmDouble>(numvertices);
     413        IssmDouble* xyz_list = NULL;
     414        IssmDouble* xyz_list_zero = NULL;
     415
     416        /* retrieve inputs and parameters */
     417        element->GetVerticesCoordinates(&xyz_list);
     418        element->GetInputListOnVertices(lsf,MaskIceLevelsetEnum);
     419
     420        /* get sign of levelset function */
     421        for(i=0;i<numvertices;i++)
     422                sign_lsf[i]=(lsf[i]>=0.?1.:-1.);
     423
     424        element->ZeroLevelsetCoordinates(&xyz_list_zero, xyz_list, MaskIceLevelsetEnum);
     425        for(d=0;d<dim;d++){
     426                s0[d]=xyz_list_zero[0+d];
     427                s1[d]=xyz_list_zero[3+d];
     428        }
     429
     430        /* get signed_distance of vertices to zero levelset straight */
     431        for(i=0;i<numvertices;i++){
     432                for(d=0;d<dim;d++)
     433                        v[d]=xyz_list[3*i+d];
     434                dist=GetDistanceToStraight(&v[0],&s0[0],&s1[0]);
     435                signed_dist[i]=sign_lsf[i]*dist;
     436        }
     437       
     438        /* insert signed_distance into vec_signed_dist, if computed distance is smaller */
     439        for(i=0;i<numvertices;i++){
     440                vec_signed_dist->GetValue(&lsf_old, element->vertices[i]->Sid());
     441                if(xIsNan<IssmDouble>(lsf_old) || fabs(signed_dist[i])<fabs(lsf_old))
     442                        vec_signed_dist->SetValue(element->vertices[i]->Sid(),signed_dist[i],INS_VAL);
     443        }
     444
     445        xDelete<IssmDouble>(lsf);
     446        xDelete<IssmDouble>(sign_lsf);
     447        xDelete<IssmDouble>(signed_dist);
     448        xDelete<IssmDouble>(xyz_list);
     449        xDelete<IssmDouble>(xyz_list_zero);
     450}/*}}}*/
     451IssmDouble LevelsetAnalysis::GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1){/*{{{*/
     452        // returns distance d of point q to straight going through points s0, s1
     453        // d=|a x b|/|b|
     454        // with a=q-s0, b=s1-s0
     455       
     456        /* Intermediaries */
     457        const int dim=2;
     458        int i;
     459        IssmDouble a[dim], b[dim];
     460        IssmDouble norm_b;
     461
     462        for(i=0;i<dim;i++){
     463                a[i]=q[i]-s0[i];
     464                b[i]=s1[i]-s0[i];
     465        }
     466       
     467        norm_b=0.;
     468        for(i=0;i<dim;i++)
     469                norm_b+=b[i]*b[i];
     470        norm_b=sqrt(norm_b);
     471        _assert_(norm_b>0.);
     472
     473        return fabs(a[0]*b[1]-a[1]*b[0])/norm_b;
     474}/*}}}*/
     475

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

    r17230 r17364  
    1111class LevelsetAnalysis: public Analysis{
    1212       
    13  public:
     13public:
    1414        /*Model processing*/
    1515        int  DofsPerNode(int** doflist,int meshtype,int approximation);
    … …  
    2727        ElementVector* CreatePVector(Element* element);
    2828        void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    29                 void InputUpdateFromSolution(IssmDouble* solution,Element* element);
    30                 void UpdateConstraints(FemModel* femmodel);
     29        void InputUpdateFromSolution(IssmDouble* solution,Element* element);
     30        void UpdateConstraints(FemModel* femmodel);
    3131        void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    3232        void GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    33 
     33        void SetDistanceOnIntersectedElements(FemModel* femmodel);
     34        void SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_signed_dist, Element* element);
     35        IssmDouble GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1);
    3436};
    3537#endif

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r17362 r17364  
    463463                                name==LevelsetfunctionSlopeYEnum ||
    464464                                name==LevelsetfunctionPicardEnum ||
     465                                name==MasstransportCalvingrateEnum ||
    465466                                name==GradientEnum ||
    466467                                name==OldGradientEnum  ||

  • issm/trunk-jpl/src/c/cores/transient_core.cpp

    r17274 r17364  
    9393        }
    9494
     95        if(islevelset){
     96                /* set distance on elements intersected by zero level set */
     97                LevelsetAnalysis* lsfanalysis = new LevelsetAnalysis();
     98                lsfanalysis->SetDistanceOnIntersectedElements(femmodel);
     99                delete lsfanalysis;
     100        }
     101       
    95102        while(time < finaltime - (yts*DBL_EPSILON)){ //make sure we run up to finaltime.
    96103
    … …  
    125132                if(islevelset){
    126133                        if(VerboseSolution()) _printf0_("   computing movement of ice boundaries\n");
    127                        
    128                         /* get slope of lsf for computation of normal on ice domain*/
     134                        /* smoothen slope of lsf for computation of normal on ice domain*/
    129135                        levelsetfunctionslope_core(femmodel);
    130136
    131                         /* extrapolate along normal */
     137                        /* extrapolate velocities onto domain with no ice */
    132138                        Analysis* extanalysis = new ExtrapolationAnalysis();
    133139                        const int nvars=2;
    … …  
    139145                        delete extanalysis;     
    140146
     147                        /* solve level set equation */
    141148                        analysis = new LevelsetAnalysis();
    142149                        analysis->Core(femmodel);
    … …  
    149156                        int outputs[1] = {IceMaskNodeActivationEnum};
    150157                        femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1);
    151 
    152158                }
    153159

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

    r17353 r17364  
    217217        MasstransportPenaltyFactorEnum,
    218218        MasstransportSpcthicknessEnum,
     219        MasstransportCalvingrateEnum,
    219220        MasstransportStabilizationEnum,
    220221        MasstransportVertexPairingEnum,

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

    r17353 r17364  
    225225                case MasstransportPenaltyFactorEnum : return "MasstransportPenaltyFactor";
    226226                case MasstransportSpcthicknessEnum : return "MasstransportSpcthickness";
     227                case MasstransportCalvingrateEnum : return "MasstransportCalvingrate";
    227228                case MasstransportStabilizationEnum : return "MasstransportStabilization";
    228229                case MasstransportVertexPairingEnum : return "MasstransportVertexPairing";

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

    r17353 r17364  
    228228              else if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum;
    229229              else if (strcmp(name,"MasstransportSpcthickness")==0) return MasstransportSpcthicknessEnum;
     230              else if (strcmp(name,"MasstransportCalvingrate")==0) return MasstransportCalvingrateEnum;
    230231              else if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum;
    231232              else if (strcmp(name,"MasstransportVertexPairing")==0) return MasstransportVertexPairingEnum;
    … …  
    259260              else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum;
    260261              else if (strcmp(name,"Surface")==0) return SurfaceEnum;
    261               else if (strcmp(name,"ThermalIsenthalpy")==0) return ThermalIsenthalpyEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"ThermalIsdynamicbasalspc")==0) return ThermalIsdynamicbasalspcEnum;
     265              if (strcmp(name,"ThermalIsenthalpy")==0) return ThermalIsenthalpyEnum;
     266              else if (strcmp(name,"ThermalIsdynamicbasalspc")==0) return ThermalIsdynamicbasalspcEnum;
    266267              else if (strcmp(name,"ThermalMaxiter")==0) return ThermalMaxiterEnum;
    267268              else if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum;
    … …  
    382383              else if (strcmp(name,"Results")==0) return ResultsEnum;
    383384              else if (strcmp(name,"GenericParam")==0) return GenericParamEnum;
    384               else if (strcmp(name,"AdolcParam")==0) return AdolcParamEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
     388              if (strcmp(name,"AdolcParam")==0) return AdolcParamEnum;
     389              else if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
    389390              else if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
    390391              else if (strcmp(name,"Contour")==0) return ContourEnum;
    … …  
    505506              else if (strcmp(name,"VzMesh")==0) return VzMeshEnum;
    506507              else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum;
    507               else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum;
     511              if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum;
     512              else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum;
    512513              else if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum;
    513514              else if (strcmp(name,"ThicknessAcrossGradient")==0) return ThicknessAcrossGradientEnum;
    … …  
    628629              else if (strcmp(name,"Sset")==0) return SsetEnum;
    629630              else if (strcmp(name,"Verbose")==0) return VerboseEnum;
    630               else if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
     634              if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum;
     635              else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
    635636              else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum;
    636637              else if (strcmp(name,"XY")==0) return XYEnum;

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

    r16764 r17364  
    77        properties (SetAccess=public)
    88                 spcthickness           = NaN;
     9                 calvingrate            = NaN;
    910                 isfreesurface          = 0;
    1011                 min_thickness          = 0;
    … …  
    8889
    8990                        md = checkfield(md,'fieldname','masstransport.spcthickness','forcing',1);
     91                        if(ismember(LevelsetAnalysisEnum(), analyses) & md.transient.islevelset)
     92                                md = checkfield(md,'fieldname','masstransport.calvingrate','NaN',1,'size',[md.mesh.numberofvertices 1],'>=',0);
     93                        end
     94
    9095                        md = checkfield(md,'fieldname','masstransport.isfreesurface','values',[0 1]);
    9196                        md = checkfield(md,'fieldname','masstransport.hydrostatic_adjustment','values',{'Absolute' 'Incremental'});
    … …  
    98103                        disp(sprintf('   Masstransport solution parameters:'));
    99104                        fielddisplay(obj,'spcthickness','thickness constraints (NaN means no constraint) [m]');
     105                        fielddisplay(obj,'calvingrate','calving rate at given location [m/a]');
    100106                        fielddisplay(obj,'isfreesurface','do we use free surfaces (FS only) are mass conservation');
    101107                        fielddisplay(obj,'min_thickness','minimum ice thickness allowed [m]');
    … …  
    110116                end % }}}
    111117                function marshall(obj,md,fid) % {{{
     118
     119                        yts=365.*24.*3600.;
     120
    112121                        WriteData(fid,'object',obj,'fieldname','spcthickness','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);
     122                        WriteData(fid,'object',obj,'fieldname','calvingrate','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1,'scale',1./yts);
    113123                        WriteData(fid,'object',obj,'fieldname','isfreesurface','format','Boolean');
    114124                        WriteData(fid,'object',obj,'fieldname','min_thickness','format','Double');

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

    r16764 r17364  
    1515        def __init__(self): # {{{
    1616                self.spcthickness           = float('NaN')
     17                self.calvingrate            = float('NaN')
    1718                self.isfreesurface          = 0
    1819                self.min_thickness          = 0
    … …  
    3031                string='   Masstransport solution parameters:'
    3132                string="%s\n%s"%(string,fielddisplay(self,'spcthickness','thickness constraints (NaN means no constraint) [m]'))
     33                string="%s\n%s"%(string,fielddisplay(self,'calvingrate','calving rate at given location [m/a]'))
    3234                string="%s\n%s"%(string,fielddisplay(self,'isfreesurface','do we use free surfaces (FS only) are mass conservation'))
    3335                string="%s\n%s"%(string,fielddisplay(self,'min_thickness','minimum ice thickness allowed [m]'))
    … …  
    6870
    6971                md = checkfield(md,'fieldname','masstransport.spcthickness','forcing',1)
     72                if LevelsetAnalysisEnum() in analyses and md.transient.islevelset:
     73                        md = checkfield(md,'fieldname','masstransport.calvingrate','NaN',1,'size',[md.mesh.numberofvertices],'>=',0)
    7074                md = checkfield(md,'fieldname','masstransport.isfreesurface','values',[0,1])
    7175                md = checkfield(md,'fieldname','masstransport.hydrostatic_adjustment','values',['Absolute','Incremental'])
    … …  
    7781        # }}}
    7882        def marshall(self,md,fid):    # {{{
     83
     84                yts=365.*24.*3600.
     85
    7986                WriteData(fid,'object',self,'fieldname','spcthickness','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1)
     87                WriteData(fid,'object',self,'fieldname','calvingrate','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1,'scale',1./yts)
    8088                WriteData(fid,'object',self,'fieldname','isfreesurface','format','Boolean')
    8189                WriteData(fid,'object',self,'fieldname','min_thickness','format','Double')

  • issm/trunk-jpl/src/m/enum/EnumDefinitions.py

    r17352 r17364  
    217217def MasstransportPenaltyFactorEnum(): return StringToEnum("MasstransportPenaltyFactor")[0]
    218218def MasstransportSpcthicknessEnum(): return StringToEnum("MasstransportSpcthickness")[0]
     219def MasstransportCalvingrateEnum(): return StringToEnum("MasstransportCalvingrate")[0]
    219220def MasstransportStabilizationEnum(): return StringToEnum("MasstransportStabilization")[0]
    220221def MasstransportVertexPairingEnum(): return StringToEnum("MasstransportVertexPairing")[0]

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

    r17251 r17364  
    3535R=0.2*(xmax-xmin)
    3636md.mask.ice_levelset=D-R
     37md.masstransport.calvingrate=0.*numpy.ones((md.mesh.numberofvertices,1))
    3738
    3839md=solve(md,TransientSolutionEnum())
Note: See TracChangeset for help on using the changeset viewer.