Commit 169b9d9a authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Added more features to FETI-DP, splitted into some more files.

parent ecba737f
...@@ -254,6 +254,8 @@ if(ENABLE_PARALLEL_DOMAIN) ...@@ -254,6 +254,8 @@ if(ENABLE_PARALLEL_DOMAIN)
${SOURCE_DIR}/parallel/PetscSolver.cc ${SOURCE_DIR}/parallel/PetscSolver.cc
${SOURCE_DIR}/parallel/PetscProblemStat.cc ${SOURCE_DIR}/parallel/PetscProblemStat.cc
${SOURCE_DIR}/parallel/PetscSolverFeti.cc ${SOURCE_DIR}/parallel/PetscSolverFeti.cc
${SOURCE_DIR}/parallel/PetscSolverFetiOperators.cc
${SOURCE_DIR}/parallel/PetscSolverFetiTimings.cc
${SOURCE_DIR}/parallel/PetscSolverGlobalMatrix.cc ${SOURCE_DIR}/parallel/PetscSolverGlobalMatrix.cc
${SOURCE_DIR}/parallel/PetscSolverGlobalBlockMatrix.cc ${SOURCE_DIR}/parallel/PetscSolverGlobalBlockMatrix.cc
${SOURCE_DIR}/parallel/PetscSolverSchur.cc) ${SOURCE_DIR}/parallel/PetscSolverSchur.cc)
......
...@@ -43,6 +43,7 @@ ...@@ -43,6 +43,7 @@
#include "FiniteElemSpace.h" #include "FiniteElemSpace.h"
#include "SurfaceQuadrature.h" #include "SurfaceQuadrature.h"
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
#include <mpi.h>
#include "parallel/ParallelDofMapping.h" #include "parallel/ParallelDofMapping.h"
#endif #endif
......
...@@ -67,10 +67,8 @@ namespace AMDiS { ...@@ -67,10 +67,8 @@ namespace AMDiS {
errorEstimate(0.0) errorEstimate(0.0)
{} {}
/** \brief /// Refinement of parent to child1 and child2.
* Refinement of parent to child1 and child2. /// @return true: must this ElementData, else not allowed to delete it
* @return true: must this ElementData, else not allowed to delete it
*/
bool refineElementData(Element* parent, bool refineElementData(Element* parent,
Element* child1, Element* child1,
Element* child2, Element* child2,
......
...@@ -12,6 +12,7 @@ ...@@ -12,6 +12,7 @@
#include "parallel/DofComm.h" #include "parallel/DofComm.h"
#include "parallel/InteriorBoundary.h" #include "parallel/InteriorBoundary.h"
#include "parallel/MeshLevelData.h"
#include "FiniteElemSpace.h" #include "FiniteElemSpace.h"
namespace AMDiS { namespace AMDiS {
......
...@@ -23,6 +23,7 @@ ...@@ -23,6 +23,7 @@
#ifndef AMDIS_DOF_COMM_H #ifndef AMDIS_DOF_COMM_H
#define AMDIS_DOF_COMM_H #define AMDIS_DOF_COMM_H
#include <mpi.h>
#include <map> #include <map>
#include "parallel/ParallelTypes.h" #include "parallel/ParallelTypes.h"
#include "FiniteElemSpace.h" #include "FiniteElemSpace.h"
......
...@@ -12,6 +12,7 @@ ...@@ -12,6 +12,7 @@
#include "VertexVector.h" #include "VertexVector.h"
#include "parallel/ElementObjectDatabase.h" #include "parallel/ElementObjectDatabase.h"
#include "parallel/MeshLevelData.h"
namespace AMDiS { namespace AMDiS {
......
...@@ -28,7 +28,6 @@ ...@@ -28,7 +28,6 @@
#include "AMDiS_fwd.h" #include "AMDiS_fwd.h"
#include "BoundaryObject.h" #include "BoundaryObject.h"
#include "parallel/MeshLevelData.h"
#include "parallel/ParallelTypes.h" #include "parallel/ParallelTypes.h"
namespace AMDiS { namespace AMDiS {
......
...@@ -20,15 +20,16 @@ ...@@ -20,15 +20,16 @@
/** \file MeshLevelData.h */ /** \file MeshLevelData.h */
#ifndef AMDIS_MESH_LEVEL_DATA_H
#define AMDIS_MESH_LEVEL_DATA_H
#include <iostream> #include <iostream>
#include <set> #include <set>
#include <vector> #include <vector>
#include <mpi.h> #include <mpi.h>
#include "Global.h" #include "Global.h"
#ifndef AMDIS_MESH_LEVEL_DATA_H
#define AMDIS_MESH_LEVEL_DATA_H
namespace AMDiS { namespace AMDiS {
using namespace std; using namespace std;
......
...@@ -10,10 +10,9 @@ ...@@ -10,10 +10,9 @@
// See also license.opensource.txt in the distribution. // See also license.opensource.txt in the distribution.
#include <mpi.h>
#include "MpiHelper.h" #include "MpiHelper.h"
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
namespace AMDiS { namespace AMDiS {
namespace mpi { namespace mpi {
...@@ -54,7 +53,10 @@ namespace AMDiS { ...@@ -54,7 +53,10 @@ namespace AMDiS {
MPI::COMM_WORLD.Allreduce(&valCopy, &value, 1, MPI_INT, MPI_MAX); MPI::COMM_WORLD.Allreduce(&valCopy, &value, 1, MPI_INT, MPI_MAX);
} }
void startRand()
{
srand(time(NULL) * (MPI::COMM_WORLD.Get_rank() + 1));
}
} }
} }
#endif
...@@ -23,12 +23,10 @@ ...@@ -23,12 +23,10 @@
#ifndef AMDIS_MPIHELPER_H #ifndef AMDIS_MPIHELPER_H
#define AMDIS_MPIHELPER_H #define AMDIS_MPIHELPER_H
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS #include <mpi.h>
#include "Global.h"
#include <time.h> #include <time.h>
#include <stdlib.h> #include <stdlib.h>
#include <mpi.h> #include "Global.h"
namespace AMDiS { namespace AMDiS {
...@@ -63,10 +61,7 @@ namespace AMDiS { ...@@ -63,10 +61,7 @@ namespace AMDiS {
WARNING("Unknown type for globalMax. Can not determine maximal value of all processors!\n"); WARNING("Unknown type for globalMax. Can not determine maximal value of all processors!\n");
} }
inline void startRand() void startRand();
{
srand(time(NULL) * (MPI::COMM_WORLD.Get_rank() + 1));
}
/** \brief /** \brief
* In many situations a rank computes a number of local DOFs. Then all * In many situations a rank computes a number of local DOFs. Then all
...@@ -97,5 +92,3 @@ namespace AMDiS { ...@@ -97,5 +92,3 @@ namespace AMDiS {
} }
#endif #endif
#endif
...@@ -191,13 +191,8 @@ namespace AMDiS { ...@@ -191,13 +191,8 @@ namespace AMDiS {
nOverallCoarseRows, nOverallCoarseRows, nOverallCoarseRows, nOverallCoarseRows,
0, nnz[i + 1][i + 1].dnnz, 0, nnz[i + 1][i + 1].onnz, 0, nnz[i + 1][i + 1].dnnz, 0, nnz[i + 1][i + 1].onnz,
&mat[i + 1][i + 1]); &mat[i + 1][i + 1]);
MSG("REMOVE THIS!\n"); cMap->createVec(vecSol[i + 1]);
MatSetOption(mat[i + 1][i + 1], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); cMap->createVec(vecRhs[i + 1]);
VecCreateMPI(mpiCommGlobal, nRankCoarseRows, nOverallCoarseRows,
&vecSol[i + 1]);
VecCreateMPI(mpiCommGlobal, nRankCoarseRows, nOverallCoarseRows,
&vecRhs[i + 1]);
} }
for (int i = 0; i < nCoarseMap + 1; i++) { for (int i = 0; i < nCoarseMap + 1; i++) {
...@@ -215,18 +210,12 @@ namespace AMDiS { ...@@ -215,18 +210,12 @@ namespace AMDiS {
nRowsRankMat, nColsRankMat, nRowsRankMat, nColsRankMat,
nRowsOverallMat, nColsOverallMat, nRowsOverallMat, nColsOverallMat,
0, nnz[i][j].dnnz, 0, nnz[i][j].onnz, 0, nnz[i][j].dnnz, 0, nnz[i][j].onnz,
&mat[i][j]); &mat[i][j]);
MSG("REMOVE THIS!\n");
MatSetOption(mat[i][j], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
MatCreateAIJ(mpiCommGlobal, MatCreateAIJ(mpiCommGlobal,
nColsRankMat, nRowsRankMat, nColsRankMat, nRowsRankMat,
nColsOverallMat, nRowsOverallMat, nColsOverallMat, nRowsOverallMat,
0, nnz[j][i].dnnz, 0, nnz[j][i].onnz, 0, nnz[j][i].dnnz, 0, nnz[j][i].onnz,
&mat[j][i]); &mat[j][i]);
MSG("REMOVE THIS!\n");
MatSetOption(mat[j][i], MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
} }
} }
} }
......
...@@ -11,6 +11,7 @@ ...@@ -11,6 +11,7 @@
#include "parallel/ParallelDofMapping.h" #include "parallel/ParallelDofMapping.h"
#include "parallel/MeshLevelData.h"
#include "parallel/StdMpi.h" #include "parallel/StdMpi.h"
namespace AMDiS { namespace AMDiS {
...@@ -37,6 +38,17 @@ namespace AMDiS { ...@@ -37,6 +38,17 @@ namespace AMDiS {
} }
FeSpaceDofMap::FeSpaceDofMap(MeshLevelData* ld)
: levelData(ld),
dofComm(NULL),
feSpace(NULL),
needGlobalMapping(false),
isNonLocal(false)
{
clear();
}
void FeSpaceDofMap::clear() void FeSpaceDofMap::clear()
{ {
dofMap.clear(); dofMap.clear();
......
...@@ -20,20 +20,22 @@ ...@@ -20,20 +20,22 @@
/** \file FeSpaceMapping.h */ /** \file FeSpaceMapping.h */
#ifndef AMDIS_FE_SPACE_MAPPING_H
#define AMDIS_FE_SPACE_MAPPING_H
#include <mpi.h>
#include <vector> #include <vector>
#include <map> #include <map>
#include <set> #include <set>
#include <petsc.h>
#include <petscis.h>
#include "AMDiS_fwd.h"
#include "parallel/DofComm.h" #include "parallel/DofComm.h"
#include "parallel/MeshLevelData.h"
#include "parallel/MpiHelper.h" #include "parallel/MpiHelper.h"
#include "parallel/ParallelTypes.h" #include "parallel/ParallelTypes.h"
#include "parallel/StdMpi.h" #include "parallel/StdMpi.h"
#include <petscis.h>
#ifndef AMDIS_FE_SPACE_MAPPING_H
#define AMDIS_FE_SPACE_MAPPING_H
namespace AMDiS { namespace AMDiS {
using namespace std; using namespace std;
...@@ -110,15 +112,7 @@ namespace AMDiS { ...@@ -110,15 +112,7 @@ namespace AMDiS {
} }
/// This is the only valid constructur to be used. /// This is the only valid constructur to be used.
FeSpaceDofMap(MeshLevelData* ld) FeSpaceDofMap(MeshLevelData* ld);
: levelData(ld),
dofComm(NULL),
feSpace(NULL),
needGlobalMapping(false),
isNonLocal(false)
{
clear();
}
/// Clears all data of the mapping. /// Clears all data of the mapping.
void clear(); void clear();
...@@ -483,6 +477,21 @@ namespace AMDiS { ...@@ -483,6 +477,21 @@ namespace AMDiS {
int firstComponent, int firstComponent,
int nComponents); int nComponents);
/// Create a parallel distributed PETSc vector based on this mapping.
inline void createVec(Vec &vec)
{
VecCreateMPI(mpiComm, getRankDofs(), getOverallDofs(), &vec);
}
/// Create a parallel distributed PETsc vector based on this mapping but
/// with a different (larger) global size. This is used in multi-level
/// method to embed a local vector into a subdomain spaned by several
/// ranks.
inline void createVec(Vec &vec, int nGlobalRows)
{
VecCreateMPI(mpiComm, getRankDofs(), nGlobalRows, &vec);
}
protected: protected:
/// Insert a new FE space DOF mapping for a given FE space. /// Insert a new FE space DOF mapping for a given FE space.
void addFeSpace(const FiniteElemSpace* feSpace); void addFeSpace(const FiniteElemSpace* feSpace);
......
...@@ -26,6 +26,10 @@ ...@@ -26,6 +26,10 @@
#include <set> #include <set>
#include <map> #include <map>
#include <mpi.h> #include <mpi.h>
#include <petsc.h>
#include <petscsys.h>
#include <petscao.h>
#include <petscksp.h>
#include "AMDiS_fwd.h" #include "AMDiS_fwd.h"
#include "Global.h" #include "Global.h"
...@@ -33,10 +37,6 @@ ...@@ -33,10 +37,6 @@
#include "DOFMatrix.h" #include "DOFMatrix.h"
#include "parallel/MeshDistributor.h" #include "parallel/MeshDistributor.h"
#include "parallel/ParallelCoarseSpaceMatVec.h" #include "parallel/ParallelCoarseSpaceMatVec.h"
#include <petsc.h>
#include <petscsys.h>
#include <petscao.h>
#include <petscksp.h>
namespace AMDiS { namespace AMDiS {
......
...@@ -14,6 +14,8 @@ ...@@ -14,6 +14,8 @@
#include "MatrixVector.h" #include "MatrixVector.h"
#include "parallel/PetscSolverFeti.h" #include "parallel/PetscSolverFeti.h"
#include "parallel/PetscSolverFetiStructs.h" #include "parallel/PetscSolverFetiStructs.h"
#include "parallel/PetscSolverFetiOperators.h"
#include "parallel/PetscSolverFetiTimings.h"
#include "parallel/StdMpi.h" #include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h" #include "parallel/MpiHelper.h"
#include "parallel/PetscSolverGlobalMatrix.h" #include "parallel/PetscSolverGlobalMatrix.h"
...@@ -23,201 +25,6 @@ namespace AMDiS { ...@@ -23,201 +25,6 @@ namespace AMDiS {
using namespace std; using namespace std;
double FetiTimings::fetiSolve = 0.0;
double FetiTimings::fetiSolve01 = 0.0;
double FetiTimings::fetiSolve02 = 0.0;
double FetiTimings::fetiPreconditioner = 0.0;
// y = mat * x
int petscMultMatSchurPrimal(Mat mat, Vec x, Vec y)
{
// S_PiPi = K_PiPi - K_PiB inv(K_BB) K_BPi
void *ctx;
MatShellGetContext(mat, &ctx);
SchurPrimalData* data = static_cast<SchurPrimalData*>(ctx);
MatMult(data->subSolver->getMatInteriorCoarse(), x, data->tmp_vec_b);
data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
MatMult(data->subSolver->getMatCoarseInterior(), data->tmp_vec_b,
data->tmp_vec_primal);
MatMult(data->subSolver->getMatCoarse(), x, y);
VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal);
return 0;
}
// y = mat * x
int petscMultMatFeti(Mat mat, Vec x, Vec y)
{
FUNCNAME("petscMultMatFeti()");
// F = L inv(K_BB) trans(L) + L inv(K_BB) K_BPi inv(S_PiPi) K_PiB inv(K_BB) trans(L)
// => F = L [I + inv(K_BB) K_BPi inv(S_PiPi) K_PiB] inv(K_BB) trans(L)
double wtime = MPI::Wtime();
void *ctx;
MatShellGetContext(mat, &ctx);
FetiData* data = static_cast<FetiData*>(ctx);
MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b);
double wtime01 = MPI::Wtime();
data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);
MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange);
MatMult(data->subSolver->getMatCoarseInterior(),
data->tmp_vec_b, data->tmp_vec_primal);
wtime01 = MPI::Wtime();
KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
FetiTimings::fetiSolve02 += (MPI::Wtime() - wtime01);
MatMult(data->subSolver->getMatInteriorCoarse(),
data->tmp_vec_primal, data->tmp_vec_b);
wtime01 = MPI::Wtime();
data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
FetiTimings::fetiSolve01 += (MPI::Wtime() - wtime01);
MatMult(*(data->mat_lagrange), data->tmp_vec_b, y);
VecAXPBY(y, 1.0, 1.0, data->tmp_vec_lagrange);
FetiTimings::fetiSolve += (MPI::Wtime() - wtime);
return 0;
}
// y = PC * x
PetscErrorCode petscApplyFetiDirichletPrecon(PC pc, Vec x, Vec y)
{
double wtime = MPI::Wtime();
// Get data for the preconditioner
void *ctx;
PCShellGetContext(pc, &ctx);
FetiDirichletPreconData* data = static_cast<FetiDirichletPreconData*>(ctx);
// Multiply with scaled Lagrange constraint matrix.
MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);
// === Restriction of the B nodes to the boundary nodes. ===
int nLocalB;
int nLocalDuals;
VecGetLocalSize(data->tmp_vec_b, &nLocalB);
VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
PetscScalar *local_b, *local_duals;
VecGetArray(data->tmp_vec_b, &local_b);
VecGetArray(data->tmp_vec_duals0, &local_duals);
for (map<int, int>::iterator it = data->localToDualMap.begin();
it != data->localToDualMap.end(); ++it)
local_duals[it->second] = local_b[it->first];
VecRestoreArray(data->tmp_vec_b, &local_b);
VecRestoreArray(data->tmp_vec_duals0, &local_duals);
// === K_DD - K_DI inv(K_II) K_ID ===
MatMult(*(data->mat_duals_duals), data->tmp_vec_duals0, data->tmp_vec_duals1);
MatMult(*(data->mat_interior_duals), data->tmp_vec_duals0, data->tmp_vec_interior);
KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
MatMult(*(data->mat_duals_interior), data->tmp_vec_interior, data->tmp_vec_duals0);
VecAXPBY(data->tmp_vec_duals0, 1.0, -1.0, data->tmp_vec_duals1);
// === Prolongation from local dual nodes to B nodes.
VecGetArray(data->tmp_vec_b, &local_b);
VecGetArray(data->tmp_vec_duals0, &local_duals);
for (map<int, int>::iterator it = data->localToDualMap.begin();
it != data->localToDualMap.end(); ++it)
local_b[it->first] = local_duals[it->second];
VecRestoreArray(data->tmp_vec_b, &local_b);
VecRestoreArray(data->tmp_vec_duals0, &local_duals);
// Multiply with scaled Lagrange constraint matrix.
MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);
FetiTimings::fetiPreconditioner += (MPI::Wtime() - wtime);
return 0;
}
// y = PC * x
PetscErrorCode petscApplyFetiLumpedPrecon(PC pc, Vec x, Vec y)
{
// Get data for the preconditioner
void *ctx;
PCShellGetContext(pc, &ctx);
FetiLumpedPreconData* data = static_cast<FetiLumpedPreconData*>(ctx);
// Multiply with scaled Lagrange constraint matrix.
MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);
// === Restriction of the B nodes to the boundary nodes. ===
int nLocalB;
int nLocalDuals;
VecGetLocalSize(data->tmp_vec_b, &nLocalB);
VecGetLocalSize(data->tmp_vec_duals0, &nLocalDuals);
PetscScalar *local_b, *local_duals;
VecGetArray(data->tmp_vec_b, &local_b);
VecGetArray(data->tmp_vec_duals0, &local_duals);
for (map<int, int>::iterator it = data->localToDualMap.begin();
it != data->localToDualMap.end(); ++it)
local_duals[it->second] = local_b[it->first];
VecRestoreArray(data->tmp_vec_b, &local_b);
VecRestoreArray(data->tmp_vec_duals0, &local_duals);
// === K_DD ===
MatMult(*(data->mat_duals_duals), data->tmp_vec_duals0, data->tmp_vec_duals1);
// === Prolongation from local dual nodes to B nodes.
VecGetArray(data->tmp_vec_b, &local_b);
VecGetArray(data->tmp_vec_duals1, &local_duals);
for (map<int, int>::iterator it = data->localToDualMap.begin();