Commit c14d2129 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Fixed eps issue in final timestep calculation.

parent 2620c115
...@@ -511,24 +511,24 @@ namespace AMDiS { ...@@ -511,24 +511,24 @@ namespace AMDiS {
return timestep; return timestep;
} }
inline void setLastProcessedTimestep(double t){ inline void setLastProcessedTimestep(double t)
lastProcessedTimestep=t; {
lastProcessedTimestep = t;
} }
inline double getLastProcessedTimestep(){ inline double getLastProcessedTimestep()
{
return lastProcessedTimestep; return lastProcessedTimestep;
} }
/** \brief /// Returns true, if the end time is reached and no more timestep
* Returns true, if the end time is reached and no more timestep /// computations must be done.
* computations must be done.
*/
inline bool reachedEndTime() inline bool reachedEndTime()
{ {
if (nTimesteps > 0) if (nTimesteps > 0)
return !(timestepNumber < nTimesteps); return !(timestepNumber < nTimesteps);
return !(time < endTime - DBL_TOL); return !(fabs(time - endTime) > DBL_TOL);
} }
......
...@@ -140,12 +140,10 @@ namespace AMDiS { ...@@ -140,12 +140,10 @@ namespace AMDiS {
/// Dimension dependent geo index of the current position in traverse. /// Dimension dependent geo index of the current position in traverse.
GeoIndex posIndex; GeoIndex posIndex;
/** \brief /// Number of DOFs at the current traverse position. Examples: independent of
* Number of DOFs at the current traverse position. Examples: independent of /// dimension and degree of basis functions there is only one DOF per vertex.
* dimension and degree of basis functions there is only one DOF per vertex. /// But in 2D and with 3rd degree lagrange basis functions there are two
* But in 2D and with 3rd degree lagrange basis functions there are two /// DOFs per edge.
* DOFs per edge.
*/
int nDofs; int nDofs;
/// Displacement of DOF indices. Used if more than one DOF admin is defined /// Displacement of DOF indices. Used if more than one DOF admin is defined
......
...@@ -154,7 +154,8 @@ namespace AMDiS { ...@@ -154,7 +154,8 @@ namespace AMDiS {
GeoIndex position, GeoIndex position,
int positionIndex) const; int positionIndex) const;
/// Calculates the number of DOFs needed for Lagrange of the given dim and degree. /// Calculates the number of DOFs needed for Lagrange of the given dim
/// and degree.
static int getNumberOfDofs(int dim, int degree); static int getNumberOfDofs(int dim, int degree);
private: private:
......
...@@ -1007,12 +1007,9 @@ namespace AMDiS { ...@@ -1007,12 +1007,9 @@ namespace AMDiS {
Mat qT, jTqT; Mat qT, jTqT;
MatTranspose(mat_augmented_lagrange, MAT_INITIAL_MATRIX, &qT); MatTranspose(mat_augmented_lagrange, MAT_INITIAL_MATRIX, &qT);
// Mat jT; // Mat jT;
MSG("START COMPUTING MAT TRANS\n");
// MatTranspose(mat_lagrange, MAT_INITIAL_MATRIX, &jT); // MatTranspose(mat_lagrange, MAT_INITIAL_MATRIX, &jT);
MSG("DONE\n");
MatTransposeMatMult(mat_lagrange, qT, MAT_INITIAL_MATRIX, PETSC_DEFAULT, MatTransposeMatMult(mat_lagrange, qT, MAT_INITIAL_MATRIX, PETSC_DEFAULT,
&jTqT); &jTqT);
MSG("DONE WITH THIS!\n");
petsc_helper::blockMatMatSolve(subdomain->getSolver(), jTqT, matTmp); petsc_helper::blockMatMatSolve(subdomain->getSolver(), jTqT, matTmp);
MatDestroy(&qT); MatDestroy(&qT);
MatDestroy(&jTqT); MatDestroy(&jTqT);
...@@ -1041,7 +1038,6 @@ namespace AMDiS { ...@@ -1041,7 +1038,6 @@ namespace AMDiS {
MatDestroy(&mat10); MatDestroy(&mat10);
MatDestroy(&mat11); MatDestroy(&mat11);
MatDestroy(&matTmp); MatDestroy(&matTmp);
MSG("FINISHED!\n");
} else { } else {
Mat tmp; Mat tmp;
......
...@@ -252,33 +252,12 @@ namespace AMDiS { ...@@ -252,33 +252,12 @@ namespace AMDiS {
double wtime = MPI::Wtime(); double wtime = MPI::Wtime();
{
MatNullSpace nullSpace;
MatGetNullSpace(mat, &nullSpace);
PetscBool hasConst;
PetscInt nVec;
const Vec *vecs;
// MatNullSpaceGetVecs(nullSpace, &hasConst, &nVec, &vecs);
MatNullSpaceRemove(nullSpace, x, PETSC_NULL);
}
Vec x_interface, x_lagrange, y_interface, y_lagrange; Vec x_interface, x_lagrange, y_interface, y_lagrange;
VecNestGetSubVec(x, 0, &x_interface); VecNestGetSubVec(x, 0, &x_interface);
VecNestGetSubVec(x, 1, &x_lagrange); VecNestGetSubVec(x, 1, &x_lagrange);
VecNestGetSubVec(y, 0, &y_interface); VecNestGetSubVec(y, 0, &y_interface);
VecNestGetSubVec(y, 1, &y_lagrange); VecNestGetSubVec(y, 1, &y_lagrange);
{
int n;
VecGetSize(x_interface, &n);
double sum;
VecSum(x_interface, &sum);
sum = -sum / static_cast<int>(n);
MSG("xbegin = %e\n", sum);
}
void *ctx; void *ctx;
MatShellGetContext(mat, &ctx); MatShellGetContext(mat, &ctx);
FetiData* data = static_cast<FetiData*>(ctx); FetiData* data = static_cast<FetiData*>(ctx);
...@@ -337,16 +316,6 @@ namespace AMDiS { ...@@ -337,16 +316,6 @@ namespace AMDiS {
// y_interface = A_{\Gamma B} tmp_vec_b0 // y_interface = A_{\Gamma B} tmp_vec_b0
MatMult(data->subSolver->getMatCoarseInterior(1), data->tmp_vec_b0, y_interface); MatMult(data->subSolver->getMatCoarseInterior(1), data->tmp_vec_b0, y_interface);
{
int n;
VecGetSize(y_interface, &n);
double sum;
VecSum(y_interface, &sum);
sum = -sum / static_cast<int>(n);
MSG("yend = %e\n", sum);
}
// tmp_vec_primal0 = S_{\Pi\Pi}^{-1} (v_{\Pi} - A_{\Pi B} A_{BB}^{-1} v_{B}) // tmp_vec_primal0 = S_{\Pi\Pi}^{-1} (v_{\Pi} - A_{\Pi B} A_{BB}^{-1} v_{B})
// tmp_vec_interface = A_{\Gamma \Pi} tmp_vec_primal0 // tmp_vec_interface = A_{\Gamma \Pi} tmp_vec_primal0
MatMult(data->subSolver->getMatCoarse(1, 0), data->tmp_vec_primal0, data->tmp_vec_interface); MatMult(data->subSolver->getMatCoarse(1, 0), data->tmp_vec_primal0, data->tmp_vec_interface);
...@@ -358,6 +327,7 @@ namespace AMDiS { ...@@ -358,6 +327,7 @@ namespace AMDiS {
FetiTimings::fetiSolve += (MPI::Wtime() - wtime); FetiTimings::fetiSolve += (MPI::Wtime() - wtime);
return 0; return 0;
} }
...@@ -497,8 +467,8 @@ namespace AMDiS { ...@@ -497,8 +467,8 @@ namespace AMDiS {
VecNestGetSubVec(yvec, 0, &y_interface); VecNestGetSubVec(yvec, 0, &y_interface);
VecNestGetSubVec(yvec, 1, &y_lagrange); VecNestGetSubVec(yvec, 1, &y_lagrange);
VecCopy(x_interface, y_interface); //VecCopy(x_interface, y_interface);
// KSPSolve(data->ksp_mass, x_interface, y_interface); KSPSolve(data->ksp_mass, x_interface, y_interface);
MatMultTranspose(*(data->mat_lagrange_scaled), x_lagrange, data->tmp_vec_b0); MatMultTranspose(*(data->mat_lagrange_scaled), x_lagrange, data->tmp_vec_b0);
......
...@@ -157,7 +157,7 @@ namespace AMDiS { ...@@ -157,7 +157,7 @@ namespace AMDiS {
if (dirichletRows.count(*cursor)) { if (dirichletRows.count(*cursor)) {
if ((!isRowCoarse && !(*interiorMap)[rowComponent].isRankDof(*cursor)) || if ((!isRowCoarse && !(*interiorMap)[rowComponent].isRankDof(*cursor)) ||
(isRowCoarse && !(*rowCoarseSpace)[rowComponent].isRankDof(*cursor))) (isRowCoarse && !(*rowCoarseSpace)[rowComponent].isRankDof(*cursor)))
continue; continue;
} }
cols.clear(); cols.clear();
...@@ -763,11 +763,16 @@ namespace AMDiS { ...@@ -763,11 +763,16 @@ namespace AMDiS {
if (rankOnly && !(*interiorMap)[rowComp].isRankDof(dof)) if (rankOnly && !(*interiorMap)[rowComp].isRankDof(dof))
continue; continue;
bool isCoarseDof = isCoarseSpace(rowComp, dof);
// Dirichlet rows can be set only be the owner ranks. // Dirichlet rows can be set only be the owner ranks.
if (dirichletValues.count(dof) && !((*interiorMap)[rowComp].isRankDof(dof))) if (dirichletValues.count(dof)) {
continue; if ((!isCoarseDof && !((*interiorMap)[rowComp].isRankDof(dof))) ||
(isCoarseDof && !((*rowCoarseSpace)[rowComp].isRankDof(dof))))
continue;
}
if (isCoarseSpace(rowComp, dof)) { if (isCoarseDof) {
TEST_EXIT_DBG(vecCoarse != PETSC_NULL)("Should not happen!\n"); TEST_EXIT_DBG(vecCoarse != PETSC_NULL)("Should not happen!\n");
int index = rowCoarseSpace->getMatIndex(rowComp, dof); int index = rowCoarseSpace->getMatIndex(rowComp, dof);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment