Commit 41836391 authored by vinciguerra_a's avatar vinciguerra_a
Browse files

Merge branch '6-redesign-solver-class' into 'master'

Resolve "Redesign solver class"

Closes #6

See merge request OPAL/Libraries/ippl-solvers!7
parents 2699e696 a80504a1
......@@ -2,8 +2,9 @@ set (_SRCS
)
set (_HDRS
SolverAlgorithm.h
Electrostatics.h
Solver.h
SolverParams.h
)
include_DIRECTORIES (
......
//
// Class Electrostatics
// Base class for solvers for electrostatics problems
//
// Copyright (c) 2021 Alessandro Vinciguerra, ETH Zürich, Zurich, Switzerland
// All rights reserved
//
// This file is part of IPPL.
//
// IPPL 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.
//
// You should have received a copy of the GNU General Public License
// along with IPPL. If not, see <https://www.gnu.org/licenses/>.
//
#ifndef IPPL_ELECTROSTATICS_H
#define IPPL_ELECTROSTATICS_H
#include "Solver/Solver.h"
namespace ippl {
template <typename Tlhs, typename Trhs, unsigned Dim,
class M=UniformCartesian<double, Dim>,
class C=typename M::DefaultCentering >
class Electrostatics : public Solver<Tlhs, Trhs, Dim, M, C>
{
public:
using grad_type = Field<Vector<Tlhs, Dim>, Dim, M, C>;
using lhs_type = typename Solver<Tlhs, Trhs, Dim, M, C>::lhs_type;
using rhs_type = typename Solver<Tlhs, Trhs, Dim, M, C>::rhs_type;
/*!
* Represents the types of fields that should
* be output by the solver
*/
enum OutputType {
SOL = 0,
GRAD,
SOL_AND_GRAD
};
/*!
* Default constructor for electrostatic solvers;
* desired output type defaults to solution only
*/
Electrostatics()
: Solver<Tlhs, Trhs, Dim, M, C>()
, grad_mp(nullptr)
{
this->params_m.add("output_type", SOL);
}
Electrostatics(lhs_type& lhs, rhs_type& rhs)
: Solver<Tlhs, Trhs, Dim, M, C>(lhs, rhs)
, grad_mp(nullptr)
{
this->params_m.add("output_type", SOL);
}
/*!
* Set the field in which the gradient of the computed potential
* should be stored
* @param grad Reference to field in which to store the gradient
*/
void setGradient(grad_type& grad) { grad_mp = &grad; }
/*!
* Solve the electrostatics problem described by
* -laplace(lhs) = rhs
*/
virtual void solve() = 0;
virtual ~Electrostatics() { }
protected:
grad_type* grad_mp;
};
}
#endif
......@@ -15,49 +15,86 @@
// You should have received a copy of the GNU General Public License
// along with IPPL. If not, see <https://www.gnu.org/licenses/>.
//
#ifndef IPPL_SOLVER_H
#define IPPL_SOLVER_H
#include "SolverParams.h"
#include <functional>
#include "Utility/ParameterList.h"
#include "Field/Field.h"
namespace ippl {
// Expands to a lambda that acts as a wrapper for a differential operator
// fun: the function for which to create the wrapper, such as ippl::laplace
// type: the argument type, which should match the LHS type for the solver
#define IPPL_SOLVER_OPERATOR_WRAPPER(fun, type) \
[] (type arg) { \
return fun(arg); \
}
template <class L, class expr_L, class R>
template <typename Tlhs, typename Trhs, unsigned Dim,
class M=UniformCartesian<double, Dim>,
class C=typename M::DefaultCentering >
class Solver
{
public:
using lhs_type = L;
using rhs_type = R;
using operator_type = std::function<expr_L(L)>;
using lhs_type = Field<Tlhs, Dim, M, C>;
using rhs_type = Field<Trhs, Dim, M, C>;
/*!
* Default constructor
*/
Solver() { }
/*!
* Convenience constructor with LHS and RHS parameters
* @param lhs The LHS for the problem to solve
* @param rhs The RHS for the problem to solve
*/
Solver(lhs_type& lhs, rhs_type& rhs) {
setLhs(lhs);
setRhs(rhs);
}
void setParameters(const SolverParams& params) {
params_m = params;
/*!
* Update one of the solver's parameters
* @param key The parameter key
* @param value The new value
* @throw IpplException Fails if there is no existing parameter with the given key
*/
template <typename T>
void updateParameter(const std::string& key, const T& value) {
params_m.update<T>(key, value);
}
void setLhs(lhs_type* lhs) { lhs_m = lhs; }
void setRhs(rhs_type* rhs) { rhs_m = rhs; }
/*!
* Updates all solver parameters based on values in another parameter set
* @param params Parameter list with updated values
* @throw IpplException Fails if the provided parameter list includes keys not already present
*/
void updateParameters(const ParameterList& params) {
params_m.update(params);
}
void setOperator(operator_type op) {
op_m = std::move(op);
/*!
* Merges another parameter set into the solver's parameters, overwriting
* existing parameters in case of conflict
* @param params Parameter list with desired values
*/
void mergeParameters(const ParameterList& params) {
params_m.merge(params);
}
virtual void solve() = 0;
/*!
* Set the problem LHS
* @param lhs Reference to problem LHS field
*/
void setLhs(lhs_type& lhs) { lhs_mp = &lhs; }
/*!
* Set the problem RHS
* @param rhs Reference to problem RHS field
*/
void setRhs(rhs_type& rhs) { rhs_mp = &rhs; }
protected:
SolverParams params_m;
ParameterList params_m;
rhs_type* rhs_mp;
rhs_type* rhs_m;
lhs_type* lhs_m;
operator_type op_m;
lhs_type* lhs_mp;
};
}
......
//
// Class SolverAlgorithm
// Base class for solver algorithms
//
// Copyright (c) 2021 Alessandro Vinciguerra, ETH Zürich, Zurich, Switzerland
// All rights reserved
//
// This file is part of IPPL.
//
// IPPL 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.
//
// You should have received a copy of the GNU General Public License
// along with IPPL. If not, see <https://www.gnu.org/licenses/>.
//
#ifndef IPPL_SOLVER_ALGORITHM_H
#define IPPL_SOLVER_ALGORITHM_H
#include "Utility/ParameterList.h"
#include <functional>
namespace ippl {
template <typename Tlhs, typename Trhs, unsigned Dim,
class M=UniformCartesian<double, Dim>,
class C=typename M::DefaultCentering >
class SolverAlgorithm
{
public:
using lhs_type = Field<Tlhs, Dim, M, C>;
using rhs_type = Field<Trhs, Dim, M, C>;
/*!
* Solve the problem described by Op(lhs) = rhs, where Op is an unspecified
* differential operator (handled by derived classes)
* @param lhs The problem's LHS
* @param rhs The problem's RHS
* @param params A set of parameters for the solver algorithm
*/
virtual void operator()(lhs_type& lhs, rhs_type& rhs, const ParameterList& params) = 0;
};
}
#endif
//
// Class SolverParams
// Solver related parameters (e.g. tolerance in case of iterative solvers).
// Example:
// ippl::SolverParams params;
// params.add<double>("tolerance", 1.0e-8);
// params.get<double>("tolerance");
//
//
// Copyright (c) 2021, Matthias Frey, University of St Andrews, St Andrews, Scotland
// All rights reserved
//
// This file is part of IPPL.
//
// IPPL 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.
//
// You should have received a copy of the GNU General Public License
// along with IPPL. If not, see <https://www.gnu.org/licenses/>.
//
#ifndef IPPL_SOLVERPARAMS_H
#define IPPL_SOLVERPARAMS_H
#include <iostream>
#include <iomanip>
#include <map>
#include <string>
#include <variant>
#include "Utility/IpplException.h"
namespace ippl {
class SolverParams {
public:
using variant_t = std::variant<double,
float,
bool,
std::string,
unsigned int,
int>;
template <typename T>
void add(const std::string& key, const T& value) {
params_m[key] = value;
}
template <typename T>
T get(const std::string& key) {
// In C++20 this can be replaced with std::map::contains
auto search = params_m.find(key);
if (search == params_m.end()) {
throw IpplException("SolverParams::get()",
"Parameter '" + key + "' not contained.");
}
return std::get<T>(params_m[key]);
}
friend
std::ostream& operator<<(std::ostream& os, const SolverParams& sp) {
for (const auto& [key, value] : sp.params_m) {
std::visit([&](auto&& arg){
os << std::left << std::setw(20) << key
<< " " << arg << '\n';
}, value);
}
return os;
}
private:
std::map<std::string, variant_t> params_m;
};
}
#endif
......@@ -2,15 +2,22 @@
#include <iostream>
#include <typeinfo>
#include <string>
#include "Solver.h"
#include "Electrostatics.h"
template <class L, class expr_L, class R>
class LaplacianSolver : public ippl::Solver<L, expr_L, R> {
constexpr unsigned int dim = 3;
class TestSolver : public ippl::Electrostatics<double, double, dim>
{
public:
void solve() override {
*this->lhs_m = this->op_m(*this->rhs_m);
*rhs_mp = *lhs_mp + *rhs_mp;
if (params_m.get<int>("output_type") > 0) {
*grad_mp = ippl::grad(*lhs_mp);
}
}
};
......@@ -18,8 +25,6 @@ int main(int argc, char *argv[]) {
Ippl ippl(argc,argv);
constexpr unsigned int dim = 3;
int pt = 4;
ippl::Index I(pt);
......@@ -34,54 +39,32 @@ int main(int argc, char *argv[]) {
//Unit box
double dx = 1.0 / double(pt);
ippl::Vector<double, 3> hx = {dx, dx, dx};
ippl::Vector<double, 3> origin = {0, 0, 0};
ippl::UniformCartesian<double, 3> mesh(owned, hx, origin);
double pi = acos(-1.0);
ippl::Vector<double, dim> hx = {dx, dx, dx};
ippl::Vector<double, dim> origin = {0, 0, 0};
ippl::UniformCartesian<double, dim> mesh(owned, hx, origin);
typedef ippl::Field<double, dim> field_type;
field_type field, Lap;
typedef ippl::detail::meta_laplace<field_type> expr_laplace;
field.initialize(mesh, layout);
Lap.initialize(mesh,layout);
typedef ippl::Field<double, dim> Field_t;
typename Field_t::view_type& view = field.getView();
Kokkos::parallel_for("Assign lfield",
Kokkos::MDRangePolicy<Kokkos::Rank<3>>({0, 0, 0},
{view.extent(0),view.extent(1),view.extent(2)}),
KOKKOS_LAMBDA(const int i, const int j, const int k){
double x = (i + 0.5) * hx[0] + origin[0];
double y = (j + 0.5) * hx[1] + origin[1];
double z = (k + 0.5) * hx[2] + origin[2];
//view(i, j, k) = 3.0 * x + 4.0 * y + 5.0 * z;
view(i, j, k) = sin(pi * x) * cos(pi * y) * exp(z);
field_type lhs(mesh, layout), rhs(mesh, layout);
});
typedef ippl::Field<ippl::Vector<double, dim>, dim> vfield_type;
vfield_type grad(mesh, layout);
lhs = 1.0;
rhs = 2.0;
Lap = 0.0;
TestSolver tsolver;
LaplacianSolver<field_type, expr_laplace, field_type> lapsolver;
tsolver.setLhs(lhs);
lapsolver.setRhs(&field);
lapsolver.setLhs(&Lap);
tsolver.setRhs(rhs);
lapsolver.setOperator(IPPL_SOLVER_OPERATOR_WRAPPER(ippl::laplace, field_type));
tsolver.setGradient(grad);
lapsolver.solve();
tsolver.solve();
Lap.write();
rhs.write();
grad.write();
return 0;
}
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