Changeset 19386

Timestamp:

06/03/15 06:25:50 (10 years ago)

Author:

cborstad

Message:

CHG: working to implement flux corrected transport for damage evolution

Location:

issm/trunk-jpl/src

Files:

• c/analyses/DamageEvolutionAnalysis.cpp (modified) (3 diffs)
• c/analyses/DamageEvolutionAnalysis.h (modified) (1 diff)
• c/analyses/StressbalanceAnalysis.cpp (modified) (1 diff)
• c/cores/damage_core.cpp (modified) (3 diffs)
• c/solutionsequences/solutionsequence_fct.cpp (modified) (8 diffs)
• m/classes/damage.m (modified) (2 diffs)
• m/classes/damage.py (modified) (2 diffs)
• m/mech/mechanicalproperties.m (modified) (1 diff)
• m/mech/mechanicalproperties.py (modified) (1 diff)

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

@@ -115 +115 @@
 
                 /* no source term; damage handled in stress balance */
-                f[i]=0;
+                f[i]=0.;
         }
 
@@ -632 +632 @@
 }/*}}}*/
 void           DamageEvolutionAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
-           _error_("not implemented yet");
+        element->GetSolutionFromInputsOneDof(solution,DamageDbarEnum);
 }/*}}}*/
 void           DamageEvolutionAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
@@ -686 +686 @@
         return;
 }/*}}}*/
+
+/*Flux Correction Transport*/
+ElementMatrix* DamageEvolutionAnalysis::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;
+}/*}}}*/
+ElementMatrix* DamageEvolutionAnalysis::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           DamageEvolutionAnalysis::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 = xDynamicCast<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;
+        if(pKfs){
+                *pKfs=Kfs;
+        }
+        else{
+                delete Kfs;
+        }
+}/*}}}*/
+void           DamageEvolutionAnalysis::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 = xDynamicCast<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;
+}/*}}}*/
+void           DamageEvolutionAnalysis::MassMatrix(Matrix<IssmDouble>** pMff,FemModel* femmodel){/*{{{*/
+
+        /*Initialize Mass matrix*/
+        Matrix<IssmDouble> *Mff = NULL;
+        AllocateSystemMatricesx(&Mff,NULL,NULL,NULL,femmodel);
+
+        /*Create and assemble matrix*/
+        for(int i=0;i<femmodel->elements->Size();i++){
+                Element*       element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                ElementMatrix* MLe     = this->CreateMassMatrix(element);
+                if(MLe){
+                        MLe->AddToGlobal(Mff);
+                }
+                delete MLe;
+        }
+        Mff->Assemble();
+
+        /*Assign output pointer*/
+        *pMff=Mff;
+}/*}}}*/

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

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

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

@@ -876 +876 @@
         iomodel->Constant(&materials_type,MaterialsEnum);
         if(materials_type==MatdamageiceEnum){
-                parameters->AddObject(iomodel->CopyConstantObject(DamageC1Enum));
+                parameters->AddObject(iomodel->CopyConstantObject(DamageLawEnum));
                 parameters->AddObject(iomodel->CopyConstantObject(DamageKappaEnum));
                 parameters->AddObject(iomodel->CopyConstantObject(DamageStressThresholdEnum));

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

@@ -1 +1 @@
 /* 
- * \brief: damgage_core.cpp: core for the damage solution
+ * \brief: damage_core.cpp: core for the damage solution
  */ 
 
@@ -15 +15 @@
         bool   save_results;
         bool   dakota_analysis     = false;
-        int    solution_type;
+        int    solution_type,stabilization;
         int    numoutputs          = 0; 
         char   **requested_outputs = NULL;
@@ -25 +25 @@
         femmodel->parameters->FindParam(&numoutputs,DamageEvolutionNumRequestedOutputsEnum);
         if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,DamageEvolutionRequestedOutputsEnum);
+        femmodel->parameters->FindParam(&stabilization,DamageStabilizationEnum);
 
         if(VerboseSolution()) _printf0_("   computing damage\n");
         femmodel->SetCurrentConfiguration(DamageEvolutionAnalysisEnum);
-        solutionsequence_linear(femmodel);
+        if(stabilization==4){
+                solutionsequence_fct(femmodel);
+        }
+        else{
+                solutionsequence_linear(femmodel);
+        }
 
         if(save_results){

issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp

@@ -1 @@
-/*!\file: solutionsequence_linear.cpp
- * \brief: numerical core of linear solutions
+/*!\file: solutionsequence_fct.cpp
+ * \brief: numerical core of flux corrected transport solution
  */ 
 
@@ -7 +7 @@
 #include "../shared/shared.h"
 #include "../modules/modules.h"
+#include "../analyses/analyses.h"
 
 #ifdef _HAVE_PETSC_
@@ -119 +120 @@
 }/*}}}*/
 void CreateRHS(Vec* pRHS,Mat K,Mat D,Vec Ml,Vec u,IssmDouble theta,IssmDouble deltat,IssmDouble dmax,FemModel* femmodel,int configuration_type){/*{{{*/
-        /*Create Left Hand side of Lower order solution
+        /*Create Right Hand side of Lower order solution
          *
          * RHS = [ML + (1 − theta) deltaT L^n] u^n
@@ -166 +167 @@
 }/*}}}*/
 void RichardsonUdot(Vec* pudot,Vec u,Vec Ml,Mat K,Mat Mc){/*{{{*/
-        /*Use Richardson's formulato get udot using 5 steps and an initial guess of 0
+        /*Use Richardson's formula to get udot using 5 steps and an initial guess of 0
          *
          * udot_new = udot_old + Ml^-1 (K^(n+1) u - Mc udot_old)
@@ -336 +337 @@
 }/*}}}*/
 #endif
-void solutionsequence_fct(FemModel* femmodel){
+void solutionsequence_fct(FemModel* femmodel){/*{{{*/
 
         /*intermediary: */
         IssmDouble           theta,deltat,dmax;
-        int                  configuration_type;
+        int                  configuration_type,analysis_type;
         Vector<IssmDouble>*  Ml = NULL;
         Matrix<IssmDouble>*  K  = NULL;
@@ -346 +347 @@
         Vector<IssmDouble>*  ug = NULL;
         Vector<IssmDouble>*  uf = NULL;
-
-        /*Create analysis*/
-        MasstransportAnalysis* analysis = new MasstransportAnalysis();
+        MasstransportAnalysis* manalysis = NULL;
+        DamageEvolutionAnalysis* danalysis = NULL;
 
         /*Recover parameters: */
@@ -356 +356 @@
         theta = 0.5;
 
-        /*Create lumped mass matrix*/
-        analysis->LumpedMassMatrix(&Ml,femmodel);
-        analysis->MassMatrix(&Mc,femmodel);
-        analysis->FctKMatrix(&K,NULL,femmodel);
-        delete analysis;
+        /*Create analysis*/
+        femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+        switch(analysis_type){
+                case MasstransportAnalysisEnum:
+                        manalysis = new MasstransportAnalysis();
+                        manalysis->LumpedMassMatrix(&Ml,femmodel);
+                        manalysis->MassMatrix(&Mc,femmodel);
+                        manalysis->FctKMatrix(&K,NULL,femmodel);
+                        break;
+                case DamageEvolutionAnalysisEnum:
+                        danalysis = new DamageEvolutionAnalysis();
+                        danalysis->LumpedMassMatrix(&Ml,femmodel);
+                        danalysis->MassMatrix(&Mc,femmodel);
+                        danalysis->FctKMatrix(&K,NULL,femmodel);
+                        break;
+                default: _error_("analysis type " << EnumToStringx(analysis_type) << " not supported for FCT\n");
+        }
+        delete manalysis;
+        delete danalysis;
 
         #ifdef _HAVE_PETSC_
@@ -442 +456 @@
         _error_("PETSc needs to be installed");
         #endif
-}
+}/*}}}*/

issm/trunk-jpl/src/m/classes/damage.m

@@ -142 +142 @@
                                 md = checkfield(md,'fieldname','damage.spcdamage','timeseries',1);
                                 md = checkfield(md,'fieldname','damage.max_damage','<',1,'>=',0);
-                                md = checkfield(md,'fieldname','damage.stabilization','numel',[1],'values',[0 1 2]);
+                                md = checkfield(md,'fieldname','damage.stabilization','numel',[1],'values',[0 1 2 4]);
                                 md = checkfield(md,'fieldname','damage.maxiter','>=0',0);
                                 md = checkfield(md,'fieldname','damage.elementinterp','values',{'P1','P2'});
@@ -178 +178 @@
                                 fielddisplay(self,'max_damage','maximum possible damage (0<=max_damage<1)');
                                 
-                                fielddisplay(self,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG');
+                                fielddisplay(self,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG (not working), 4: flux corrected transport');
                                 fielddisplay(self,'maxiter','maximum number of non linear iterations');
                                 fielddisplay(self,'elementinterp','interpolation scheme for finite elements {''P1'',''P2''}');

issm/trunk-jpl/src/m/classes/damage.py

@@ -56 +56 @@
                         s+="%s\n" % fielddisplay(self,"max_damage","maximum possible damage (0<=max_damage<1)")
 
-                        s+="%s\n" % fielddisplay(self,"stabilization","0: no, 1: artificial_diffusivity, 2: SUPG")
+                        s+="%s\n" % fielddisplay(self,"stabilization","0: no, 1: artificial_diffusivity, 2: SUPG (not working), 4: Flux corrected transport")
                         s+="%s\n" % fielddisplay(self,"maxiter","maximum number of non linear iterations")
                         s+="%s\n" %     fielddisplay(self,"elementinterp","interpolation scheme for finite elements [''P1'',''P2'']")
@@ -126 +126 @@
                         md = checkfield(md,'fieldname','damage.law','numel',[1],'values',[0,1,2,3])
                         md = checkfield(md,'fieldname','damage.spcdamage','timeseries',1)
-                        md = checkfield(md,'fieldname','damage.stabilization','numel',[1],'values',[0,1,2])
+                        md = checkfield(md,'fieldname','damage.stabilization','numel',[1],'values',[0,1,2,4])
                         md = checkfield(md,'fieldname','damage.maxiter','>=0',0)
                         md = checkfield(md,'fieldname','damage.elementinterp','values',['P1','P2'])

issm/trunk-jpl/src/m/mech/mechanicalproperties.m

@@ -62 +62 @@
         nu(location)=B_bar(location);
 elseif isa(md.materials,'matdamageice')
+        % FIXME this needs to use md.damage.D or damage from a solution
         Zinv=md.materials.rheology_Z(index)*summation/3;
         location=find(second_inv~=0);

issm/trunk-jpl/src/m/mech/mechanicalproperties.py

@@ -77 +77 @@
         elif 'matdamageice' in md.materials.__module__:
                 print 'computing damage-dependent properties!'
+                # FIXME might be using damage from a solution, not the initial md.damage.D
                 Zinv=npy.dot(1-md.damage.D[index-1],summation/3.).reshape(-1,)
                 location=npy.nonzero(second_inv)