[cig-commits] r5709 - in long/3D/Gale/trunk/src/StgFEM: .
SLE/LinearAlgebra SLE/LinearAlgebra/src
SLE/LinearAlgebra/tests SLE/LinearAlgebra/tests/data
SLE/ProvidedSystems/AdvectionDiffusion/src
SLE/ProvidedSystems/Energy/src SLE/ProvidedSystems/StokesFlow/src
walter at geodynamics.org
walter at geodynamics.org
Mon Jan 8 23:50:53 PST 2007
Author: walter
Date: 2007-01-08 23:50:52 -0800 (Mon, 08 Jan 2007)
New Revision: 5709
Added:
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/tests/
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/tests/data/
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/tests/data/wallVC.xml
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/tests/makefile
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/tests/testSROpGenerator.c
Modified:
long/3D/Gale/trunk/src/StgFEM/
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Finalise.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Init.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MGOpGenerator.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Matrix.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MatrixSolver.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MatrixSolver.h
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MultigridSolver.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MultigridSolver.h
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrix.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrix.h
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrixSolver.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrixSolver.h
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScVector.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/SROpGenerator.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/SROpGenerator.h
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Vector.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/types.h
long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/AdvectionDiffusion/src/Multicorrector.c
long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/Energy/src/Energy_SLE_Solver.c
long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_PenaltySolver.c
long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_UzawaSolver.c
Log:
r927 at earth (orig r696): LukeHodkinson | 2007-01-08 22:32:51 -0800
* Adding multigrid as a sub-class of MatrixSolver. This
changes its interface to the more familiar component
style.
* Adding a test for simple regular inter-grid operators
for use with multigrid.
Property changes on: long/3D/Gale/trunk/src/StgFEM
___________________________________________________________________
Name: svk:merge
- 38867592-cf10-0410-9e16-a142ea72ac34:/cig:880
db209038-57f2-0310-97fa-b160e0ae9d04:/branches/decomp3d:695
db209038-57f2-0310-97fa-b160e0ae9d04:/trunk:669
+ 38867592-cf10-0410-9e16-a142ea72ac34:/cig:880
db209038-57f2-0310-97fa-b160e0ae9d04:/branches/decomp3d:696
db209038-57f2-0310-97fa-b160e0ae9d04:/trunk:669
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Finalise.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Finalise.c 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Finalise.c 2007-01-09 07:50:52 UTC (rev 5709)
@@ -33,7 +33,7 @@
#include <mpi.h>
#include "StGermain/StGermain.h"
-
+#include "Discretisation/Discretisation.h"
#include "LinearAlgebra.h"
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Init.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Init.c 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Init.c 2007-01-09 07:50:52 UTC (rev 5709)
@@ -33,7 +33,7 @@
#include <mpi.h>
#include "StGermain/StGermain.h"
-
+#include "Discretisation/Discretisation.h"
#include "LinearAlgebra.h"
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MGOpGenerator.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MGOpGenerator.c 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MGOpGenerator.c 2007-01-09 07:50:52 UTC (rev 5709)
@@ -35,6 +35,7 @@
#include <mpi.h>
#include "StGermain/StGermain.h"
+#include "Discretisation/Discretisation.h"
#include "LinearAlgebra.h"
@@ -55,6 +56,7 @@
/* Virtual info */
self->setMultigridSolverFunc = setMultigridSolverFunc;
+ self->hasExpiredFunc = hasExpiredFunc;
self->generateFunc = generateFunc;
/* MGOpGenerator info */
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Matrix.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Matrix.c 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Matrix.c 2007-01-09 07:50:52 UTC (rev 5709)
@@ -35,6 +35,7 @@
#include <mpi.h>
#include "StGermain/StGermain.h"
+#include "Discretisation/Discretisation.h"
#include "LinearAlgebra.h"
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MatrixSolver.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MatrixSolver.c 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MatrixSolver.c 2007-01-09 07:50:52 UTC (rev 5709)
@@ -35,6 +35,7 @@
#include <mpi.h>
#include "StGermain/StGermain.h"
+#include "Discretisation/Discretisation.h"
#include "LinearAlgebra.h"
@@ -67,6 +68,7 @@
self->getSolveStatusFunc = getSolveStatusFunc;
self->getIterationsFunc = getIterationsFunc;
self->getMaxIterationsFunc = getMaxIterationsFunc;
+ self->getResidualNormFunc = getResidualNormFunc;
/* MatrixSolver info */
_MatrixSolver_Init( self );
@@ -76,6 +78,15 @@
void _MatrixSolver_Init( MatrixSolver* self ) {
assert( self && Stg_CheckType( self, MatrixSolver ) );
+
+ self->comm = MPI_COMM_WORLD;
+ self->matrix = NULL;
+ self->inversion = NULL;
+ self->residual = NULL;
+ self->expiredResidual = True;
+
+ self->curRHS = NULL;
+ self->curSolution = NULL;
}
@@ -151,6 +162,7 @@
self->curRHS = rhs;
self->curSolution = solution;
+ self->expiredResidual = True;
}
@@ -179,6 +191,11 @@
assert( self && Stg_CheckType( self, MatrixSolver ) );
+ if( !self->residual || self->expiredResidual ) {
+ MatrixSolver_CalcResidual( self );
+ self->expiredResidual = False;
+ }
+
return self->residual;
}
@@ -186,3 +203,15 @@
/*----------------------------------------------------------------------------------------------------------------------------------
** Private Functions
*/
+
+void MatrixSolver_CalcResidual( MatrixSolver* self ) {
+ assert( self && Stg_CheckType( self, MatrixSolver ) );
+ assert( self->curSolution && self->curRHS && self->matrix );
+
+ if( !self->residual ) {
+ Vector_Duplicate( self->curSolution, (void**)&self->residual );
+ Vector_SetLocalSize( self->residual, Vector_GetLocalSize( self->curSolution ) );
+ }
+ Matrix_Multiply( self->matrix, self->curSolution, self->residual );
+ Vector_ScaleAdd( self->residual, -1.0, self->curRHS );
+}
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MatrixSolver.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MatrixSolver.h 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MatrixSolver.h 2007-01-09 07:50:52 UTC (rev 5709)
@@ -58,6 +58,7 @@
typedef MatrixSolver_Status (MatrixSolver_GetSolveStatusFunc)( void* matrixSolver );
typedef unsigned (MatrixSolver_GetIterationsFunc)( void* matrixSolver );
typedef unsigned (MatrixSolver_GetMaxIterationsFunc)( void* matrixSolver );
+ typedef double (MatrixSolver_GetResidualNormFunc)( void* matrixSolver );
/** MatrixSolver class contents */
#define __MatrixSolver \
@@ -78,12 +79,14 @@
MatrixSolver_GetSolveStatusFunc* getSolveStatusFunc; \
MatrixSolver_GetIterationsFunc* getIterationsFunc; \
MatrixSolver_GetMaxIterationsFunc* getMaxIterationsFunc; \
+ MatrixSolver_GetResidualNormFunc* getResidualNormFunc; \
\
/* MatrixSolver info */ \
MPI_Comm comm; \
Matrix* matrix; \
Matrix* inversion; \
Vector* residual; \
+ Bool expiredResidual; \
\
Vector* curRHS; \
Vector* curSolution;
@@ -108,7 +111,8 @@
\
MatrixSolver_GetSolveStatusFunc* getSolveStatusFunc, \
MatrixSolver_GetIterationsFunc* getIterationsFunc, \
- MatrixSolver_GetMaxIterationsFunc* getMaxIterationsFunc
+ MatrixSolver_GetMaxIterationsFunc* getMaxIterationsFunc, \
+ MatrixSolver_GetResidualNormFunc* getResidualNormFunc
#define MATRIXSOLVER_PASSARGS \
STG_COMPONENT_PASSARGS, \
@@ -124,7 +128,8 @@
\
getSolveStatusFunc, \
getIterationsFunc, \
- getMaxIterationsFunc
+ getMaxIterationsFunc, \
+ getResidualNormFunc
MatrixSolver* _MatrixSolver_New( MATRIXSOLVER_DEFARGS );
void _MatrixSolver_Init( MatrixSolver* self );
@@ -178,6 +183,9 @@
#define MatrixSolver_GetMaxIterations( self ) \
VirtualCall( self, getMaxIterationsFunc, self )
+ #define MatrixSolver_GetResidualNorm( self ) \
+ VirtualCall( self, getResidualNormFunc, self )
+
/*--------------------------------------------------------------------------------------------------------------------------
** Public functions
*/
@@ -190,4 +198,6 @@
** Private Member functions
*/
+ void MatrixSolver_CalcResidual( MatrixSolver* self );
+
#endif /* __StgFEM_SLE_LinearAlgebra_MatrixSolver_h__ */
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MultigridSolver.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MultigridSolver.c 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MultigridSolver.c 2007-01-09 07:50:52 UTC (rev 5709)
@@ -35,6 +35,7 @@
#include <mpi.h>
#include "StGermain/StGermain.h"
+#include "Discretisation/Discretisation.h"
#include "LinearAlgebra.h"
@@ -70,7 +71,8 @@
MultigridSolver_Setup,
MultigridSolver_GetSolveStatus,
MultigridSolver_GetIterations,
- MultigridSolver_GetMaxIterations );
+ MultigridSolver_GetMaxIterations,
+ MultigridSolver_GetResidualNorm );
}
MultigridSolver* _MultigridSolver_New( MULTIGRIDSOLVER_DEFARGS ) {
@@ -96,6 +98,10 @@
self->opGen = NULL;
self->ready = False;
+ self->curIt = 0;
+ self->maxIts = 500;
+ self->useInitial = False;
+
#ifdef HAVE_PETSC
self->outerSolver = (MatrixSolver*)PETScMatrixSolver_New( "" );
PETScMatrixSolver_SetKSPType( self->outerSolver, PETScMatrixSolver_KSPType_Richardson );
@@ -139,9 +145,18 @@
void _MultigridSolver_Construct( void* matrixSolver, Stg_ComponentFactory* cf, void* data ) {
MultigridSolver* self = (MultigridSolver*)matrixSolver;
+ unsigned nLevels;
assert( self && Stg_CheckType( self, MultigridSolver ) );
assert( cf );
+
+ _MatrixSolver_Construct( self, cf, data );
+
+ nLevels = Stg_ComponentFactory_GetUnsignedInt( cf, self->name, "levels", 1 );
+ MultigridSolver_SetLevels( self, nLevels );
+ self->opGen = Stg_ComponentFactory_ConstructByKey( cf, self->name, "opGenerator", MGOpGenerator,
+ True, data );
+ MGOpGenerator_SetMultigridSolver( self->opGen, self );
}
void _MultigridSolver_Build( void* matrixSolver, void* data ) {
@@ -204,7 +219,6 @@
MultigridSolver* self = (MultigridSolver*)matrixSolver;
Vector* rhs = (Vector*)_rhs;
Vector* solution = (Vector*)_solution;
- unsigned curIt = 0;
assert( self && Stg_CheckType( self, MultigridSolver ) );
assert( rhs && Stg_CheckType( rhs, Vector ) );
@@ -215,7 +229,7 @@
MultigridSolver_Setup( self, rhs, solution );
self->curIt = 0;
- while( MatrixSolver_GetSolveStatus( self ) == MatrixSolver_Status_Iterating && curIt++ < self->maxIts ) {
+ while( MatrixSolver_GetSolveStatus( self ) == MatrixSolver_Status_Iterating && self->curIt < self->maxIts ) {
MultigridSolver_LevelCycle( self, self->nLevels - 1, rhs, solution );
self->curIt++;
}
@@ -232,7 +246,7 @@
assert( rhs && Stg_CheckType( rhs, Vector ) );
assert( solution && Stg_CheckType( solution, Vector ) );
- MatrixSolver_Setup( self, rhs, solution );
+ _MatrixSolver_Setup( self, rhs, solution );
if( self->opGen )
rebuildOps = MGOpGenerator_HasExpired( self->opGen );
@@ -258,9 +272,13 @@
Matrix* mat;
mat = MatrixSolver_GetMatrix( level->downSolver );
- if( !mat )
+ if( !mat ) {
Matrix_Duplicate( self->matrix, (void**)&mat );
- MultigridSolver_RestrictMatrix( self, level, mat );
+ Matrix_AssemblyBegin( mat );
+ Matrix_AssemblyEnd( mat );
+ }
+ MultigridSolver_RestrictMatrix( self, self->levels + l_i + 1, mat );
+
MatrixSolver_SetMatrix( level->downSolver, mat );
if( level->upSolver != level->downSolver )
MatrixSolver_SetMatrix( level->upSolver, mat );
@@ -271,7 +289,7 @@
MatrixSolver_SetMatrix( level->upSolver, self->matrix );
}
- MatrixSolver_SetUseInitialSolution( level->downSolver, True );
+ MatrixSolver_SetUseInitialSolution( level->downSolver, False );
if( level->upSolver != level->downSolver )
MatrixSolver_SetUseInitialSolution( level->upSolver, True );
}
@@ -302,6 +320,8 @@
MatrixSolver_Solve( self->outerSolver, self->curRHS, self->curSolution );
outerStatus = MatrixSolver_GetSolveStatus( self->outerSolver );
if( outerStatus == MatrixSolver_Status_ConvergedIterations )
+ return MatrixSolver_Status_ConvergedIterations;
+ else if( outerStatus == MatrixSolver_Status_DivergedIterations )
return MatrixSolver_Status_Iterating;
else
return outerStatus;
@@ -323,13 +343,22 @@
return self->maxIts;
}
+double MultigridSolver_GetResidualNorm( void* matrixSolver ) {
+ MultigridSolver* self = (MultigridSolver*)matrixSolver;
+ assert( self && Stg_CheckType( self, MultigridSolver ) );
+
+ return MatrixSolver_GetResidualNorm( self->outerSolver );
+}
+
+
/*--------------------------------------------------------------------------------------------------------------------------
** Public Functions
*/
void MultigridSolver_SetLevels( void* matrixSolver, unsigned nLevels ) {
MultigridSolver* self = (MultigridSolver*)matrixSolver;
+ unsigned nProcs;
unsigned l_i;
assert( self && Stg_CheckType( self, MultigridSolver ) );
@@ -337,26 +366,47 @@
MultigridSolver_DestructLevels( self );
self->nLevels = nLevels;
self->levels = AllocArray( MultigridSolver_Level, nLevels );
+ MPI_Comm_size( self->comm, (int*)&nProcs );
for( l_i = 0; l_i < nLevels; l_i++ ) {
MultigridSolver_Level* level = self->levels + l_i;
level->nDownIts = 1;
- level->nUpIts = 1;
level->nCycles = 1;
if( l_i == 0 ) {
+ level->nUpIts = 0;
+
/* Create direct solve. */
- level->downSolver = NULL;
+#ifdef HAVE_PETSC
+ level->downSolver = (MatrixSolver*)PETScMatrixSolver_New( "" );
+ PETScMatrixSolver_SetKSPType( level->downSolver,
+ PETScMatrixSolver_KSPType_PreOnly );
+ if( nProcs == 1 ) {
+ PETScMatrixSolver_SetPCType( level->downSolver,
+ PETScMatrixSolver_PCType_LU );
+ }
+ else {
+ PETScMatrixSolver_SetPCType( level->downSolver,
+ PETScMatrixSolver_PCType_RedundantLU );
+ }
+#else
+ self->downSolver = NULL;
+#endif
}
- else if( l_i == 1 ) {
+ else {
+ level->nUpIts = 1;
+
/* Create smoother. */
+#ifdef HAVE_PETSC
+ level->downSolver = (MatrixSolver*)PETScMatrixSolver_New( "" );
+ PETScMatrixSolver_SetKSPType( level->downSolver,
+ PETScMatrixSolver_KSPType_Richardson );
+ PETScMatrixSolver_SetPCType( level->downSolver,
+ PETScMatrixSolver_PCType_SOR );
+#else
level->downSolver = NULL;
+#endif
}
- else {
- /* Copy previous smoother. */
- level->downSolver = self->levels[l_i - 1].downSolver;
- Stg_Class_AddRef( level->downSolver );
- }
level->upSolver = level->downSolver;
Stg_Class_AddRef( level->upSolver );
@@ -375,7 +425,8 @@
if( self->levels[level].R != R )
self->ready = False;
- Stg_Class_RemoveRef( self->levels[level].R );
+ if( self->levels[level].R )
+ Stg_Class_RemoveRef( self->levels[level].R );
self->levels[level].R = R;
if( R )
Stg_Class_AddRef( R );
@@ -391,7 +442,8 @@
if( self->levels[level].P != P )
self->ready = False;
- Stg_Class_RemoveRef( self->levels[level].P );
+ if( self->levels[level].P )
+ Stg_Class_RemoveRef( self->levels[level].P );
self->levels[level].P = P;
if( P )
Stg_Class_AddRef( P );
@@ -521,13 +573,21 @@
MultigridSolver_SetLevelDownSolver( self, 0, solver );
}
+unsigned MultigridSolver_GetNumLevels( void* matrixSolver ) {
+ MultigridSolver* self = (MultigridSolver*)matrixSolver;
+ assert( self && Stg_CheckType( self, MultigridSolver ) );
+
+ return self->nLevels;
+}
+
+
/*----------------------------------------------------------------------------------------------------------------------------------
** Private Functions
*/
void MultigridSolver_RestrictMatrix( MultigridSolver* self, MultigridSolver_Level* level, Matrix* dstMatrix ) {
- Matrix *srcMat, *R;
+ Matrix *srcMat;
unsigned nSrcRows, nSrcCols;
unsigned nRRows, nRCols;
@@ -535,13 +595,13 @@
assert( level );
assert( dstMatrix );
- srcMat = MatrixSolver_GetMatrix( level->downSolver ); assert( srcMat );
+ srcMat = MatrixSolver_GetMatrix( level->downSolver );
Matrix_GetGlobalSize( srcMat, &nSrcRows, &nSrcCols );
- Matrix_GetGlobalSize( R, &nRRows, &nRCols );
+ Matrix_GetGlobalSize( level->R, &nRRows, &nRCols );
if( nRCols == nSrcRows )
- Matrix_PAPt( dstMatrix, R, srcMat );
+ Matrix_PAPt( srcMat, level->R, dstMatrix );
else
- Matrix_PtAP( dstMatrix, R, srcMat );
+ Matrix_PtAP( srcMat, level->R, dstMatrix );
}
void MultigridSolver_LevelCycle( MultigridSolver* self, unsigned levelInd, Vector* rhs, Vector* solution ) {
@@ -552,18 +612,31 @@
level = self->levels + levelInd;
- if( level->nDownIts )
+ if( level->nDownIts ) {
+ if( levelInd == self->nLevels - 1 && self->curIt > 0 )
+ MatrixSolver_SetUseInitialSolution( level->downSolver, True );
MatrixSolver_Solve( level->downSolver, rhs, solution );
+ }
if( levelInd > 0 ) {
unsigned c_i;
- Matrix_Multiply( level->R, MatrixSolver_GetResidual( level->downSolver ), level->nextRHS );
+ if( level->R == level->P )
+ Matrix_TransposeMultiply( level->R, MatrixSolver_GetResidual( level->downSolver ), level->nextRHS );
+ else
+ Matrix_Multiply( level->R, MatrixSolver_GetResidual( level->downSolver ), level->nextRHS );
for( c_i = 0; c_i < level->nCycles; c_i++ )
- MultigridSolver_LevelCycle( self, levelInd - 1, level->nextSol, level->nextRHS );
+ MultigridSolver_LevelCycle( self, levelInd - 1, level->nextRHS, level->nextSol );
+ Matrix_MultiplyAdd( level->P, level->nextSol, solution, solution );
}
- Matrix_MultiplyAdd( level->P, level->nextSol, solution, solution );
+ if( level->nUpIts ) {
+ if( level->upSolver == level->downSolver )
+ MatrixSolver_SetUseInitialSolution( level->upSolver, True );
+ MatrixSolver_Solve( level->upSolver, rhs, solution );
+ if( level->upSolver == level->downSolver )
+ MatrixSolver_SetUseInitialSolution( level->upSolver, False );
+ }
}
void MultigridSolver_Destruct( MultigridSolver* self ) {
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MultigridSolver.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MultigridSolver.h 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/MultigridSolver.h 2007-01-09 07:50:52 UTC (rev 5709)
@@ -118,6 +118,7 @@
MatrixSolver_Status MultigridSolver_GetSolveStatus( void* matrixSolver );
unsigned MultigridSolver_GetIterations( void* matrixSolver );
unsigned MultigridSolver_GetMaxIterations( void* matrixSolver );
+ double MultigridSolver_GetResidualNorm( void* matrixSolver );
/*--------------------------------------------------------------------------------------------------------------------------
** Public functions
@@ -138,6 +139,8 @@
void MultigridSolver_SetAllSolver( void* matrixSolver, MatrixSolver* solver );
void MultigridSolver_SetCoarseSolver( void* matrixSolver, MatrixSolver* solver );
+ unsigned MultigridSolver_GetNumLevels( void* matrixSolver );
+
/*--------------------------------------------------------------------------------------------------------------------------
** Private Member functions
*/
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrix.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrix.c 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrix.c 2007-01-09 07:50:52 UTC (rev 5709)
@@ -35,6 +35,7 @@
#include <mpi.h>
#include "StGermain/StGermain.h"
+#include "Discretisation/Discretisation.h"
#include "LinearAlgebra.h"
@@ -108,9 +109,12 @@
}
void _PETScMatrix_Init( PETScMatrix* self ) {
+ PetscErrorCode ec;
+
assert( self && Stg_CheckType( self, PETScMatrix ) );
- self->petscMat = PETSC_NULL;
+ ec = MatCreate( self->comm, &self->petscMat );
+ CheckPETScError( ec );
}
@@ -429,15 +433,18 @@
PETScMatrix* P = (PETScMatrix*)_P;
PETScMatrix* dstMatrix = (PETScMatrix*)_dstMatrix;
PetscErrorCode ec;
+#if 0
double fillRatio;
double nRowsA, nRowsC;
double nzA, nzP, nzC;
MatInfo mInfo;
+#endif
assert( self && Stg_CheckType( self, PETScMatrix ) );
assert( P && Stg_CheckType( P, PETScMatrix ) );
assert( dstMatrix && Stg_CheckType( dstMatrix, PETScMatrix ) );
+#if 0
MatGetInfo( self->petscMat, MAT_GLOBAL_SUM, &mInfo );
nRowsA = mInfo.rows_global;
nzA = mInfo.nz_used;
@@ -448,11 +455,12 @@
nzC = (nRowsC / nRowsA) * nzA;
fillRatio = nzC / (nzA + nzP);
+#endif
if( dstMatrix->petscMat != PETSC_NULL )
ec = MatDestroy( dstMatrix->petscMat );
CheckPETScError( ec );
- ec = MatPtAP( self->petscMat, P->petscMat, MAT_INITIAL_MATRIX, (PetscScalar)fillRatio, &dstMatrix->petscMat );
+ ec = MatPtAP( self->petscMat, P->petscMat, MAT_INITIAL_MATRIX, 1.0, &dstMatrix->petscMat );
CheckPETScError( ec );
}
@@ -461,15 +469,18 @@
PETScMatrix* matrix0 = (PETScMatrix*)_matrix0;
PETScMatrix* dstMatrix = (PETScMatrix*)_dstMatrix;
PetscErrorCode ec;
+#if 0
double fillRatio;
double nRowsA, nRowsC;
double nzA, nzP, nzC;
MatInfo mInfo;
+#endif
assert( self && Stg_CheckType( self, PETScMatrix ) );
assert( matrix0 && Stg_CheckType( matrix0, PETScMatrix ) );
assert( dstMatrix && Stg_CheckType( dstMatrix, PETScMatrix ) );
+#if 0
MatGetInfo( self->petscMat, MAT_GLOBAL_SUM, &mInfo );
nRowsA = mInfo.rows_global;
nzA = mInfo.nz_used;
@@ -480,11 +491,12 @@
nzC = (nRowsC / nRowsA) * nzA;
fillRatio = nzC / (nzA + nzP);
+#endif
if( dstMatrix->petscMat != PETSC_NULL )
ec = MatDestroy( dstMatrix->petscMat );
CheckPETScError( ec );
- ec = MatPtAP( self->petscMat, matrix0->petscMat, MAT_INITIAL_MATRIX, (PetscScalar)fillRatio, &dstMatrix->petscMat );
+ ec = MatPtAP( self->petscMat, matrix0->petscMat, MAT_INITIAL_MATRIX, 1.0, &dstMatrix->petscMat );
CheckPETScError( ec );
}
@@ -672,7 +684,17 @@
CheckPETScError( ec );
}
+void PETScMatrix_Draw( void* matrix ) {
+ PETScMatrix* self = (PETScMatrix*)matrix;
+ PetscErrorCode ec;
+ assert( self && Stg_CheckType( self, PETScMatrix ) );
+
+ ec = MatView( self->petscMat, PETSC_VIEWER_DRAW_WORLD );
+ CheckPETScError( ec );
+}
+
+
/*----------------------------------------------------------------------------------------------------------------------------------
** Private Functions
*/
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrix.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrix.h 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrix.h 2007-01-09 07:50:52 UTC (rev 5709)
@@ -129,6 +129,7 @@
void PETScMatrix_SetNonZeroStructure( void* matrix, unsigned nNonZeros,
unsigned* diagonalNonZeroIndices, unsigned* offDiagonalNonZeroIndices);
+ void PETScMatrix_Draw( void* matrix );
/*--------------------------------------------------------------------------------------------------------------------------
** Private Member functions
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrixSolver.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrixSolver.c 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrixSolver.c 2007-01-09 07:50:52 UTC (rev 5709)
@@ -35,6 +35,7 @@
#include <mpi.h>
#include "StGermain/StGermain.h"
+#include "Discretisation/Discretisation.h"
#include "LinearAlgebra.h"
@@ -70,7 +71,8 @@
_MatrixSolver_Setup,
PETScMatrixSolver_GetSolveStatus,
PETScMatrixSolver_GetIterations,
- PETScMatrixSolver_GetMaxIterations );
+ PETScMatrixSolver_GetMaxIterations,
+ PETScMatrixSolver_GetResidualNorm );
}
PETScMatrixSolver* _PETScMatrixSolver_New( PETSCMATRIXSOLVER_DEFARGS ) {
@@ -213,6 +215,7 @@
assert( solution && Stg_CheckType( solution, PETScVector ) );
assert( rhs && Stg_CheckType( rhs, PETScVector ) );
+ MatrixSolver_Setup( self, rhs, solution );
KSPSolve( self->ksp, rhs->petscVec, solution->petscVec );
}
@@ -249,7 +252,18 @@
return (unsigned)nIts;
}
+double PETScMatrixSolver_GetResidualNorm( void* matrixSolver ) {
+ PETScMatrixSolver* self = (PETScMatrixSolver*)matrixSolver;
+ PetscScalar norm;
+ assert( self && Stg_CheckType( self, PETScMatrixSolver ) );
+
+ KSPGetResidualNorm( self->ksp, &norm );
+
+ return (double)norm;
+}
+
+
/*--------------------------------------------------------------------------------------------------------------------------
** Public Functions
*/
@@ -281,28 +295,44 @@
assert( self && Stg_CheckType( self, PETScMatrixSolver ) );
KSPGetPC( self->ksp, &pc );
- switch( type ) {
- case PETScMatrixSolver_PCType_Jacobi:
+ if( type == PETScMatrixSolver_PCType_Jacobi ) {
PCSetType( pc, PCJACOBI );
- break;
- case PETScMatrixSolver_PCType_BlockJacobi:
+ }
+ else if( type == PETScMatrixSolver_PCType_BlockJacobi ) {
PCSetType( pc, PCBJACOBI );
- break;
- case PETScMatrixSolver_PCType_SOR:
+ }
+ else if( type == PETScMatrixSolver_PCType_SOR ) {
PCSetType( pc, PCSOR );
- break;
- case PETScMatrixSolver_PCType_LU:
+ }
+#if 0
+ else if( type == PETScMatrixSolver_PCType_ParallelSOR ) {
+ KSP* subKSPs;
+ PC subPC;
+
+ PCSetType( pc, BJACOBI );
+ PCSetUp( pc );
+ PCBJacobiGetSubKSP( pc, PETSC_NULL, PETSC_NULL, &subKSPs );
+ KSPGetPC( subKSPs[0], &subPC );
+ PCSetType( subPC, PCSOR );
+ }
+#endif
+ else if( type == PETScMatrixSolver_PCType_LU ) {
PCSetType( pc, PCLU );
- break;
- case PETScMatrixSolver_PCType_ILU:
+ }
+ else if( type == PETScMatrixSolver_PCType_RedundantLU ) {
+ PCSetType( pc, PCREDUNDANT );
+ PCRedundantGetPC( pc, &pc );
+ PCSetType( pc, PCLU );
+ }
+ else if( type == PETScMatrixSolver_PCType_ILU ) {
PCSetType( pc, PCILU );
- break;
- case PETScMatrixSolver_PCType_None:
+ }
+ else if( type == PETScMatrixSolver_PCType_None ) {
PCSetType( pc, PCNONE );
- break;
- default:
- assert( 0 );
}
+ else {
+ assert( 0 ); /* TODO */
+ }
}
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrixSolver.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrixSolver.h 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScMatrixSolver.h 2007-01-09 07:50:52 UTC (rev 5709)
@@ -96,6 +96,7 @@
MatrixSolver_Status PETScMatrixSolver_GetSolveStatus( void* matrixSolver );
unsigned PETScMatrixSolver_GetIterations( void* matrixSolver );
unsigned PETScMatrixSolver_GetMaxIterations( void* matrixSolver );
+ double PETScMatrixSolver_GetResidualNorm( void* matrixSolver );
/*--------------------------------------------------------------------------------------------------------------------------
** Public functions
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScVector.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScVector.c 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/PETScVector.c 2007-01-09 07:50:52 UTC (rev 5709)
@@ -35,6 +35,7 @@
#include <mpi.h>
#include "StGermain/StGermain.h"
+#include "Discretisation/Discretisation.h"
#include "LinearAlgebra.h"
@@ -431,7 +432,8 @@
*dstVector = (void*)self->_defaultConstructor( "" );
ec = VecDestroy( ((PETScVector*)*dstVector)->petscVec );
CheckPETScError( ec );
- ec = VecDuplicate( self->petscVec, &((PETScVector*)*dstVector)->petscVec );
+
+ ec = VecCreate( self->comm, &((PETScVector*)*dstVector)->petscVec );
CheckPETScError( ec );
}
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/SROpGenerator.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/SROpGenerator.c 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/SROpGenerator.c 2007-01-09 07:50:52 UTC (rev 5709)
@@ -35,6 +35,7 @@
#include <mpi.h>
#include "StGermain/StGermain.h"
+#include "Discretisation/Discretisation.h"
#include "LinearAlgebra.h"
@@ -82,6 +83,13 @@
void _SROpGenerator_Init( SROpGenerator* self ) {
assert( self && Stg_CheckType( self, SROpGenerator ) );
+
+ self->fineVar = NULL;
+ self->fineEqNum = NULL;
+ self->meshes = NULL;
+ self->topMaps = NULL;
+ self->eqNums = NULL;
+ self->nLocalEqNums = NULL;
}
@@ -113,10 +121,14 @@
}
void _SROpGenerator_Construct( void* srOpGenerator, Stg_ComponentFactory* cf, void* data ) {
- SROpGenerator* self = (SROpGenerator*)srOpGenerator;
+ SROpGenerator* self = (SROpGenerator*)srOpGenerator;
+ FeVariable* var;
assert( self && Stg_CheckType( self, SROpGenerator ) );
- assert( cf );
+
+ var = Stg_ComponentFactory_ConstructByKey( cf, self->name, "fineVariable", FeVariable,
+ True, data );
+ SROpGenerator_SetFineVariable( self, var );
}
void _SROpGenerator_Build( void* srOpGenerator, void* data ) {
@@ -143,28 +155,11 @@
SROpGenerator* self = (SROpGenerator*)srOpGenerator;
assert( self && Stg_CheckType( self, SROpGenerator ) );
-#if 0
- assert( ExtensionManager_GetHandle( self->feMesh->info, "vertexGrid" ) != -1 );
- eqNum = self->eqNum;
- feMesh = eqNum->feMesh;
- dofLayout = eqNum->dofLayout;
- nDims = Mesh_GetDimSize( feMesh );
- nLocalNodes = FeMesh_GetNodeLocalSize( feMesh );
- nodeGrid = *(Grid**)ExtensionManager_Get( self->feMesh->info, self->feMesh,
- ExtensionManager_GetHandle( self->feMesh->info, "vertexGrid" ) );
- globalInds = AllocArray( unsigned, nDims );
-
- for( n_i = 0; n_i < nLocalNodes; n_i++ ) {
- insist( FeMesh_NodeDomainToGlobal( feMesh, n_i, &globalNode ) );
- Grid_Lift( nodeGrid, globalNode, globalInds );
-
- for( l_i = nLevels - 1; l_i < nLevels; l_i-- )
- SROpGenerator_ProcessLevel( self, l_i, globalNode, globalInds );
- }
-
- FreeArray( globalInds );
-#endif
+ self->fineEqNum = self->fineVar->eqNum;
+ SROpGenerator_GenMeshes( self );
+ SROpGenerator_GenOps( self );
+ SROpGenerator_DestructMeshes( self );
}
@@ -172,35 +167,440 @@
** Public Functions
*/
+void SROpGenerator_SetFineVariable( void* srOpGenerator, void* _variable ) {
+ SROpGenerator* self = (SROpGenerator*)srOpGenerator;
+ FeVariable* variable = (FeVariable*)_variable;
+ assert( self && Stg_CheckType( self, SROpGenerator ) );
+ assert( !variable || Stg_CheckType( variable, FeVariable ) );
+
+ SROpGenerator_DestructMeshes( self );
+ self->fineVar = variable;
+ self->fineEqNum = NULL;
+}
+
+
/*----------------------------------------------------------------------------------------------------------------------------------
** Private Functions
*/
-#if 0
-void SROpGenerator_ProcessLevel( SROpGenerator* self, unsigned level, unsigned* ) {
- SROpGenerator* self = (SROpGenerator*)srOpGenerator;
+void SROpGenerator_GenMeshes( SROpGenerator* self ) {
+ unsigned nLevels;
+ unsigned l_i;
assert( self && Stg_CheckType( self, SROpGenerator ) );
- assert( ExtensionManager_GetHandle( self->feMesh->info, "vertexGrid" ) != -1 );
- eqNum = self->eqNum;
- feMesh = eqNum->feMesh;
- dofLayout = eqNum->dofLayout;
- nDims = Mesh_GetDimSize( feMesh );
- nLocalNodes = FeMesh_GetNodeLocalSize( feMesh );
- nodeGrid = *(Grid**)ExtensionManager_Get( self->feMesh->info, self->feMesh,
- ExtensionManager_GetHandle( self->feMesh->info, "vertexGrid" ) );
- globalInds = AllocArray( unsigned, nDims );
+ nLevels = MultigridSolver_GetNumLevels( self->solver );
+ self->meshes = AllocArray( Mesh*, nLevels );
+ memset( self->meshes, 0, nLevels * sizeof(Mesh*) );
+ self->topMaps = AllocArray( unsigned*, nLevels );
+ memset( self->topMaps, 0, nLevels * sizeof(unsigned*) );
+ self->eqNums = AllocArray( unsigned**, nLevels );
+ memset( self->eqNums, 0, nLevels * sizeof(unsigned**) );
+ self->nLocalEqNums = AllocArray( unsigned, nLevels );
+ memset( self->nLocalEqNums, 0, nLevels * sizeof(unsigned) );
+ self->eqNumBases = AllocArray( unsigned, nLevels );
+ memset( self->eqNumBases, 0, nLevels * sizeof(unsigned) );
+ self->meshes[nLevels - 1] = (Mesh*)self->fineEqNum->feMesh;
+ self->eqNumBases[nLevels - 1] = self->fineEqNum->firstOwnedEqNum;
+ for( l_i = nLevels - 2; l_i < nLevels; l_i-- ) {
+ SROpGenerator_GenLevelMesh( self, l_i );
+ SROpGenerator_GenLevelTopMap( self, l_i );
+ SROpGenerator_GenLevelEqNums( self, l_i );
+ }
+}
+
+void SROpGenerator_GenLevelMesh( SROpGenerator* self, unsigned level ) {
+ Stream* errorStream = Journal_Register( ErrorStream_Type, "SROpGenerator::GenLevelMesh" );
+ Mesh *fMesh, *cMesh;
+ CartesianGenerator *fGen, *cGen;
+ unsigned nDims;
+ unsigned* cSize;
+ unsigned d_i;
+
+ assert( self && Stg_CheckType( self, SROpGenerator ) );
+ assert( self->meshes );
+ assert( level < MultigridSolver_GetNumLevels( self->solver ) );
+
+ fMesh = self->meshes[level + 1];
+ nDims = Mesh_GetDimSize( fMesh );
+ fGen = (CartesianGenerator*)fMesh->generator;
+ Journal_Firewall( fGen && !strcmp( fGen->type, CartesianGenerator_Type ),
+ errorStream,
+ "\n" \
+ "****************************************************************\n" \
+ "* Error: Simple regular multigrid operator generation requires *\n" \
+ "* a fine mesh that has been generated with a *\n" \
+ "* cartesian generator. *\n" \
+ "****************************************************************\n" \
+ "\n" );
+
+ cGen = CartesianGenerator_New( "" );
+ CartesianGenerator_SetDimSize( cGen, nDims );
+ cSize = AllocArray( unsigned, nDims );
+ for( d_i = 0; d_i < nDims; d_i++ )
+ cSize[d_i] = fGen->elGrid->sizes[d_i] / 2;
+ CartesianGenerator_SetTopologyParams( cGen, cSize, fGen->maxDecompDims, fGen->minDecomp, fGen->maxDecomp );
+ CartesianGenerator_SetGeometryParams( cGen, fGen->crdMin, fGen->crdMax );
+ CartesianGenerator_SetShadowDepth( cGen, 0 );
+ FreeArray( cSize );
+
+ cMesh = (Mesh*)FeMesh_New( "" );
+ Mesh_SetGenerator( cMesh, cGen );
+ FeMesh_SetElementFamily( cMesh, ((FeMesh*)fMesh)->feElFamily );
+ Build( cMesh, NULL, False );
+ Initialise( cMesh, NULL, False );
+ self->meshes[level] = cMesh;
+}
+
+void SROpGenerator_GenLevelTopMap( SROpGenerator* self, unsigned level ) {
+ Stream* errorStream = Journal_Register( ErrorStream_Type, "SROpGenerator::GenLevelTopMap" );
+ Mesh *fMesh, *cMesh;
+ unsigned nDomainNodes;
+ unsigned nLevels;
+ unsigned nDims;
+ unsigned nearest;
+ double *cVert, *fVert;
+ unsigned d_i, n_i;
+
+ assert( self && Stg_CheckType( self, SROpGenerator ) );
+ assert( self->meshes );
+ assert( level < MultigridSolver_GetNumLevels( self->solver ) );
+
+ fMesh = self->meshes[level + 1];
+ cMesh = self->meshes[level];
+ nDims = Mesh_GetDimSize( fMesh );
+ nLevels = MultigridSolver_GetNumLevels( self->solver );
+ nDomainNodes = Mesh_GetDomainSize( cMesh, MT_VERTEX );
+ self->topMaps[level] = ReallocArray( self->topMaps[level], unsigned, nDomainNodes );
+ for( n_i = 0; n_i < nDomainNodes; n_i++ ) {
+ cVert = Mesh_GetVertex( cMesh, n_i );
+ nearest = Mesh_NearestVertex( fMesh, cVert );
+ fVert = Mesh_GetVertex( fMesh, nearest );
+ for( d_i = 0; d_i < nDims; d_i++ ) {
+ if( !Num_Approx( cVert[d_i], fVert[d_i] ) )
+ break;
+ }
+
+ Journal_Firewall( d_i == nDims,
+ errorStream,
+ "\n" \
+ "*****************************************************************\n" \
+ "* Error: The finest grid could not be sub-sampled to all coarse *\n" \
+ "* levels. This is due to the size of the finest grid. *\n" \
+ "*****************************************************************\n" \
+ "\n" );
+
+ if( level < nLevels - 2 )
+ self->topMaps[level][n_i] = self->topMaps[level + 1][nearest];
+ else
+ self->topMaps[level][n_i] = nearest;
+ }
+}
+
+void SROpGenerator_GenLevelEqNums( SROpGenerator* self, unsigned level ) {
+ Mesh* cMesh;
+ unsigned* nNodalDofs;
+ unsigned nDomainNodes, nLocalNodes;
+ DofLayout* dofLayout;
+ unsigned** dstArray;
+ unsigned curEqNum;
+ unsigned maxDofs;
+ unsigned topNode;
+ unsigned base, subTotal;
+ MPI_Comm comm;
+ unsigned nProcs, rank;
+ MPI_Status status;
+ unsigned* tuples;
+ Decomp_Sync* sync;
+ Decomp_Sync_Array* array;
+ unsigned n_i, dof_i;
+
+ assert( self && Stg_CheckType( self, SROpGenerator ) );
+ assert( self->meshes && self->topMaps && self->eqNums );
+ assert( level < MultigridSolver_GetNumLevels( self->solver ) - 1 );
+
+ cMesh = self->meshes[level];
+ nDomainNodes = Mesh_GetDomainSize( cMesh, MT_VERTEX );
+ nLocalNodes = Mesh_GetLocalSize( cMesh, MT_VERTEX );
+ dofLayout = self->fineEqNum->dofLayout;
+ comm = CommTopology_GetComm( Mesh_GetCommTopology( cMesh, MT_VERTEX ) );
+ MPI_Comm_size( comm, (int*)&nProcs );
+ MPI_Comm_rank( comm, (int*)&rank );
+
+ /* Allocate for destination array. */
+ nNodalDofs = AllocArray( unsigned, nDomainNodes );
+ for( n_i = 0; n_i < nDomainNodes; n_i++ )
+ nNodalDofs[n_i] = dofLayout->dofCounts[self->topMaps[level][n_i]];
+ dstArray = AllocComplex2D( unsigned, nDomainNodes, nNodalDofs );
+
+ /* Build initial destination array and store max dofs. */
+ curEqNum = 0;
+ maxDofs = 0;
for( n_i = 0; n_i < nLocalNodes; n_i++ ) {
- insist( FeMesh_NodeDomainToGlobal( feMesh, n_i, &globalNode ) );
- Grid_Lift( nodeGrid, globalNode, globalInds );
+ if( nNodalDofs[n_i] > maxDofs )
+ maxDofs = nNodalDofs[n_i];
- for( l_i = nLevels - 1; l_i < nLevels; l_i-- )
- SROpGenerator_ProcessLevel( self, l_i, globalNode, globalInds );
+ topNode = self->topMaps[level][n_i];
+ for( dof_i = 0; dof_i < nNodalDofs[n_i]; dof_i++ ) {
+ if( self->fineEqNum->destinationArray[topNode][dof_i] != (unsigned)-1 )
+ dstArray[n_i][dof_i] = curEqNum++;
+ else
+ dstArray[n_i][dof_i] = (unsigned)-1;
+ }
}
- FreeArray( globalInds );
+ /* Order the equation numbers based on processor rank; cascade counts forward. */
+ base = 0;
+ subTotal = curEqNum;
+ if( rank > 0 ) {
+ insist( !MPI_Recv( &base, 1, MPI_UNSIGNED, rank - 1, 6669, comm, &status ) );
+ subTotal = base + curEqNum;
+ }
+ if( rank < nProcs - 1 )
+ insist( !MPI_Send( &subTotal, 1, MPI_UNSIGNED, rank + 1, 6669, comm ) );
+
+ /* Modify existing destination array and dump to a tuple array. */
+ tuples = AllocArray( unsigned, nDomainNodes * maxDofs );
+ for( n_i = 0; n_i < nLocalNodes; n_i++ ) {
+ for( dof_i = 0; dof_i < nNodalDofs[n_i]; dof_i++ ) {
+ if( dstArray[n_i][dof_i] != (unsigned)-1 )
+ dstArray[n_i][dof_i] += base;
+ }
+ }
+
+ /* Update all other procs. */
+ sync = Mesh_GetSync( cMesh, MT_VERTEX );
+ array = Decomp_Sync_Array_New();
+ Decomp_Sync_Array_SetSync( array, sync );
+ Decomp_Sync_Array_SetMemory( array, tuples, tuples + nLocalNodes * maxDofs,
+ maxDofs * sizeof(unsigned), maxDofs * sizeof(unsigned),
+ maxDofs * sizeof(unsigned) );
+ Decomp_Sync_Array_Sync( array );
+ FreeObject( array );
+
+ /* Update destination array's domain indices. */
+ for( n_i = nLocalNodes; n_i < nDomainNodes; n_i++ ) {
+ topNode = self->topMaps[level][n_i];
+ for( dof_i = 0; dof_i < nNodalDofs[n_i]; dof_i++ ) {
+ if( self->fineEqNum->destinationArray[topNode][dof_i] != (unsigned)-1 )
+ dstArray[n_i][dof_i] = tuples[n_i * maxDofs + dof_i];
+ else
+ dstArray[n_i][dof_i] = -1;
+ }
+ }
+
+ /* Destroy arrays. */
+ FreeArray( nNodalDofs );
+ FreeArray( tuples );
+
+ self->eqNums[level] = dstArray;
+ self->nLocalEqNums[level] = curEqNum;
+ self->eqNumBases[level] = base;
}
-#endif
+
+void SROpGenerator_GenOps( SROpGenerator* self ) {
+ unsigned nLevels;
+ Matrix *fineMat, *P;
+ unsigned nRows, nCols;
+ unsigned l_i;
+
+ assert( self && Stg_CheckType( self, SROpGenerator ) );
+
+ fineMat = MatrixSolver_GetMatrix( self->solver );
+ nLevels = MultigridSolver_GetNumLevels( self->solver );
+
+ for( l_i = nLevels - 1; l_i > 0; l_i-- ) {
+ nRows = self->eqNums[l_i] ? self->nLocalEqNums[l_i] :
+ self->fineEqNum->localEqNumsOwnedCount;
+ nCols = self->nLocalEqNums[l_i - 1];
+
+ Matrix_Duplicate( fineMat, (void**)&P );
+ Matrix_SetComm( P, fineMat->comm );
+ Matrix_SetLocalSize( P, nRows, nCols );
+ if( !strcmp( P->type, PETScMatrix_Type ) ) {
+ unsigned *nDiagNonZeros, *nOffDiagNonZeros;
+
+ SROpGenerator_CalcOpNonZeros( self, l_i, &nDiagNonZeros, &nOffDiagNonZeros );
+ PETScMatrix_SetNonZeroStructure( P, 0, nDiagNonZeros, nOffDiagNonZeros );
+ FreeArray( nDiagNonZeros );
+ FreeArray( nOffDiagNonZeros );
+ }
+
+ SROpGenerator_GenLevelOp( self, l_i, P );
+ MultigridSolver_SetRestriction( self->solver, l_i, P );
+ MultigridSolver_SetProlongation( self->solver, l_i, P );
+ }
+}
+
+void SROpGenerator_GenLevelOp( SROpGenerator* self, unsigned level, Matrix* P ) {
+ Mesh *fMesh, *cMesh;
+ unsigned nDims;
+ unsigned nLocalFineNodes;
+ DofLayout* dofLayout;
+ unsigned ind;
+ unsigned nInc, *inc;
+ unsigned maxInc;
+ unsigned fTopNode, cTopNode;
+ unsigned fEqNum, cEqNum;
+ double *localCoord, *basis;
+ unsigned n_i, dof_i, dof_j, inc_i, e_i;
+
+ assert( self && Stg_CheckType( self, SROpGenerator ) );
+ assert( self->meshes );
+ assert( level < MultigridSolver_GetNumLevels( self->solver ) );
+ assert( P );
+
+ fMesh = self->meshes[level];
+ cMesh = self->meshes[level - 1];
+ nDims = Mesh_GetDimSize( fMesh );
+ nLocalFineNodes = Mesh_GetLocalSize( fMesh, MT_VERTEX );
+ dofLayout = self->fineEqNum->dofLayout;
+ localCoord = AllocArray( double, nDims );
+
+ maxInc = 0;
+ for( e_i = 0; e_i < Mesh_GetDomainSize( cMesh, nDims ); e_i++ ) {
+ nInc = Mesh_GetIncidenceSize( cMesh, nDims, e_i, MT_VERTEX );
+ if( nInc > maxInc )
+ maxInc = nInc;
+ }
+ basis = AllocArray( double, maxInc );
+
+ for( n_i = 0; n_i < nLocalFineNodes; n_i++ ) {
+ if( self->topMaps[level] )
+ fTopNode = self->topMaps[level][n_i];
+ else
+ fTopNode = n_i;
+
+ for( dof_i = 0; dof_i < dofLayout->dofCounts[fTopNode]; dof_i++ ) {
+ if( self->eqNums[level] )
+ fEqNum = self->eqNums[level][n_i][dof_i];
+ else
+ fEqNum = self->fineEqNum->destinationArray[fTopNode][dof_i];
+
+ if( fEqNum == (unsigned)-1 )
+ continue;
+
+ insist( Mesh_SearchElements( cMesh, Mesh_GetVertex( fMesh, n_i ), &ind ) );
+ FeMesh_CoordGlobalToLocal( cMesh, ind, Mesh_GetVertex( fMesh, n_i ), localCoord );
+ FeMesh_EvalBasis( cMesh, ind, localCoord, basis );
+ Mesh_GetIncidence( cMesh, nDims, ind, MT_VERTEX, &nInc, &inc );
+ for( inc_i = 0; inc_i < nInc; inc_i++ ) {
+ cTopNode = self->topMaps[level - 1][inc[inc_i]];
+ for( dof_j = 0; dof_j < dofLayout->dofCounts[cTopNode]; dof_j++ ) {
+ cEqNum = self->eqNums[level - 1][inc[inc_i]][dof_j];
+ if( cEqNum != (unsigned)-1 && !Num_Approx( basis[inc_i], 0.0 ) )
+ Matrix_InsertEntries( P, 1, &fEqNum, 1, &cEqNum, basis + inc_i );
+ }
+ }
+ }
+ }
+
+ FreeArray( localCoord );
+ FreeArray( basis );
+
+ Matrix_AssemblyBegin( P );
+ Matrix_AssemblyEnd( P );
+}
+
+void SROpGenerator_CalcOpNonZeros( SROpGenerator* self, unsigned level,
+ unsigned** nDiagNonZeros, unsigned** nOffDiagNonZeros )
+{
+ Mesh *fMesh, *cMesh;
+ unsigned nLocalFineNodes;
+ DofLayout* dofLayout;
+ unsigned dim, ind;
+ unsigned nInc, *inc;
+ unsigned fTopNode, cTopNode;
+ unsigned fEqNum, cEqNum;
+ unsigned nLocalEqNums;
+ unsigned n_i, dof_i, dof_j, inc_i;
+
+ assert( self && Stg_CheckType( self, SROpGenerator ) );
+ assert( self->meshes );
+ assert( level < MultigridSolver_GetNumLevels( self->solver ) );
+
+ fMesh = self->meshes[level];
+ cMesh = self->meshes[level - 1];
+ nLocalFineNodes = Mesh_GetLocalSize( fMesh, MT_VERTEX );
+ dofLayout = self->fineEqNum->dofLayout;
+ if( self->eqNums[level] )
+ nLocalEqNums = self->nLocalEqNums[level];
+ else
+ nLocalEqNums = self->fineEqNum->localEqNumsOwnedCount;
+
+ *nDiagNonZeros = AllocArray( unsigned, nLocalEqNums );
+ memset( *nDiagNonZeros, 0, nLocalEqNums * sizeof(unsigned) );
+ *nOffDiagNonZeros = AllocArray( unsigned, nLocalEqNums );
+ memset( *nOffDiagNonZeros, 0, nLocalEqNums * sizeof(unsigned) );
+
+ for( n_i = 0; n_i < nLocalFineNodes; n_i++ ) {
+ if( self->topMaps[level] )
+ fTopNode = self->topMaps[level][n_i];
+ else
+ fTopNode = n_i;
+
+ for( dof_i = 0; dof_i < dofLayout->dofCounts[fTopNode]; dof_i++ ) {
+ if( self->eqNums[level] )
+ fEqNum = self->eqNums[level][n_i][dof_i];
+ else
+ fEqNum = self->fineEqNum->destinationArray[fTopNode][dof_i];
+
+ if( fEqNum == (unsigned)-1 )
+ continue;
+
+ insist( Mesh_Search( cMesh, Mesh_GetVertex( fMesh, n_i ), &dim, &ind ) );
+ if( dim == MT_VERTEX ) {
+ cTopNode = self->topMaps[level - 1][ind];
+ for( dof_j = 0; dof_j < dofLayout->dofCounts[cTopNode]; dof_j++ ) {
+ cEqNum = self->eqNums[level - 1][ind][dof_j];
+ if( cEqNum == (unsigned)-1 )
+ continue;
+
+ if( cEqNum - self->eqNumBases[level - 1] < nLocalEqNums )
+ (*nDiagNonZeros)[fEqNum - self->eqNumBases[level]]++;
+ else
+ (*nOffDiagNonZeros)[fEqNum - self->eqNumBases[level]]++;
+ }
+ }
+ else {
+ Mesh_GetIncidence( cMesh, dim, ind, MT_VERTEX, &nInc, &inc );
+ for( inc_i = 0; inc_i < nInc; inc_i++ ) {
+ cTopNode = self->topMaps[level - 1][inc[inc_i]];
+ for( dof_j = 0; dof_j < dofLayout->dofCounts[cTopNode]; dof_j++ ) {
+ cEqNum = self->eqNums[level - 1][inc[inc_i]][dof_j];
+ if( cEqNum == (unsigned)-1 )
+ continue;
+
+ if( cEqNum - self->eqNumBases[level - 1] < nLocalEqNums )
+ (*nDiagNonZeros)[fEqNum - self->eqNumBases[level]]++;
+ else
+ (*nOffDiagNonZeros)[fEqNum - self->eqNumBases[level]]++;
+ }
+ }
+ }
+ }
+ }
+}
+
+void SROpGenerator_DestructMeshes( SROpGenerator* self ) {
+ unsigned nLevels;
+ unsigned l_i;
+
+ assert( self && Stg_CheckType( self, SROpGenerator ) );
+
+ if( self->meshes ) {
+ nLevels = MultigridSolver_GetNumLevels( self->solver );
+ for( l_i = 0; l_i < nLevels - 1; l_i++ ) {
+ FreeObject( self->meshes[l_i] );
+ FreeArray( self->topMaps[l_i] );
+ FreeArray( self->eqNums[l_i] );
+ }
+ KillArray( self->meshes );
+ KillArray( self->topMaps );
+ KillArray( self->eqNums );
+ KillArray( self->nLocalEqNums );
+ KillArray( self->eqNumBases );
+ }
+}
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/SROpGenerator.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/SROpGenerator.h 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/SROpGenerator.h 2007-01-09 07:50:52 UTC (rev 5709)
@@ -47,13 +47,20 @@
/** Virtual function types */
/** SROpGenerator class contents */
- #define __SROpGenerator \
- /* General info */ \
- __MGOpGenerator \
- \
- /* Virtual info */ \
- \
- /* SROpGenerator info */
+ #define __SROpGenerator \
+ /* General info */ \
+ __MGOpGenerator \
+ \
+ /* Virtual info */ \
+ \
+ /* SROpGenerator info */ \
+ FeVariable* fineVar; \
+ FeEquationNumber* fineEqNum; \
+ Mesh** meshes; \
+ unsigned** topMaps; \
+ unsigned*** eqNums; \
+ unsigned* nLocalEqNums; \
+ unsigned* eqNumBases;
struct SROpGenerator { __SROpGenerator };
@@ -90,8 +97,22 @@
** Public functions
*/
+ void SROpGenerator_SetFineVariable( void* srOpGenerator, void* variable );
+
/*--------------------------------------------------------------------------------------------------------------------------
** Private Member functions
*/
+ void SROpGenerator_GenMeshes( SROpGenerator* self );
+ void SROpGenerator_GenLevelMesh( SROpGenerator* self, unsigned level );
+ void SROpGenerator_GenLevelTopMap( SROpGenerator* self, unsigned level );
+ void SROpGenerator_GenLevelEqNums( SROpGenerator* self, unsigned level );
+
+ void SROpGenerator_GenOps( SROpGenerator* self );
+ void SROpGenerator_GenLevelOp( SROpGenerator* self, unsigned level, Matrix* P );
+ void SROpGenerator_CalcOpNonZeros( SROpGenerator* self, unsigned level,
+ unsigned** nDiagNonZeros, unsigned** nOffDiagNonZeros );
+
+ void SROpGenerator_DestructMeshes( SROpGenerator* self );
+
#endif /* __StgFEM_SLE_LinearAlgebra_SROpGenerator_h__ */
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Vector.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Vector.c 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Vector.c 2007-01-09 07:50:52 UTC (rev 5709)
@@ -35,6 +35,7 @@
#include <mpi.h>
#include "StGermain/StGermain.h"
+#include "Discretisation/Discretisation.h"
#include "LinearAlgebra.h"
@@ -159,8 +160,10 @@
unsigned entry_i;
assert( self && Stg_CheckType( self, Vector ) );
- assert( stream );
+ if( !stream )
+ stream = Journal_Register( Info_Type, "tmp" );
+
size = Vector_GetLocalSize( self );
Vector_GetArray( self, &array );
Journal_Printf( stream, "%s = [\n", self->name );
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/types.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/types.h 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/types.h 2007-01-09 07:50:52 UTC (rev 5709)
@@ -85,7 +85,9 @@
PETScMatrixSolver_PCType_Jacobi,
PETScMatrixSolver_PCType_BlockJacobi,
PETScMatrixSolver_PCType_SOR,
+ PETScMatrixSolver_PCType_ParallelSOR,
PETScMatrixSolver_PCType_LU,
+ PETScMatrixSolver_PCType_RedundantLU,
PETScMatrixSolver_PCType_ILU,
PETScMatrixSolver_PCType_None
} PETScMatrixSolver_PCType;
Added: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/tests/data/wallVC.xml
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/tests/data/wallVC.xml 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/tests/data/wallVC.xml 2007-01-09 07:50:52 UTC (rev 5709)
@@ -0,0 +1,24 @@
+<?xml version="1.0"?>
+<!DOCTYPE StGermainData SYSTEM "stgermain.dtd">
+<!-- A StGermain input file -->
+<!-- DTD to validate against -->
+<StGermainData xmlns="http://www.vpac.org/StGermain/XML_IO_Handler/Jun2003">
+
+ <struct name="wallVC">
+ <param name="type"> WallVC </param>
+ <param name="wall"> bottom </param>
+ <list name="variables">
+ <struct>
+ <param name="name"> one </param>
+ <param name="type"> double </param>
+ <param name="value" type="double"> 1 </param>
+ </struct>
+ <struct>
+ <param name="name"> two </param>
+ <param name="type"> double </param>
+ <param name="value" type="double"> 3 </param>
+ </struct>
+ </list>
+ </struct>
+
+</StGermainData>
Added: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/tests/makefile
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/tests/makefile 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/tests/makefile 2007-01-09 07:50:52 UTC (rev 5709)
@@ -0,0 +1,63 @@
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+##
+## Copyright (C), 2003-2006, Victorian Partnership for Advanced Computing (VPAC) Ltd, 110 Victoria Street,
+## Melbourne, 3053, Australia.
+##
+## Primary Contributing Organisations:
+## Victorian Partnership for Advanced Computing Ltd, Computational Software Development - http://csd.vpac.org
+## Australian Computational Earth Systems Simulator - http://www.access.edu.au
+## Monash Cluster Computing - http://www.mcc.monash.edu.au
+## Computational Infrastructure for Geodynamics - http://www.geodynamics.org
+##
+## Contributors:
+## Patrick D. Sunter, Software Engineer, VPAC. (pds at vpac.org)
+## Robert Turnbull, Research Assistant, Monash University. (robert.turnbull at sci.monash.edu.au)
+## Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+## David May, PhD Student, Monash University (david.may at sci.monash.edu.au)
+## Louis Moresi, Associate Professor, Monash University. (louis.moresi at sci.monash.edu.au)
+## Luke J. Hodkinson, Computational Engineer, VPAC. (lhodkins at vpac.org)
+## Alan H. Lo, Computational Engineer, VPAC. (alan at vpac.org)
+## Raquibul Hassan, Computational Engineer, VPAC. (raq at vpac.org)
+## Julian Giordani, Research Assistant, Monash University. (julian.giordani at sci.monash.edu.au)
+## Vincent Lemiale, Postdoctoral Fellow, Monash University. (vincent.lemiale at sci.monash.edu.au)
+##
+## This library is free software; you can redistribute it and/or
+## modify it under the terms of the GNU Lesser General Public
+## License as published by the Free Software Foundation; either
+## version 2.1 of the License, or (at your option) any later version.
+##
+## This library is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+## Lesser General Public License for more details.
+##
+## You should have received a copy of the GNU Lesser General Public
+## License along with this library; if not, write to the Free Software
+## Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+##
+## $Id: makefile 676 2006-12-22 01:44:03Z LukeHodkinson $
+##
+##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+#Finds the Absolute path to the Project Root directory
+SHELL := /bin/bash
+PROJ_ROOT := $(shell until test -r ./Makefile.system ; do cd .. ; done ; echo `pwd`)
+include ${PROJ_ROOT}/Makefile.system
+
+NAME = $(shell basename `pwd | sed s/tests//g`)
+
+tests = lib${PROJECT}_SLE_LinearAlgebra
+checks = $(wildcard *.sh)
+SRCS = $(wildcard *.c)
+
+EXTERNAL_INCLUDES = -I${INC_DIR}/${PROJECT}
+EXTERNAL_LIBS= -L${LIB_DIR} -l${PROJECT}_SLE_LinearAlgebra -l${PROJECT}_Discretisation
+
+packages = STGERMAIN MPI XML MATH DL HYPRE
+
+ifeq (true,$(shell if test -r $(PETSC_INCDIR)/petsc.h; then echo true; fi ))
+ USER_CFLAGS += -DHAVE_PETSC
+ packages += PETSC
+endif
+
+include ${PROJ_ROOT}/Makefile.vmake
Added: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/tests/testSROpGenerator.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/tests/testSROpGenerator.c 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/tests/testSROpGenerator.c 2007-01-09 07:50:52 UTC (rev 5709)
@@ -0,0 +1,190 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003, Victorian Partnership for Advanced Computing (VPAC) Ltd, 110 Victoria Street, Melbourne, 3053, Australia.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Patrick D. Sunter, Software Engineer, VPAC. (pds at vpac.org)
+** Luke J. Hodkinson, Computational Engineer, VPAC. (lhodkins at vpac.org)
+** Siew-Ching Tan, Software Engineer, VPAC. (siew at vpac.org)
+** Alan H. Lo, Computational Engineer, VPAC. (alan at vpac.org)
+** Raquibul Hassan, Computational Engineer, VPAC. (raq at vpac.org)
+**
+** This library is free software; you can redistribute it and/or
+** modify it under the terms of the GNU Lesser General Public
+** License as published by the Free Software Foundation; either
+** version 2.1 of the License, or (at your option) any later version.
+**
+** This library is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+** Lesser General Public License for more details.
+**
+** You should have received a copy of the GNU Lesser General Public
+** License along with this library; if not, write to the Free Software
+** Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+**
+** $Id: testSROpGenerator.c 678 2006-12-22 01:45:53Z LukeHodkinson $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <assert.h>
+#include <mpi.h>
+
+#include "StGermain/StGermain.h"
+#include "Discretisation/Discretisation.h"
+#include "SLE/LinearAlgebra/LinearAlgebra.h"
+
+
+FeEquationNumber* buildEqNumWithBCs( unsigned nProcs ) {
+ CartesianGenerator* gen;
+ FeMesh* feMesh;
+ DofLayout* dofs;
+ FeEquationNumber* eqNum;
+ VariableCondition* bcs;
+ Variable_Register* varReg;
+ ConditionFunction_Register* cfReg;
+ Variable* vars[2];
+ unsigned maxDecomp[3] = {0, 1, 1};
+ unsigned sizes[3];
+ double minCrd[3];
+ double maxCrd[3];
+ SizeT dataOffs = 1;
+ Variable_DataType dataType = Variable_DataType_Double;
+ unsigned nDataTypes = 1;
+ char* dataNames = "nothing";
+ static SizeT structSize = sizeof(double);
+ static unsigned arraySize;
+ static void* arrayPtrs[2];
+ Dictionary* dict;
+ XML_IO_Handler* ioHandler;
+
+ sizes[1] = sizes[2] = sizes[0] = nProcs * 8;
+ minCrd[0] = minCrd[1] = minCrd[2] = 0.0;
+ maxCrd[0] = maxCrd[1] = maxCrd[2] = (double)nProcs;
+
+ gen = CartesianGenerator_New( "" );
+ CartesianGenerator_SetDimSize( gen, 3 );
+ CartesianGenerator_SetTopologyParams( gen, sizes, 0, NULL, maxDecomp );
+ CartesianGenerator_SetGeometryParams( gen, minCrd, maxCrd );
+ CartesianGenerator_SetShadowDepth( gen, 0 );
+
+ feMesh = FeMesh_New( "" );
+ Mesh_SetGenerator( feMesh, gen );
+ FeMesh_SetElementFamily( feMesh, "linear" );
+ Build( feMesh, NULL, False );
+
+ varReg = Variable_Register_New();
+ cfReg = ConditionFunction_Register_New();
+
+ arraySize = Mesh_GetDomainSize( feMesh, MT_VERTEX );
+ arrayPtrs[0] = Memory_Alloc_Array_Unnamed( double, arraySize );
+ arrayPtrs[1] = Memory_Alloc_Array_Unnamed( double, arraySize );
+ vars[0] = Variable_New( "one", 1, &dataOffs, &dataType, &nDataTypes, &dataNames,
+ &structSize, &arraySize, arrayPtrs, varReg );
+ vars[1] = Variable_New( "two", 1, &dataOffs, &dataType, &nDataTypes, &dataNames,
+ &structSize, &arraySize, arrayPtrs + 1, varReg );
+
+ dofs = DofLayout_New( "", varReg, 0, feMesh );
+ dofs->nBaseVariables = 2;
+ dofs->baseVariables = Memory_Alloc_Array_Unnamed( Variable*, 2 );
+ dofs->baseVariables[0] = vars[0];
+ dofs->baseVariables[1] = vars[1];
+ Build( dofs, NULL, False );
+ Initialise( dofs, NULL, False );
+
+ ioHandler = XML_IO_Handler_New();
+ dict = Dictionary_New();
+ IO_Handler_ReadAllFromFile( ioHandler, "data/wallVC.xml", dict );
+ bcs = (VariableCondition*)WallVC_New( "", "wallVC", varReg, cfReg, dict, feMesh );
+ Build( bcs, NULL, False );
+ Initialise( bcs, NULL, False );
+
+ eqNum = FeEquationNumber_New( "", feMesh, dofs, bcs, NULL );
+ Build( eqNum, NULL, False );
+ Initialise( eqNum, NULL, False );
+
+ return eqNum;
+}
+
+
+Bool testOps( unsigned rank, unsigned nProcs, unsigned watch ) {
+ Bool result = True;
+ FeEquationNumber* eqNum;
+ MultigridSolver* solver;
+ SROpGenerator* opGen;
+ Matrix* mat;
+
+ mat = (Matrix*)PETScMatrix_New( "" );
+ Matrix_SetLocalSize( mat, nProcs, nProcs );
+ Matrix_AssemblyBegin( mat );
+ Matrix_AssemblyEnd( mat );
+
+ eqNum = buildEqNumWithBCs( nProcs );
+ solver = MultigridSolver_New( "" );
+ MultigridSolver_SetLevels( solver, 4 );
+
+ MatrixSolver_SetMatrix( solver, mat );
+ opGen = SROpGenerator_New( "" );
+ MGOpGenerator_SetMultigridSolver( opGen, solver );
+ /*SROpGenerator_SetVariable( opGen, eqNum ); TODO */
+ Build( opGen, NULL, False );
+ SROpGenerator_Generate( opGen );
+
+ if( rank == watch ) {
+ }
+
+done:
+ FreeObject( opGen );
+ FreeObject( mat );
+ FreeObject( solver );
+ FreeObject( eqNum );
+
+ return result;
+}
+
+
+#define nTests 1
+
+TestSuite_Test tests[nTests] = {{"test operator generation", testOps, 1}};
+
+
+int main( int argc, char* argv[] ) {
+ TestSuite* suite;
+
+ /* Initialise MPI, get world info. */
+ MPI_Init( &argc, &argv );
+
+ /* Initialise StGermain. */
+ StGermain_Init( &argc, &argv );
+ StgFEM_Discretisation_Init( &argc, &argv );
+ StgFEM_SLE_LinearAlgebra_Init( &argc, &argv );
+
+ /* Disable all journaling. */
+ Journal_Enable_TypedStream( Debug_Type, False );
+ Journal_Enable_TypedStream( Info_Type, False );
+
+ /* Create the test suite. */
+ suite = TestSuite_New();
+ TestSuite_SetProcToWatch( suite, (argc >= 2) ? atoi( argv[1] ) : 0 );
+ TestSuite_SetTests( suite, nTests, tests );
+
+ /* Run the tests. */
+ TestSuite_Run( suite );
+
+ /* Destroy test suites. */
+ FreeObject( suite );
+
+ /* Finalise StGermain. */
+ BaseContainer_Finalise();
+ BaseIO_Finalise();
+ BaseFoundation_Finalise();
+
+ /* Close off MPI */
+ MPI_Finalize();
+
+ return MPI_SUCCESS;
+}
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/AdvectionDiffusion/src/Multicorrector.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/AdvectionDiffusion/src/Multicorrector.c 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/AdvectionDiffusion/src/Multicorrector.c 2007-01-09 07:50:52 UTC (rev 5709)
@@ -232,6 +232,7 @@
/* Allocate Memory For Corrector Step */
Vector_Duplicate( sle->phiVector->vector, (void**)&deltaPhiDot );
+ Vector_SetLocalSize( deltaPhiDot, Vector_GetLocalSize( sle->phiVector->vector ) );
/* Multi-corrector Steps */
for ( iteration_I = 0 ; iteration_I < self->multiCorrectorIterations ; iteration_I++ ) {
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/Energy/src/Energy_SLE_Solver.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/Energy/src/Energy_SLE_Solver.c 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/Energy/src/Energy_SLE_Solver.c 2007-01-09 07:50:52 UTC (rev 5709)
@@ -129,6 +129,7 @@
void _Energy_SLE_Solver_Init( Energy_SLE_Solver* self ) {
self->isConstructed = True;
self->residual = NULL;
+ self->matrixSolver = NULL;
}
void Energy_SLE_Solver_InitAll( Energy_SLE_Solver* solver, Bool useStatSolve, int statReps ) {
Energy_SLE_Solver* self = (Energy_SLE_Solver*)solver;
@@ -163,7 +164,7 @@
void* _Energy_SLE_Solver_Copy( void* standardSleSolver, void* dest, Bool deep, Name nameExt, PtrMap* ptrMap ) {
Energy_SLE_Solver* self = (Energy_SLE_Solver*)standardSleSolver;
Energy_SLE_Solver* newEnergySleSolver;
-
+
newEnergySleSolver = _SLE_Solver_Copy( self, dest, deep, nameExt, ptrMap );
newEnergySleSolver->matrixSolver = self->matrixSolver;
@@ -172,11 +173,22 @@
}
void _Energy_SLE_Solver_Construct( void* sleSolver, Stg_ComponentFactory* cf, void* data ) {
- Energy_SLE_Solver *self = (Energy_SLE_Solver*)sleSolver;
+ Energy_SLE_Solver *self = (Energy_SLE_Solver*)sleSolver;
+ assert( self && Stg_CheckType( self, Energy_SLE_Solver ) );
+ assert( cf && Stg_CheckType( cf, Stg_ComponentFactory ) );
+
_SLE_Solver_Construct( self, cf, data );
- _Energy_SLE_Solver_Init( self );
+ self->matrixSolver = Stg_ComponentFactory_ConstructByKey( cf, self->name, "matrixSolver", MatrixSolver,
+ False, data );
+ if( !self->matrixSolver ) {
+#ifdef HAVE_PETSC
+ self->matrixSolver = (MatrixSolver*)PETScMatrixSolver_New( "" );
+#else
+ self->matrixSolver = NULL;
+#endif
+ }
}
/* Build */
@@ -191,16 +203,12 @@
Build( stiffMat, standardSLE, False );
-#ifdef HAVE_PETSC
- self->matrixSolver = (MatrixSolver*)PETScMatrixSolver_New( "" );
-#else
-#error No linear algebra package found.
-#endif
- MatrixSolver_SetComm( self->matrixSolver, sle->comm );
+ if( self->matrixSolver ) {
+ MatrixSolver_SetComm( self->matrixSolver, sle->comm );
+ Build( self->matrixSolver, NULL, False );
+ }
Stream_UnIndentBranch( StgFEM_SLE_ProvidedSystems_Energy_Debug );
-
- Vector_Duplicate( ((ForceVector**)sle->forceVectors->data)[0]->vector, (void**)&self->residual );
}
@@ -210,14 +218,15 @@
/* Initialise parent. */
_SLE_Solver_Initialise( self, sle );
+
+ if( self->matrixSolver )
+ Initialise( self->matrixSolver, NULL, False );
}
void _Energy_SLE_Solver_Execute( void* sleSolver, void* data ) {
-
}
void _Energy_SLE_Solver_Destroy( void* sleSolver, void* data ) {
-
}
void _Energy_SLE_Solver_SolverSetup( void* sleSolver, void* standardSLE ) {
@@ -235,6 +244,10 @@
if( self->maxIterations > 0 ) {
MatrixSolver_SetMaxIterations( self->matrixSolver, self->maxIterations );
}
+
+ Vector_Duplicate( ((ForceVector**)sle->forceVectors->data)[0]->vector, (void**)&self->residual );
+ Vector_SetLocalSize( self->residual,
+ Vector_GetLocalSize( ((ForceVector**)sle->forceVectors->data)[0]->vector ) );
}
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_PenaltySolver.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_PenaltySolver.c 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_PenaltySolver.c 2007-01-09 07:50:52 UTC (rev 5709)
@@ -242,8 +242,11 @@
Journal_DPrintf( self->debug, "In %s():\n", __func__ );
Vector_Duplicate( hVec, (void**)&hTempVec );
+ Vector_SetLocalSize( hTempVec, Vector_GetLocalSize( hVec ) );
Vector_Duplicate( fVec, (void**)&fTempVec );
+ Vector_SetLocalSize( fTempVec, Vector_GetLocalSize( fVec ) );
Vector_Duplicate( pVec, (void**)&diagC );
+ Vector_SetLocalSize( diagC, Vector_GetLocalSize( pVec ) );
if( sle->dStiffMat == NULL ) {
Journal_DPrintf( self->debug, "Div matrix == NULL : Problem is assumed to be symmetric. ie Div = GTrans \n");
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_UzawaSolver.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_UzawaSolver.c 2007-01-09 07:50:48 UTC (rev 5708)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_UzawaSolver.c 2007-01-09 07:50:52 UTC (rev 5709)
@@ -247,11 +247,16 @@
Journal_DPrintfL( self->debug, 2, "Allocate the auxillary vectors pTemp, r, s, fTemp and vStar.\n" );
Vector_Duplicate( sle->pSolnVec->vector, (void**)&self->pTempVec );
+ Vector_SetLocalSize( self->pTempVec, Vector_GetLocalSize( sle->pSolnVec->vector ) );
Vector_Duplicate( sle->pSolnVec->vector, (void**)&self->rVec );
+ Vector_SetLocalSize( self->rVec, Vector_GetLocalSize( sle->pSolnVec->vector ) );
Vector_Duplicate( sle->pSolnVec->vector, (void**)&self->sVec );
+ Vector_SetLocalSize( self->sVec, Vector_GetLocalSize( sle->pSolnVec->vector ) );
Vector_Duplicate( sle->fForceVec->vector, (void**)&self->fTempVec );
+ Vector_SetLocalSize( self->fTempVec, Vector_GetLocalSize( sle->fForceVec->vector ) );
Vector_Duplicate( sle->fForceVec->vector, (void**)&self->vStarVec );
+ Vector_SetLocalSize( self->vStarVec, Vector_GetLocalSize( sle->fForceVec->vector ) );
Stream_UnIndentBranch( StgFEM_Debug );
}
More information about the cig-commits
mailing list