Changeset 26894


Ignore:
Timestamp:
02/18/22 02:55:33 (3 years ago)
Author:
rueckamp
Message:

CHG: Added SUPG stabilization and min thickness consraint for Free Surface Analysis

Location:
issm/trunk-jpl/src/c
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r26825 r26894  
    7373        }
    7474
    75         iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum);
     75        iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum);
    7676        iomodel->FetchDataToInput(inputs,elements,"md.initialization.sealevel",SealevelEnum,0);
    7777        iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
    … …  
    9696                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.floatingice_melting_rate",BasalforcingsFloatingiceMeltingRateEnum);
    9797                        if(isstochastic){
    98             iomodel->FetchDataToInput(inputs,elements,"md.stochasticforcing.default_id",StochasticForcingDefaultIdEnum);
    99             iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.floatingice_melting_rate",BaselineBasalforcingsFloatingiceMeltingRateEnum);
    100          }
     98                                iomodel->FetchDataToInput(inputs,elements,"md.stochasticforcing.default_id",StochasticForcingDefaultIdEnum);
     99                                iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.floatingice_melting_rate",BaselineBasalforcingsFloatingiceMeltingRateEnum);
     100                        }
    101101                        break;
    102102                case LinearFloatingMeltRateEnum:
    … …  
    113113                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.upperwater_elevation",BasalforcingsSpatialUpperwaterElevationEnum);
    114114                        if(isstochastic){
    115             iomodel->FetchDataToInput(inputs,elements,"md.stochasticforcing.default_id",StochasticForcingDefaultIdEnum);
    116             iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_melting_rate",BaselineBasalforcingsSpatialDeepwaterMeltingRateEnum);
    117          }
     115                                iomodel->FetchDataToInput(inputs,elements,"md.stochasticforcing.default_id",StochasticForcingDefaultIdEnum);
     116                                iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.deepwater_melting_rate",BaselineBasalforcingsSpatialDeepwaterMeltingRateEnum);
     117                        }
    118118                        break;
    119119                case BasalforcingsPicoEnum:
    … …  
    148148}/*}}}*/
    149149ElementMatrix* FreeSurfaceBaseAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
    150 _error_("Not implemented");
     150        _error_("Not implemented");
    151151}/*}}}*/
    152152ElementMatrix* FreeSurfaceBaseAnalysis::CreateKMatrix(Element* element){/*{{{*/
    … …  
    157157        IssmDouble *xyz_list  = NULL;
    158158        IssmDouble  Jdet,D_scalar,dt,h;
    159         IssmDouble  vel,vx,vy;
     159        IssmDouble  vel,vx,vy,tau;
    160160
    161161        /*Get basal element*/
    … …  
    230230                }
    231231
    232                 if(stabilization==2){
    233                         /*Streamline upwinding*/
    234                         if(dim==1){
    235                          vel=fabs(vx)+1.e-8;
    236                          D[0] = h/(2.*vel)*vx*vx;
    237                         }
    238                         else{
    239                          vel=sqrt(vx*vx+vy*vy)+1.e-8;
    240                          D[0*dim+0]=h/(2*vel)*vx*vx;
    241                          D[1*dim+0]=h/(2*vel)*vy*vx;
    242                          D[0*dim+1]=h/(2*vel)*vx*vy;
    243                          D[1*dim+1]=h/(2*vel)*vy*vy;
    244                         }
    245                 }
    246                 else if(stabilization==1){
     232                if(stabilization==1){
    247233                        /*SSA*/
    248234                        if(dim==1){
    … …  
    255241                                D[0*dim+0]=h/2.0*fabs(vx);
    256242                                D[1*dim+1]=h/2.0*fabs(vy);
     243                        }
     244                }
     245                else if(stabilization==2){
     246                        /*Streamline upwinding*/
     247                        if(dim==1){
     248                                vel=fabs(vx)+1.e-8;
     249                                D[0] = h/(2.*vel)*vx*vx;
     250                        }
     251                        else{
     252                                vel=sqrt(vx*vx+vy*vy)+1.e-8;
     253                                D[0*dim+0]=h/(2*vel)*vx*vx;
     254                                D[1*dim+0]=h/(2*vel)*vy*vx;
     255                                D[0*dim+1]=h/(2*vel)*vx*vy;
     256                                D[1*dim+1]=h/(2*vel)*vy*vy;
     257                        }
     258
     259                }
     260                else if(stabilization==5){
     261                        /*SUPG*/
     262                        if(dim==1){
     263                                vx_input->GetInputAverage(&vx);
     264                                tau=h/(2.*fabs(vx)+1e-10);
     265                        }
     266                        else{
     267                                vx_input->GetInputAverage(&vx);
     268                                vy_input->GetInputAverage(&vy);
     269                                tau=1*h/(2.*pow(vx*vx+vy*vy,0.5)+1e-10);
    257270                        }
    258271                }
    … …  
    263276                                        for(int j=0;j<numnodes;j++){
    264277                                                Ke->values[i*numnodes+j] += (
    265                                                                         dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) +
    266                                                                         dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j])
    267                                                                         );
     278                                                                dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) +
     279                                                                dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j])
     280                                                                );
    268281                                        }
    269282                                }
    … …  
    271284                        else{
    272285                                for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += dbasis[0*numnodes+i]*D[0]*dbasis[0*numnodes+j];
     286                        }
     287                }
     288                else if(stabilization==5){
     289                        D_scalar=gauss->weight*Jdet*dt;
     290                        if(dim==2){
     291                                for(int i=0;i<numnodes;i++){
     292                                        for(int j=0;j<numnodes;j++){
     293                                                Ke->values[i*numnodes+j]+=tau*D_scalar*
     294                                                        (vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i])*
     295                                                        (vx*dbasis[0*numnodes+j]+vy*dbasis[1*numnodes+j]);
     296                                        }
     297                                }
     298                        }
     299                        else{
     300                                for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j]+=tau*D_scalar*(vx*dbasis[0*numnodes+i])*(vx*dbasis[0*numnodes+j]);
    273301                        }
    274302                }
    … …  
    285313}/*}}}*/
    286314ElementVector* FreeSurfaceBaseAnalysis::CreatePVector(Element* element){/*{{{*/
     315
    287316        /*Intermediaries*/
    288         int         domaintype,dim;
     317        int         domaintype,dim,stabilization;
    289318        IssmDouble  Jdet,dt;
    290         IssmDouble  gmb,fmb,mb,bed,phi,vz;
     319        IssmDouble  gmb,fmb,mb,bed,vx,vy,vz,tau;
    291320        Element*    basalelement = NULL;
    292321        IssmDouble *xyz_list  = NULL;
    … …  
    314343        /*Fetch number of nodes and dof for this finite element*/
    315344        int numnodes = basalelement->GetNumberOfNodes();
     345        int         melt_style,point1;
     346        IssmDouble  fraction1,fraction2;
     347        bool        mainlyfloating;
    316348
    317349        /*Initialize Element vector and other vectors*/
    318350        ElementVector* pe    = basalelement->NewElementVector();
    319351        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     352        IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
     353        IssmDouble  gllevelset,phi=1.;
    320354
    321355        /*Retrieve all inputs and parameters*/
    … …  
    326360        Input* fmb_input           = basalelement->GetInput(BasalforcingsFloatingiceMeltingRateEnum);  _assert_(fmb_input);
    327361        Input* base_input          = basalelement->GetInput(BaseEnum);                                 _assert_(base_input);
    328         Input* vz_input      = NULL;
     362        Input* gllevelset_input = basalelement->GetInput(MaskOceanLevelsetEnum);              _assert_(gllevelset_input);
     363        Input* vz_input = NULL;
     364        Input* vx_input = NULL;
     365        Input* vy_input = NULL;
    329366        switch(dim){
    330                 case 1: vz_input = basalelement->GetInput(VyEnum); _assert_(vz_input); break;
    331                 case 2: vz_input = basalelement->GetInput(VzEnum); _assert_(vz_input); break;
     367                case 1:
     368                        vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
     369                        vz_input=basalelement->GetInput(VyEnum) ; _assert_(vz_input);
     370                        break;
     371                case 2:
     372                        vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
     373                        vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input);
     374                        vz_input=basalelement->GetInput(VzEnum); _assert_(vz_input);
     375                        break;
    332376                default: _error_("not implemented");
    333377        }
     378        IssmDouble h = basalelement->CharacteristicLength();
     379
     380        /*Recover portion of element that is grounded*/
     381        basalelement->GetVerticesCoordinates(&xyz_list);
     382        phi=basalelement->GetGroundedPortion(xyz_list);
     383        Gauss*      gauss     = NULL;
     384        if(melt_style==SubelementMelt2Enum){
     385                basalelement->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
     386                gauss = basalelement->NewGauss(point1,fraction1,fraction2,3);
     387        }
     388        else{
     389                gauss = basalelement->NewGauss(3);   
     390        }
    334391
    335392        /* Start  looping on the number of gaussian points: */
    336         Gauss* gauss=basalelement->NewGauss(2);
    337393        while(gauss->next()){
    338394
    … …  
    340396                basalelement->NodalFunctions(basis,gauss);
    341397
    342                 vz_input->GetInputValue(&vz,gauss);
     398                vx_input->GetInputValue(&vx,gauss); 
     399                vy_input->GetInputValue(&vy,gauss); 
     400                vz_input->GetInputValue(&vz,gauss); 
    343401                gmb_input->GetInputValue(&gmb,gauss);
    344402                fmb_input->GetInputValue(&fmb,gauss);
    345403                base_input->GetInputValue(&bed,gauss);
    346404                groundedice_input->GetInputValue(&phi,gauss);
    347                 if(phi>0) mb=gmb;
    348                 else mb=fmb;
     405                gllevelset_input->GetInputValue(&gllevelset,gauss);
     406                if(melt_style==SubelementMelt1Enum){
     407                        //if (phi>0.999999999) mb=gmb;
     408                        //else mb=(1-phi)*fmb+phi*gmb; // phi is the fraction of grounded ice so (1-phi) is floating
     409                        if(phi>0) mb=gmb;
     410                        else mb=fmb;
     411                }
     412                else if(melt_style==SubelementMelt2Enum){
     413                        if(gllevelset>0.) mb=gmb;
     414                        else mb=fmb;
     415                }
     416                else if(melt_style==NoMeltOnPartiallyFloatingEnum){
     417                        if (phi<0.00000001) mb=fmb; 
     418                        else mb=gmb;
     419                }
     420                else if(melt_style==FullMeltOnPartiallyFloatingEnum){
     421                        if (phi<0.99999999) mb=fmb; 
     422                        else mb=gmb;
     423                }
     424                else  _error_("melt interpolation "<<EnumToStringx(melt_style)<<" not implemented yet");
    349425
    350426                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(bed+dt*(mb) + dt*vz)*basis[i];
     427
     428                if(stabilization==5){
     429                        /*SUPG*/
     430                        basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     431                        if(dim==1){
     432                                vx_input->GetInputAverage(&vx);
     433                                tau=h/(2.*fabs(vx)+1e-10);
     434                                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(dt*mb+dt*vz)*tau*(vx*dbasis[0*numnodes+i]);
     435                        }
     436                        else{
     437                                vx_input->GetInputAverage(&vx);
     438                                vy_input->GetInputAverage(&vy);
     439                                tau=1*h/(2.*pow(vx*vx+vy*vy,0.5)+1e-10);
     440                                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(bed*0.+dt*mb+dt*vz)*tau*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
     441                        }
     442                }
    351443        }
    352444
    … …  
    360452}/*}}}*/
    361453void           FreeSurfaceBaseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    362            _error_("not implemented yet");
     454        _error_("not implemented yet");
    363455}/*}}}*/
    364456void           FreeSurfaceBaseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/

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

    r26090 r26894  
    8080        iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
    8181        iomodel->FetchDataToInput(inputs,elements,"md.initialization.vx",VxEnum);
     82        iomodel->FetchDataToInput(inputs,elements,"md.initialization.vy",VyEnum);
    8283        if(iomodel->domaintype!=Domain2DhorizontalEnum){
    8384                iomodel->FetchDataToInput(inputs,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
    … …  
    111112}/*}}}*/
    112113ElementMatrix* FreeSurfaceTopAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
    113 _error_("Not implemented");
     114        _error_("Not implemented");
    114115}/*}}}*/
    115116ElementMatrix* FreeSurfaceTopAnalysis::CreateKMatrix(Element* element){/*{{{*/
    … …  
    120121        IssmDouble *xyz_list  = NULL;
    121122        IssmDouble  Jdet,D_scalar,dt,h;
    122         IssmDouble  vel,vx,vy;
     123        IssmDouble  vel,vx,vy,tau;
    123124
    124125        /*Get top element*/
    … …  
    195196                }
    196197
    197                 if(stabilization==2){
     198                if(stabilization==1){
     199                        /*artifical diffusion*/
     200                        if(dim==1){
     201                                vx_input->GetInputAverage(&vx);
     202                                D[0]=h/2.*fabs(vx);
     203                        }
     204                        else{
     205                                vx_input->GetInputAverage(&vx);
     206                                vy_input->GetInputAverage(&vy);
     207
     208                                D[0*dim+0]=h/2.0*fabs(vx);
     209                                D[1*dim+1]=h/2.0*fabs(vy);
     210                        }
     211                }
     212                else if(stabilization==2){
    198213                        /*Streamline upwinding*/
    199214                        if(dim==1){
    … …  
    209224                        }
    210225                }
    211                 else if(stabilization==1){
    212                         /*SSA*/
     226                else if(stabilization==5){
     227                        /*SUPG*/
    213228                        if(dim==1){
    214229                                vx_input->GetInputAverage(&vx);
    215                                 D[0]=h/2.*fabs(vx);
     230                                tau=h/(2.*fabs(vx)+1e-10);
    216231                        }
    217232                        else{
    218233                                vx_input->GetInputAverage(&vx);
    219234                                vy_input->GetInputAverage(&vy);
    220 
    221                                 D[0*dim+0]=h/2.0*fabs(vx);
    222                                 D[1*dim+1]=h/2.0*fabs(vy);
     235                                tau=1*h/(2.*pow(vx*vx+vy*vy,0.5)+1e-10);
    223236                        }
    224237                }
    … …  
    229242                                        for(int j=0;j<numnodes;j++){
    230243                                                Ke->values[i*numnodes+j] += (
    231                                                                         dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) +
    232                                                                         dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j])
    233                                                                         );
     244                                                                dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) +
     245                                                                dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j])
     246                                                                );
    234247                                        }
    235248                                }
    … …  
    237250                        else{
    238251                                for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += dbasis[0*numnodes+i]*D[0]*dbasis[0*numnodes+j];
     252                        }
     253                }
     254                else if(stabilization==5){
     255                        D_scalar=gauss->weight*Jdet*dt;
     256                        if(dim==2){
     257                                for(int i=0;i<numnodes;i++){
     258                                        for(int j=0;j<numnodes;j++){
     259                                                Ke->values[i*numnodes+j]+=tau*D_scalar*
     260                                                        (vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i])*
     261                                                        (vx*dbasis[0*numnodes+j]+vy*dbasis[1*numnodes+j]);
     262                                        }
     263                                }
     264                        }
     265                        else{
     266                                for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j]+=tau*D_scalar*(vx*dbasis[0*numnodes+i])*(vx*dbasis[0*numnodes+j]);
    239267                        }
    240268                }
    … …  
    251279}/*}}}*/
    252280ElementVector* FreeSurfaceTopAnalysis::CreatePVector(Element* element){/*{{{*/
     281
    253282        /*Intermediaries*/
    254         int         domaintype,dim;
     283        int         domaintype,dim,stabilization;
    255284        IssmDouble  Jdet,dt;
    256         IssmDouble  ms,surface,vz;
     285        IssmDouble  ms,surface,vx,vy,vz,tau;
    257286        Element*    topelement = NULL;
    258287        IssmDouble *xyz_list  = NULL;
    … …  
    284313        ElementVector* pe    = topelement->NewElementVector();
    285314        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     315        IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
    286316
    287317        /*Retrieve all inputs and parameters*/
    … …  
    291321        Input* surface_input = topelement->GetInput(SurfaceEnum);                     _assert_(surface_input);
    292322        Input* vz_input      = NULL;
     323        Input* vx_input = NULL;
     324        Input* vy_input = NULL;
    293325        switch(dim){
    294                 case 1: vz_input = topelement->GetInput(VyEnum); _assert_(vz_input); break;
    295                 case 2: vz_input = topelement->GetInput(VzEnum); _assert_(vz_input); break;
     326                case 1:
     327                        vx_input=topelement->GetInput(VxEnum); _assert_(vx_input);
     328                        vz_input = topelement->GetInput(VyEnum) ; _assert_(vz_input);
     329                        break;
     330                case 2:
     331                        vx_input=topelement->GetInput(VxEnum); _assert_(vx_input);
     332                        vy_input = topelement->GetInput(VyEnum); _assert_(vy_input);
     333                        vz_input = topelement->GetInput(VzEnum); _assert_(vz_input);
     334                        break;
    296335                default: _error_("not implemented");
    297336        }
     337        IssmDouble h = topelement->CharacteristicLength();
    298338
    299339        /*Initialize mb_correction to 0, do not forget!:*/
    … …  
    312352        }
    313353
     354        if(stabilization==5){
     355                /*SUPG*/
     356                topelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     357                if(dim==1){
     358                        vx_input->GetInputAverage(&vx);
     359                        tau=h/(2.*fabs(vx)+1e-10);
     360                        for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(dt*ms+dt*vz)*tau*(vx*dbasis[0*numnodes+i]);
     361                }
     362                else{
     363                        vx_input->GetInputAverage(&vx);
     364                        vy_input->GetInputAverage(&vy);
     365                        tau=h/(2.*pow(vx*vx+vy*vy,0.5)+1e-10);
     366                        for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(dt*ms+dt*vz)*tau*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
     367                }
     368        }
     369
    314370        /*Clean up and return*/
    315371        xDelete<IssmDouble>(xyz_list);
    316372        xDelete<IssmDouble>(basis);
     373        xDelete<IssmDouble>(dbasis);
    317374        delete gauss;
    318375        if(topelement->IsSpawnedElement()){topelement->DeleteMaterials(); delete topelement;};
    … …  
    321378}/*}}}*/
    322379void           FreeSurfaceTopAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    323            _error_("not implemented yet");
     380        _error_("not implemented yet");
    324381}/*}}}*/
    325382void           FreeSurfaceTopAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element*  element,int control_type,int control_interp,int control_index){/*{{{*/
    … …  
    329386
    330387        element->InputUpdateFromSolutionOneDof(solution,SurfaceEnum);
     388
     389        /*Now, we need to do some "processing"*/
     390        int numvertices = element->GetNumberOfVertices();
     391        int        migration_style;
     392
     393        IssmDouble* surface = xNew<IssmDouble>(numvertices);
     394        IssmDouble* newsurface = xNew<IssmDouble>(numvertices);
     395        IssmDouble* thickness = xNew<IssmDouble>(numvertices);
     396        IssmDouble* base = xNew<IssmDouble>(numvertices);
     397        IssmDouble* bed = xNew<IssmDouble>(numvertices);
     398        IssmDouble* phi = xNew<IssmDouble>(numvertices);
     399        IssmDouble* sealevel = xNew<IssmDouble>(numvertices);
     400
     401        IssmDouble minthickness = element->FindParam(MasstransportMinThicknessEnum);
     402        IssmDouble rho_ice      = element->FindParam(MaterialsRhoIceEnum);
     403        IssmDouble rho_water    = element->FindParam(MaterialsRhoSeawaterEnum);
     404
     405        bool isgroundingline;
     406        element->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
     407        element->FindParam(&migration_style,GroundinglineMigrationEnum);
     408        if(isgroundingline) element->GetInputListOnVertices(&bed[0],BedEnum);
     409
     410        element->GetInputListOnVertices(&base[0],BaseEnum);
     411        element->GetInputListOnVertices(&surface[0],SurfaceEnum);
     412        element->GetInputListOnVertices(&phi[0],MaskOceanLevelsetEnum);
     413        element->GetInputListOnVertices(&sealevel[0],SealevelEnum);
     414
     415        for(int i=0;i<numvertices;i++){
     416                newsurface[i]=surface[i];
     417                thickness[i]=surface[i]-base[i];
     418                /*Check solution*/
     419                if(xIsNan<IssmDouble>(thickness[i])) _error_("NaN found in solution vector");
     420                if(xIsInf<IssmDouble>(thickness[i])) _error_("Inf found in solution vector");
     421
     422                /* check for thickness<minthickness */
     423                if(thickness[i]<minthickness){
     424                        thickness[i]=minthickness;
     425                        if(phi[i]>0.){
     426                                if(base[i]<=bed[i]) base[i] = bed[i];
     427                                newsurface[i] = base[i]+minthickness;
     428                        }else{
     429                                // assume floatation condition
     430                                newsurface[i] = (1.-rho_ice/rho_water)*minthickness;
     431                                base[i] = -rho_ice/rho_water*minthickness;
     432                        }
     433                }
     434
     435                /* update thickness */
     436                thickness[i]=newsurface[i]-base[i];
     437                /* some checks */
     438                if(thickness[i]<0.) _error_("thickness<0");
     439                if(newsurface[i]<base[i]) _error_("surface<base");
     440        }
     441
     442        /* update inputs */
     443        element->AddInput(BaseEnum,base,element->GetElementType());
     444        element->AddInput(SurfaceEnum,newsurface,element->GetElementType());
     445        element->AddInput(ThicknessEnum,thickness,element->GetElementType());
     446
     447        /* Free resources */
     448        xDelete<IssmDouble>(newsurface);
     449        xDelete<IssmDouble>(surface);
     450        xDelete<IssmDouble>(thickness);
     451        xDelete<IssmDouble>(base);
     452        xDelete<IssmDouble>(bed);
     453        xDelete<IssmDouble>(phi);
     454        xDelete<IssmDouble>(sealevel);
    331455}/*}}}*/
    332456void           FreeSurfaceTopAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/

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

    r26468 r26894  
    4949                solutionsequence_linear(femmodel);
    5050                femmodel->parameters->SetParam(SurfaceEnum,InputToExtrudeEnum);
     51                extrudefromtop_core(femmodel);
     52                femmodel->parameters->SetParam(ThicknessEnum,InputToExtrudeEnum);
     53                extrudefromtop_core(femmodel);
     54                femmodel->parameters->SetParam(BaseEnum,InputToExtrudeEnum);
    5155                extrudefromtop_core(femmodel);
    5256        }
Note: See TracChangeset for help on using the changeset viewer.