Changeset 24679

Timestamp:

04/01/20 00:38:51 (5 years ago)

Author:

yfischler

Message:

NEW: reusing allocated matrices in solutionsequence_nonlinear

Location:

issm/trunk-jpl/src/c

Files:

• modules/SystemMatricesx/SystemMatricesx.cpp (modified) (2 diffs)
• modules/SystemMatricesx/SystemMatricesx.h (modified) (1 diff)
• solutionsequences/solutionsequence_nonlinear.cpp (modified) (4 diffs)
• toolkits/issm/IssmAbsMat.h (modified) (1 diff)
• toolkits/issm/IssmMat.h (modified) (1 diff)
• toolkits/objects/Matrix.h (modified) (1 diff)
• toolkits/petsc/objects/PetscMat.cpp (modified) (1 diff)
• toolkits/petsc/objects/PetscMat.h (modified) (1 diff)

issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp

@@ -8 +8 @@
 #include "../AllocateSystemMatricesx/AllocateSystemMatricesx.h"
 
-void SystemMatricesx(Matrix<IssmDouble>** pKff, Matrix<IssmDouble>** pKfs, Vector<IssmDouble>** ppf, Vector<IssmDouble>** pdf, IssmDouble* pkmax,FemModel* femmodel){
+void SystemMatricesx(Matrix<IssmDouble>** pKff, Matrix<IssmDouble>** pKfs, Vector<IssmDouble>** ppf, Vector<IssmDouble>** pdf, IssmDouble* pkmax,FemModel* femmodel, bool isAllocated){
 
         /*intermediary: */
@@ -64 +64 @@
 
         /*Allocate stiffness matrices and load vector*/
-        AllocateSystemMatricesx(&Kff,&Kfs,&df,&pf,femmodel);
+        if(isAllocated) {
+                Kff  = *pKff;
+                Kfs  = *pKfs;
+                pf   = *ppf;
+                df   = *pdf;
+        }
+        else {
+                AllocateSystemMatricesx(&Kff,&Kfs,&df,&pf,femmodel);
+        }
 
         /*Display size*/

issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.h

@@ -9 +9 @@
 
 /* local prototypes: */
-void SystemMatricesx(Matrix<IssmDouble>** pKff, Matrix<IssmDouble>** pKfs, Vector<IssmDouble>** ppf, Vector<IssmDouble>** pdf, IssmDouble* pkmax,FemModel* femmodel);
+void SystemMatricesx(Matrix<IssmDouble>** pKff, Matrix<IssmDouble>** pKfs, Vector<IssmDouble>** ppf, Vector<IssmDouble>** pdf, IssmDouble* pkmax,FemModel* femmodel, bool isAllocated=false);
 
 #endif  /* _SYSTEMMATRICESX_H */

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

@@ -58 +58 @@
         InputUpdateFromSolutionx(femmodel,ug);
 
+        // allocate the matrices once and reuce them per iteration
+        AllocateSystemMatricesx(&Kff,&Kfs,&df,&pf,femmodel);
+
         for(;;){
 
@@ -64 +67 @@
                 delete ug;
 
-                SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
+                SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel, true);
                 CreateNodalConstraintsx(&ys,femmodel->nodes);
-                Reduceloadx(pf, Kfs, ys); delete Kfs;
+                Reduceloadx(pf, Kfs, ys);
                 femmodel->profiler->Start(SOLVER);
                 Solverx(&uf, Kff, pf, old_uf, df, femmodel->parameters);
@@ -73 +76 @@
                 Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete ys;
 
-                convergence(&converged,Kff,pf,uf,old_uf,eps_res,eps_rel,eps_abs); delete Kff; delete pf; delete df;
+                convergence(&converged,Kff,pf,uf,old_uf,eps_res,eps_rel,eps_abs);
                 InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
                 InputUpdateFromSolutionx(femmodel,ug);
@@ -102 +105 @@
                         break;
                 }
+
+                // Set the matrix entries to zero if we do an other iteration
+                Kff->SetZero();
+                Kfs->SetZero();
+                df->Set(0);
+                pf->Set(0);
         }
+
+        // delete matrices after the iteration loop
+        delete Kff; delete pf; delete df;
+        delete Kfs;
 
         if(VerboseConvergence()) _printf0_("\n   total number of iterations: " << count << "\n");

issm/trunk-jpl/src/c/toolkits/issm/IssmAbsMat.h

@@ -43 +43 @@
                 virtual void SetValues(int m,int* idxm,int n,int* idxn,doubletype* values,InsMode mode)=0;
                 virtual void Convert(MatrixType type)=0;
+                virtual void SetZero(void){};
                 #ifndef _HAVE_WRAPPERS_
                 virtual IssmAbsVec<IssmDouble>* Solve(IssmAbsVec<IssmDouble>* pf, Parameters* parameters)=0;

issm/trunk-jpl/src/c/toolkits/issm/IssmMat.h

@@ -239 +239 @@
                         matrix->convert(type);
                 }/*}}}*/
+                void SetZero(void){/*{{{*/
+                        matrix->SetZero();
+                }/*}}}*/
                 #ifndef _HAVE_WRAPPERS_
                 IssmVec<doubletype>* Solve(IssmVec<doubletype>* pf, Parameters* parameters){ /*{{{*/

issm/trunk-jpl/src/c/toolkits/objects/Matrix.h

@@ -307 +307 @@
                 }
                 /*}}}*/
-
+                /*
+                * sets all values to 0 but keeps the structure of a sparse matrix
+                */
+                void SetZero(void) {/*{{{*/
+                        if(type==PetscMatType){
+                                #ifdef _HAVE_PETSC_
+                                this->pmatrix->SetZero();
+                                #endif
+                        }
+                        else{
+                                this->imatrix->SetZero();
+                        }
+                }
+                /*}}}*/
 };
 

issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.cpp

@@ -193 +193 @@
 }
 /*}}}*/
+void PetscMat::SetZero(void){/*{{{*/
+        MatZeroEntries(this->matrix);
+}
+/*}}}*/

issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.h

@@ -52 +52 @@
                 void SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode);
                 void Convert(MatrixType type);
+                void SetZero(void);
 };