Ignore:
Timestamp:
11/01/19 12:01:57 (5 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 24310

Location:
issm/trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/src

  • issm/trunk/src/c

    • Property svn:ignore

      •  

        old new  
        55*.obj
        66*.exe
         7*.mod
        78appscan.*
        89ouncemake*
        … …  
        1415probe.results
        1516stXXXX*
        16 *.deps
        17 *.dirstamp
         17.deps
         18.dirstamp
        1819.libs
        1920issm
        … …  
        2122issm_slr
        2223issm_ocean
        23 lnb_param.mod
        24 lovenb_sub.mod
        25 model.mod
        26 util.mod
         24issm_dakota

  • issm/trunk/src/c/analyses/MasstransportAnalysis.cpp

    r23394 r24313  
    55#include "../modules/modules.h"
    66
     7#define FINITEELEMENT P1Enum
     8
    79/*Model processing*/
    810void MasstransportAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
    … …  
    1416        /*Do not add constraints in DG,  they are weakly imposed*/
    1517        if(stabilization!=3){
    16                 IoModelToConstraintsx(constraints,iomodel,"md.masstransport.spcthickness",MasstransportAnalysisEnum,P1Enum);
     18                IoModelToConstraintsx(constraints,iomodel,"md.masstransport.spcthickness",MasstransportAnalysisEnum,FINITEELEMENT);
    1719        }
    1820
    … …  
    2527
    2628        /*Intermediaries*/
    27         int element;
    2829        int penpair_ids[2];
    2930        int count=0;
    … …  
    3536
    3637        /*Loads only in DG*/
    37         if (stabilization==3){
     38        if(stabilization==3){
    3839
    3940                /*Get faces and elements*/
    … …  
    4546
    4647                        /*Get left and right elements*/
    47                         element=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
     48                        int element=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
    4849
    4950                        /*Now, if this element is not in the partition, pass: */
    5051                        if(!iomodel->my_elements[element]) continue;
    51 
    52                         /* Add load */
    53                         loads->AddObject(new Numericalflux(iomodel->loadcounter+i+1,i,i,iomodel,MasstransportAnalysisEnum));
     52                        loads->AddObject(new Numericalflux(i+1,i,i,iomodel));
    5453                }
    5554
    … …  
    7776
    7877                        /*Get node ids*/
    79                         penpair_ids[0]=iomodel->nodecounter+reCast<int>(vertex_pairing[2*i+0]);
    80                         penpair_ids[1]=iomodel->nodecounter+reCast<int>(vertex_pairing[2*i+1]);
     78                        penpair_ids[0]=reCast<int>(vertex_pairing[2*i+0]);
     79                        penpair_ids[1]=reCast<int>(vertex_pairing[2*i+1]);
    8180
    8281                        /*Create Load*/
    83                         loads->AddObject(new Penpair(
    84                                                         iomodel->loadcounter+count+1,
    85                                                         &penpair_ids[0],
    86                                                         MasstransportAnalysisEnum));
     82                        loads->AddObject(new Penpair(count+1,&penpair_ids[0]));
    8783                        count++;
    8884                }
    … …  
    9288        iomodel->DeleteData(vertex_pairing,"md.masstransport.vertex_pairing");
    9389        iomodel->DeleteData(nodeonbase,"md.mesh.vertexonbase");
    94 
    95 }/*}}}*/
    96 void MasstransportAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
     90}/*}}}*/
     91void MasstransportAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/
    9792
    9893        /*Fetch parameters: */
    … …  
    106101        if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
    107102        if(stabilization!=3){
    108                 ::CreateNodes(nodes,iomodel,MasstransportAnalysisEnum,P1Enum);
     103                ::CreateNodes(nodes,iomodel,MasstransportAnalysisEnum,FINITEELEMENT,isamr);
    109104        }
    110105        else{
    111                 ::CreateNodes(nodes,iomodel,MasstransportAnalysisEnum,P1DGEnum);
     106                ::CreateNodes(nodes,iomodel,MasstransportAnalysisEnum,P1DGEnum,isamr);
    112107        }
    113108        iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
    … …  
    134129
    135130        /*Finite element type*/
    136         finiteelement = P1Enum;
     131        finiteelement = FINITEELEMENT;
    137132        if(stabilization==3){
    138133                finiteelement = P1DGEnum;
     134                //finiteelement = P0DGEnum;
    139135        }
    140136
    … …  
    158154        iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum);
    159155        iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum);
    160 
     156        if(isgroundingline)     iomodel->FetchDataToInput(elements,"md.geometry.bed",BedEnum);
    161157        /*Initialize cumdeltalthickness input*/
    162158        InputUpdateFromConstantx(elements,0.,SealevelriseCumDeltathicknessEnum);
     159        /*Initialize ThicknessResidual input*/
     160        InputUpdateFromConstantx(elements,0.,ThicknessResidualEnum);
    163161
    164162        /*Get what we need for ocean-induced basal melting*/
    … …  
    183181                        iomodel->FetchDataToInput(elements,"md.basalforcings.basin_id",BasalforcingsPicoBasinIdEnum);
    184182                        break;
     183                case BasalforcingsIsmip6Enum:{
     184                        iomodel->FetchDataToInput(elements,"md.basalforcings.basin_id",BasalforcingsIsmip6BasinIdEnum);
     185                        iomodel->FetchDataToInput(elements,"md.basalforcings.melt_anomaly",BasalforcingsIsmip6MeltAnomalyEnum,0.);
     186                        IssmDouble** array3d = NULL; int* Ms = NULL; int* Ns = NULL; int K;
     187                        iomodel->FetchData(&array3d,&Ms,&Ns,&K,"md.basalforcings.tf");
     188                        if(!array3d) _error_("md.basalforcings.tf not found in binary file");
     189                        for(int i=0;i<elements->Size();i++){
     190                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
     191                                if(iomodel->domaintype!=Domain2DhorizontalEnum && !element->IsOnBase()) continue;
     192                                for(int kk=0;kk<K;kk++){
     193                                        element->DatasetInputAdd(BasalforcingsIsmip6TfEnum,array3d[kk],iomodel,Ms[kk],Ns[kk],1,BasalforcingsIsmip6TfEnum,7,kk);
     194                                }
     195                        }
     196                        xDelete<int>(Ms); xDelete<int>(Ns);
     197                        for(int i=0;i<K;i++) xDelete<IssmDouble>(array3d[i]);
     198                        xDelete<IssmDouble*>(array3d);
     199                        }
     200                        break;
     201                case BeckmannGoosseFloatingMeltRateEnum:
     202                        iomodel->FetchDataToInput(elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum);
     203                        iomodel->FetchDataToInput(elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum);
     204                        break;
    185205                default:
    186206                        _error_("Basal forcing model "<<EnumToStringx(basalforcing_model)<<" not supported yet");
    … …  
    246266                        Ke = CreateKMatrixCG(basalelement);
    247267                        break;
     268                case P0DGEnum:
    248269                case P1DGEnum:
    249270                        Ke = CreateKMatrixDG(basalelement);
    … …  
    268289        IssmDouble Jdet,D_scalar,dt,h;
    269290        IssmDouble vel,vx,vy,dvxdx,dvydy;
     291        IssmDouble xi,tau;
    270292        IssmDouble dvx[2],dvy[2];
    271293        IssmDouble* xyz_list = NULL;
    … …  
    286308        ElementMatrix* Ke     = element->NewElementMatrix();
    287309        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     310        IssmDouble*             dbasis = xNew<IssmDouble>(dim*numnodes);
    288311        IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
    289312        IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
    … …  
    349372                switch(stabilization){
    350373                        case 0:
    351                                 /*Nothing to be onde*/
     374                                /*Nothing to be done*/
    352375                                break;
    353376                        case 1:
    … …  
    359382                                break;
    360383                        case 2:
     384                                /*Streamline upwinding*/
     385                                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     386                                vxaverage_input->GetInputAverage(&vx);
    361387                                if(dim==1){
    362388                                        vel=fabs(vx)+1.e-8;
    363                                         D[0]=h/(2*vel)*vx*vx;
    364389                                }
    365390                                else{
    366                                         /*Streamline upwinding*/
     391                                        vyaverage_input->GetInputAverage(&vy);
    367392                                        vel=sqrt(vx*vx+vy*vy)+1.e-8;
    368                                         D[0*dim+0]=h/(2*vel)*vx*vx;
    369                                         D[1*dim+0]=h/(2*vel)*vy*vx;
    370                                         D[0*dim+1]=h/(2*vel)*vx*vy;
    371                                         D[1*dim+1]=h/(2*vel)*vy*vy;
    372393                                }
     394                                tau=h/(2*vel);
     395                                break;
     396                        case 5:
     397                                /*SUPG*/
     398                                if(dim!=2) _error_("Stabilization "<<stabilization<<" not supported yet for dim != 2");
     399                                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     400                                vxaverage_input->GetInputAverage(&vx);
     401                                vyaverage_input->GetInputAverage(&vy);
     402                                vel=sqrt(vx*vx+vy*vy)+1.e-8;
     403                                //xi=0.3130;
     404                                xi=1;
     405                                tau=xi*h/(2*vel);
    373406                                break;
    374407                        default:
    375408                                _error_("Stabilization "<<stabilization<<" not supported yet");
    376409                }
    377                 if(stabilization==1 || stabilization==2){
     410                if(stabilization==1){
     411                        /*SSA*/
    378412                        if(dim==1) D[0]=D_scalar*D[0];
    379413                        else{
    … …  
    389423                                                &Ke->values[0],1);
    390424                }
     425                if(stabilization==2){
     426                        /*Streamline upwind*/
     427                        for(int i=0;i<numnodes;i++){
     428                                for(int j=0;j<numnodes;j++){
     429                                        Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i])*(vx*dbasis[0*numnodes+j]+vy*dbasis[1*numnodes+j]);
     430                                }
     431                        }
     432                }
     433                if(stabilization==5){/*{{{*/
     434                         /*Mass matrix - part 2*/
     435                        for(int i=0;i<numnodes;i++){
     436                                for(int j=0;j<numnodes;j++){
     437                                        Ke->values[i*numnodes+j]+=gauss->weight*Jdet*tau*basis[j]*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);   
     438                                }
     439                        }
     440                        /*Mass matrix - part 3*/
     441                        for(int i=0;i<numnodes;i++){
     442                                for(int j=0;j<numnodes;j++){
     443                                        Ke->values[i*numnodes+j]+=gauss->weight*Jdet*tau*basis[j]*(basis[i]*dvxdx+basis[i]*dvydy);     
     444                                }
     445                        }
     446                       
     447                        /*Advection matrix - part 2, A*/
     448                        for(int i=0;i<numnodes;i++){
     449            for(int j=0;j<numnodes;j++){
     450               Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(vx*dbasis[0*numnodes+j]+vy*dbasis[1*numnodes+j])*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
     451            }
     452         }
     453                        /*Advection matrix - part 3, A*/
     454                        for(int i=0;i<numnodes;i++){
     455            for(int j=0;j<numnodes;j++){
     456                                        Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(vx*dbasis[0*numnodes+j]+vy*dbasis[1*numnodes+j])*(basis[i]*dvxdx+basis[i]*dvydy);
     457                                }
     458         }
     459
     460                        /*Advection matrix - part 2, B*/
     461                        for(int i=0;i<numnodes;i++){
     462            for(int j=0;j<numnodes;j++){
     463                                        Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(basis[j]*dvxdx+basis[j]*dvydy)*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
     464                                }
     465         }
     466                        /*Advection matrix - part 3, B*/
     467                        for(int i=0;i<numnodes;i++){
     468            for(int j=0;j<numnodes;j++){
     469                                        Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(basis[j]*dvxdx+basis[j]*dvydy)*(basis[i]*dvxdx+basis[i]*dvydy);   
     470                                }
     471                        }
     472                }/*}}}*/
    391473        }
    392474
    … …  
    394476        xDelete<IssmDouble>(xyz_list);
    395477        xDelete<IssmDouble>(basis);
     478        xDelete<IssmDouble>(dbasis);
    396479        xDelete<IssmDouble>(B);
    397480        xDelete<IssmDouble>(Bprime);
    … …  
    481564                        pe = CreatePVectorCG(basalelement);
    482565                        break;
     566                case P0DGEnum:
    483567                case P1DGEnum:
    484568                        pe = CreatePVectorDG(basalelement);
    … …  
    499583
    500584        /*Intermediaries */
     585        int                     stabilization,dim,domaintype;
    501586        int         melt_style,point1;
    502587        bool        mainlyfloating;
    … …  
    504589        IssmDouble  Jdet,dt;
    505590        IssmDouble  ms,mb,gmb,fmb,thickness;
     591        IssmDouble  vx,vy,vel,dvxdx,dvydy,xi,h,tau;
     592        IssmDouble  dvx[2],dvy[2];
    506593        IssmDouble  gllevelset,phi=1.;
    507594        IssmDouble* xyz_list = NULL;
    508595        Gauss*      gauss     = NULL;
    509596
     597        /*Get problem dimension*/
     598        element->FindParam(&domaintype,DomainTypeEnum);
     599        switch(domaintype){
     600                case Domain2DverticalEnum:   dim = 1; break;
     601                case Domain2DhorizontalEnum: dim = 2; break;
     602                case Domain3DEnum:           dim = 2; break;
     603                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     604        }
     605         
    510606        /*Fetch number of nodes and dof for this finite element*/
    511607        int numnodes = element->GetNumberOfNodes();
    … …  
    514610        ElementVector* pe    = element->NewElementVector();
    515611        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     612        IssmDouble*    dbasis= xNew<IssmDouble>(dim*numnodes);
    516613
    517614        /*Retrieve all inputs and parameters*/
    … …  
    519616        element->FindParam(&melt_style,GroundinglineMeltInterpolationEnum);
    520617        element->FindParam(&dt,TimesteppingTimeStepEnum);
    521         Input* gmb_input           = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);  _assert_(gmb_input);
    522         Input* fmb_input           = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum);  _assert_(fmb_input);
    523         Input* gllevelset_input   = element->GetInput(MaskGroundediceLevelsetEnum);              _assert_(gllevelset_input);
    524         Input* ms_input            = element->GetInput(SmbMassBalanceEnum);                       _assert_(ms_input);
    525         Input* thickness_input     = element->GetInput(ThicknessEnum);                            _assert_(thickness_input);
    526 
     618        element->FindParam(&stabilization,MasstransportStabilizationEnum);
     619        Input* gmb_input        = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);  _assert_(gmb_input);
     620        Input* fmb_input        = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum);  _assert_(fmb_input);
     621        Input* gllevelset_input = element->GetInput(MaskGroundediceLevelsetEnum);              _assert_(gllevelset_input);
     622        Input* ms_input         = element->GetInput(SmbMassBalanceEnum);                       _assert_(ms_input);
     623        Input* thickness_input  = element->GetInput(ThicknessEnum);                            _assert_(thickness_input);
     624        Input* vxaverage_input  = element->GetInput(VxAverageEnum);                                                                             _assert_(vxaverage_input);
     625        Input* vyaverage_input  = element->GetInput(VyAverageEnum);                                                                             _assert_(vyaverage_input);
     626
     627        h=element->CharacteristicLength();
     628       
    527629        /*Recover portion of element that is grounded*/
    528630        phi=element->GetGroundedPortion(xyz_list);
    529631        if(melt_style==SubelementMelt2Enum){
    530                 gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
    531632                element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
    532633           gauss = element->NewGauss(point1,fraction1,fraction2,3);
    … …  
    563664                else if(melt_style==FullMeltOnPartiallyFloatingEnum){
    564665                        if (phi<0.99999999) mb=fmb;
    565                         else mb=fmb;
     666                        else mb=gmb;
    566667                }
    567668                else  _error_("melt interpolation "<<EnumToStringx(melt_style)<<" not implemented yet");
    568669
    569670                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*basis[i];
     671       
     672                if(stabilization==5){ //SUPG
     673                        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     674                        vxaverage_input->GetInputAverage(&vx);
     675                        vyaverage_input->GetInputAverage(&vy);
     676                        vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
     677                        vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
     678                        vel=sqrt(vx*vx+vy*vy)+1.e-8;
     679                        dvxdx=dvx[0];
     680                        dvydy=dvy[1];
     681                        //xi=0.3130;
     682                        xi=1;
     683                        tau=xi*h/(2*vel);
     684                       
     685                        /*Force vector - part 2*/
     686                        for(int i=0;i<numnodes;i++){
     687                                pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*(tau*vx*dbasis[0*numnodes+i]+tau*vy*dbasis[1*numnodes+i]);
     688                        }
     689                        /*Force vector - part 3*/
     690                        for(int i=0;i<numnodes;i++){
     691                                pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*(tau*basis[i]*dvxdx+tau*basis[i]*dvydy);
     692                        }
     693                }
     694       
    570695        }
    571696
    … …  
    573698        xDelete<IssmDouble>(xyz_list);
    574699        xDelete<IssmDouble>(basis);
     700        xDelete<IssmDouble>(dbasis);
    575701        delete gauss;
    576702        return pe;
    … …  
    582708
    583709        /*Intermediaries */
     710        int melt_style, point1;
     711        bool mainlyfloating;
     712        IssmDouble  fraction1,fraction2,gllevelset;
    584713        IssmDouble  Jdet,dt;
    585714        IssmDouble  ms,mb,gmb,fmb,thickness,phi=1.;
    … …  
    596725        element->GetVerticesCoordinates(&xyz_list);
    597726        element->FindParam(&dt,TimesteppingTimeStepEnum);
    598         Input* gmb_input           = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(gmb_input);
    599         Input* fmb_input           = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fmb_input);
    600         Input* ms_input            = element->GetInput(SmbMassBalanceEnum);          _assert_(ms_input);
    601         Input* gllevelset_input   = element->GetInput(MaskGroundediceLevelsetEnum);             _assert_(gllevelset_input);
    602         Input* thickness_input     = element->GetInput(ThicknessEnum);                           _assert_(thickness_input);
     727        element->FindParam(&melt_style,GroundinglineMeltInterpolationEnum);
     728        Input* gmb_input        = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(gmb_input);
     729        Input* fmb_input        = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fmb_input);
     730        Input* ms_input         = element->GetInput(SmbMassBalanceEnum);                      _assert_(ms_input);
     731        Input* gllevelset_input = element->GetInput(MaskGroundediceLevelsetEnum);             _assert_(gllevelset_input);
     732        Input* thickness_input  = element->GetInput(ThicknessEnum);                           _assert_(thickness_input);
     733
     734   /*Recover portion of element that is grounded*/
     735   Gauss* gauss=NULL;
     736   phi=element->GetGroundedPortion(xyz_list);
     737   if(melt_style==SubelementMelt2Enum){
     738      element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
     739      gauss = element->NewGauss(point1,fraction1,fraction2,3);
     740   }
     741   else{
     742      gauss = element->NewGauss(3);
     743   }
    603744
    604745        /* Start  looping on the number of gaussian points: */
    605         Gauss* gauss=element->NewGauss(2);
    606746        for(int ig=gauss->begin();ig<gauss->end();ig++){
    607747                gauss->GaussPoint(ig);
    … …  
    613753                gmb_input->GetInputValue(&gmb,gauss);
    614754                fmb_input->GetInputValue(&fmb,gauss);
    615                 gllevelset_input->GetInputValue(&phi,gauss);
     755                gllevelset_input->GetInputValue(&gllevelset,gauss);
    616756                thickness_input->GetInputValue(&thickness,gauss);
    617757
    618                 if(phi>0.) mb=gmb;
    619                 else mb=fmb;
     758                if(melt_style==SubelementMelt1Enum){
     759         if (phi>0.999999999) mb=gmb;
     760         else mb=(1-phi)*fmb+phi*gmb; // phi is the fraction of grounded ice so (1-phi) is floating
     761      }
     762      else if(melt_style==SubelementMelt2Enum){
     763         if(gllevelset>0.) mb=gmb;
     764         else mb=fmb;
     765      }
     766      else if(melt_style==NoMeltOnPartiallyFloatingEnum){
     767         if (phi<0.00000001) mb=fmb;
     768         else mb=gmb;
     769      }
     770      else if(melt_style==FullMeltOnPartiallyFloatingEnum){
     771         if (phi<0.99999999) mb=fmb;
     772         else mb=gmb;
     773      }
     774      else  _error_("melt interpolation "<<EnumToStringx(melt_style)<<" not implemented yet");
    620775
    621776                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*basis[i];
    … …  
    693848void           MasstransportAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    694849
    695         int        i,hydroadjustment,domaintype;
    696         int*       doflist=NULL;
    697         IssmDouble rho_ice,rho_water,minthickness;
    698         Element*   basalelement=NULL;
    699         bool        isgroundingline=0;
    700 
    701         element->FindParam(&domaintype,DomainTypeEnum);
    702         if(domaintype!=Domain2DhorizontalEnum){
    703                 if(!element->IsOnBase()) return;
    704                 basalelement=element->SpawnBasalElement();
    705         }
    706         else{
    707                 basalelement = element;
    708         }
    709 
    710         /*Fetch number of nodes for this finite element*/
    711         int numnodes = basalelement->GetNumberOfNodes();
     850        /*Only update if on base*/
     851        if(!element->IsOnBase()) return;
    712852
    713853        /*Fetch dof list and allocate solution vector*/
    714         basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    715         IssmDouble* newthickness   = xNew<IssmDouble>(numnodes);
    716         IssmDouble* cumdeltathickness = xNew<IssmDouble>(numnodes);
    717         IssmDouble* deltathickness    = xNew<IssmDouble>(numnodes);
    718         IssmDouble* newbase        = xNew<IssmDouble>(numnodes);
    719         IssmDouble* bed            = xNew<IssmDouble>(numnodes);
    720         IssmDouble* newsurface     = xNew<IssmDouble>(numnodes);
    721         IssmDouble* oldthickness   = xNew<IssmDouble>(numnodes);
    722         IssmDouble* oldbase        = xNew<IssmDouble>(numnodes);
    723         IssmDouble* oldsurface     = xNew<IssmDouble>(numnodes);
    724         IssmDouble* phi            = xNew<IssmDouble>(numnodes);
    725         IssmDouble* sealevel       = xNew<IssmDouble>(numnodes);
     854        int *doflist = NULL;
     855        element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
     856
     857        int numnodes = element->GetNumberOfNodes();
     858        IssmDouble* newthickness = xNew<IssmDouble>(numnodes);
     859        IssmDouble* thicknessresidual = xNew<IssmDouble>(numnodes);
    726860
    727861        /*Use the dof list to index into the solution vector: */
    728         basalelement->FindParam(&minthickness,MasstransportMinThicknessEnum);
    729         for(i=0;i<numnodes;i++){
     862        IssmDouble minthickness = element->FindParam(MasstransportMinThicknessEnum);
     863        for(int i=0;i<numnodes;i++){
    730864                newthickness[i]=solution[doflist[i]];
     865                thicknessresidual[i]=0.;       
     866                /*Check solution*/
    731867                if(xIsNan<IssmDouble>(newthickness[i])) _error_("NaN found in solution vector");
    732868                if(xIsInf<IssmDouble>(newthickness[i])) _error_("Inf found in solution vector");
    733                 /*Constrain thickness to be at least 1m*/
    734                 if(newthickness[i]<minthickness) newthickness[i]=minthickness;
    735         }
     869                if(newthickness[i]<minthickness){
     870                        thicknessresidual[i]=minthickness-newthickness[i];
     871                        newthickness[i]=minthickness;
     872                }
     873        }
     874        element->AddBasalInput(ThicknessEnum,newthickness,element->GetElementType());
     875        element->AddBasalInput(ThicknessResidualEnum,thicknessresidual,element->GetElementType());
     876       
     877        xDelete<int>(doflist);
     878        xDelete<IssmDouble>(newthickness);
     879        xDelete<IssmDouble>(thicknessresidual);
     880
     881        /*Update bed and surface accordingly*/
     882
     883        /*Get basal element*/
     884        int domaintype; element->FindParam(&domaintype,DomainTypeEnum);
     885        Element* basalelement=element;
     886        if(domaintype!=Domain2DhorizontalEnum) basalelement = element->SpawnBasalElement();
     887
     888        /*Fetch number of nodes and dof for this finite element*/
     889        int numvertices = basalelement->GetNumberOfVertices();
     890
     891        /*Now, we need to do some "processing"*/
     892        newthickness  = xNew<IssmDouble>(numvertices);
     893        IssmDouble* oldthickness      = xNew<IssmDouble>(numvertices);
     894        IssmDouble* cumdeltathickness = xNew<IssmDouble>(numvertices);
     895        IssmDouble* deltathickness    = xNew<IssmDouble>(numvertices);
     896        IssmDouble* newbase           = xNew<IssmDouble>(numvertices);
     897        IssmDouble* bed               = xNew<IssmDouble>(numvertices);
     898        IssmDouble* newsurface        = xNew<IssmDouble>(numvertices);
     899        IssmDouble* oldbase           = xNew<IssmDouble>(numvertices);
     900        IssmDouble* oldsurface        = xNew<IssmDouble>(numvertices);
     901        IssmDouble* phi               = xNew<IssmDouble>(numvertices);
     902        IssmDouble* sealevel          = xNew<IssmDouble>(numvertices);
    736903
    737904        /*Get previous base, thickness, surfac and current sealevel and bed:*/
    738         basalelement->GetInputListOnNodes(&oldbase[0],BaseEnum);
    739         basalelement->GetInputListOnNodes(&oldsurface[0],SurfaceEnum);
    740         basalelement->GetInputListOnNodes(&oldthickness[0],ThicknessEnum);
    741         basalelement->GetInputListOnNodes(&phi[0],MaskGroundediceLevelsetEnum);
    742         basalelement->GetInputListOnNodes(&sealevel[0],SealevelEnum);
    743         basalelement->GetInputListOnNodes(&cumdeltathickness[0],SealevelriseCumDeltathicknessEnum);
    744 
     905        basalelement->GetInputListOnVertices(&newthickness[0],ThicknessEnum);
     906        basalelement->GetInputListOnVertices(&oldthickness[0],ThicknessOldEnum);
     907        basalelement->GetInputListOnVertices(&oldbase[0],BaseEnum);
     908        basalelement->GetInputListOnVertices(&oldsurface[0],SurfaceEnum);
     909        basalelement->GetInputListOnVertices(&phi[0],MaskGroundediceLevelsetEnum);
     910        basalelement->GetInputListOnVertices(&sealevel[0],SealevelEnum);
     911        basalelement->GetInputListOnVertices(&cumdeltathickness[0],SealevelriseCumDeltathicknessEnum);
     912
     913        /*Do we do grounding line migration?*/
     914        bool isgroundingline;
    745915        element->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
    746         if(isgroundingline) basalelement->GetInputListOnNodes(&bed[0],BedEnum);
     916        if(isgroundingline) basalelement->GetInputListOnVertices(&bed[0],BedEnum);
    747917
    748918        /*What is the delta thickness forcing the sea-level rise core: cumulated over time, hence the +=:*/
    749         for(i=0;i<numnodes;i++){
    750                 cumdeltathickness[i]+=newthickness[i]-oldthickness[i];
    751                 deltathickness[i]=newthickness[i]-oldthickness[i];
     919        for(int i=0;i<numvertices;i++){
     920                cumdeltathickness[i] += newthickness[i]-oldthickness[i];
     921                deltathickness[i]     = newthickness[i]-oldthickness[i];
    752922        }
    753923
    754924        /*Find MasstransportHydrostaticAdjustment to figure out how to update the geometry:*/
     925        int hydroadjustment;
    755926        basalelement->FindParam(&hydroadjustment,MasstransportHydrostaticAdjustmentEnum);
    756         rho_ice   = basalelement->GetMaterialParameter(MaterialsRhoIceEnum);
    757         rho_water = basalelement->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    758 
    759         for(i=0;i<numnodes;i++) {
    760                 if (phi[i]>0.){ //this is an ice sheet: just add thickness to base.
    761                         /*Update! actually, the bed has moved too!:*/
     927        IssmDouble rho_ice   = basalelement->FindParam(MaterialsRhoIceEnum);
     928        IssmDouble rho_water = basalelement->FindParam(MaterialsRhoSeawaterEnum);
     929
     930        for(int i=0;i<numvertices;i++) {
     931                if (phi[i]>0.){ //this is grounded ice: just add thickness to base.
    762932                        if(isgroundingline){
    763933                                newsurface[i] = bed[i]+newthickness[i]; //surface = bed + newthickness
    764                                 newbase[i]     = bed[i];                 //new base at new bed
     934                                newbase[i]    = bed[i];                 //new base at new bed
    765935                        }
    766936                        else{
    767937                                 newsurface[i] = oldbase[i]+newthickness[i]; //surface = oldbase + newthickness
    768                                  newbase[i]     = oldbase[i];                 //same base: do nothing
     938                                 newbase[i]    = oldbase[i];                 //same base: do nothing
    769939                        }
    770940                }
    771941                else{ //this is an ice shelf: hydrostatic equilibrium*/
    772942                        if(hydroadjustment==AbsoluteEnum){
    773                                 newsurface[i] = newthickness[i]*(1-rho_ice/rho_water)+sealevel[i];
    774                                 newbase[i]     = newthickness[i]*(-rho_ice/rho_water)+sealevel[i];
     943                                newsurface[i] = newthickness[i]*(1.-rho_ice/rho_water)+sealevel[i];
     944                                newbase[i]    = newthickness[i]*(-rho_ice/rho_water)+sealevel[i];
    775945                        }
    776946                        else if(hydroadjustment==IncrementalEnum){
    777947                                newsurface[i] = oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i])+sealevel[i]; //surface = oldsurface + (1-di) * dH
    778                                 newbase[i]     = oldbase[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i])+sealevel[i]; //base               = oldbed + di * dH
     948                                newbase[i]    = oldbase[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i])+sealevel[i]; //base               = oldbed + di * dH
    779949                        }
    780950                        else _error_("Hydrostatic adjustment " << hydroadjustment << " (" << EnumToStringx(hydroadjustment) << ") not supported yet");
    … …  
    783953
    784954        /*Add input to the element: */
    785         element->AddBasalInput(ThicknessEnum,newthickness,P1Enum);
    786955        element->AddBasalInput(SurfaceEnum,newsurface,P1Enum);
    787956        element->AddBasalInput(BaseEnum,newbase,P1Enum);
    … …  
    801970        xDelete<IssmDouble>(sealevel);
    802971        xDelete<IssmDouble>(bed);
    803 
    804972        xDelete<int>(doflist);
    805973        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    … …  
    9341102void           MasstransportAnalysis::LumpedMassMatrix(Vector<IssmDouble>** pMlff,FemModel* femmodel){/*{{{*/
    9351103
    936         /*Intermediaries*/
    937         int  configuration_type;
    938 
    9391104        /*Initialize Lumped mass matrix (actually we just save its diagonal)*/
    940         femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
    941         int fsize      = femmodel->nodes->NumberOfDofs(configuration_type,FsetEnum);
    942         int flocalsize = femmodel->nodes->NumberOfDofsLocal(configuration_type,FsetEnum);
     1105        int fsize      = femmodel->nodes->NumberOfDofs(FsetEnum);
     1106        int flocalsize = femmodel->nodes->NumberOfDofsLocal(FsetEnum);
    9431107        Vector<IssmDouble>* Mlff = new Vector<IssmDouble>(flocalsize,fsize);
    9441108
Note: See TracChangeset for help on using the changeset viewer.