Changeset 18344

Timestamp:

08/07/14 16:21:11 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added Low order scheme of Flux Transport Correction ALgorithm

Location:

issm/trunk-jpl/src/c

Files:

• . (modified) (1 prop)
• Makefile.am (modified) (1 diff)
• analyses/MasstransportAnalysis.cpp (modified) (4 diffs)
• analyses/MasstransportAnalysis.h (modified) (1 diff)
• classes/Constraints/Constraint.h (modified) (1 diff)
• classes/Constraints/SpcDynamic.h (modified) (1 diff)
• classes/Constraints/SpcStatic.h (modified) (1 diff)
• classes/Constraints/SpcTransient.cpp (modified) (5 diffs)
• classes/Constraints/SpcTransient.h (modified) (2 diffs)
• classes/matrix/ElementMatrix.cpp (modified) (1 diff)
• classes/matrix/ElementMatrix.h (modified) (1 diff)
• cores/masstransport_core.cpp (modified) (3 diffs)
• solutionsequences/solutionsequence_fct.cpp (added)
• solutionsequences/solutionsequences.h (modified) (1 diff)

issm/trunk-jpl/src/c

@@ -5 +5 @@
 *.obj
 *.exe
+issm
+kriging
 appscan.*
 ouncemake*

@@ -5 +5 @@
 *.obj
 *.exe
+issm
+kriging
 appscan.*
 ouncemake*

issm/trunk-jpl/src/c/Makefile.am

@@ -354 +354 @@
                                         ./solutionsequences/solutionsequence_nonlinear.cpp\
                                         ./solutionsequences/solutionsequence_newton.cpp\
+                                        ./solutionsequences/solutionsequence_fct.cpp\
                                         ./solutionsequences/convergence.cpp\
                                         ./classes/Options/Options.h\

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

@@ -70 +70 @@
         if(stabilization==3){
                 iomodel->FetchDataToInput(elements,MasstransportSpcthicknessEnum); //for DG, we need the spc in the element
+        }
+        if(stabilization==4){
+                iomodel->FetchDataToInput(elements,MasstransportSpcthicknessEnum); //for FCT, we need the spc in the element (penlaties)
         }
 
@@ -148 +151 @@
         iomodel->Constant(&stabilization,MasstransportStabilizationEnum);
 
-        /*Do not add constraints in DG, they are weakly imposed*/
+        /*Do not add constraints in DG,  they are weakly imposed*/
         if(stabilization!=3){
                 IoModelToConstraintsx(constraints,iomodel,MasstransportSpcthicknessEnum,MasstransportAnalysisEnum,P1Enum);
         }
+        /*Do not add constraints in FCT, they are imposed using penalties*/
 }/*}}}*/
 void MasstransportAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
@@ -675 +679 @@
 }/*}}}*/
 void MasstransportAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
-           _error_("not implemented yet");
+        element->GetSolutionFromInputsOneDof(solution,ThicknessEnum);
 }/*}}}*/
 void MasstransportAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
@@ -773 +777 @@
         return;
 }/*}}}*/
+
+/*Flux Correction Transport*/
+void           MasstransportAnalysis::LumpedMassMatrix(Vector<IssmDouble>** pMlff,FemModel* femmodel){/*{{{*/
+
+        /*Intermediaries*/
+        int  configuration_type;
+
+        /*Initialize Lumped mass matrix (actually we just save its diagonal)*/
+        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+        int fsize      = femmodel->nodes->NumberOfDofs(configuration_type,FsetEnum);
+        int flocalsize = femmodel->nodes->NumberOfDofsLocal(configuration_type,FsetEnum);
+        Vector<IssmDouble>* Mlff = new Vector<IssmDouble>(flocalsize,fsize);
+
+        /*Create and assemble matrix*/
+        for(int i=0;i<femmodel->elements->Size();i++){
+                Element*       element = dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                ElementMatrix* MLe     = this->CreateMassMatrix(element);
+                if(MLe){
+                        MLe->Lump();
+                        MLe->AddDiagonalToGlobal(Mlff);
+                }
+                delete MLe;
+        }
+        Mlff->Assemble();
+
+        /*Assign output pointer*/
+        *pMlff=Mlff;
+}/*}}}*/
+ElementMatrix* MasstransportAnalysis::CreateMassMatrix(Element* element){/*{{{*/
+
+        /* Check if ice in element */
+        if(!element->IsIceInElement()) return NULL;
+
+        /*Intermediaries*/
+        IssmDouble  D,Jdet;
+        IssmDouble* xyz_list = NULL;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Initialize Element vector and other vectors*/
+        ElementMatrix* Me     = element->NewElementMatrix();
+        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinates(&xyz_list);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGauss(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                element->NodalFunctions(basis,gauss);
+
+                D=gauss->weight*Jdet;
+                TripleMultiply(basis,1,numnodes,1,
+                                        &D,1,1,0,
+                                        basis,1,numnodes,0,
+                                        &Me->values[0],1);
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(basis);
+        delete gauss;
+        return Me;
+}/*}}}*/
+void           MasstransportAnalysis::FctKMatrix(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,FemModel* femmodel){/*{{{*/
+
+        /*Output*/
+        Matrix<IssmDouble>* Kff = NULL;
+        Matrix<IssmDouble>* Kfs = NULL;
+
+        /*Initialize Jacobian Matrix*/
+        AllocateSystemMatricesx(&Kff,&Kfs,NULL,NULL,femmodel);
+
+        /*Create and assemble matrix*/
+        for(int i=0;i<femmodel->elements->Size();i++){
+                Element*       element = dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                ElementMatrix* Ke     = this->CreateFctKMatrix(element);
+                if(Ke) Ke->AddToGlobal(Kff,Kfs);
+                delete Ke;
+        }
+        Kff->Assemble();
+        Kfs->Assemble();
+
+        /*Assign output pointer*/
+        *pKff=Kff;
+        *pKfs=Kfs;
+}/*}}}*/
+ElementMatrix* MasstransportAnalysis::CreateFctKMatrix(Element* element){/*{{{*/
+
+        /* Check if ice in element */
+        if(!element->IsIceInElement()) return NULL;
+
+        /*Intermediaries */
+        IssmDouble Jdet;
+        IssmDouble vx,vy;
+        IssmDouble* xyz_list = NULL;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+        int dim      = 2;
+
+        /*Initialize Element vector and other vectors*/
+        ElementMatrix* Ke     = element->NewElementMatrix();
+        IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
+        IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
+        IssmDouble*    D      = xNewZeroInit<IssmDouble>(dim*dim);
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinates(&xyz_list);
+        Input* vxaverage_input=element->GetInput(VxEnum); _assert_(vxaverage_input);
+        Input* vyaverage_input=element->GetInput(VyEnum); _assert_(vyaverage_input);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGauss(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                GetB(B,element,dim,xyz_list,gauss);
+                GetBprime(Bprime,element,dim,xyz_list,gauss);
+                vxaverage_input->GetInputValue(&vx,gauss);
+                vyaverage_input->GetInputValue(&vy,gauss);
+
+                D[0*dim+0] = -gauss->weight*vx*Jdet;
+                D[1*dim+1] = -gauss->weight*vy*Jdet;
+
+                TripleMultiply(B,dim,numnodes,1,
+                                        D,dim,dim,0,
+                                        Bprime,dim,numnodes,0,
+                                        &Ke->values[0],1);
+
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(B);
+        xDelete<IssmDouble>(Bprime);
+        xDelete<IssmDouble>(D);
+        delete gauss;
+        return Ke;
+}/*}}}*/

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

@@ -36 +36 @@
                 void InputUpdateFromSolution(IssmDouble* solution,Element* element);
                 void UpdateConstraints(FemModel* femmodel);
+
+                /*FCT*/
+                void           LumpedMassMatrix(Vector<IssmDouble>** pMLff,FemModel* femmodel);
+                void           FctKMatrix(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,FemModel* femmodel);
+                ElementMatrix* CreateMassMatrix(Element* element);
+                ElementMatrix* CreateFctKMatrix(Element* element);
 };
 #endif

issm/trunk-jpl/src/c/classes/Constraints/Constraint.h

@@ -21 +21 @@
                 virtual      ~Constraint(){};
                 virtual void ConstrainNode(Nodes* nodes,Parameters* parameters)=0;
+                virtual void PenaltyDofAndValue(int* dof,IssmDouble* value,Nodes* nodes,Parameters* parameters)=0;
                 virtual bool InAnalysis(int analysis_type)=0;
 

issm/trunk-jpl/src/c/classes/Constraints/SpcDynamic.h

@@ -38 +38 @@
                 void   ConstrainNode(Nodes* nodes,Parameters* parameters);
                 bool   InAnalysis(int analysis_type);
+                void   PenaltyDofAndValue(int* dof,IssmDouble* value,Nodes* nodes,Parameters* parameters){_error_("not implemented yet");};
                 /*}}}*/
                 /*SpcDynamic management:{{{ */

issm/trunk-jpl/src/c/classes/Constraints/SpcStatic.h

@@ -37 +37 @@
                 void   ConstrainNode(Nodes* nodes,Parameters* parameters);
                 bool   InAnalysis(int analysis_type);
+                void   PenaltyDofAndValue(int* dof,IssmDouble* value,Nodes* nodes,Parameters* parameters){_error_("not implemented yet");};
                 /*}}}*/
                 /*SpcStatic management:{{{ */

issm/trunk-jpl/src/c/classes/Constraints/SpcTransient.cpp

@@ -15 +15 @@
 /*SpcTransient constructors and destructor*/
 SpcTransient::SpcTransient(){/*{{{*/
-        sid=-1;
-        nodeid=-1;
-        dof=-1;
-        values=NULL;
-        times=NULL;
-        nsteps=-1;
-        analysis_type=-1;
+        penalty       = false;
+        sid           = -1;
+        nodeid        = -1;
+        dof           = -1;
+        values        = NULL;
+        times         = NULL;
+        nsteps        = -1;
+        analysis_type = -1;
         return;
 }
@@ -27 +28 @@
 SpcTransient::SpcTransient(int spc_sid,int spc_nodeid, int spc_dof,int spc_nsteps, IssmDouble* spc_times, IssmDouble* spc_values,int spc_analysis_type){/*{{{*/
 
-        sid=spc_sid;
-        nodeid=spc_nodeid;
-        dof=spc_dof;
-        nsteps=spc_nsteps;
+        penalty= false;
+        sid    = spc_sid;
+        nodeid = spc_nodeid;
+        dof    = spc_dof;
+        nsteps = spc_nsteps;
         if(spc_nsteps){
-                values=xNew<IssmDouble>(spc_nsteps);
-                times=xNew<IssmDouble>(spc_nsteps);
+                values = xNew<IssmDouble>(spc_nsteps);
+                times  = xNew<IssmDouble>(spc_nsteps);
                 xMemCpy<IssmDouble>(values,spc_values,nsteps);
                 xMemCpy<IssmDouble>(times,spc_times,nsteps);
@@ -69 +71 @@
 }               
 /*}}}*/
-int    SpcTransient::Id(void){ return sid; }/*{{{*/
+int  SpcTransient::Id(void){/*{{{*/
+        return sid;
+}
 /*}}}*/
 int SpcTransient::ObjectEnum(void){/*{{{*/
@@ -101 +105 @@
         node=(Node*)nodes->GetObjectById(NULL,nodeid);
 
-        if(node){ //in case the spc is dealing with a node on another cpu
+        if(!this->penalty && node){ //in case the spc is dealing with a node on another cpu
 
                 /*Retrieve time in parameters: */
@@ -137 +141 @@
 }
 /*}}}*/
+void SpcTransient::PenaltyDofAndValue(int* pdof,IssmDouble* pvalue,Nodes* nodes,Parameters* parameters){/*{{{*/
+
+        if(!this->penalty) _error_("cannot return dof and value for non penalty constraint");
+
+        Node       *node  = NULL;
+        IssmDouble  time  = 0.;
+        int         i,gdof;
+        IssmDouble  alpha = -1.;
+        IssmDouble  value;
+        bool        found = false;
+
+        /*Chase through nodes and find the node to which this SpcTransient applys: */
+        node=(Node*)nodes->GetObjectById(NULL,nodeid);
+
+        if(node){ //in case the spc is dealing with a node on another cpu
+
+                /*Retrieve time in parameters: */
+                parameters->FindParam(&time,TimeEnum);
+
+                /*Now, go fetch value for this time: */
+                if (time<=times[0]){
+                        value=values[0];
+                        found=true;
+                }
+                else if (time>=times[nsteps-1]){
+                        value=values[nsteps-1];
+                        found=true;
+                }
+                else{
+                        for(i=0;i<nsteps-1;i++){
+                                if (times[i]<=time && time<times[i+1]){
+                                        alpha=(time-times[i])/(times[i+1]-times[i]);
+                                        value=(1-alpha)*values[i]+alpha*values[i+1];
+                                        found=true;
+                                        break;
+                                }
+                        }
+                }
+                if(!found)_error_("could not find time segment for constraint");
+
+                /*Get gdof */
+                gdof = node->GetDof(dof,GsetEnum);
+                if(xIsNan<IssmDouble>(value)){
+                        gdof = -1;
+                }
+        }
+        else{
+                value = NAN;
+                gdof = -1;
+        }
+
+        /*Assign output pointers*/
+        *pdof   = gdof;
+        *pvalue = value;
+}
+/*}}}*/
 
 /*SpcTransient functions*/

issm/trunk-jpl/src/c/classes/Constraints/SpcTransient.h

@@ -14 +14 @@
 
         private: 
-                int     sid; /*! id, to track it*/
-                int     nodeid; /*!node id*/
-                int dof; /*!component*/
-                IssmDouble* values; /*different values in time*/
-                IssmDouble* times; /*different time steps*/
-                int nsteps; /*number of time steps*/
-                int analysis_type;
+                bool        penalty;         /*Is this a penalty constraint */
+                int         sid;             /* id, to track it             */
+                int         nodeid;          /*node id                      */
+                int         dof;             /*component                    */
+                IssmDouble *values;          /*different values in time     */
+                IssmDouble *times;           /*different time steps         */
+                int         nsteps;          /*number of time steps         */
+                int         analysis_type;
 
         public:
@@ -39 +40 @@
                 void   ConstrainNode(Nodes* nodes,Parameters* parameters);
                 bool   InAnalysis(int analysis_type);
+                void   PenaltyDofAndValue(int* dof,IssmDouble* value,Nodes* nodes,Parameters* parameters);
                 /*}}}*/
                 /*SpcTransient management:{{{ */

issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.cpp

@@ -322 +322 @@
                 }
 
+        }
+        else{
+                _error_("non dofsymmetrical matrix AddToGlobal routine not support yet!");
+        }
+
+}
+/*}}}*/
+void ElementMatrix::AddDiagonalToGlobal(Vector<IssmDouble>* pf){/*{{{*/
+
+        IssmDouble* localvalues=NULL;
+
+        /*Check that pf is not NULL*/
+        _assert_(pf); 
+
+        /*In debugging mode, check consistency (no NaN, and values not too big)*/
+        this->CheckConsistency();
+
+        if(this->dofsymmetrical){
+                /*only use row dofs to add values into global matrices: */
+
+                if(this->row_fsize){
+                        /*first, retrieve values that are in the f-set from the g-set values matrix: */
+                        localvalues=xNew<IssmDouble>(this->row_fsize);
+                        for(int i=0;i<this->row_fsize;i++){
+                                localvalues[i] = this->values[this->ncols*this->row_flocaldoflist[i]+ this->row_flocaldoflist[i]];
+                        }
+
+                        /*add local values into global  matrix, using the fglobaldoflist: */
+                        pf->SetValues(this->row_fsize,this->row_fglobaldoflist,localvalues,ADD_VAL);
+
+                        /*Free ressources:*/
+                        xDelete<IssmDouble>(localvalues);
+                }
         }
         else{

issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.h

@@ -60 +60 @@
                 void AddToGlobal(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs);
                 void AddToGlobal(Matrix<IssmDouble>* Jff);
+                void AddDiagonalToGlobal(Vector<IssmDouble>* pf);
                 void Echo(void);
                 void CheckConsistency(void);

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

@@ -16 +16 @@
         bool   save_results;
         bool   isFS,isfreesurface,dakota_analysis;
-        int    solution_type;
+        int    solution_type,stabilization;
         char** requested_outputs = NULL;
 
@@ -29 +29 @@
         femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
         femmodel->parameters->FindParam(&numoutputs,MasstransportNumRequestedOutputsEnum);
+        femmodel->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);
         if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,MasstransportRequestedOutputsEnum);
 
@@ -49 +50 @@
         else{
                 if(VerboseSolution()) _printf0_("   call computational core\n");
-                solutionsequence_linear(femmodel);
+                if(stabilization==4){
+                        solutionsequence_fct(femmodel);
+                }
+                else{
+                        solutionsequence_linear(femmodel);
+                }
                 femmodel->parameters->SetParam(ThicknessEnum,InputToExtrudeEnum);
                 extrudefrombase_core(femmodel);

issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h

@@ -16 +16 @@
 void solutionsequence_nonlinear(FemModel* femmodel,bool conserve_loads);
 void solutionsequence_newton(FemModel* femmodel);
+void solutionsequence_fct(FemModel* femmodel);
 void solutionsequence_FScoupling_nonlinear(FemModel* femmodel,bool conserve_loads);
 void solutionsequence_linear(FemModel* femmodel);

@@ -5 +5 @@
 *.obj
 *.exe
+issm
+kriging
 appscan.*
 ouncemake*