Unverified Commit 9938b2a7 authored by greole's avatar greole Committed by GitHub
Browse files

Unify lduLdu Base (#23)

* implement a unified lduLduBase

* change location of lduMatrix based solvers

* Wrap rhs and initial guesses in vectors

* prepare init_initial_guess_vector

* unify solve_impl_ member functions

* prepare independent preconditioner generation

* add basic implementation for reusable preconditioner

* Add separate preconditioner class

* Add IOPreconditioner files

* remove missing folder

* Fix formating, minor fixes

* fix overflow issue for sorting

* Use csr matrix as default matrix format (#25)

* basic test if csr has any impact

* use csr matrix as matrix type
parent f99a4fc2
/*---------------------------------------------------------------------------*\
License
This file is part of OGL.
OGL is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::GKOCG
Author: Gregor Olenik <go@hpsim.de>
SourceFiles
GKOCG.C
\*---------------------------------------------------------------------------*/
#ifndef GKO_LDUBase_H
#define GKO_LDUBase_H
#include "../common/common.H"
#include "LduMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam {
/*---------------------------------------------------------------------------*\
Class GKOCG Declaration
\*---------------------------------------------------------------------------*/
template <class Type, class LDUType, class SolverFactory>
class GKOLduBaseSolver
: public lduLduBase<LduMatrix<Type, LDUType, LDUType>, SolverFactory> {
public:
GKOLduBaseSolver(const word &fieldName,
const LduMatrix<Type, LDUType, LDUType> &matrix,
const dictionary &solverControls)
: lduLduBase<LduMatrix<Type, LDUType, LDUType>, SolverFactory>(
fieldName, matrix, solverControls)
{}
virtual SolverPerformance<Type> solve_impl(Field<Type> &psi) const
{
return this->solve_impl_(psi);
}
};
} // End namespace Foam
#endif
......@@ -47,25 +47,7 @@ namespace Foam {
template <class SolverFactory>
class GKOlduBaseSolver : public HostMatrix<lduMatrix>,
public SolverFactory,
public StoppingCriterion,
public IOSortingIdxHandler,
public IOGKOMatrixHandler {
public:
// Some shortcuts
using IndexType = label;
using vec = gko::matrix::Dense<scalar>;
using idx_vec = gko::matrix::Dense<label>;
using mtx = gko::matrix::Coo<scalar>;
using val_array = gko::Array<scalar>;
using idx_array = gko::Array<label>;
using cg = gko::solver::Cg<scalar>;
// for now a copy of the OF matrix is stored
// to keep values contiguous in memory
class GKOlduBaseSolver : public lduLduBase<lduMatrix, SolverFactory> {
public:
//- Construct from matrix components and solver controls
GKOlduBaseSolver(const word &fieldName, const lduMatrix &matrix,
......@@ -73,124 +55,15 @@ public:
const FieldField<Field, scalar> &interfaceIntCoeffs,
const lduInterfaceFieldPtrsList &interfaces,
const dictionary &solverControls)
: HostMatrix<lduMatrix>(fieldName, matrix, interfaceBouCoeffs,
interfaceIntCoeffs, interfaces, solverControls),
SolverFactory(solverControls),
StoppingCriterion(solverControls),
IOSortingIdxHandler(
matrix.mesh().thisDb(), nElems(),
solverControls.lookupOrDefault<Switch>("sort", true)),
IOGKOMatrixHandler(matrix.mesh().thisDb(), solverControls, fieldName)
{
// if sys_matrix is not stored updating is neccesary
// initially
bool stored = get_sys_matrix_stored();
if (!stored) {
SIMPLE_TIME(init_host_sparsity_pattern,
init_host_sparsity_pattern();)
SIMPLE_TIME(update_host_mtx, update_host_matrix_data();)
} else {
// if sys_matrix is stored updating is only neccesary
// when requested explictly
if (get_update_sys_matrix()) {
// since sparsity pattern should already be stored
// on device no init_host_sparsity_pattern call is needed
SIMPLE_TIME(exp_update_host_mtx, update_host_matrix_data();)
}
}
// TODO move compute_sorting_idxs in the if clause above to have all
// init calls together
// after updating the host matrix the host matrix needs to be sorted
if (!get_is_sorted()) {
SIMPLE_TIME(compute_sorting_idxs,
compute_sorting_idxs(row_idxs(), col_idxs(), nCells());)
}
// TODO move compute_sorting_idxs in the if clause above to have all
// init calls together
if (!stored && get_sort()) {
SIMPLE_TIME(sort_host_mtx_sparsity_pattern,
sort_host_matrix_sparsity_pattern(get_sorting_idxs());)
SIMPLE_TIME(sort_host_mtx_data,
sort_host_matrix_data(get_sorting_idxs());)
}
init_device_matrix(matrix.mesh().thisDb(), values(), col_idxs(),
row_idxs(), nElems(), nCells(),
!get_update_sys_matrix());
};
: lduLduBase<lduMatrix, SolverFactory>(
fieldName, matrix, interfaceBouCoeffs, interfaceIntCoeffs,
interfaces, solverControls){};
virtual solverPerformance solve_impl(word typeName, scalarField &psi,
const scalarField &source,
const direction cmpt = 0) const
{
std::shared_ptr<gko::Executor> device_exec =
this->get_device_executor();
// --- Setup class containing solver performance data
solverPerformance solverPerf(
lduMatrix::preconditioner::getName(controlDict_) + typeName,
fieldName_);
scalarField pA(nCells());
scalarField wA(nCells());
matrix_.Amul(wA, psi, interfaceBouCoeffs_, interfaces_, cmpt);
scalar norm_factor = this->normFactor(psi, source, wA, pA);
auto source_view = val_array::view(ref_exec(), nCells(),
const_cast<scalar *>(&source[0]));
auto b =
vec::create(ref_exec(), gko::dim<2>(nCells(), 1), source_view, 1);
std::cout << "solve_impl init_initial_guess" << std::endl;
init_initial_guess(psi, matrix().mesh().thisDb(), nCells());
std::cout << "solve_impl init_initial_guess done" << std::endl;
auto x = this->get_initial_guess();
std::cout << "solve_impl get_initial_guess done" << std::endl;
std::vector<std::shared_ptr<const gko::stop::CriterionFactory>>
criterion_vec{};
criterion_vec.push_back(
build_stopping_criterion(device_exec, norm_factor));
// Generate solver
auto solver_gen = this->create_solver(device_exec, criterion_vec);
// Instantiate a ResidualLogger logger.
auto logger = std::make_shared<IterationLogger>(device_exec);
// Add the previously created logger to the solver factory. The
// logger will be automatically propagated to all solvers created
// from this factory.
solver_gen->add_logger(logger);
// auto b_clone = gko::clone(b);
std::shared_ptr<mtx> gkomatrix = get_gkomatrix();
if (get_export())
export_system(fieldName(), gko::lend(gkomatrix), gko::lend(x),
gko::lend(b));
solverPerf.initialResidual() = 0.0;
auto solver = solver_gen->generate(gko::share(gkomatrix));
// Solve system
SIMPLE_TIME(solve, solver->apply(gko::lend(b), gko::lend(x));)
this->copy_result_back(psi, nCells());
// b_clone->copy_from(gko::lend(b));
solverPerf.finalResidual() = 0.0;
solverPerf.nIterations() = logger->get_iters();
return solverPerf;
return this->solve_impl_(typeName, psi, source, cmpt);
};
};
......
......@@ -28,7 +28,7 @@ SourceFiles
#ifndef GKOBiCGStab_H
#define GKOBiCGStab_H
#include "GKOlduBase.H"
#include "../BaseWrapper/lduBase/GKOlduBase.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -38,6 +38,8 @@ namespace Foam {
Class GKOBiCGStab Declaration
\*---------------------------------------------------------------------------*/
class GKOBiCGStabFactory {
using bj = gko::preconditioner::Jacobi<>;
private:
// executor where Ginkgo will perform the computation
......@@ -56,9 +58,11 @@ public:
create_solver(
std::shared_ptr<gko::Executor> exec,
std::vector<std::shared_ptr<const gko::stop::CriterionFactory>>
criterion_vec) const
criterion_vec,
std::shared_ptr<const bj> precond) const
{
if (preconditioner_ == "BJ") return create_BJ(exec, criterion_vec);
if (precond != NULL)
return create_precond(exec, criterion_vec, precond);
return create_default(exec, criterion_vec);
};
......@@ -76,15 +80,15 @@ public:
std::unique_ptr<gko::solver::Bicgstab<double>::Factory,
std::default_delete<gko::solver::Bicgstab<double>::Factory>>
create_BJ(std::shared_ptr<gko::Executor> exec,
std::vector<std::shared_ptr<const gko::stop::CriterionFactory>>
criterion_vec) const
create_precond(
std::shared_ptr<gko::Executor> exec,
std::vector<std::shared_ptr<const gko::stop::CriterionFactory>>
criterion_vec,
std::shared_ptr<const bj> precond) const
{
using bj = gko::preconditioner::Jacobi<>;
return gko::solver::Bicgstab<scalar>::build()
.with_criteria(criterion_vec)
.with_preconditioner(
bj::build().with_max_block_size(blockSize_).on(exec))
.with_generated_preconditioner(precond)
.on(exec);
};
};
......
......@@ -28,13 +28,15 @@ SourceFiles
#ifndef GKOCG_H
#define GKOCG_H
#include "GKOlduBase.H"
#include "../BaseWrapper/lduBase/GKOlduBase.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam {
class GKOCGFactory {
using bj = gko::preconditioner::Jacobi<>;
private:
// executor where Ginkgo will perform the computation
......@@ -53,9 +55,11 @@ public:
create_solver(
std::shared_ptr<gko::Executor> exec,
std::vector<std::shared_ptr<const gko::stop::CriterionFactory>>
criterion_vec) const
criterion_vec,
std::shared_ptr<const bj> precond) const
{
if (preconditioner_ == "BJ") return create_BJ(exec, criterion_vec);
if (precond != NULL)
return create_precond(exec, criterion_vec, precond);
return create_default(exec, criterion_vec);
};
......@@ -73,15 +77,16 @@ public:
std::unique_ptr<gko::solver::Cg<double>::Factory,
std::default_delete<gko::solver::Cg<double>::Factory>>
create_BJ(std::shared_ptr<gko::Executor> exec,
std::vector<std::shared_ptr<const gko::stop::CriterionFactory>>
criterion_vec) const
create_precond(
std::shared_ptr<gko::Executor> exec,
std::vector<std::shared_ptr<const gko::stop::CriterionFactory>>
criterion_vec,
std::shared_ptr<const bj> precond) const
{
using bj = gko::preconditioner::Jacobi<>;
std::cout << " CG create precond " << std::endl;
return gko::solver::Cg<scalar>::build()
.with_criteria(criterion_vec)
.with_preconditioner(
bj::build().with_max_block_size(blockSize_).on(exec))
.with_generated_preconditioner(precond)
.on(exec);
};
};
......
......@@ -97,11 +97,12 @@ IOHandler/IOPtr/IOPtr.C
IOHandler/IOSortingIdxHandler/IOSortingIdxHandler.C
IOHandler/IOExecutorHandler/IOExecutorHandler.C
IOHandler/IOGKOMatrixHandler/IOGKOMatrixHandler.C
lduMatrix/GKOlduBase/GKOlduBase.C
lduMatrix/GKOCG/GKOCG.C
lduMatrix/GKOIR/GKOIR.C
lduMatrix/GKOBiCGStab/GKOBiCGStab.C
LduMatrix/GKOLduBase/GKOLduBase.C
IOHandler/IOPreconditioner/IOPreconditioner.C
BaseWrapper/lduBase/GKOlduBase.C
BaseWrapper/LduBase/GKOLduBase.C
CG/GKOCG.C
BiCGStab/GKOBiCGStab.C
IR/GKOIR.C
LduMatrix/GKOACG/GKOACG.C
PUBLIC
common/common.H
......@@ -109,11 +110,12 @@ IOHandler/IOPtr/IOPtr.H
IOHandler/IOSortingIdxHandler/IOSortingIdxHandler.H
IOHandler/IOExecutorHandler/IOExecutorHandler.H
IOHandler/IOGKOMatrixHandler/IOGKOMatrixHandler.H
lduMatrix/GKOlduBase/GKOlduBase.H
lduMatrix/GKOCG/GKOCG.H
lduMatrix/GKOIR/GKOIR.H
lduMatrix/GKOBiCGStab/GKOBiCGStab.H
LduMatrix/GKOLduBase/GKOLduBase.H
IOHandler/IOPreconditioner/IOPreconditioner.H
BaseWrapper/lduBase/GKOlduBase.H
BaseWrapper/LduBase/GKOLduBase.H
CG/GKOCG.H
IR/GKOIR.H
BiCGStab/GKOBiCGStab.H
LduMatrix/GKOACG/GKOACG.H
)
......@@ -125,12 +127,13 @@ target_include_directories(OGL
$ENV{FOAM_SRC}/OpenFOAM/lnInclude
$ENV{FOAM_SRC}/OSspecific/POSIX/lnInclude
common/
BaseWrapper/lduBase/
BaseWrapper/LduBase/
IOHandler/IOPtr
IOHandler/IOExecutorHandler/
IOHandler/IOSortingIdxHandler/
IOHandler/IOGKOMatrixHandler/
lduMatrix/GKOlduBase
LduMatrix/GKOLduBase
IOHandler/IOPreconditioner/
)
......
......@@ -36,10 +36,6 @@ SourceFiles
namespace Foam {
using vec = gko::matrix::Dense<scalar>;
using mtx = gko::matrix::Coo<scalar>;
class IOExecutorHandler {
private:
const word device_executor_name_;
......
......@@ -37,7 +37,7 @@ void IOGKOMatrixHandler::init_device_matrix(
std::shared_ptr<gko::Executor> device_exec = get_device_executor();
if (sys_matrix_stored_) {
gkomatrix_ptr_ = &db.lookupObjectRef<GKOCOOIOPtr>(sys_matrix_name_);
gkomatrix_ptr_ = &db.lookupObjectRef<GKOCSRIOPtr>(sys_matrix_name_);
return;
}
......@@ -66,16 +66,22 @@ void IOGKOMatrixHandler::init_device_matrix(
}
// if system matrix is not stored create it and set shared pointer
auto coo_mtx = gko::share(
coo_mtx::create(device_exec, gko::dim<2>(nCells, nCells),
val_array::view(gko::ReferenceExecutor::create(),
nElems, &values_host[0]),
*col_idx.get(), *row_idx.get()));
auto gkomatrix =
gko::share(mtx::create(device_exec, gko::dim<2>(nCells, nCells),
val_array::view(gko::ReferenceExecutor::create(),
nElems, &values_host[0]),
*col_idx.get(), *row_idx.get()));
gko::share(mtx::create(device_exec, gko::dim<2>(nCells, nCells)));
coo_mtx->convert_to(gkomatrix.get());
// if updating system matrix is not needed store ptr in obj registry
if (store) {
const fileName path = sys_matrix_name_;
gkomatrix_ptr_ = new GKOCOOIOPtr(IOobject(path, db), gkomatrix);
gkomatrix_ptr_ = new GKOCSRIOPtr(IOobject(path, db), gkomatrix);
} else {
gkomatrix_ptr_ = NULL;
gkomatrix_ = gkomatrix;
......@@ -91,7 +97,7 @@ void IOGKOMatrixHandler::init_device_matrix(
defineTemplateTypeNameWithName(GKOIDXIOPtr, "IDXIOPtr");
defineTemplateTypeNameWithName(GKOCOOIOPtr, "COOIOPtr");
defineTemplateTypeNameWithName(GKOCSRIOPtr, "CSRIOPtr");
defineTemplateTypeNameWithName(GKOVECIOPtr, "VECIOPtr");
} // namespace Foam
......@@ -35,7 +35,8 @@ SourceFiles
namespace Foam {
using mtx = gko::matrix::Coo<scalar>;
using coo_mtx = gko::matrix::Coo<scalar>;
using mtx = gko::matrix::Csr<scalar>;
using vec = gko::matrix::Dense<scalar>;
using val_array = gko::Array<scalar>;
using idx_array = gko::Array<label>;
......@@ -64,7 +65,7 @@ private:
mutable std::shared_ptr<mtx> gkomatrix_ = NULL;
mutable GKOCOOIOPtr *gkomatrix_ptr_ = NULL;
mutable GKOCSRIOPtr *gkomatrix_ptr_ = NULL;
mutable GKOIDXIOPtr *io_col_idxs_ptr_ = NULL;
......@@ -72,7 +73,7 @@ private:
mutable std::shared_ptr<vec> init_guess_ = NULL;
mutable GKOVECIOPtr *io_init_guess_ptr_ = NULL;
mutable std::vector<GKOVECIOPtr *> io_init_guess_ptrs_ = {};
public:
IOGKOMatrixHandler(const objectRegistry &db, const dictionary &controlDict,
......@@ -111,51 +112,81 @@ public:
return gkomatrix_ptr_->get_ptr();
};
void init_initial_guess(const scalarField &psi, const objectRegistry &db,
const label nCells) const
void init_initial_guess_vector(const Field<vector> &psi,
const objectRegistry &db,
const label nCells) const
{
std::vector<word> postFixes{"", "y", "z"};
for (int i = 0; i < 3; i++) {
init_initial_guess(&psi[0][i], db, nCells, postFixes[i]);
}
}
void init_initial_guess_vector(const Field<scalar> &psi,
const objectRegistry &db,
const label nCells) const
{
init_initial_guess(&psi[0], db, nCells, "");
}
void init_initial_guess(const scalar *psi, const objectRegistry &db,
const label nCells, const word postFix) const
{
std::shared_ptr<gko::Executor> device_exec = get_device_executor();
if (init_guess_vector_stored_ && !update_init_guess_vector_) {
std::cout << " init_initial_guess do nothing" << std::endl;
io_init_guess_ptr_ =
&db.lookupObjectRef<GKOVECIOPtr>(init_guess_vector_name_);
io_init_guess_ptrs_.push_back(&db.lookupObjectRef<GKOVECIOPtr>(
init_guess_vector_name_ + postFix));
return;
}
std::cout << " init_initial_guess create psi_view" << std::endl;
auto psi_view = val_array::view(gko::ReferenceExecutor::create(),
nCells, const_cast<scalar *>(&psi[0]));
nCells, const_cast<scalar *>(psi));
std::cout << " init_initial_guess create x" << std::endl;
auto x = gko::share(
vec::create(device_exec, gko::dim<2>(nCells, 1), psi_view, 1));
const fileName path_init_guess = init_guess_vector_name_;
std::cout << " init_initial_guess create io_init_guess_ptr_"
<< std::endl;
io_init_guess_ptr_ = new GKOVECIOPtr(IOobject(path_init_guess, db), x);
const fileName path_init_guess = init_guess_vector_name_ + postFix;
io_init_guess_ptrs_.push_back(
new GKOVECIOPtr(IOobject(path_init_guess, db), x));
}
std::shared_ptr<vec> get_initial_guess() const
void get_initial_guess(std::vector<std::shared_ptr<vec>> &init_guess) const
{
std::cout << "get_initial_guess" << std::endl;
if (io_init_guess_ptr_ == NULL) {
std::cout << "is Null" << std::endl;
return init_guess_;
if (io_init_guess_ptrs_.size() == 0) {
init_guess.push_back(init_guess_);
}
for (int i = 0; i < io_init_guess_ptrs_.size(); i++) {
init_guess.push_back(io_init_guess_ptrs_[i]->get_ptr());
}
}
template <class OFField>
void copy_result_back(const OFField &psi, const label nCells) const
{
std::vector<std::shared_ptr<vec>> device_xs_ptr{};
get_initial_guess(device_xs_ptr);
for (int i = 0; i < 3; i++) {
auto host_x = vec::create(ref_exec(), gko::dim<2>(nCells, 1));
host_x->copy_from(gko::lend(device_xs_ptr[i]));
auto host_x_view =
val_array::view(ref_exec(), nCells, host_x->get_values());
auto psi_view = val_array::view(ref_exec(), nCells,
const_cast<scalar *>(&psi[0][i]));
psi_view = host_x_view;
}
return io_init_guess_ptr_->get_ptr();
}
void copy_result_back(const scalarField &psi, const label nCells) const
{
std::cout << "copy back" << std::endl;
auto device_x = vec::create(ref_exec(), gko::dim<2>(nCells, 1));
std::cout << "copy back 2" << std::endl;
device_x->copy_from(gko::lend(get_initial_guess()));
std::cout << "create x_view back" << std::endl;
std::vector<std::shared_ptr<vec>> device_xs_ptr{};
get_initial_guess(device_xs_ptr);
device_x->copy_from(gko::lend(device_xs_ptr[0]));
auto x_view =
val_array::view(ref_exec(), nCells, device_x->get_values());
// for (label i = 0; i < nCells; i++) {
......@@ -163,12 +194,9 @@ public:
// }
// move frome device
std::cout << "create psi_view " << std::endl;
auto psi_view =
val_array::view(ref_exec(), nCells, const_cast<scalar *>(&psi[0]));
std::cout << "copy x_view to psi_view " << std::endl;
psi_view = x_view;
}
......
/*---------------------------------------------------------------------------*\
License
This file is part of OGL.
OGL is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OGL. If not, see <http://www.gnu.org/licenses/>.
Class