diff --git a/MathLib/LinAlg/LinAlg.cpp b/MathLib/LinAlg/LinAlg.cpp index 1f804d46c6e..927fe14b792 100644 --- a/MathLib/LinAlg/LinAlg.cpp +++ b/MathLib/LinAlg/LinAlg.cpp @@ -22,6 +22,11 @@ namespace MathLib { namespace LinAlg // Vector +void setLocalAccessibleVector(PETScVector const& x) +{ + x.setLocalAccessibleVector(); +} + void set(PETScVector& x, double const a) { VecSet(x.getRawVector(), a); @@ -177,6 +182,10 @@ namespace MathLib { namespace LinAlg // Vector +void setLocalAccessibleVector(EigenVector const& /*x*/) +{ +} + void set(EigenVector& x, double const a) { x.getRawVector().setConstant(a); diff --git a/MathLib/LinAlg/LinAlg.h b/MathLib/LinAlg/LinAlg.h index 06388524bce..8577adf1d4e 100644 --- a/MathLib/LinAlg/LinAlg.h +++ b/MathLib/LinAlg/LinAlg.h @@ -148,6 +148,10 @@ namespace LinAlg // Vector +/// Set local accessible vector in order to get entries. +/// Call this before call operator[] or get(...) of x. +void setLocalAccessibleVector(PETScVector const& x); + void set(PETScVector& x, double const a); void copy(PETScVector const& x, PETScVector& y); @@ -206,6 +210,16 @@ namespace LinAlg // Vector +/** + Set local accessible vector in order to get entries. Call this + before call operator[] or get(...) of x. + The function only has computation if DDC is enabled, + e.g. parallel computing. Up to now, Eigen vector is not used for + global vectors in parallel computing. Therefore this function is + empty in the current status. +*/ +void setLocalAccessibleVector(EigenVector const& x); + void set(EigenVector& x, double const a); void copy(EigenVector const& x, EigenVector& y); diff --git a/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp b/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp index e4fa6ee8052..a563e4e47c9 100644 --- a/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp +++ b/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp @@ -83,7 +83,7 @@ bool PETScLinearSolver::solve(PETScMatrix& A, PETScVector &b, PETScVector &x) KSPSetOperators(_solver, A.getRawMatrix(), A.getRawMatrix(), DIFFERENT_NONZERO_PATTERN); #endif - KSPSolve(_solver, b.getData(), x.getData()); + KSPSolve(_solver, b.getRawVector(), x.getRawVector()); KSPConvergedReason reason; KSPGetConvergedReason(_solver, &reason); diff --git a/MathLib/LinAlg/PETSc/PETScVector.cpp b/MathLib/LinAlg/PETSc/PETScVector.cpp index e5ea7875ce0..12e1a2b10e1 100644 --- a/MathLib/LinAlg/PETSc/PETScVector.cpp +++ b/MathLib/LinAlg/PETSc/PETScVector.cpp @@ -28,10 +28,13 @@ namespace MathLib { PETScVector::PETScVector(const PetscInt vec_size, const bool is_global_size) { - if( is_global_size ) { + if( is_global_size ) + { VecCreate(PETSC_COMM_WORLD, &_v); VecSetSizes(_v, PETSC_DECIDE, vec_size); - } else { + } + else + { // Fix size partitioning // the size can be associated to specific memory allocation of a matrix VecCreateMPI(PETSC_COMM_WORLD, vec_size, PETSC_DECIDE, &_v); @@ -43,8 +46,8 @@ PETScVector::PETScVector(const PetscInt vec_size, const bool is_global_size) PETScVector::PETScVector(const PetscInt vec_size, const std::vector& ghost_ids, const bool is_global_size) - : _size_ghosts{static_cast(ghost_ids.size())} - , _has_ghost_id{true} + : _size_ghosts {static_cast(ghost_ids.size())}, + _has_ghost_id {true} { PetscInt nghosts = static_cast( ghost_ids.size() ); if ( is_global_size ) @@ -61,6 +64,10 @@ PETScVector::PETScVector(const PetscInt vec_size, } config(); + + for (PetscInt i=0; i& u) const { #ifdef TEST_MEM_PETSC @@ -138,10 +146,12 @@ void PETScVector::getGlobalVector(PetscScalar u[]) PetscMemoryGetCurrentUsage(&mem1); #endif + assert(u.size() == _size); + PetscScalar *xp = nullptr; VecGetArray(_v, &xp); - gatherLocalVectors(xp, u); + gatherLocalVectors(xp, u.data()); //This following line may be needed late on // for a communication load balance: @@ -156,15 +166,77 @@ void PETScVector::getGlobalVector(PetscScalar u[]) #endif } -void PETScVector::copyValues(std::vector& u) const +void PETScVector::setLocalAccessibleVector() const +{ + if (_entry_array.empty()) + { + const PetscInt array_size + = _global_ids2local_ids_ghost.size() > 0 ? + _size_loc + _size_ghosts: _size; + _entry_array.resize(array_size); + } + + if ( !_global_ids2local_ids_ghost.empty() ) + { + PetscScalar* loc_x = getLocalVector(); + std::copy_n(loc_x, _size_loc + _size_ghosts, + _entry_array.begin()); + restoreArray(loc_x); + } + else + getGlobalVector(_entry_array); +} + +void PETScVector::copyValues(std::vector& u) const { assert(u.size() == (std::size_t) (getLocalSize() + getGhostSize())); - double* loc_x = getLocalVector(); + PetscScalar* loc_x = getLocalVector(); std::copy_n(loc_x, getLocalSize() + getGhostSize(), u.begin()); restoreArray(loc_x); } +PetscScalar PETScVector::get(const PetscInt idx) const +{ + if (_global_ids2local_ids_ghost.size() > 0) + { + return _entry_array[getLocalIndex(idx)]; + } + + // Ghost entries, and its original index is 0. + const PetscInt id_p = (idx == -_size) ? 0 : std::abs(idx); + return _entry_array[id_p]; +} + + +std::vector PETScVector::get(std::vector const& indices) const +{ + std::vector local_x(indices.size()); + // If VecGetValues can get values from different processors, + // use VecGetValues(_v, indices.size(), indices.data(), + // local_x.data()); + + if (_global_ids2local_ids_ghost.size() > 0) + { + for (std::size_t i=0; i= 0) // non-ghost entry. + return global_index - _start_rank; + + // A special case for a ghost location with global index equal to + // the size of the local vector: + PetscInt real_global_index = + (-global_index == _size) ? 0 : -global_index; + + return _global_ids2local_ids_ghost.at(real_global_index); +} + void PETScVector::restoreArray(PetscScalar* array) const { if (_has_ghost_id) @@ -221,6 +306,7 @@ void PETScVector::shallowCopy(const PETScVector &v) _size_loc = v._size_loc; _size_ghosts = v._size_ghosts; _has_ghost_id = v._has_ghost_id; + _global_ids2local_ids_ghost = v._global_ids2local_ids_ghost; VecSetOption(_v, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE); } diff --git a/MathLib/LinAlg/PETSc/PETScVector.h b/MathLib/LinAlg/PETSc/PETScVector.h index cba4692313e..bb5193109d5 100644 --- a/MathLib/LinAlg/PETSc/PETScVector.h +++ b/MathLib/LinAlg/PETSc/PETScVector.h @@ -16,8 +16,10 @@ #pragma once +#include #include #include + #include namespace MathLib @@ -154,59 +156,41 @@ class PETScVector VecSetValues(_v, e_idxs.size(), &e_idxs[0], &sub_vec[0], INSERT_VALUES); } - //! Get several entries - std::vector get(std::vector const& indices) const - { - std::vector local_x(indices.size()); + // TODO preliminary + void setZero() { VecSet(_v, 0.0); } - VecGetValues(_v, indices.size(), indices.data(), local_x.data()); + /// Disallow moving. + /// \todo This operator should be implemented properly when doing a + /// general cleanup of all matrix and vector classes. + PETScVector& operator = (PETScVector &&) = delete; - return local_x; - } + /// Set local accessible vector in order to get entries. + /// Call this before call operator[] or get(...). + void setLocalAccessibleVector() const; - // TODO preliminary - double operator[] (PetscInt idx) const + /// Get several entries. setLocalAccessibleVector() must be + /// called beforehand. + std::vector get(std::vector const& indices) const; + + /// Get the value of an entry by [] operator. + /// setLocalAccessibleVector() must be called beforehand. + PetscScalar operator[] (PetscInt idx) const { - double value; - VecGetValues(_v, 1, &idx, &value); - return value; + return get(idx); } /*! Get global vector \param u Array to store the global vector. Memory allocation is needed in advance */ - void getGlobalVector(PetscScalar u[]); + void getGlobalVector(std::vector& u) const; - /*! - Copy local entries including ghost ones to an array - \param u Preallocated vector for the values of local entries. + /* Get an entry value. This is an expensive operation, + and it only get local value. Use it for only test purpose + Get the value of an entry by [] operator. + setLocalAccessibleVector() must be called beforehand. */ - void copyValues(std::vector& u) const; - - /// Get an entry value. This is an expensive operation, - /// and it only get local value. Use it for only test purpose - PetscScalar get(const PetscInt idx) const - { - PetscScalar x; - VecGetValues(_v, 1, &idx, &x); - return x; - } - - // TODO eliminate in favour of getRawVector() - /// Get PETsc vector. Use it only for test purpose - const PETSc_Vec &getData() const - { - return _v; - } - - // TODO preliminary - void setZero() { VecSet(_v, 0.0); } - - /// Disallow moving. - /// \todo This operator should be implemented properly when doing a - /// general cleanup of all matrix and vector classes. - PETScVector& operator = (PETScVector &&) = delete; + PetscScalar get(const PetscInt idx) const; //! Exposes the underlying PETSc vector. PETSc_Vec getRawVector() { return _v; } @@ -217,8 +201,13 @@ class PETScVector * This method is dangerous insofar as you can do arbitrary things also * with a const PETSc vector! */ - PETSc_Vec getRawVector() const {return _v; } + PETSc_Vec getRawVector() const { return _v; } + /*! + Copy local entries including ghost ones to an array + \param u Preallocated vector for the values of local entries. + */ + void copyValues(std::vector& u) const; /*! View the global vector for test purpose. Do not use it for output a big vector. \param file_name File name for output @@ -253,7 +242,7 @@ class PETScVector PETSc_Vec _v = nullptr; /// Local vector, which is only for the case that _v is created - /// with ghost entries. + /// with ghost entries. mutable PETSc_Vec _v_loc = nullptr; /// Starting index in a rank @@ -271,19 +260,34 @@ class PETScVector /// Flag to indicate whether the vector is created with ghost entry indices bool _has_ghost_id = false; + /*! + \brief Array containing the entries of the vector. + If the vector is created without given ghost IDs, the array + contains all entries of the global vector, _v. Otherwise it + only contains the entries of the global vector owned by the + current rank. + */ + mutable std::vector _entry_array; + + /// Map global indices of ghost enrties to local indices + mutable std::map _global_ids2local_ids_ghost; + /*! \brief Collect local vectors \param local_array Local array \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 */ PetscScalar* getLocalVector() const; + /// Get local index by a global index + PetscInt getLocalIndex(const PetscInt global_index) const; + /*! Restore array after finish access local array \param array Pointer to the local array fetched by VecGetArray diff --git a/NumLib/DOF/MeshComponentMap.cpp b/NumLib/DOF/MeshComponentMap.cpp index 19a4b01c701..a365bffd7bf 100644 --- a/NumLib/DOF/MeshComponentMap.cpp +++ b/NumLib/DOF/MeshComponentMap.cpp @@ -12,6 +12,7 @@ #include "MeshComponentMap.h" +#include "BaseLib/Error.h" #include "MeshLib/MeshSubsets.h" #ifdef USE_PETSC @@ -70,17 +71,18 @@ MeshComponentMap::MeshComponentMap( for (std::size_t j = 0; j < mesh_subset.getNumberOfNodes(); j++) { GlobalIndexType global_id = 0; - if (order == ComponentOrder::BY_LOCATION) + if (order != ComponentOrder::BY_LOCATION) { - global_id = static_cast( - components.size() * mesh.getGlobalNodeID(j) + comp_id); - } - else - { - // _num_global_dof is used as the global index offset - global_id = static_cast( - _num_global_dof + mesh.getGlobalNodeID(j)); + // Deactivated since this case is not suitable to + // arrange non ghost entries of a partition within + // a rank in the parallel computing. + OGS_FATAL("Global index in the system of equations" + " can only be numbered by the oder type" + " of ComponentOrder::BY_LOCATION"); } + global_id = static_cast( + components.size() * mesh.getGlobalNodeID(j) + + comp_id); const bool is_ghost = mesh.isGhostNode(mesh.getNode(j)->getID()); if (is_ghost) @@ -294,16 +296,19 @@ GlobalIndexType MeshComponentMap::getLocalIndex( // A special case for a ghost location with global index equal to the size // of the local vector: - if (-global_index == static_cast(_num_global_dof)) - return 0; + GlobalIndexType const real_global_index = + (-global_index == static_cast(_num_global_dof)) + ? 0 : -global_index; // TODO Find in ghost indices is O(n^2/2) for n being the length of // _ghosts_indices. Providing an inverted table would be faster. auto const ghost_index_it = std::find(_ghosts_indices.begin(), - _ghosts_indices.end(), -global_index); + _ghosts_indices.end(), + real_global_index); if (ghost_index_it == _ghosts_indices.end()) { - OGS_FATAL("index %d not found in ghost_indices", -global_index); + OGS_FATAL("index %d not found in ghost_indices", + real_global_index); } // Using std::distance on a std::vector is O(1). As long as _ghost_indices diff --git a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp index 9bc318cd9d8..5c5a15784ef 100644 --- a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp +++ b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp @@ -27,13 +27,14 @@ LocalLinearLeastSquaresExtrapolator::LocalLinearLeastSquaresExtrapolator( MathLib::MatrixSpecifications(dof_table.dofSizeWithoutGhosts(), dof_table.dofSizeWithoutGhosts(), &dof_table.getGhostIndices(), - nullptr))) + nullptr))), #ifndef USE_PETSC - , _residuals(dof_table.size()) + _residuals(dof_table.size()), #else - , _residuals(dof_table.size(), false) + _residuals(dof_table.dofSizeWithoutGhosts(), + dof_table.getGhostIndices(), false), #endif - , _local_to_global(dof_table) + _local_to_global(dof_table) { /* Note in case the following assertion fails: * If you copied the extrapolation code, for your processes from diff --git a/ProcessLib/GlobalVectorFromNamedFunction.cpp b/ProcessLib/GlobalVectorFromNamedFunction.cpp index 6affaad855d..82020a384e3 100644 --- a/ProcessLib/GlobalVectorFromNamedFunction.cpp +++ b/ProcessLib/GlobalVectorFromNamedFunction.cpp @@ -8,6 +8,7 @@ */ #include "GlobalVectorFromNamedFunction.h" +#include "MathLib/LinAlg/FinalizeVectorAssembly.h" #include "MathLib/LinAlg/MatrixVectorTraits.h" #include "NumLib/DOF/DOFTableUtil.h" @@ -53,10 +54,10 @@ GlobalVector const& GlobalVectorFromNamedFunction::call( _context.index = node_id; auto const value = _function_caller.call(args); - // TODO Problems with PETSc? (local vs. global index) result->set(node_id, value); } + MathLib::finalizeVectorAssembly(*result); return *result; } diff --git a/ProcessLib/HT/Tests.cmake b/ProcessLib/HT/Tests.cmake index 2618492602d..c84d4ad8eba 100644 --- a/ProcessLib/HT/Tests.cmake +++ b/ProcessLib/HT/Tests.cmake @@ -11,3 +11,24 @@ AddTest( ThermalConvection_const_viscosity_expected.vtu ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000.000000.vtu temperature T VIS ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000.000000.vtu ) + +# MPI tests +AddTest( + NAME Parallel_LARGE_2D_ThermalConvection_constviscosity + PATH Parabolic/HT/ConstViscosity + EXECUTABLE_ARGS square_5500x5500.prj + WRAPPER mpirun + WRAPPER_ARGS -np 4 + TESTER vtkdiff + REQUIREMENTS OGS_USE_MPI + ABSTOL 1e-15 RELTOL 1e-14 + DIFF_DATA + ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_0.vtu ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_0.vtu p p + ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_1.vtu ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_1.vtu p p + ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_2.vtu ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_2.vtu p p + ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_3.vtu ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_3.vtu p p + ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_0.vtu ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_0.vtu T T + ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_1.vtu ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_1.vtu T T + ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_2.vtu ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_2.vtu T T + ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_3.vtu ConstViscosityThermalConvection_pcs_0_ts_149_t_50000000000_000000_3.vtu T T +) diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp index 5bf2b92329e..25e265394bb 100644 --- a/ProcessLib/Process.cpp +++ b/ProcessLib/Process.cpp @@ -122,6 +122,8 @@ void Process::assemble(const double t, GlobalVector const& x, GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b, StaggeredCouplingTerm const& coupling_term) { + MathLib::LinAlg::setLocalAccessibleVector(x); + assembleConcreteProcess(t, x, M, K, b, coupling_term); _boundary_conditions.applyNaturalBC(t, x, K, b); @@ -134,6 +136,9 @@ void Process::assembleWithJacobian(const double t, GlobalVector const& x, GlobalVector& b, GlobalMatrix& Jac, StaggeredCouplingTerm const& coupling_term) { + MathLib::LinAlg::setLocalAccessibleVector(x); + MathLib::LinAlg::setLocalAccessibleVector(xdot); + assembleWithJacobianConcreteProcess(t, x, xdot, dxdot_dx, dx_dx, M, K, b, Jac, coupling_term); @@ -235,12 +240,14 @@ void Process::computeSparsityPattern() void Process::preTimestep(GlobalVector const& x, const double t, const double delta_t) { + MathLib::LinAlg::setLocalAccessibleVector(x); preTimestepConcreteProcess(x, t, delta_t); _boundary_conditions.preTimestep(t); } void Process::postTimestep(GlobalVector const& x) { + MathLib::LinAlg::setLocalAccessibleVector(x); postTimestepConcreteProcess(x); } @@ -248,6 +255,8 @@ void Process::computeSecondaryVariable(const double t, GlobalVector const& x, StaggeredCouplingTerm const& coupled_term) { + MathLib::LinAlg::setLocalAccessibleVector(x); + computeSecondaryVariableConcrete(t, x, coupled_term); } @@ -258,11 +267,13 @@ void Process::preIteration(const unsigned iter, const GlobalVector &x) cached_var->expire(); } + MathLib::LinAlg::setLocalAccessibleVector(x); preIterationConcreteProcess(iter, x); } NumLib::IterationResult Process::postIteration(const GlobalVector &x) { + MathLib::LinAlg::setLocalAccessibleVector(x); return postIterationConcreteProcess(x); } diff --git a/ProcessLib/StaggeredCouplingTerm.cpp b/ProcessLib/StaggeredCouplingTerm.cpp index a6c7b0da65d..7047b87f829 100644 --- a/ProcessLib/StaggeredCouplingTerm.cpp +++ b/ProcessLib/StaggeredCouplingTerm.cpp @@ -11,10 +11,40 @@ */ #include "StaggeredCouplingTerm.h" + +#include "MathLib/LinAlg/LinAlg.h" #include "Process.h" namespace ProcessLib { + +StaggeredCouplingTerm::StaggeredCouplingTerm( + std::unordered_map const& + coupled_processes_, + std::unordered_map const& coupled_xs_, + const double dt_, const bool empty_) + : coupled_processes(coupled_processes_), + coupled_xs(coupled_xs_), + dt(dt_), + empty(empty_) +{ + for (auto const& coupled_x_pair : coupled_xs) + { + auto const& coupled_x = coupled_x_pair.second; + MathLib::LinAlg::setLocalAccessibleVector(coupled_x); + } + + for (auto const& coupled_process_pair : coupled_processes) + { + auto const& coupled_pcs = coupled_process_pair.second; + auto const prevous_time_x = coupled_pcs.getPreviousTimeStepSolution(); + if (prevous_time_x) + { + MathLib::LinAlg::setLocalAccessibleVector(*prevous_time_x); + } + } +} + const StaggeredCouplingTerm createVoidStaggeredCouplingTerm() { std::unordered_map coupled_processes; diff --git a/ProcessLib/StaggeredCouplingTerm.h b/ProcessLib/StaggeredCouplingTerm.h index 8e362536a17..6cbebcc99ab 100644 --- a/ProcessLib/StaggeredCouplingTerm.h +++ b/ProcessLib/StaggeredCouplingTerm.h @@ -36,13 +36,7 @@ struct StaggeredCouplingTerm coupled_processes_, std::unordered_map const& coupled_xs_, - const double dt_, const bool empty_ = false) - : coupled_processes(coupled_processes_), - coupled_xs(coupled_xs_), - dt(dt_), - empty(empty_) - { - } + const double dt_, const bool empty_ = false); /// References to the coupled processes are distinguished by the keys of /// process types. diff --git a/ProcessLib/TES/TESProcess.cpp b/ProcessLib/TES/TESProcess.cpp index 55d9b8566ac..db778b975b8 100644 --- a/ProcessLib/TES/TESProcess.cpp +++ b/ProcessLib/TES/TESProcess.cpp @@ -285,6 +285,8 @@ NumLib::IterationResult TESProcess::postIterationConcreteProcess( std::vector local_x_cache; std::vector local_x_prev_ts_cache; + MathLib::LinAlg::setLocalAccessibleVector(*_x_previous_timestep); + auto check_variable_bounds = [&](std::size_t id, TESLocalAssemblerInterface& loc_asm) { auto const r_c_indices = NumLib::getRowColumnIndices( @@ -340,7 +342,6 @@ TESProcess::computeVapourPartialPressure( auto const x_nV = Adsorption::AdsorptionReaction::getMolarFraction( x_mV, _assembly_params.M_react, _assembly_params.M_inert); - // TODO Problems with PETSc? (local vs. global index) result_cache->set(node_id, p * x_nV); } @@ -378,7 +379,6 @@ TESProcess::computeRelativeHumidity( auto const p_S = Adsorption::AdsorptionReaction::getEquilibriumVapourPressure(T); - // TODO Problems with PETSc? (local vs. global index) result_cache->set(node_id, p * x_nV / p_S); } @@ -419,7 +419,6 @@ TESProcess::computeEquilibriumLoading( : _assembly_params.react_sys->getEquilibriumLoading( p_V, T, _assembly_params.M_react); - // TODO Problems with PETSc? (local vs. global index) result_cache->set(node_id, C_eq); } diff --git a/Tests/Data b/Tests/Data index adf88547b78..92e79f16e56 160000 --- a/Tests/Data +++ b/Tests/Data @@ -1 +1 @@ -Subproject commit adf88547b78bf4a38568c1efbab5e848ca2911e3 +Subproject commit 92e79f16e56e4051d0d3844df9854466c2cea921 diff --git a/Tests/MathLib/TestGlobalVectorInterface.cpp b/Tests/MathLib/TestGlobalVectorInterface.cpp index 360ceed95ed..51582e59f1a 100644 --- a/Tests/MathLib/TestGlobalVectorInterface.cpp +++ b/Tests/MathLib/TestGlobalVectorInterface.cpp @@ -95,6 +95,7 @@ void checkGlobalVectorInterfacePETSc() const int r0 = x.getRangeBegin(); //x.get(0) is expensive, only get local value. Use it for test purpose + x.setLocalAccessibleVector(); ASSERT_EQ(.0, x.get(r0)); set(x, 10.); @@ -103,18 +104,22 @@ void checkGlobalVectorInterfacePETSc() // Value of x is not copied to y const bool deep_copy = false; T_VECTOR y(x, deep_copy); + y.setLocalAccessibleVector(); ASSERT_EQ(0, y.get(r0)); set(y, 10.0); + y.setLocalAccessibleVector(); ASSERT_EQ(10, y.get(r0)); // y += x axpy(y, 1., x); + y.setLocalAccessibleVector(); ASSERT_EQ(20, y.get(r0)); ASSERT_EQ(80., norm2(y)); // y -= x axpy(y, -1., x); + y.setLocalAccessibleVector(); ASSERT_EQ(10, y.get(r0)); ASSERT_EQ(40., norm2(y)); @@ -130,7 +135,7 @@ void checkGlobalVectorInterfacePETSc() ASSERT_NEAR(0.0, normy-norm2(y), 1.e-10); - double x0[16]; + std::vector x0(16); double z[] = { 2.0000000000000000e+01, diff --git a/Tests/MathLib/TestLinearSolver.cpp b/Tests/MathLib/TestLinearSolver.cpp index 358cea04d44..f1f1e92afa2 100644 --- a/Tests/MathLib/TestLinearSolver.cpp +++ b/Tests/MathLib/TestLinearSolver.cpp @@ -166,8 +166,7 @@ void checkLinearSolverInterface(T_MATRIX& A, T_VECTOR& b, local_vec[1] = 2. * (mrank+1); x.set(row_pos, local_vec); - double x0[6]; - double x1[6]; + std::vector x0(6); x.getGlobalVector(x0); MathLib::LinAlg::matMult(A, x, b); @@ -194,6 +193,7 @@ void checkLinearSolverInterface(T_MATRIX& A, T_VECTOR& b, EXPECT_GT(ls.getNumberOfIterations(), 0u); + std::vector x1(6); x.getGlobalVector(x1); ASSERT_ARRAY_NEAR(x0, x1, 6, 1e-5); } 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 @@ class ODE2 final : public NumLib::ODESystem< ) override { MathLib::setMatrix(M, {1.0}); + MathLib::LinAlg::setLocalAccessibleVector(x); MathLib::setMatrix(K, {x[0]}); MathLib::setVector(b, {0.0}); } @@ -200,6 +201,7 @@ class ODE3 final : public NumLib::ODESystem< 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 @@ class ODE3 final : public NumLib::ODESystem< 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]; diff --git a/Tests/NumLib/TestODEInt.cpp b/Tests/NumLib/TestODEInt.cpp index 9e8ee9b0ce6..229f5fe7be7 100644 --- a/Tests/NumLib/TestODEInt.cpp +++ b/Tests/NumLib/TestODEInt.cpp @@ -19,6 +19,7 @@ #include "BaseLib/BuildInfo.h" #include "NumLib/NumericsConfig.h" #include "NumLib/ODESolver/ConvergenceCriterionDeltaX.h" +#include "MathLib/LinAlg/LinAlg.h" #include "Tests/TestTools.h" #include "ODEs.h" @@ -103,6 +104,9 @@ class TestOutput if (num_timesteps > 0) EXPECT_TRUE(loop.loop(t0, x0, t_end, delta_t, cb)); + for (auto& x : sol.solutions) + MathLib::LinAlg::setLocalAccessibleVector(x); + return sol; } @@ -258,7 +262,8 @@ class NumLibODEIntTyped : public ::testing::Test TestParams::tol_picard_newton); auto const t = sol_picard.ts[i]; - auto const sol_analyt = ODETraits::solution(t); + auto sol_analyt = ODETraits::solution(t); + MathLib::LinAlg::setLocalAccessibleVector(sol_analyt); EXPECT_NEAR(sol_picard.solutions[i][comp], sol_analyt[comp], TestParams::tol_analyt); diff --git a/Tests/NumLib/TimeLoopSingleODE.h b/Tests/NumLib/TimeLoopSingleODE.h index 531ef5432da..1458f420a01 100644 --- a/Tests/NumLib/TimeLoopSingleODE.h +++ b/Tests/NumLib/TimeLoopSingleODE.h @@ -12,6 +12,7 @@ #include "NumLib/DOF/GlobalMatrixProviders.h" #include "NumLib/ODESolver/TimeDiscretizedODESystem.h" #include "NumLib/ODESolver/NonlinearSolver.h" +#include "MathLib/LinAlg/LinAlg.h" #include "ProcessLib/StaggeredCouplingTerm.h" @@ -97,6 +98,7 @@ bool TimeLoopSingleODE::loop(const double t0, GlobalVector const& x0, if (time_disc.needsPreload()) { + MathLib::LinAlg::setLocalAccessibleVector(x); _nonlinear_solver->assemble(x, ProcessLib::createVoidStaggeredCouplingTerm()); time_disc.pushState(t0, x0,