From c2eaea75bb3445f206411e8e895379d887426044 Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Thu, 1 Jun 2017 17:57:00 +0200
Subject: [PATCH] [PETSc] Moved the calling of setLocalAccessibleVector to
Process members
---
MathLib/LinAlg/LinAlg.cpp | 4 +--
MathLib/LinAlg/LinAlg.h | 4 +--
MathLib/LinAlg/PETSc/PETScVector.cpp | 6 ++--
MathLib/LinAlg/PETSc/PETScVector.h | 8 ++---
NumLib/ODESolver/NonlinearSolver.cpp | 23 --------------
NumLib/ODESolver/TimeDiscretization.h | 31 -------------------
NumLib/ODESolver/TimeDiscretizedODESystem.cpp | 7 -----
ProcessLib/Process.cpp | 13 ++++++++
ProcessLib/UncoupledProcessesTimeLoop.cpp | 16 ++--------
Tests/NumLib/ODEs.h | 4 +++
10 files changed, 30 insertions(+), 86 deletions(-)
diff --git a/MathLib/LinAlg/LinAlg.cpp b/MathLib/LinAlg/LinAlg.cpp
index 3bad1aa5ed3..927fe14b792 100644
--- a/MathLib/LinAlg/LinAlg.cpp
+++ b/MathLib/LinAlg/LinAlg.cpp
@@ -22,7 +22,7 @@ namespace MathLib { namespace LinAlg
// Vector
-void setLocalAccessibleVector(PETScVector& x)
+void setLocalAccessibleVector(PETScVector const& x)
{
x.setLocalAccessibleVector();
}
@@ -182,7 +182,7 @@ namespace MathLib { namespace LinAlg
// Vector
-void setLocalAccessibleVector(EigenVector& /*x*/)
+void setLocalAccessibleVector(EigenVector const& /*x*/)
{
}
diff --git a/MathLib/LinAlg/LinAlg.h b/MathLib/LinAlg/LinAlg.h
index 43ad6a018c8..0dee45dd135 100644
--- a/MathLib/LinAlg/LinAlg.h
+++ b/MathLib/LinAlg/LinAlg.h
@@ -150,7 +150,7 @@ namespace LinAlg
/// Set local accessible vector in order to get entries.
/// Call this before call operator[] or get(...).
-void setLocalAccessibleVector(PETScVector& x);
+void setLocalAccessibleVector(PETScVector const& x);
void set(PETScVector& x, double const a);
@@ -212,7 +212,7 @@ namespace LinAlg
/// Set local accessible vector in order to get entries.
/// Call this before call operator[] or get(...). Not used for EIGEN
-void setLocalAccessibleVector(EigenVector& x);
+void setLocalAccessibleVector(EigenVector const& x);
void set(EigenVector& x, double const a);
diff --git a/MathLib/LinAlg/PETSc/PETScVector.cpp b/MathLib/LinAlg/PETSc/PETScVector.cpp
index d1d40e7ddeb..a9a37c676ae 100644
--- a/MathLib/LinAlg/PETSc/PETScVector.cpp
+++ b/MathLib/LinAlg/PETSc/PETScVector.cpp
@@ -104,7 +104,7 @@ void PETScVector::finalizeAssembly()
}
void PETScVector::gatherLocalVectors( PetscScalar local_array[],
- PetscScalar global_array[])
+ PetscScalar global_array[]) const
{
// Collect vectors from processors.
int size_rank;
@@ -130,7 +130,7 @@ void PETScVector::gatherLocalVectors( PetscScalar local_array[],
}
-void PETScVector::getGlobalVector(PetscScalar u[])
+void PETScVector::getGlobalVector(PetscScalar u[]) const
{
#ifdef TEST_MEM_PETSC
@@ -156,7 +156,7 @@ void PETScVector::getGlobalVector(PetscScalar u[])
#endif
}
-void PETScVector::setLocalAccessibleVector()
+void PETScVector::setLocalAccessibleVector() const
{
// TODO: use getLocalVector
if (!_global_v)
diff --git a/MathLib/LinAlg/PETSc/PETScVector.h b/MathLib/LinAlg/PETSc/PETScVector.h
index 67e0d0d5261..fb982912770 100644
--- a/MathLib/LinAlg/PETSc/PETScVector.h
+++ b/MathLib/LinAlg/PETSc/PETScVector.h
@@ -157,7 +157,7 @@ class PETScVector
/// Set local accessible vector in order to get entries.
/// Call this before call operator[] or get(...).
- void setLocalAccessibleVector();
+ void setLocalAccessibleVector() const;
/// Get several entries
std::vector<double> get(std::vector<IndexType> const& indices) const;
@@ -173,7 +173,7 @@ class PETScVector
Get global vector
\param u Array to store the global vector. Memory allocation is needed in advance
*/
- void getGlobalVector(PetscScalar u[]);
+ void getGlobalVector(PetscScalar u[]) const;
/*!
Copy local entries including ghost ones to an array
@@ -275,7 +275,7 @@ class PETScVector
and its associated computations can be dropped if
VecGetValues can get values from different processors.
*/
- std::unique_ptr<PetscScalar[]> _global_v = nullptr;
+ mutable std::unique_ptr<PetscScalar[]> _global_v;
/*!
\brief Collect local vectors
@@ -283,7 +283,7 @@ class PETScVector
\param global_array Global array
*/
void gatherLocalVectors(PetscScalar local_array[],
- PetscScalar global_array[]);
+ PetscScalar global_array[]) const;
/*!
Get local vector, i.e. entries in the same rank
diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index 3a0c4a24046..c8a36c347e1 100644
--- a/NumLib/ODESolver/NonlinearSolver.cpp
+++ b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -55,15 +55,10 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
BaseLib::RunTime time_iteration;
time_iteration.start();
- // The function only has computation if DDC is appied,
- // e.g. Parallel comuting.
- MathLib::LinAlg::setLocalAccessibleVector(x);
-
sys.preIteration(iteration, x);
BaseLib::RunTime time_assembly;
time_assembly.start();
-
sys.assemble(x, coupling_term);
sys.getA(A);
sys.getRhs(rhs);
@@ -93,10 +88,6 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
}
else
{
- // The function only has computation if DDC is appied,
- // e.g. Parallel comuting.
- LinAlg::setLocalAccessibleVector(x_new);
-
if (postIterationCallback)
postIterationCallback(iteration, x_new);
@@ -114,7 +105,6 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
// Thereby the failed solution can be used by the caller for
// debugging purposes.
LinAlg::copy(x_new, x);
- LinAlg::setLocalAccessibleVector(x);
break;
case IterationResult::REPEAT_ITERATION:
INFO(
@@ -146,7 +136,6 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
// Update x s.t. in the next iteration we will compute the right delta x
LinAlg::copy(x_new, x);
- LinAlg::setLocalAccessibleVector(x);
INFO("[time] Iteration #%u took %g s.", iteration,
time_iteration.elapsed());
@@ -210,15 +199,10 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
BaseLib::RunTime time_iteration;
time_iteration.start();
- // The function only has computation if DDC is appied,
- // e.g. Parallel comuting.
- LinAlg::setLocalAccessibleVector(x);
-
sys.preIteration(iteration, x);
BaseLib::RunTime time_assembly;
time_assembly.start();
-
sys.assemble(x, coupling_term);
sys.getResidual(x, res);
sys.getJacobian(J);
@@ -251,11 +235,6 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
x, _x_new_id);
LinAlg::axpy(x_new, -_alpha, minus_delta_x);
- // The function only has computation if DDC is appied,
- // e.g. Parallel comuting.
- LinAlg::setLocalAccessibleVector(x_new);
-
-
if (postIterationCallback)
postIterationCallback(iteration, x_new);
@@ -282,8 +261,6 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
// TODO could be done via swap. Note: that also requires swapping
// the ids. Same for the Picard scheme.
LinAlg::copy(x_new, x); // copy new solution to x
- LinAlg::setLocalAccessibleVector(x);
-
NumLib::GlobalVectorProvider::provider.releaseVector(
x_new);
}
diff --git a/NumLib/ODESolver/TimeDiscretization.h b/NumLib/ODESolver/TimeDiscretization.h
index eec9b707f57..0db170d920e 100644
--- a/NumLib/ODESolver/TimeDiscretization.h
+++ b/NumLib/ODESolver/TimeDiscretization.h
@@ -220,18 +220,12 @@ public:
{
_t = t0;
MathLib::LinAlg::copy(x0, _x_old);
- // The function only has computation if DDC is appied,
- // e.g. Parallel comuting.
- MathLib::LinAlg::setLocalAccessibleVector(_x_old);
}
void pushState(const double /*t*/, GlobalVector const& x,
InternalMatrixStorage const&) override
{
MathLib::LinAlg::copy(x, _x_old);
- // The function only has computation if DDC is appied,
- // e.g. Parallel comuting.
- MathLib::LinAlg::setLocalAccessibleVector(_x_old);
}
void nextTimestep(const double t, const double delta_t) override
@@ -276,18 +270,12 @@ public:
_t = t0;
_t_old = t0;
MathLib::LinAlg::copy(x0, _x_old);
- // The function only has computation if DDC is appied,
- // e.g. Parallel comuting.
- MathLib::LinAlg::setLocalAccessibleVector(_x_old);
}
void pushState(const double /*t*/, GlobalVector const& x,
InternalMatrixStorage const&) override
{
MathLib::LinAlg::copy(x, _x_old);
- // The function only has computation if DDC is appied,
- // e.g. Parallel comuting.
- MathLib::LinAlg::setLocalAccessibleVector(_x_old);
}
void nextTimestep(const double t, const double delta_t) override
@@ -356,18 +344,12 @@ public:
{
_t = t0;
MathLib::LinAlg::copy(x0, _x_old);
- // The function only has computation if DDC is appied,
- // e.g. Parallel comuting.
- MathLib::LinAlg::setLocalAccessibleVector(_x_old);
}
void pushState(const double, GlobalVector const& x,
InternalMatrixStorage const& strg) override
{
MathLib::LinAlg::copy(x, _x_old);
- // The function only has computation if DDC is appied,
- // e.g. Parallel comuting.
- MathLib::LinAlg::setLocalAccessibleVector(_x_old);
strg.pushMatrices();
}
@@ -451,12 +433,6 @@ public:
_t = t0;
_xs_old.push_back(
&NumLib::GlobalVectorProvider::provider.getVector(x0));
- for (auto x_old : _xs_old)
- {
- // The function only has computation if DDC is appied,
- // e.g. Parallel comuting.
- MathLib::LinAlg::setLocalAccessibleVector(*x_old);
- }
}
void pushState(const double, GlobalVector const& x,
@@ -477,13 +453,6 @@ public:
_offset = (_offset + 1) %
_num_steps; // treat _xs_old as a circular buffer
}
-
- for (auto x_old : _xs_old)
- {
- // The function only has computation if DDC is appied,
- // e.g. Parallel comuting.
- MathLib::LinAlg::setLocalAccessibleVector(*x_old);
- }
}
void nextTimestep(const double t, const double delta_t) override
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
index 48de599f3d6..0684b2657f4 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
@@ -10,7 +10,6 @@
#include "TimeDiscretizedODESystem.h"
#include "MathLib/LinAlg/ApplyKnownSolution.h"
-#include "MathLib/LinAlg/LinAlg.h"
#include "MathLib/LinAlg/UnifiedMatrixSetters.h"
#include "NumLib/IndexValueVector.h"
@@ -80,9 +79,6 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
auto& xdot = NumLib::GlobalVectorProvider::provider.getVector(_xdot_id);
_time_disc.getXdot(x_new_timestep, xdot);
- // The function only has computation if DDC is appied,
- // e.g. Parallel comuting.
- MathLib::LinAlg::setLocalAccessibleVector(xdot);
_M->setZero();
_K->setZero();
@@ -110,9 +106,6 @@ void TimeDiscretizedODESystem<
// fragile.
auto& xdot = NumLib::GlobalVectorProvider::provider.getVector(_xdot_id);
_time_disc.getXdot(x_new_timestep, xdot);
- // The function only has computation if DDC is appied,
- // e.g. Parallel comuting.
- MathLib::LinAlg::setLocalAccessibleVector(xdot);
_mat_trans->computeResidual(*_M, *_K, *_b, x_new_timestep, xdot, res);
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 5bf2b92329e..46a0caeb660 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -122,6 +122,10 @@ void Process::assemble(const double t, GlobalVector const& x, GlobalMatrix& M,
GlobalMatrix& K, GlobalVector& b,
StaggeredCouplingTerm const& coupling_term)
{
+ // The function only has computation if DDC is appied,
+ // e.g. Parallel comuting.
+ MathLib::LinAlg::setLocalAccessibleVector(x);
+
assembleConcreteProcess(t, x, M, K, b, coupling_term);
_boundary_conditions.applyNaturalBC(t, x, K, b);
@@ -134,6 +138,11 @@ void Process::assembleWithJacobian(const double t, GlobalVector const& x,
GlobalVector& b, GlobalMatrix& Jac,
StaggeredCouplingTerm const& coupling_term)
{
+ // The function only has computation if DDC is appied,
+ // e.g. Parallel comuting.
+ MathLib::LinAlg::setLocalAccessibleVector(x);
+ MathLib::LinAlg::setLocalAccessibleVector(xdot);
+
assembleWithJacobianConcreteProcess(t, x, xdot, dxdot_dx, dx_dx, M, K, b,
Jac, coupling_term);
@@ -248,6 +257,10 @@ void Process::computeSecondaryVariable(const double t, GlobalVector const& x,
StaggeredCouplingTerm const&
coupled_term)
{
+ // The function only has computation if DDC is appied,
+ // e.g. Parallel comuting.
+ MathLib::LinAlg::setLocalAccessibleVector(x);
+
computeSecondaryVariableConcrete(t, x, coupled_term);
}
diff --git a/ProcessLib/UncoupledProcessesTimeLoop.cpp b/ProcessLib/UncoupledProcessesTimeLoop.cpp
index f0465c12d4e..f31970e6894 100644
--- a/ProcessLib/UncoupledProcessesTimeLoop.cpp
+++ b/ProcessLib/UncoupledProcessesTimeLoop.cpp
@@ -520,7 +520,6 @@ bool UncoupledProcessesTimeLoop::setCoupledSolutions()
const std::size_t c_id =
std::distance(_per_process_data.begin(), found_item);
- MathLib::LinAlg::setLocalAccessibleVector(*_process_solutions[c_id]);
BaseLib::insertIfTypeIndexKeyUniqueElseError(
coupled_xs, coupled_process_pair.first,
*_process_solutions[c_id], "global_coupled_x");
@@ -533,7 +532,6 @@ bool UncoupledProcessesTimeLoop::setCoupledSolutions()
// Create a vector to store the solution of the last coupling iteration
auto& x_coupling0 = NumLib::GlobalVectorProvider::provider.getVector(x);
MathLib::LinAlg::copy(x, x_coupling0);
- MathLib::LinAlg::setLocalAccessibleVector(x_coupling0);
// append a solution vector of suitable size
_solutions_of_last_cpl_iteration.emplace_back(&x_coupling0);
@@ -589,7 +587,7 @@ bool UncoupledProcessesTimeLoop::loop()
for (auto& spd : _per_process_data)
{
auto& pcs = spd->process;
- auto& x0 = *_process_solutions[pcs_idx];
+ auto const& x0 = *_process_solutions[pcs_idx];
pcs.preTimestep(x0, t0, _timestepper->getTimeStep().dt());
_output->doOutput(pcs, spd->process_output, 0, t0, x0);
@@ -638,7 +636,7 @@ bool UncoupledProcessesTimeLoop::loop()
for (auto& spd : _per_process_data)
{
auto& pcs = spd->process;
- auto& x = *_process_solutions[pcs_idx];
+ auto const& x = *_process_solutions[pcs_idx];
_output->doOutputLastTimestep(pcs, spd->process_output, timestep, t,
x);
@@ -669,11 +667,6 @@ bool UncoupledProcessesTimeLoop::solveUncoupledEquationSystems(
const auto nonlinear_solver_succeeded =
solveOneTimeStepOneProcess(x, timestep_id, t, dt, *spd,
void_staggered_coupling_term, *_output);
-
- // The function only has computation if DDC is appied,
- // e.g. Parallel comuting.
- MathLib::LinAlg::setLocalAccessibleVector(x);
-
pcs.postTimestep(x);
pcs.computeSecondaryVariable(t, x, void_staggered_coupling_term);
@@ -805,11 +798,6 @@ bool UncoupledProcessesTimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
StaggeredCouplingTerm coupled_term(
spd->coupled_processes, _solutions_of_coupled_processes[pcs_idx],
0.0);
-
- // The function only has computation if DDC is appied,
- // e.g. Parallel comuting.
- MathLib::LinAlg::setLocalAccessibleVector(x);
-
pcs.computeSecondaryVariable(t, x, coupled_term);
_output->doOutput(pcs, spd->process_output, timestep_id, t, x);
diff --git a/Tests/NumLib/ODEs.h b/Tests/NumLib/ODEs.h
index 0350116e052..74515c9b33b 100644
--- a/Tests/NumLib/ODEs.h
+++ b/Tests/NumLib/ODEs.h
@@ -113,6 +113,7 @@ public:
) override
{
MathLib::setMatrix(M, {1.0});
+ MathLib::LinAlg::setLocalAccessibleVector(x);
MathLib::setMatrix(K, {x[0]});
MathLib::setVector(b, {0.0});
}
@@ -200,6 +201,7 @@ public:
ProcessLib::StaggeredCouplingTerm const& /*coupling_term*/
) override
{
+ MathLib::LinAlg::setLocalAccessibleVector(x_curr);
auto const u = x_curr[0];
auto const v = x_curr[1];
@@ -217,6 +219,8 @@ public:
ProcessLib::StaggeredCouplingTerm const&
coupling_term) override
{
+ MathLib::LinAlg::setLocalAccessibleVector(x_curr);
+ MathLib::LinAlg::setLocalAccessibleVector(xdot);
assemble(t, x_curr, M, K, b, coupling_term);
auto const u = x_curr[0];
--
GitLab