diff --git a/AMDiS/src/ProblemInstat.h b/AMDiS/src/ProblemInstat.h
index d4f5b76826564af5bf0f03543bb0c9a1ac90dc08..9528d0bdbf42ebb73849e76fc988ccb6600c1d73 100644
--- a/AMDiS/src/ProblemInstat.h
+++ b/AMDiS/src/ProblemInstat.h
@@ -64,7 +64,7 @@ namespace AMDiS {
virtual void solve(AdaptInfo* adaptInfo) {}
- virtual void solve(AdaptInfo *adaptInfo, bool fixedMatrix)
+ virtual void solve(AdaptInfo *adaptInfo, bool, bool)
{
solve(adaptInfo);
}
diff --git a/AMDiS/src/ProblemStat.cc b/AMDiS/src/ProblemStat.cc
index 79bc39b8b05cf4a211401c63e60c7cffd04ec924..91686c7bf0b4a0779f9e2d00a3e103702e02e565 100644
--- a/AMDiS/src/ProblemStat.cc
+++ b/AMDiS/src/ProblemStat.cc
@@ -540,7 +540,7 @@ namespace AMDiS {
}
- void ProblemStatSeq::solve(AdaptInfo *adaptInfo, bool fixedMatrix)
+ void ProblemStatSeq::solve(AdaptInfo *adaptInfo, bool, bool)
{
FUNCNAME("Problem::solve()");
diff --git a/AMDiS/src/ProblemStat.h b/AMDiS/src/ProblemStat.h
index 35a7230f87a03085d55b3268da2847581aad96e3..4c487044df1ebf351ebf3ebd7aef2afeb26011a5 100644
--- a/AMDiS/src/ProblemStat.h
+++ b/AMDiS/src/ProblemStat.h
@@ -67,7 +67,7 @@ namespace AMDiS {
public:
/// Constructor
ProblemStatSeq(string nameStr,
- ProblemIterationInterface *problemIteration = NULL)
+ ProblemIterationInterface *problemIteration = NULL)
: StandardProblemIteration(this),
name(nameStr),
nComponents(-1),
@@ -147,7 +147,9 @@ namespace AMDiS {
* Implementation of ProblemStatBase::solve(). Deligates the solving
* of problems system to \ref solver.
*/
- virtual void solve(AdaptInfo *adaptInfo, bool fixedMatrix = false);
+ void solve(AdaptInfo *adaptInfo,
+ bool createMatrixData = true,
+ bool storeMatrixData = false);
/** \brief
* Implementation of ProblemStatBase::estimate(). Deligates the estimation
diff --git a/AMDiS/src/ProblemStatBase.h b/AMDiS/src/ProblemStatBase.h
index ce544ea6cb9cdf3223fa97351a87113219fec0bb..8ba55350049d89bb364c650fedb91a3956d1dfa3 100644
--- a/AMDiS/src/ProblemStatBase.h
+++ b/AMDiS/src/ProblemStatBase.h
@@ -104,8 +104,27 @@ namespace AMDiS {
/// Coarsening of the mesh.
virtual Flag coarsenMesh(AdaptInfo *adaptInfo) = 0;
- /// Solves the assembled system. The result is an approximative solution.
- virtual void solve(AdaptInfo *adaptInfo, bool fixedMatrix) = 0;
+ /** \brief
+ * Solves the assembled system. The result is an approximative solution.
+ * The last two boolean arguments can be used to controll successive
+ * solutions of systems with the same matrix.
+ *
+ * \param adaptInfo Pointer to an \ref AdaptInfo object.
+ * \param createMatrixData If false, the solver assumes that all of its
+ * internal data structures for the system
+ * matrix are already created. This is the case,
+ * if we solve different systems but with the
+ * same matrix. After the first call to this
+ * function (with this parameter set to true),
+ * all other calls may set it to false.
+ * \param storeMatrixData If true, all internal data structures for the
+ * system matrix are not deleted such that they
+ * can be used for next solutions with the same
+ * system matrix.
+ */
+ virtual void solve(AdaptInfo *adaptInfo,
+ bool createMatrixData = true,
+ bool storeMatrixData = false) = 0;
/** \brief
* A posteriori error estimation of the calculated solution. Should store
diff --git a/AMDiS/src/StandardProblemIteration.cc b/AMDiS/src/StandardProblemIteration.cc
index 135cc613b8745eb2eb1e64261ef4841ad6d4c808..148d9a50d66c56a8192e358d2df1d6f00dbb44f6 100644
--- a/AMDiS/src/StandardProblemIteration.cc
+++ b/AMDiS/src/StandardProblemIteration.cc
@@ -37,10 +37,10 @@ namespace AMDiS {
Flag flag = buildAndAdapt(adaptInfo, toDo);
if (toDo.isSet(SOLVE))
- problem->solve(adaptInfo, false);
+ problem->solve(adaptInfo, true, false);
if (toDo.isSet(SOLVE_RHS))
- problem->solve(adaptInfo, true);
+ problem->solve(adaptInfo, true, false);
if (toDo.isSet(ESTIMATE))
problem->estimate(adaptInfo);
diff --git a/AMDiS/src/parallel/PetscProblemStat.cc b/AMDiS/src/parallel/PetscProblemStat.cc
index 58ff640377031189afe44f0b734794df154fde96..89eb3f7d069f393aa5148e2c683d5c2c02fa12fb 100644
--- a/AMDiS/src/parallel/PetscProblemStat.cc
+++ b/AMDiS/src/parallel/PetscProblemStat.cc
@@ -61,7 +61,9 @@ namespace AMDiS {
}
- void PetscProblemStat::solve(AdaptInfo *adaptInfo, bool fixedMatrix)
+ void PetscProblemStat::solve(AdaptInfo *adaptInfo,
+ bool createMatrixData,
+ bool storeMatrixData)
{
FUNCNAME("PetscProblemStat::solve()");
@@ -79,8 +81,13 @@ namespace AMDiS {
MSG("Overall memory usage is VM = %.1f MB RSS = %.1f MB\n", vm, rss);
#endif
- petscSolver->setMeshDistributor(meshDistributor);
- petscSolver->fillPetscMatrix(systemMatrix, rhs);
+ if (createMatrixData) {
+ petscSolver->setMeshDistributor(meshDistributor);
+ petscSolver->fillPetscMatrix(systemMatrix);
+ }
+
+ petscSolver->fillPetscRhs(rhs);
+
#if 0
processMemUsage(vm, rss);
@@ -93,6 +100,9 @@ namespace AMDiS {
petscSolver->solvePetscMatrix(*solution, adaptInfo);
+ if (!storeMatrixData)
+ petscSolver->destroyMatrixData();
+
#if 0
processMemUsage(vm, rss);
MSG("STAGE 3\n");
@@ -102,7 +112,6 @@ namespace AMDiS {
MSG("Overall memory usage is VM = %.1f MB RSS = %.1f MB\n", vm, rss);
#endif
-
INFO(info, 8)("solution of discrete system needed %.5f seconds\n",
MPI::Wtime() - wtime);
}
diff --git a/AMDiS/src/parallel/PetscProblemStat.h b/AMDiS/src/parallel/PetscProblemStat.h
index f2bb85ab844a972838142215fafa50ae1d4e92b2..9638cc7cf467aa4202c4949989acccb8c692b2fc 100644
--- a/AMDiS/src/parallel/PetscProblemStat.h
+++ b/AMDiS/src/parallel/PetscProblemStat.h
@@ -48,7 +48,9 @@ namespace AMDiS {
ProblemStatSeq *adoptProblem = NULL,
Flag adoptFlag = INIT_NOTHING);
- void solve(AdaptInfo *adaptInfo, bool fixedMatrix = false);
+ void solve(AdaptInfo *adaptInfo,
+ bool createMatrixData = true,
+ bool storeMatrixData = false);
protected:
PetscSolver *petscSolver;
diff --git a/AMDiS/src/parallel/PetscSolver.h b/AMDiS/src/parallel/PetscSolver.h
index a0d28124b0ddcb09743b9c648b09e191224b8523..925d9cd4579b2f47a3203963f4f2962b51a757b9 100644
--- a/AMDiS/src/parallel/PetscSolver.h
+++ b/AMDiS/src/parallel/PetscSolver.h
@@ -60,18 +60,26 @@ namespace AMDiS {
}
/** \brief
- * Create a PETSc matrix and PETSc vectors. The given DOF matrices are used to
- * create the nnz structure of the PETSc matrix and the values are transfered to it.
- * The given DOF vectors are used to the the values of the PETSc rhs vector.
+ * Create a PETSc matrix. The given DOF matrices are used to create the nnz
+ * structure of the PETSc matrix and the values are transfered to it.
*
* \param[in] mat
+ */
+ virtual void fillPetscMatrix(Matrix<DOFMatrix*> *mat) = 0;
+
+ /** \brief
+ * Create a PETSc vector and fills it with the rhs values of the system.
+ *
* \param[in] vec
*/
- virtual void fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec) = 0;
+ virtual void fillPetscRhs(SystemVector *vec) = 0;
/// Use PETSc to solve the linear system of equations
virtual void solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo) = 0;
+ /// Destroys all matrix data structures.
+ virtual void destroyMatrixData() = 0;
+
virtual Flag getBoundaryDofRequirement()
{
return 0;
diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc
index ab4dfe212c8e05b58f9fe43cda3fdc6e5ac81a39..222af4d66b172dd8c26e8709250267b9e7422be0 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.cc
+++ b/AMDiS/src/parallel/PetscSolverFeti.cc
@@ -595,8 +595,13 @@ namespace AMDiS {
schurPrimalData.mat_b_primal = &mat_b_primal;
schurPrimalData.ksp_b = &ksp_b;
- VecDuplicate(f_b, &(schurPrimalData.tmp_vec_b));
- VecDuplicate(f_primal, &(schurPrimalData.tmp_vec_primal));
+ VecCreateMPI(PETSC_COMM_WORLD,
+ nRankB * nComponents, nOverallB * nComponents,
+ &(schurPrimalData.tmp_vec_b));
+ VecCreateMPI(PETSC_COMM_WORLD,
+ nRankPrimals * nComponents, nOverallPrimals * nComponents,
+ &(schurPrimalData.tmp_vec_primal));
+
MatCreateShell(PETSC_COMM_WORLD,
nRankPrimals * nComponents, nRankPrimals * nComponents,
@@ -643,9 +648,15 @@ namespace AMDiS {
fetiData.ksp_b = &ksp_b;
fetiData.ksp_schur_primal = &ksp_schur_primal;
- VecDuplicate(f_b, &(fetiData.tmp_vec_b0));
- VecDuplicate(f_b, &(fetiData.tmp_vec_b1));
- VecDuplicate(f_primal, &(fetiData.tmp_vec_primal));
+ VecCreateMPI(PETSC_COMM_WORLD,
+ nRankB * nComponents, nOverallB * nComponents,
+ &(fetiData.tmp_vec_b0));
+ VecCreateMPI(PETSC_COMM_WORLD,
+ nRankB * nComponents, nOverallB * nComponents,
+ &(fetiData.tmp_vec_b1));
+ VecCreateMPI(PETSC_COMM_WORLD,
+ nRankPrimals * nComponents, nOverallPrimals * nComponents,
+ &(fetiData.tmp_vec_primal));
MatCreateShell(PETSC_COMM_WORLD,
nRankLagrange * nComponents, nRankLagrange * nComponents,
@@ -681,7 +692,9 @@ namespace AMDiS {
fetiDirichletPreconData.mat_duals_interior = &mat_duals_interior;
fetiDirichletPreconData.ksp_interior = &ksp_interior;
- VecDuplicate(f_b, &(fetiDirichletPreconData.tmp_vec_b));
+ VecCreateMPI(PETSC_COMM_WORLD,
+ nRankB * nComponents, nOverallB * nComponents,
+ &(fetiDirichletPreconData.tmp_vec_b));
MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiDirichletPreconData.tmp_vec_duals0));
MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiDirichletPreconData.tmp_vec_duals1));
MatGetVecs(mat_interior_interior, PETSC_NULL, &(fetiDirichletPreconData.tmp_vec_interior));
@@ -697,7 +710,9 @@ namespace AMDiS {
fetiLumpedPreconData.mat_lagrange_scaled = &mat_lagrange_scaled;
fetiLumpedPreconData.mat_duals_duals = &mat_duals_duals;
- VecDuplicate(f_b, &(fetiLumpedPreconData.tmp_vec_b));
+ VecCreateMPI(PETSC_COMM_WORLD,
+ nRankB * nComponents, nOverallB * nComponents,
+ &(fetiLumpedPreconData.tmp_vec_b));
MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiLumpedPreconData.tmp_vec_duals0));
MatGetVecs(mat_duals_duals, PETSC_NULL, &(fetiLumpedPreconData.tmp_vec_duals1));
@@ -851,12 +866,11 @@ namespace AMDiS {
}
- void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat,
- SystemVector *vec)
+ void PetscSolverFeti::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
{
FUNCNAME("PetscSolverFeti::fillPetscMatrix()");
- nComponents = vec->getSize();
+ nComponents = mat->getSize();
// === Create all sets and indices. ===
@@ -1156,16 +1170,40 @@ namespace AMDiS {
}
- // === Create and fill PETSc's right hand side vectors. ===
+ // === Create and fill PETSc matrix for Lagrange constraints. ===
+
+ createMatLagrange();
+
+
+ // === Create PETSc solver for the Schur complement on primal variables. ===
+
+ createSchurPrimalKsp();
+
+
+ // === Create PETSc solver for the FETI-DP operator. ===
+
+ createFetiKsp();
+
+ // === Create solver for the non primal (thus local) variables. ===
+
+ KSPCreate(PETSC_COMM_WORLD, &ksp_b);
+ KSPSetOperators(ksp_b, mat_b_b, mat_b_b, SAME_NONZERO_PATTERN);
+ KSPSetOptionsPrefix(ksp_b, "solver_b_");
+ KSPSetFromOptions(ksp_b);
+ }
+
- VecCreate(PETSC_COMM_WORLD, &f_b);
- VecSetSizes(f_b, nRankB * nComponents, nOverallB * nComponents);
- VecSetType(f_b, VECMPI);
+ void PetscSolverFeti::fillPetscRhs(SystemVector *vec)
+ {
+ FUNCNAME("PetscSolverFeti::fillPetscRhs()");
+
+ int nComponents = vec->getSize();
- VecCreate(PETSC_COMM_WORLD, &f_primal);
- VecSetSizes(f_primal, nRankPrimals * nComponents,
- nOverallPrimals * nComponents);
- VecSetType(f_primal, VECMPI);
+ VecCreateMPI(PETSC_COMM_WORLD,
+ nRankB * nComponents, nOverallB * nComponents, &f_b);
+ VecCreateMPI(PETSC_COMM_WORLD,
+ nRankPrimals * nComponents,
+ nOverallPrimals * nComponents, &f_primal);
for (int i = 0; i < nComponents; i++) {
DOFVector<double>::Iterator dofIt(vec->getDOFVector(i), USED_DOFS);
@@ -1194,21 +1232,6 @@ namespace AMDiS {
VecAssemblyBegin(f_primal);
VecAssemblyEnd(f_primal);
-
-
- // === Create and fill PETSc matrix for Lagrange constraints. ===
-
- createMatLagrange();
-
-
- // === Create PETSc solver for the Schur complement on primal variables. ===
-
- createSchurPrimalKsp();
-
-
- // === Create PETSc solver for the FETI-DP operator. ===
-
- createFetiKsp();
}
@@ -1216,11 +1239,9 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverFeti::solveFetiMatrix()");
- // Create transpose of Lagrange matrix.
- Mat mat_lagrange_transpose;
+ Mat mat_lagrange_transpose;
MatTranspose(mat_lagrange, MAT_INITIAL_MATRIX, &mat_lagrange_transpose);
-
// === Create nested matrix which will contain the overall FETI system. ===
Mat A;
@@ -1261,12 +1282,8 @@ namespace AMDiS {
// === Create rhs and solution vectors for the overall FETI system. ===
- Vec f;
- VecCreate(PETSC_COMM_WORLD, &f);
- VecSetSizes(f, nRankNest, nOverallNest);
- VecSetType(f, VECMPI);
-
- Vec b;
+ Vec f, b;
+ VecCreateMPI(PETSC_COMM_WORLD, nRankNest, nOverallNest, &f);
VecDuplicate(f, &b);
@@ -1352,14 +1369,7 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverFeti::solveReducedFetiMatrix()");
- // === Create solver for the non primal (thus local) variables. ===
-
- KSPCreate(PETSC_COMM_WORLD, &ksp_b);
- KSPSetOperators(ksp_b, mat_b_b, mat_b_b, SAME_NONZERO_PATTERN);
- KSPSetOptionsPrefix(ksp_b, "solver_b_");
- KSPSetFromOptions(ksp_b);
-
- // RHS and solution vector.
+ // RHS vector.
Vec vec_rhs;
// Some temporary vectors.
@@ -1444,8 +1454,14 @@ namespace AMDiS {
VecDestroy(&tmp_primal0);
VecDestroy(&tmp_primal1);
+ VecDestroy(&f_b);
+ VecDestroy(&f_primal);
+ }
- KSPDestroy(&ksp_b);
+
+ void PetscSolverFeti::destroyMatrixData()
+ {
+ FUNCNAME("PetscSolverFeti::destroyMatrixData()");
MatDestroy(&mat_b_b);
MatDestroy(&mat_primal_primal);
@@ -1453,14 +1469,6 @@ namespace AMDiS {
MatDestroy(&mat_primal_b);
MatDestroy(&mat_lagrange);
- VecDestroy(&f_b);
- VecDestroy(&f_primal);
-
- destroySchurPrimalKsp();
-
- destroyFetiKsp();
-
-
// === Destroy preconditioner data structures. ===
if (fetiPreconditioner != FETI_NONE)
@@ -1471,8 +1479,13 @@ namespace AMDiS {
MatDestroy(&mat_interior_duals);
MatDestroy(&mat_duals_interior);
}
- }
+ destroySchurPrimalKsp();
+
+ destroyFetiKsp();
+
+ KSPDestroy(&ksp_b);
+ }
void PetscSolverFeti::solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo)
{
@@ -1490,6 +1503,7 @@ namespace AMDiS {
}
}
+
#endif
}
diff --git a/AMDiS/src/parallel/PetscSolverFeti.h b/AMDiS/src/parallel/PetscSolverFeti.h
index 4dd72bca95dedbec2b8178d7160d94a178a613b6..6607b79cef06560eb5bf00e23da9a550096432b3 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.h
+++ b/AMDiS/src/parallel/PetscSolverFeti.h
@@ -133,13 +133,19 @@ namespace AMDiS {
public:
PetscSolverFeti();
- /// Assemble the sequentially created matrices and vectors to the
- /// global matrices and vectors required by the FETI-DP method.
- void fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec);
+ /// Assemble the sequentially created matrices to the global matrices
+ /// required by the FETI-DP method.
+ void fillPetscMatrix(Matrix<DOFMatrix*> *mat);
+
+ /// Assembles the global rhs vectors from the sequentially created ones.
+ void fillPetscRhs(SystemVector *vec);
/// Solve the system using FETI-DP method.
void solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo);
+ /// Destroys all matrix data structures.
+ void destroyMatrixData();
+
/// Returns flags to denote which information of the boundary DOFs are
/// required by the FETI-DP solver.
Flag getBoundaryDofRequirement()
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
index 5eff17d991d86f6242d7a3a483a9060ec2730446..0bf7c516ee113bb3cac828d7bcd0e3424808dc9d 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc
@@ -25,36 +25,28 @@ namespace AMDiS {
}
- void PetscSolverGlobalMatrix::fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec)
+ void PetscSolverGlobalMatrix::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
{
FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrix()");
TEST_EXIT_DBG(meshDistributor)("No mesh distributor object defined!\n");
TEST_EXIT_DBG(mat)("No DOF matrix defined!\n");
- TEST_EXIT_DBG(vec)("NO DOF vector defined!\n");
double wtime = MPI::Wtime();
int nComponents = mat->getNumRows();
int nRankRows = meshDistributor->getNumberRankDofs() * nComponents;
int nOverallRows = meshDistributor->getNumberOverallDofs() * nComponents;
- // === Create PETSc vector (rhs, solution and a temporary vector). ===
-
- VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
- VecSetSizes(petscRhsVec, nRankRows, nOverallRows);
- VecSetType(petscRhsVec, VECMPI);
-
- VecCreate(PETSC_COMM_WORLD, &petscSolVec);
- VecSetSizes(petscSolVec, nRankRows, nOverallRows);
- VecSetType(petscSolVec, VECMPI);
+ // === Create PETSc vector (solution and a temporary vector). ===
- VecCreate(PETSC_COMM_WORLD, &petscTmpVec);
- VecSetSizes(petscTmpVec, nRankRows, nOverallRows);
- VecSetType(petscTmpVec, VECMPI);
+ VecCreateMPI(PETSC_COMM_WORLD, nRankRows, nOverallRows, &petscSolVec);
+ VecCreateMPI(PETSC_COMM_WORLD, nRankRows, nOverallRows, &petscTmpVec);
int recvAllValues = 0;
- int sendValue = static_cast<int>(meshDistributor->getLastMeshChangeIndex() != lastMeshNnz);
- meshDistributor->getMpiComm().Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);
+ int sendValue =
+ static_cast<int>(meshDistributor->getLastMeshChangeIndex() != lastMeshNnz);
+ meshDistributor->getMpiComm().Allreduce(&sendValue, &recvAllValues,
+ 1, MPI_INT, MPI_SUM);
if (!d_nnz || recvAllValues != 0) {
if (d_nnz) {
@@ -102,10 +94,21 @@ namespace AMDiS {
MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(petscMatrix, MAT_FINAL_ASSEMBLY);
-#if (DEBUG != 0)
- MSG("Fill petsc matrix 3 needed %.5f seconds\n",
- TIME_USED(MPI::Wtime(), wtime));
-#endif
+ MSG("Fill petsc matrix needed %.5f seconds\n", MPI::Wtime() - wtime);
+ }
+
+
+ void PetscSolverGlobalMatrix::fillPetscRhs(SystemVector *vec)
+ {
+ FUNCNAME("PetscSolverGlobalMatrix::fillPetscRhs()");
+
+ TEST_EXIT_DBG(vec)("NO DOF vector defined!\n");
+
+ int nComponents = vec->getSize();
+ int nRankRows = meshDistributor->getNumberRankDofs() * nComponents;
+ int nOverallRows = meshDistributor->getNumberOverallDofs() * nComponents;
+
+ VecCreateMPI(PETSC_COMM_WORLD, nRankRows, nOverallRows, &petscRhsVec);
// === Transfer values from DOF vector to the PETSc vector. ===
@@ -114,12 +117,11 @@ namespace AMDiS {
VecAssemblyBegin(petscRhsVec);
VecAssemblyEnd(petscRhsVec);
-
- MSG("Fill petsc matrix needed %.5f seconds\n", MPI::Wtime() - wtime);
}
- void PetscSolverGlobalMatrix::solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo)
+ void PetscSolverGlobalMatrix::solvePetscMatrix(SystemVector &vec,
+ AdaptInfo *adaptInfo)
{
FUNCNAME("PetscSolverGlobalMatrix::solvePetscMatrix()");
diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.h b/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
index c3ef0c0632e4a32cbb7add5562333880fd82b9af..34bc9af1ba7962699ecfbcee4a9d23dfe9cb6d63 100644
--- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
+++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.h
@@ -44,10 +44,15 @@ namespace AMDiS {
Parameters::get("parallel->use zero start vector", zeroStartVector);
}
- void fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec);
+ void fillPetscMatrix(Matrix<DOFMatrix*> *mat);
+
+ void fillPetscRhs(SystemVector *vec);
void solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo);
+ void destroyMatrixData()
+ {}
+
protected:
/// Creates a new non zero pattern structure for the PETSc matrix.
void createPetscNnzStructure(Matrix<DOFMatrix*> *mat);
diff --git a/AMDiS/src/parallel/PetscSolverSchur.cc b/AMDiS/src/parallel/PetscSolverSchur.cc
index a4c5761eb7db68858565e1851f6b08748de3a111..fc39ef164a6ff41643726317586dd3d23e92d856 100644
--- a/AMDiS/src/parallel/PetscSolverSchur.cc
+++ b/AMDiS/src/parallel/PetscSolverSchur.cc
@@ -178,7 +178,7 @@ namespace AMDiS {
}
- void PetscSolverSchur::fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec)
+ void PetscSolverSchur::fillPetscMatrix(Matrix<DOFMatrix*> *mat)
{
FUNCNAME("PetscSolverSchur::fillPetscMatrix()");
@@ -248,18 +248,20 @@ namespace AMDiS {
int nRankRows = meshDistributor->getNumberRankDofs() * nComponents;
int nOverallRows = meshDistributor->getNumberOverallDofs() * nComponents;
- VecCreate(PETSC_COMM_WORLD, &petscRhsVec);
- VecSetSizes(petscRhsVec, nRankRows, nOverallRows);
- VecSetType(petscRhsVec, VECMPI);
+ VecCreateMPI(PETSC_COMM_WORLD, nRankRows, nOverallRows, &petscSolVec);
+ VecCreateMPI(PETSC_COMM_WORLD, nRankRows, nOverallRows, &petscTmpVec);
+ }
+
- VecCreate(PETSC_COMM_WORLD, &petscSolVec);
- VecSetSizes(petscSolVec, nRankRows, nOverallRows);
- VecSetType(petscSolVec, VECMPI);
+ void PetscSolverSchur::fillPetscRhs(SystemVector *vec)
+ {
+ FUNCNAME("PetscSolverSchur::fillPetscRhs()");
- VecCreate(PETSC_COMM_WORLD, &petscTmpVec);
- VecSetSizes(petscTmpVec, nRankRows, nOverallRows);
- VecSetType(petscTmpVec, VECMPI);
+ int nComponents = vec->getSize();
+ int nRankRows = meshDistributor->getNumberRankDofs() * nComponents;
+ int nOverallRows = meshDistributor->getNumberOverallDofs() * nComponents;
+ VecCreateMPI(PETSC_COMM_WORLD, nRankRows, nOverallRows, &petscRhsVec);
for (int i = 0; i < nComponents; i++)
setDofVector(petscRhsVec, vec->getDOFVector(i), nComponents, i);
@@ -269,7 +271,8 @@ namespace AMDiS {
}
- void PetscSolverSchur::solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo)
+ void PetscSolverSchur::solvePetscMatrix(SystemVector &vec,
+ AdaptInfo *adaptInfo)
{
FUNCNAME("PetscSolverSchur::solvePetscMatrix()");
diff --git a/AMDiS/src/parallel/PetscSolverSchur.h b/AMDiS/src/parallel/PetscSolverSchur.h
index 8f9034f3ec6ca4fc78f400bed24acf0c65c9b167..f6ee16e7e175743b97321daaca1683a5da837067 100644
--- a/AMDiS/src/parallel/PetscSolverSchur.h
+++ b/AMDiS/src/parallel/PetscSolverSchur.h
@@ -38,10 +38,15 @@ namespace AMDiS {
: PetscSolver()
{}
- void fillPetscMatrix(Matrix<DOFMatrix*> *mat, SystemVector *vec);
+ void fillPetscMatrix(Matrix<DOFMatrix*> *mat);
+
+ void fillPetscRhs(SystemVector *vec);
void solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo);
+ void destroyMatrixData()
+ {}
+
Flag getBoundaryDofRequirement()
{
return
diff --git a/AMDiS/src/time/RosenbrockAdaptInstationary.cc b/AMDiS/src/time/RosenbrockAdaptInstationary.cc
index 2b1477b081c999e7970fa2ad95c7b621cb024c38..262521edc7272e8eaa9db8572ff608835b7407a4 100644
--- a/AMDiS/src/time/RosenbrockAdaptInstationary.cc
+++ b/AMDiS/src/time/RosenbrockAdaptInstationary.cc
@@ -53,15 +53,15 @@ namespace AMDiS {
FUNCNAME("RosenbrockAdaptInstationary::RosenbrockAdaptInstationary()");
initConstructor(&problemStat);
}
- void RosenbrockAdaptInstationary::initConstructor(
- RosenbrockStationary *problemStat)
- {
+
+ void RosenbrockAdaptInstationary::initConstructor(RosenbrockStationary *problemStat)
+ {
std::string str;
Parameters::get(name + "->rosenbrock method", str);
RosenbrockMethodCreator *creator =
- dynamic_cast<RosenbrockMethodCreator*>(CreatorMap<RosenbrockMethod>::getCreator(str));
+ dynamic_cast<RosenbrockMethodCreator*>(CreatorMap<RosenbrockMethod>::getCreator(str));
rosenbrockMethod = creator->create();
-
+
TEST_EXIT_DBG(rosenbrockMethod)("This should not happen!\n");
Parameters::get(name + "->fix first timesteps", fixFirstTimesteps);
@@ -71,17 +71,16 @@ namespace AMDiS {
problemStat->setTauGamma(&tauGamma);
problemStat->setTau(&tau);
- }
-
+ }
- void RosenbrockAdaptInstationary::oneTimestep()
- {
+ void RosenbrockAdaptInstationary::oneTimestep()
+ {
FUNCNAME("RosenbrockAdaptInstationary::oneTimestep()");
-
+
// estimate before first adaption
if (adaptInfo->getTime() <= adaptInfo->getStartTime())
problemIteration->oneIteration(adaptInfo, ESTIMATE);
-
+
bool rejected = false;
double timeTol = adaptInfo->getTimeTolerance(0);
@@ -94,7 +93,7 @@ namespace AMDiS {
INFO(info, 6)("time = %e timestep = %e\n",
adaptInfo->getTime(), adaptInfo->getTimestep());
-
+
adaptInfo->setSpaceIteration(0);
// do the iteration
@@ -147,8 +146,8 @@ namespace AMDiS {
lastTimestepRejected = true;
}
-
-
+
+
if (errorEst < timeTol || fixFirstTimesteps
|| !(adaptInfo->getTimestep()>adaptInfo->getMinTimestep()) ) {
INFO(info, 6)("Accepted timestep at time = %e with timestep = %e\n",
diff --git a/AMDiS/src/time/RosenbrockAdaptInstationary.h b/AMDiS/src/time/RosenbrockAdaptInstationary.h
index e9fc31a9452b0849c8f07b70b8c92341120113b1..4d55c6897211c37314f48d246b88dbe5ba10aa25 100644
--- a/AMDiS/src/time/RosenbrockAdaptInstationary.h
+++ b/AMDiS/src/time/RosenbrockAdaptInstationary.h
@@ -34,11 +34,11 @@ namespace AMDiS {
class RosenbrockAdaptInstationary : public AdaptInstationary
{
public:
-/** \brief
- * Creates a AdaptInstationary object for Rosenbrock method
- * with the given name for the time
- * dependent problem problemInstat. TODO: Make obsolete!
- */
+ /** \brief
+ * Creates a AdaptInstationary object for Rosenbrock method
+ * with the given name for the time
+ * dependent problem problemInstat. TODO: Make obsolete!
+ */
RosenbrockAdaptInstationary(std::string name,
RosenbrockStationary *problemStat,
AdaptInfo *info,
@@ -46,11 +46,11 @@ namespace AMDiS {
AdaptInfo *initialInfo,
time_t initialTimestamp = 0);
-/** \brief
- * Creates a AdaptInstationary object for Rosenbrock method
- * with the given name for the time
- * dependent problem problemInstat.
- */
+ /** \brief
+ * Creates a AdaptInstationary object for Rosenbrock method
+ * with the given name for the time
+ * dependent problem problemInstat.
+ */
RosenbrockAdaptInstationary(std::string name,
RosenbrockStationary &problemStat,
AdaptInfo &info,
@@ -60,7 +60,8 @@ namespace AMDiS {
/// Runs the Rosenbrock loop until one timestep is accepted.
void oneTimestep();
-/** \brief
+
+ /** \brief
* This funciton is used only to avoid double code in both constructors. If the
* obsolte constructure, which uses pointers instead of references, will be
* removed, remove also this function.
diff --git a/AMDiS/src/time/RosenbrockStationary.cc b/AMDiS/src/time/RosenbrockStationary.cc
index 87ec8da9ec9ffb7194f570c2ef8ecb14765cd288..f6f4495c66a220c6b91154de7cc93774eb5ae928 100644
--- a/AMDiS/src/time/RosenbrockStationary.cc
+++ b/AMDiS/src/time/RosenbrockStationary.cc
@@ -11,6 +11,7 @@
#include "time/RosenbrockStationary.h"
+#include "io/VtkWriter.h"
#include "ProblemStat.h"
#include "SystemVector.h"
#include "OEMSolver.h"
@@ -32,10 +33,15 @@ namespace AMDiS {
newUn = new SystemVector(*solution);
tmp = new SystemVector(*solution);
lowSol = new SystemVector(*solution);
+
+ stageSolution->set(0.0);
+ unVec->set(0.0);
stageSolutions.resize(rm->getStages());
- for (int i = 0; i < rm->getStages(); i++)
+ for (int i = 0; i < rm->getStages(); i++) {
stageSolutions[i] = new SystemVector(*solution);
+ stageSolutions[i]->set(0.0);
+ }
}
@@ -47,8 +53,9 @@ namespace AMDiS {
TEST_EXIT(tauPtr)("No tau pointer defined in stationary problem!\n");
if (first) {
+ MeshDistributor::globalMeshDistributor->synchVector(*solution);
first = false;
- *unVec = *solution;
+ *unVec = *solution;
}
*newUn = *unVec;
@@ -62,16 +69,11 @@ namespace AMDiS {
*stageSolution += *tmp;
}
- clock_t first = clock();
-
for (unsigned int j = 0; j < boundaries.size(); j++) {
boundaries[j].vec->interpol(boundaries[j].fct);
*(boundaries[j].vec) -= *(stageSolution->getDOFVector(boundaries[j].row));
}
- MSG("Solution %.5f seconds\n", TIME_USED(first, clock()));
-
-
timeRhsVec->set(0.0);
for (int j = 0; j < i; j++) {
*tmp = *(stageSolutions[j]);
@@ -79,16 +81,15 @@ namespace AMDiS {
*timeRhsVec += *tmp;
}
- ProblemStatSeq::buildAfterCoarsen(adaptInfo, flag, (i == 0), asmVector);
- // solver->setMultipleRhs(i != 0);
- ProblemStatSeq::solve(adaptInfo);
-
+ ProblemStat::buildAfterCoarsen(adaptInfo, flag, (i == 0), asmVector);
+ ProblemStat::solve(adaptInfo, i == 0, i + 1 < rm->getStages());
+
*(stageSolutions[i]) = *solution;
*tmp = *solution;
*tmp *= rm->getM1(i);
*newUn += *tmp;
-
+
*tmp = *solution;
*tmp *= rm->getM2(i);
*lowSol += *tmp;
@@ -97,11 +98,11 @@ namespace AMDiS {
for (int i = 0; i < nComponents; i++) {
(*(lowSol->getDOFVector(i))) -= (*(newUn->getDOFVector(i)));
adaptInfo->setTimeEstSum(lowSol->getDOFVector(i)->l2norm(), i);
- }
+ }
}
- void RosenbrockStationary::solve(AdaptInfo *adaptInfo, bool fixedMatrix)
+ void RosenbrockStationary::solve(AdaptInfo *adaptInfo, bool, bool)
{}
@@ -113,7 +114,7 @@ namespace AMDiS {
TEST_EXIT(op.getUhOld() == NULL)("UhOld is not allowed to be set!\n");
op.setUhOld(stageSolution->getDOFVector(col));
- ProblemStatSeq::addVectorOperator(op, row);
+ ProblemStat::addVectorOperator(op, row);
}
@@ -125,7 +126,7 @@ namespace AMDiS {
TEST_EXIT(factor == NULL)("Not yet implemented!\n");
TEST_EXIT(estFactor == NULL)("Not yet implemented!\n");
- ProblemStatSeq::addMatrixOperator(op, row, col, &minusOne, &minusOne);
+ ProblemStat::addMatrixOperator(op, row, col, &minusOne, &minusOne);
}
@@ -137,11 +138,11 @@ namespace AMDiS {
Operator *op = new Operator(componentSpaces[row], componentSpaces[col]);
op->addZeroOrderTerm(new Simple_ZOT);
- ProblemStatSeq::addMatrixOperator(op, row, col, tauGamma, tauGamma);
+ ProblemStat::addMatrixOperator(op, row, col, tauGamma, tauGamma);
Operator *opRhs = new Operator(componentSpaces[row]);
opRhs->addZeroOrderTerm(new VecAtQP_ZOT(timeRhsVec->getDOFVector(col), new IdFunc()));
- ProblemStatSeq::addVectorOperator(opRhs, row, &minusOne, &minusOne);
+ ProblemStat::addVectorOperator(opRhs, row, &minusOne, &minusOne);
}
@@ -154,7 +155,7 @@ namespace AMDiS {
RosenbrockBoundary bound = {fct, vec, row, col};
boundaries.push_back(bound);
- ProblemStatSeq::addDirichletBC(type, row, col, vec);
+ ProblemStat::addDirichletBC(type, row, col, vec);
}
}
diff --git a/AMDiS/src/time/RosenbrockStationary.h b/AMDiS/src/time/RosenbrockStationary.h
index a4f545da39f86abc6c954efb50615ffa5e292a46..e5e3e09f19568bac9292cb96103c6b075463ea10 100644
--- a/AMDiS/src/time/RosenbrockStationary.h
+++ b/AMDiS/src/time/RosenbrockStationary.h
@@ -28,6 +28,10 @@
#include "SystemVector.h"
#include "time/RosenbrockMethod.h"
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+#include "parallel/PetscProblemStat.h"
+#endif
+
namespace AMDiS {
struct RosenbrockBoundary
@@ -42,7 +46,7 @@ namespace AMDiS {
};
- class RosenbrockStationary : public ProblemStatSeq
+ class RosenbrockStationary : public ProblemStat
{
public:
class IdFunc : public AbstractFunction<double, double>
@@ -60,7 +64,7 @@ namespace AMDiS {
RosenbrockStationary(std::string name)
- : ProblemStatSeq(name),
+ : ProblemStat(name),
first(true),
minusOne(-1.0),
tauPtr(NULL),
@@ -82,7 +86,7 @@ namespace AMDiS {
void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
bool asmMatrix, bool asmVector);
- void solve(AdaptInfo *adaptInfo, bool fixedMatrix);
+ void solve(AdaptInfo *adaptInfo, bool, bool);
void setRosenbrockMethod(RosenbrockMethod *method)
{