Changeset 24730

Timestamp:

04/23/20 11:14:03 (5 years ago)

Author:

tsantos

Message:

CHG: inserting PVector in the FCT method

Location:

issm/trunk-jpl/src/c/analyses

Files:

• MasstransportAnalysis.cpp (modified) (5 diffs)
• MasstransportAnalysis.h (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

@@ -632 +632 @@
 //              _error_("S");
 //      }
-
+        
         h=element->CharacteristicLength();
         
@@ -936 +936 @@
         /*Intermediaries */
         IssmDouble Jdet,D_scalar;
-        IssmDouble vx,vy;
+        IssmDouble vx,vy,dvxdx,dvydy;
         IssmDouble* xyz_list = NULL;
+        IssmDouble dvx[2],dvy[2];
 
         /*Fetch number of nodes and dof for this finite element*/
@@ -964 +965 @@
                 vxaverage_input->GetInputValue(&vx,gauss);
                 vyaverage_input->GetInputValue(&vy,gauss);
+                vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
+                vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
+                dvxdx=dvx[0];
+                dvydy=dvy[1];
 
                 D_scalar = gauss->weight*Jdet;
                 for(int i=0;i<numnodes;i++){
                         for(int j=0;j<numnodes;j++){
+                                /*\phi_i \phi_j \nabla\cdot v*/
+                                Ke->values[i*numnodes+j] += -D_scalar*basis[i]*basis[j]*(dvxdx+dvydy);
+                                /*\phi_i v\cdot\nabla\phi_j*/
                                 Ke->values[i*numnodes+j] += -D_scalar*(vx*dbasis[0*numnodes+j]*basis[i] + vy*dbasis[1*numnodes+j]*basis[i]);
                         }
@@ -1021 +1029 @@
         return Me;
 }/*}}}*/
+ElementVector* MasstransportAnalysis::CreateFctPVector(Element* element){/*{{{*/
+
+        /* Check if ice in element */
+        if(!element->IsIceInElement()) return NULL;
+
+        /*Intermediaries */
+        int                     stabilization,dim,domaintype;
+        int         melt_style,point1;
+        bool        mainlyfloating;
+        IssmDouble  fraction1,fraction2;
+        IssmDouble  Jdet;
+        IssmDouble  ms,mb,gmb,fmb;
+        IssmDouble  vx,vy,dvxdx,dvydy;
+        IssmDouble  dvx[2],dvy[2];
+        IssmDouble  gllevelset,phi=1.;
+        IssmDouble* xyz_list = NULL;
+        Gauss*      gauss     = NULL;
+
+        /*Get problem dimension*/
+        element->FindParam(&domaintype,DomainTypeEnum);
+        switch(domaintype){
+                case Domain2DverticalEnum:   dim = 1; break;
+                case Domain2DhorizontalEnum: dim = 2; break;
+                case Domain3DEnum:           dim = 2; break;
+                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+        }
+          
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Initialize Element vector and other vectors*/
+        ElementVector* pe    = element->NewElementVector();
+        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+        IssmDouble*    dbasis= xNew<IssmDouble>(dim*numnodes);
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinates(&xyz_list);
+        element->FindParam(&melt_style,GroundinglineMeltInterpolationEnum);
+        element->FindParam(&stabilization,MasstransportStabilizationEnum);
+        Input2* gmb_input        = element->GetInput2(BasalforcingsGroundediceMeltingRateEnum);  _assert_(gmb_input);
+        Input2* fmb_input        = element->GetInput2(BasalforcingsFloatingiceMeltingRateEnum);  _assert_(fmb_input);
+        Input2* gllevelset_input = element->GetInput2(MaskGroundediceLevelsetEnum);              _assert_(gllevelset_input);
+        Input2* ms_input         = element->GetInput2(SmbMassBalanceEnum);                       _assert_(ms_input);
+        Input2* vxaverage_input  = element->GetInput2(VxAverageEnum);                                                                           _assert_(vxaverage_input);
+        Input2* vyaverage_input  = element->GetInput2(VyAverageEnum);                                                                           _assert_(vyaverage_input);
+        
+        /*Recover portion of element that is grounded*/
+        phi=element->GetGroundedPortion(xyz_list);
+        if(melt_style==SubelementMelt2Enum){
+                element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
+           gauss = element->NewGauss(point1,fraction1,fraction2,3);
+        }
+        else{
+                gauss = element->NewGauss(3);
+        }
+
+        /* Start  looping on the number of gaussian points: */
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                element->NodalFunctions(basis,gauss);
+
+                ms_input->GetInputValue(&ms,gauss);
+                gmb_input->GetInputValue(&gmb,gauss);
+                fmb_input->GetInputValue(&fmb,gauss);
+                gllevelset_input->GetInputValue(&gllevelset,gauss);
+
+                if(melt_style==SubelementMelt1Enum){
+                        if (phi>0.999999999) mb=gmb;
+                        else mb=(1-phi)*fmb+phi*gmb; // phi is the fraction of grounded ice so (1-phi) is floating
+                }
+                else if(melt_style==SubelementMelt2Enum){
+                        if(gllevelset>0.) mb=gmb;
+                        else mb=fmb;
+                }
+                else if(melt_style==NoMeltOnPartiallyFloatingEnum){
+                        if (phi<0.00000001) mb=fmb;
+                        else mb=gmb;
+                }
+                else if(melt_style==FullMeltOnPartiallyFloatingEnum){
+                        if (phi<0.99999999) mb=fmb;
+                        else mb=gmb;
+                }
+                else  _error_("melt interpolation "<<EnumToStringx(melt_style)<<" not implemented yet");
+
+                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(ms-mb)*basis[i];
+        
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(basis);
+        xDelete<IssmDouble>(dbasis);
+        delete gauss;
+        return pe;
+}/*}}}*/
 void           MasstransportAnalysis::FctKMatrix(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,FemModel* femmodel){/*{{{*/
 
@@ -1049 +1154 @@
         }
 }/*}}}*/
+void           MasstransportAnalysis::FctPVector(Vector<IssmDouble>** ppf,FemModel* femmodel){/*{{{*/
+
+        /*Output*/
+        Vector<IssmDouble>* pf = NULL;
+
+        /*Initialize P vector*/
+        AllocateSystemMatricesx(NULL,NULL,NULL,&pf,femmodel);
+
+        /*Create and assemble matrix*/
+        for(int i=0;i<femmodel->elements->Size();i++){
+                Element*       element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                ElementVector* pe      = this->CreateFctPVector(element);
+                if(pe) pe->AddToGlobal(pf);
+                delete pe;
+        }
+        pf->Assemble();
+
+        /*Assign output pointer*/
+        *ppf=pf;
+}/*}}}*/
 void           MasstransportAnalysis::MassMatrix(Matrix<IssmDouble>** pMff,FemModel* femmodel){/*{{{*/
 

issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.h

@@ -38 +38 @@
                 ElementMatrix* CreateFctKMatrix(Element* element);
                 ElementMatrix* CreateMassMatrix(Element* element);
+                ElementVector* CreateFctPVector(Element* element);
                 void           FctKMatrix(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,FemModel* femmodel);
+                void                            FctPVector(Vector<IssmDouble>** ppf,FemModel* femmodel);
                 void           MassMatrix(Matrix<IssmDouble>** pMff,FemModel* femmodel);
 };