diff --git a/language_bindings/python/bind_py_fftengine.cc b/language_bindings/python/bind_py_fftengine.cc
index 2cfb82a..fa29e07 100644
--- a/language_bindings/python/bind_py_fftengine.cc
+++ b/language_bindings/python/bind_py_fftengine.cc
@@ -1,106 +1,105 @@
/**
* @file bind_py_fftengine.cc
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 17 Jan 2018
*
* @brief Python bindings for the FFT engines
*
* Copyright © 2018 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include "fft/fftw_engine.hh"
#ifdef WITH_FFTWMPI
#include "fft/fftwmpi_engine.hh"
#endif
#ifdef WITH_PFFT
#include "fft/pfft_engine.hh"
#endif
#include "bind_py_declarations.hh"
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/eigen.h>
using namespace muSpectre;
namespace py=pybind11;
using namespace pybind11::literals;
template <Dim_t dim, class Engine>
void add_engine_helper(py::module & mod, std::string name) {
using Ccoord = Ccoord_t<dim>;
- using Rcoord = Rcoord_t<dim>;
using ArrayXXc = Eigen::Array<Complex, Eigen::Dynamic,
Eigen::Dynamic>;
py::class_<Engine>(mod, name.c_str())
#ifdef WITH_MPI
- .def(py::init([](Ccoord res, Rcoord lengths, size_t comm) {
- return new Engine(res, lengths,
- std::move(Communicator(MPI_Comm(comm))));
+ .def(py::init([](Ccoord res, size_t comm) {
+ return new Engine(res, std::move(Communicator(MPI_Comm(comm))));
}),
"resolutions"_a,
- "lengths"_a,
"communicator"_a=size_t(MPI_COMM_SELF))
#else
- .def(py::init<Ccoord, Rcoord>())
+ .def(py::init<Ccoord>())
#endif
.def("fft",
[](Engine & eng, py::EigenDRef<Eigen::ArrayXXd> v) {
using Coll_t = typename Engine::GFieldCollection_t;
using Field_t = typename Engine::Field_t;
Coll_t coll{};
- coll.initialise(eng.get_resolutions(), eng.get_locations());
+ coll.initialise(eng.get_subdomain_resolutions(),
+ eng.get_subdomain_locations());
Field_t & temp{make_field<Field_t>("temp_field", coll)};
temp.eigen() = v;
return ArrayXXc{eng.fft(temp).eigen()};
},
"array"_a)
.def("ifft",
[](Engine & eng,
py::EigenDRef<ArrayXXc> v) {
using Coll_t = typename Engine::GFieldCollection_t;
using Field_t = typename Engine::Field_t;
Coll_t coll{};
- coll.initialise(eng.get_resolutions(), eng.get_locations());
+ coll.initialise(eng.get_subdomain_resolutions(),
+ eng.get_subdomain_locations());
Field_t & temp{make_field<Field_t>("temp_field", coll)};
eng.get_work_space().eigen() = v;
eng.ifft(temp);
return Eigen::ArrayXXd{temp.eigen()};
},
"array"_a)
.def("initialise", &Engine::initialise,
"flags"_a=FFT_PlanFlags::estimate)
.def("normalisation", &Engine::normalisation);
}
void add_fft_engines(py::module & mod) {
auto fft{mod.def_submodule("fft")};
fft.doc() = "bindings for µSpectre's fft engines";
add_engine_helper< twoD, FFTWEngine< twoD, twoD>>(fft, "FFTW_2d");
add_engine_helper<threeD, FFTWEngine<threeD, threeD>>(fft, "FFTW_3d");
#ifdef WITH_FFTWMPI
add_engine_helper< twoD, FFTWMPIEngine< twoD, twoD>>(fft, "FFTWMPI_2d");
add_engine_helper<threeD, FFTWMPIEngine<threeD, threeD>>(fft, "FFTWMPI_3d");
#endif
#ifdef WITH_PFFT
add_engine_helper< twoD, PFFTEngine< twoD, twoD>>(fft, "PFFT_2d");
add_engine_helper<threeD, PFFTEngine<threeD, threeD>>(fft, "PFFT_3d");
#endif
add_projections(fft);
}
diff --git a/language_bindings/python/bind_py_projections.cc b/language_bindings/python/bind_py_projections.cc
index e1c4d38..1dc0ab8 100644
--- a/language_bindings/python/bind_py_projections.cc
+++ b/language_bindings/python/bind_py_projections.cc
@@ -1,139 +1,140 @@
/**
* @file bind_py_projections.cc
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 18 Jan 2018
*
* @brief Python bindings for the Projection operators
*
* Copyright © 2018 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include "fft/projection_small_strain.hh"
#include "fft/projection_finite_strain.hh"
#include "fft/projection_finite_strain_fast.hh"
#include "fft/fftw_engine.hh"
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/eigen.h>
#include <sstream>
#include <memory>
using namespace muSpectre;
namespace py=pybind11;
using namespace pybind11::literals;
/**
* "Trampoline" class for handling the pure virtual methods, see
* [http://pybind11.readthedocs.io/en/stable/advanced/classes.html#overriding-virtual-functions-in-python]
* for details
*/
template <Dim_t DimS, Dim_t DimM>
class PyProjectionBase: public ProjectionBase<DimS, DimM> {
public:
//! base class
using Parent = ProjectionBase<DimS, DimM>;
//! field type on which projection is applied
using Field_t = typename Parent::Field_t;
void apply_projection(Field_t & field) override {
PYBIND11_OVERLOAD_PURE
(void,
Parent,
apply_projection,
field
);
}
Eigen::Map<Eigen::ArrayXXd> get_operator() override {
PYBIND11_OVERLOAD_PURE
(Eigen::Map<Eigen::ArrayXXd>,
Parent,
get_operator
);
}
};
template <class Proj, Dim_t DimS, Dim_t DimM=DimS>
void add_proj_helper(py::module & mod, std::string name_start) {
using Ccoord = Ccoord_t<DimS>;
using Rcoord = Rcoord_t<DimS>;
using Engine = FFTWEngine<DimS, DimM>;
using Field_t = typename Proj::Field_t;
static_assert(DimS == DimM,
"currently only for DimS==DimM");
std::stringstream name{};
name << name_start << '_' << DimS << 'd';
py::class_<Proj>(mod, name.str().c_str())
.def(py::init([](Ccoord res, Rcoord lengths) {
- auto engine = std::make_unique<Engine>(res, lengths);
- return Proj(std::move(engine));
+ auto engine = std::make_unique<Engine>(res);
+ return Proj(std::move(engine), lengths);
}))
.def("initialise", &Proj::initialise,
"flags"_a=FFT_PlanFlags::estimate,
"initialises the fft engine (plan the transform)")
.def("apply_projection",
[](Proj & proj, py::EigenDRef<Eigen::ArrayXXd> v){
typename Engine::GFieldCollection_t coll{};
- coll.initialise(proj.get_resolutions(), proj.get_locations());
+ coll.initialise(proj.get_subdomain_resolutions(),
+ proj.get_subdomain_locations());
Field_t & temp{make_field<Field_t>("temp_field", coll)};
temp.eigen() = v;
proj.apply_projection(temp);
return Eigen::ArrayXXd{temp.eigen()};
})
.def("get_operator", &Proj::get_operator)
.def("get_formulation", &Proj::get_formulation,
"return a Formulation enum indicating whether the projection is small"
" or finite strain");
}
void add_proj_dispatcher(py::module & mod) {
add_proj_helper<
ProjectionSmallStrain< twoD, twoD>,
twoD>(mod, "ProjectionSmallStrain");
add_proj_helper<
ProjectionSmallStrain<threeD, threeD>,
threeD>(mod, "ProjectionSmallStrain");
add_proj_helper<
ProjectionFiniteStrain< twoD, twoD>,
twoD>(mod, "ProjectionFiniteStrain");
add_proj_helper<
ProjectionFiniteStrain<threeD, threeD>,
threeD>(mod, "ProjectionFiniteStrain");
add_proj_helper<
ProjectionFiniteStrainFast< twoD, twoD>,
twoD>(mod, "ProjectionFiniteStrainFast");
add_proj_helper<
ProjectionFiniteStrainFast<threeD, threeD>,
threeD>(mod, "ProjectionFiniteStrainFast");
}
void add_projections(py::module & mod) {
add_proj_dispatcher(mod);
}
diff --git a/language_bindings/python/muSpectre/fft.py b/language_bindings/python/muSpectre/fft.py
index 3fdc5f1..698c399 100644
--- a/language_bindings/python/muSpectre/fft.py
+++ b/language_bindings/python/muSpectre/fft.py
@@ -1,95 +1,97 @@
#
# @file fft.py
#
# @author Lars Pastewka <lars.pastewka@imtek.uni-freiburg.de>
#
# @date 27 Mar 2018
#
# @brief Wrapper for muSpectre's FFT engines
#
# Copyright © 2018 Till Junge
#
# µSpectre 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, or (at
# your option) any later version.
#
# µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
# Free Software Foundation, Inc., 59 Temple Place - Suite 330,
# Boston, MA 02111-1307, USA.
#
try:
from mpi4py import MPI
except ImportError:
MPI = None
import _muSpectre
_factories = {'fftw': ('FFTW_2d', 'FFTW_3d', False),
'fftwmpi': ('FFTWMPI_2d', 'FFTWMPI_3d', True),
'pfft': ('PFFT_2d', 'PFFT_3d', True),
'p3dfft': ('P3DFFT_2d', 'P3DFFT_3d', True)}
-def FFT(resolutions, lengths, fft='fftw', communicator=None):
+def FFT(resolutions, fft='fftw', communicator=None):
"""
Instantiate a muSpectre FFT class.
Parameters
----------
resolutions: list
Grid resolutions in the Cartesian directions.
- lengths: list
- Physical size of the cell in the Cartesian directions.
fft: string
FFT engine to use. Options are 'fftw', 'fftwmpi', 'pfft' and 'p3dfft'.
Default is 'fftw'.
communicator: mpi4py communicator
mpi4py communicator object passed to parallel FFT engines. Note that
the default 'fftw' engine does not support parallel execution.
Returns
-------
cell: object
Return a muSpectre Cell object.
"""
- if len(resolutions) != len(lengths):
- raise ValueError("'resolutions' and 'lengths' must have identical "
- "lengths.")
try:
factory_name_2d, factory_name_3d, is_parallel = _factories[fft]
except KeyError:
raise KeyError("Unknown FFT engine '{}'.".format(fft))
if len(resolutions) == 2:
factory_name = factory_name_2d
elif len(resolutions) == 3:
factory_name = factory_name_3d
else:
raise ValueError('{}-d transforms are not supported'
.format(len(resolutions)))
try:
factory = _muSpectre.fft.__dict__[factory_name]
except KeyError:
raise KeyError("FFT engine '{}' has not been compiled into the "
"muSpectre library.".format(fft))
if is_parallel:
if MPI is None:
raise RuntimeError('Parallel solver requested but mpi4py could'
' not be imported.')
if communicator is None:
+#<<<<<<< HEAD
+# communicator = MPI.COMM_SELF
+# return factory(resolutions, lengths,
+# MPI._handleof(communicator))
+#=======
+# communicator = mpi4py.MPI.COMM_SELF
+# return factory(resolutions, mpi4py.MPI._handleof(communicator))
+#>>>>>>> 6acade5140d1e3e1d35ec7b90b7edb79976d0a08
communicator = MPI.COMM_SELF
- return factory(resolutions, lengths,
- MPI._handleof(communicator))
+ return factory(resolutions, MPI._handleof(communicator))
else:
if communicator is not None:
raise ValueError("FFT engine '{}' does not support parallel "
"execution.".format(fft))
- return factory(resolutions, lengths)
+ return factory(resolutions)
diff --git a/src/cell/cell_base.cc b/src/cell/cell_base.cc
index f8b5b5a..502ae4c 100644
--- a/src/cell/cell_base.cc
+++ b/src/cell/cell_base.cc
@@ -1,292 +1,295 @@
/**
* @file cell_base.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 01 Nov 2017
*
* @brief Implementation for cell base class
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include "cell/cell_base.hh"
#include "common/ccoord_operations.hh"
#include "common/iterators.hh"
#include "common/tensor_algebra.hh"
#include <sstream>
#include <algorithm>
namespace muSpectre {
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
CellBase<DimS, DimM>::CellBase(Projection_ptr projection_)
- :resolutions{projection_->get_resolutions()},
- locations{projection_->get_locations()},
+ :subdomain_resolutions{projection_->get_subdomain_resolutions()},
+ subdomain_locations{projection_->get_subdomain_locations()},
domain_resolutions{projection_->get_domain_resolutions()},
- pixels(resolutions, locations),
- lengths{projection_->get_lengths()},
+ pixels(subdomain_resolutions, subdomain_locations),
+ domain_lengths{projection_->get_domain_lengths()},
fields{std::make_unique<FieldCollection_t>()},
F{make_field<StrainField_t>("Gradient", *this->fields)},
P{make_field<StressField_t>("Piola-Kirchhoff-1", *this->fields)},
projection{std::move(projection_)},
form{projection->get_formulation()}
{ }
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
typename CellBase<DimS, DimM>::Material_t &
CellBase<DimS, DimM>::add_material(Material_ptr mat) {
this->materials.push_back(std::move(mat));
return *this->materials.back();
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
typename CellBase<DimS, DimM>::FullResponse_t
CellBase<DimS, DimM>::evaluate_stress_tangent(StrainField_t & grad) {
if (this->initialised == false) {
this->initialise();
}
//! High level compatibility checks
if (grad.size() != this->F.size()) {
throw std::runtime_error("Size mismatch");
}
constexpr bool create_tangent{true};
this->get_tangent(create_tangent);
for (auto & mat: this->materials) {
mat->compute_stresses_tangent(grad, this->P, this->K.value(),
this->form);
}
return std::tie(this->P, this->K.value());
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
typename CellBase<DimS, DimM>::StressField_t &
CellBase<DimS, DimM>::directional_stiffness(const TangentField_t &K,
const StrainField_t &delF,
StressField_t &delP) {
for (auto && tup:
akantu::zip(K.get_map(), delF.get_map(), delP.get_map())){
auto & k = std::get<0>(tup);
auto & df = std::get<1>(tup);
auto & dp = std::get<2>(tup);
dp = Matrices::tensmult(k, df);
}
return this->project(delP);
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
typename CellBase<DimS, DimM>::SolvVectorOut
CellBase<DimS, DimM>::directional_stiffness_vec(const SolvVectorIn &delF) {
if (!this->K) {
throw std::runtime_error
("corrently only implemented for cases where a stiffness matrix "
"exists");
}
if (delF.size() != this->nb_dof()) {
std::stringstream err{};
- err << "input should be of size ndof = ¶(" << this->resolutions <<") × "
- << DimS << "² = "<< this->nb_dof() << " but I got " << delF.size();
+ err << "input should be of size ndof = ¶(" << this->subdomain_resolutions
+ << ") × " << DimS << "² = "<< this->nb_dof() << " but I got "
+ << delF.size();
throw std::runtime_error(err.str());
}
const std::string out_name{"temp output for directional stiffness"};
const std::string in_name{"temp input for directional stiffness"};
auto & out_tempref = this->get_managed_field(out_name);
auto & in_tempref = this->get_managed_field(in_name);
SolvVectorOut(in_tempref.data(), this->nb_dof()) = delF;
this->directional_stiffness(this->K.value(), in_tempref, out_tempref);
return SolvVectorOut(out_tempref.data(), this->nb_dof());
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
Eigen::ArrayXXd
CellBase<DimS, DimM>::
directional_stiffness_with_copy
(Eigen::Ref<Eigen::ArrayXXd> delF) {
if (!this->K) {
throw std::runtime_error
("corrently only implemented for cases where a stiffness matrix "
"exists");
}
const std::string out_name{"temp output for directional stiffness"};
const std::string in_name{"temp input for directional stiffness"};
auto & out_tempref = this->get_managed_field(out_name);
auto & in_tempref = this->get_managed_field(in_name);
in_tempref.eigen() = delF;
this->directional_stiffness(this->K.value(), in_tempref, out_tempref);
return out_tempref.eigen();
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
typename CellBase<DimS, DimM>::StressField_t &
CellBase<DimS, DimM>::project(StressField_t &field) {
this->projection->apply_projection(field);
return field;
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
typename CellBase<DimS, DimM>::StrainField_t &
CellBase<DimS, DimM>::get_strain() {
if (this->initialised == false) {
this->initialise();
}
return this->F;
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
const typename CellBase<DimS, DimM>::StressField_t &
CellBase<DimS, DimM>::get_stress() const {
return this->P;
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
const typename CellBase<DimS, DimM>::TangentField_t &
CellBase<DimS, DimM>::get_tangent(bool create) {
if (!this->K) {
if (create) {
this->K = make_field<TangentField_t>("Tangent Stiffness", *this->fields);
} else {
throw std::runtime_error
("K does not exist");
}
}
return this->K.value();
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
typename CellBase<DimS, DimM>::StrainField_t &
CellBase<DimS, DimM>::get_managed_field(std::string unique_name) {
if (!this->fields->check_field_exists(unique_name)) {
return make_field<StressField_t>(unique_name, *this->fields);
} else {
return static_cast<StressField_t&>(this->fields->at(unique_name));
}
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void CellBase<DimS, DimM>::initialise(FFT_PlanFlags flags) {
// check that all pixels have been assigned exactly one material
this->check_material_coverage();
// resize all global fields (strain, stress, etc)
- this->fields->initialise(this->resolutions, this->locations);
+ this->fields->initialise(this->subdomain_resolutions, this->subdomain_locations);
// initialise the projection and compute the fft plan
this->projection->initialise(flags);
this->initialised = true;
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void CellBase<DimS, DimM>::initialise_materials(bool stiffness) {
for (auto && mat: this->materials) {
mat->initialise(stiffness);
}
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void CellBase<DimS, DimM>::save_history_variables() {
for (auto && mat: this->materials) {
mat->save_history_variables();
}
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
typename CellBase<DimS, DimM>::iterator
CellBase<DimS, DimM>::begin() {
return this->pixels.begin();
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
typename CellBase<DimS, DimM>::iterator
CellBase<DimS, DimM>::end() {
return this->pixels.end();
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
CellAdaptor<CellBase<DimS, DimM>>
CellBase<DimS, DimM>::get_adaptor() {
return Adaptor(*this);
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void CellBase<DimS, DimM>::check_material_coverage() {
- auto nb_pixels = CcoordOps::get_size(this->resolutions);
+ auto nb_pixels = CcoordOps::get_size(this->subdomain_resolutions);
std::vector<MaterialBase<DimS, DimM>*> assignments(nb_pixels, nullptr);
for (auto & mat: this->materials) {
for (auto & pixel: *mat) {
- auto index = CcoordOps::get_index(this->resolutions, this->locations,
+ auto index = CcoordOps::get_index(this->subdomain_resolutions,
+ this->subdomain_locations,
pixel);
auto& assignment{assignments.at(index)};
if (assignment != nullptr) {
std::stringstream err{};
err << "Pixel " << pixel << "is already assigned to material '"
<< assignment->get_name()
<< "' and cannot be reassigned to material '" << mat->get_name();
throw std::runtime_error(err.str());
} else {
assignments[index] = mat.get();
}
}
}
// find and identify unassigned pixels
std::vector<Ccoord> unassigned_pixels;
for (size_t i = 0; i < assignments.size(); ++i) {
if (assignments[i] == nullptr) {
unassigned_pixels.push_back(
- CcoordOps::get_ccoord(this->resolutions, this->locations, i));
+ CcoordOps::get_ccoord(this->subdomain_resolutions,
+ this->subdomain_locations, i));
}
}
if (unassigned_pixels.size() != 0) {
std::stringstream err {};
err << "The following pixels have were not assigned a material: ";
for (auto & pixel: unassigned_pixels) {
err << pixel << ", ";
}
err << "and that cannot be handled";
throw std::runtime_error(err.str());
}
}
template class CellBase<twoD, twoD>;
template class CellBase<threeD, threeD>;
} // muSpectre
diff --git a/src/cell/cell_base.hh b/src/cell/cell_base.hh
index aa906cd..018ae27 100644
--- a/src/cell/cell_base.hh
+++ b/src/cell/cell_base.hh
@@ -1,311 +1,313 @@
/**
* @file cell_base.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 01 Nov 2017
*
* @brief Base class representing a unit cell cell with single
* projection operator
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef CELL_BASE_H
#define CELL_BASE_H
#include "common/common.hh"
#include "common/ccoord_operations.hh"
#include "common/field.hh"
#include "common/utilities.hh"
#include "materials/material_base.hh"
#include "fft/projection_base.hh"
#include "cell/cell_traits.hh"
#include <vector>
#include <memory>
#include <tuple>
#include <functional>
namespace muSpectre {
template <class Cell>
class CellAdaptor;
//! DimS spatial dimension (dimension of problem
//! DimM material_dimension (dimension of constitutive law)
template <Dim_t DimS, Dim_t DimM=DimS>
class CellBase
{
public:
using Ccoord = Ccoord_t<DimS>; //!< cell coordinates type
using Rcoord = Rcoord_t<DimS>; //!< physical coordinates type
//! global field collection
using FieldCollection_t = GlobalFieldCollection<DimS>;
//! the collection is handled in a `std::unique_ptr`
using Collection_ptr = std::unique_ptr<FieldCollection_t>;
//! polymorphic base material type
using Material_t = MaterialBase<DimS, DimM>;
//! materials handled through `std::unique_ptr`s
using Material_ptr = std::unique_ptr<Material_t>;
//! polymorphic base projection type
using Projection_t = ProjectionBase<DimS, DimM>;
//! projections handled through `std::unique_ptr`s
using Projection_ptr = std::unique_ptr<Projection_t>;
//! expected type for strain fields
using StrainField_t =
TensorField<FieldCollection_t, Real, secondOrder, DimM>;
//! expected type for stress fields
using StressField_t =
TensorField<FieldCollection_t, Real, secondOrder, DimM>;
//! expected type for tangent stiffness fields
using TangentField_t =
TensorField<FieldCollection_t, Real, fourthOrder, DimM>;
//! combined stress and tangent field
using FullResponse_t =
std::tuple<const StressField_t&, const TangentField_t&>;
//! iterator type over all cell pixel's
using iterator = typename CcoordOps::Pixels<DimS>::iterator;
//! input vector for solvers
using SolvVectorIn = Eigen::Ref<const Eigen::VectorXd>;
//! output vector for solvers
using SolvVectorOut = Eigen::Map<Eigen::VectorXd>;
//! sparse matrix emulation
using Adaptor = CellAdaptor<CellBase>;
//! Default constructor
CellBase() = delete;
//! constructor using sizes and resolution
CellBase(Projection_ptr projection);
//! Copy constructor
CellBase(const CellBase &other) = delete;
//! Move constructor
CellBase(CellBase &&other) = default;
//! Destructor
virtual ~CellBase() = default;
//! Copy assignment operator
CellBase& operator=(const CellBase &other) = delete;
//! Move assignment operator
CellBase& operator=(CellBase &&other) = default;
/**
* Materials can only be moved. This is to assure exclusive
* ownership of any material by this cell
*/
Material_t & add_material(Material_ptr mat);
/**
* evaluates all materials
*/
FullResponse_t evaluate_stress_tangent(StrainField_t & F);
/**
* evaluate directional stiffness (i.e. G:K:δF or G:K:δε)
*/
StressField_t & directional_stiffness(const TangentField_t & K,
const StrainField_t & delF,
StressField_t & delP);
/**
* vectorized version for eigen solvers, no copy, but only works
* when fields have ArrayStore=false
*/
SolvVectorOut directional_stiffness_vec(const SolvVectorIn & delF);
/**
* Evaluate directional stiffness into a temporary array and
* return a copy. This is a costly and wasteful interface to
* directional_stiffness and should only be used for debugging or
* in the python interface
*/
Eigen::ArrayXXd directional_stiffness_with_copy
(Eigen::Ref<Eigen::ArrayXXd> delF);
/**
* Convenience function circumventing the neeed to use the
* underlying projection
*/
StressField_t & project(StressField_t & field);
//! returns a ref to the cell's strain field
StrainField_t & get_strain();
//! returns a ref to the cell's stress field
const StressField_t & get_stress() const;
//! returns a ref to the cell's tangent stiffness field
const TangentField_t & get_tangent(bool create = false);
//! returns a ref to a temporary field managed by the cell
StrainField_t & get_managed_field(std::string unique_name);
/**
* general initialisation; initialises the projection and
* fft_engine (i.e. infrastructure) but not the materials. These
* need to be initialised separately
*/
void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate);
/**
* initialise materials (including resetting any history variables)
*/
void initialise_materials(bool stiffness=false);
/**
* for materials with state variables, these typically need to be
* saved/updated an the end of each load increment, this function
* calls this update for each material in the cell
*/
void save_history_variables();
iterator begin(); //!< iterator to the first pixel
iterator end(); //!< iterator past the last pixel
//! number of pixels in the cell
size_t size() const {return pixels.size();}
//! return the subdomain resolutions of the cell
- const Ccoord & get_resolutions() const {return this->resolutions;}
+ const Ccoord & get_subdomain_resolutions() const {
+ return this->subdomain_resolutions;}
//! return the subdomain locations of the cell
- const Ccoord & get_locations() const {return this->locations;}
+ const Ccoord & get_subdomain_locations() const {
+ return this->subdomain_locations;}
//! return the domain resolutions of the cell
const Ccoord & get_domain_resolutions() const {
return this->domain_resolutions;}
//! return the sizes of the cell
- const Rcoord & get_lengths() const {return this->lengths;}
+ const Rcoord & get_domain_lengths() const {return this->domain_lengths;}
/**
* formulation is hard set by the choice of the projection class
*/
const Formulation & get_formulation() const {
return this->projection->get_formulation();}
//! for handling double initialisations right
bool is_initialised() const {return this->initialised;}
/**
* get a reference to the projection object. should only be
* required for debugging
*/
Eigen::Map<Eigen::ArrayXXd> get_projection() {
return this->projection->get_operator();}
//! returns the spatial size
constexpr static Dim_t get_sdim() {return DimS;};
//! return a sparse matrix adaptor to the cell
Adaptor get_adaptor();
//! returns the number of degrees of freedom in the cell
Dim_t nb_dof() const {return this->size()*ipow(DimS, 2);};
//! return the communicator object
const Communicator & get_communicator() const {
return this->projection->get_communicator();
}
protected:
//! make sure that every pixel is assigned to one and only one material
void check_material_coverage();
- const Ccoord & resolutions; //!< the cell's subdomain resolutions
- const Ccoord & locations; //!< the cell's subdomain resolutions
+ const Ccoord & subdomain_resolutions; //!< the cell's subdomain resolutions
+ const Ccoord & subdomain_locations; //!< the cell's subdomain resolutions
const Ccoord & domain_resolutions; //!< the cell's domain resolutions
CcoordOps::Pixels<DimS> pixels; //!< helper to iterate over the pixels
- const Rcoord & lengths; //!< the cell's lengths
+ const Rcoord & domain_lengths; //!< the cell's lengths
Collection_ptr fields; //!< handle for the global fields of the cell
StrainField_t & F; //!< ref to strain field
StressField_t & P; //!< ref to stress field
//! Tangent field might not even be required; so this is an
//! optional ref_wrapper instead of a ref
optional<std::reference_wrapper<TangentField_t>> K{};
//! container of the materials present in the cell
std::vector<Material_ptr> materials{};
Projection_ptr projection; //!< handle for the projection operator
bool initialised{false}; //!< to handle double initialisation right
const Formulation form; //!< formulation for solution
private:
};
/**
* lightweight resource handle wrapping a `muSpectre::CellBase` or
* a subclass thereof into `Eigen::EigenBase`, so it can be
* interpreted as a sparse matrix by Eigen solvers
*/
template <class Cell>
class CellAdaptor: public Eigen::EigenBase<CellAdaptor<Cell>> {
public:
using Scalar = double; //!< sparse matrix traits
using RealScalar = double; //!< sparse matrix traits
using StorageIndex = int; //!< sparse matrix traits
enum {
ColsAtCompileTime = Eigen::Dynamic,
MaxColsAtCompileTime = Eigen::Dynamic,
RowsAtCompileTime = Eigen::Dynamic,
MaxRowsAtCompileTime = Eigen::Dynamic,
IsRowMajor = false
};
//! constructor
CellAdaptor(Cell & cell):cell{cell}{}
//!returns the number of logical rows
Eigen::Index rows() const {return this->cell.nb_dof();}
//!returns the number of logical columns
Eigen::Index cols() const {return this->rows();}
//! implementation of the evaluation
template<typename Rhs>
Eigen::Product<CellAdaptor,Rhs,Eigen::AliasFreeProduct>
operator*(const Eigen::MatrixBase<Rhs>& x) const {
return Eigen::Product<CellAdaptor,Rhs,Eigen::AliasFreeProduct>
(*this, x.derived());
}
Cell & cell; //!< ref to the cell
};
} // muSpectre
namespace Eigen {
namespace internal {
//! Implementation of `muSpectre::CellAdaptor` * `Eigen::DenseVector` through a
//! specialization of `Eigen::internal::generic_product_impl`:
template<typename Rhs, class CellAdaptor>
struct generic_product_impl<CellAdaptor, Rhs, SparseShape, DenseShape, GemvProduct> // GEMV stands for matrix-vector
: generic_product_impl_base<CellAdaptor,Rhs,generic_product_impl<CellAdaptor,Rhs> >
{
//! undocumented
typedef typename Product<CellAdaptor,Rhs>::Scalar Scalar;
//! undocumented
template<typename Dest>
static void scaleAndAddTo(Dest& dst, const CellAdaptor& lhs, const Rhs& rhs, const Scalar& /*alpha*/)
{
// This method should implement "dst += alpha * lhs * rhs" inplace,
// however, for iterative solvers, alpha is always equal to 1, so let's not bother about it.
// Here we could simply call dst.noalias() += lhs.my_matrix() * rhs,
dst.noalias() += const_cast<CellAdaptor&>(lhs).cell.directional_stiffness_vec(rhs);
}
};
}
}
#endif /* CELL_BASE_H */
diff --git a/src/cell/cell_factory.hh b/src/cell/cell_factory.hh
index 8b15cf6..c6b491e 100644
--- a/src/cell/cell_factory.hh
+++ b/src/cell/cell_factory.hh
@@ -1,159 +1,159 @@
/**
* @file cell_factory.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 15 Dec 2017
*
* @brief Cell factories to help create cells with ease
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef CELL_FACTORY_H
#define CELL_FACTORY_H
#include "common/common.hh"
#include "common/ccoord_operations.hh"
#include "cell/cell_base.hh"
#include "fft/projection_finite_strain_fast.hh"
#include "fft/projection_small_strain.hh"
#include "fft/fftw_engine.hh"
#ifdef WITH_MPI
#include "common/communicator.hh"
#include "fft/fftwmpi_engine.hh"
#endif
#include <memory>
namespace muSpectre {
/**
* Create a unique ptr to a Projection operator (with appropriate
* FFT_engine) to be used in a cell constructor
*/
template <Dim_t DimS, Dim_t DimM,
typename FFTEngine=FFTWEngine<DimS, DimM>>
inline
std::unique_ptr<ProjectionBase<DimS, DimM>>
cell_input(Ccoord_t<DimS> resolutions,
Rcoord_t<DimS> lengths,
Formulation form) {
- auto fft_ptr{std::make_unique<FFTEngine>(resolutions, lengths)};
+ auto fft_ptr{std::make_unique<FFTEngine>(resolutions)};
switch (form)
{
case Formulation::finite_strain: {
using Projection = ProjectionFiniteStrainFast<DimS, DimM>;
- return std::make_unique<Projection>(std::move(fft_ptr));
+ return std::make_unique<Projection>(std::move(fft_ptr), lengths);
break;
}
case Formulation::small_strain: {
using Projection = ProjectionSmallStrain<DimS, DimM>;
- return std::make_unique<Projection>(std::move(fft_ptr));
+ return std::make_unique<Projection>(std::move(fft_ptr), lengths);
break;
}
default: {
throw std::runtime_error("unknow formulation");
break;
}
}
}
/**
* convenience function to create a cell (avoids having to build
* and move the chain of unique_ptrs
*/
template <size_t DimS, size_t DimM=DimS,
typename Cell=CellBase<DimS, DimM>,
typename FFTEngine=FFTWEngine<DimS, DimM>>
inline
Cell make_cell(Ccoord_t<DimS> resolutions,
Rcoord_t<DimS> lengths,
Formulation form) {
auto && input = cell_input<DimS, DimM, FFTEngine>(resolutions, lengths,
form);
auto cell{Cell{std::move(input)}};
return cell;
}
#ifdef WITH_MPI
/**
* Create a unique ptr to a parallel Projection operator (with appropriate
* FFT_engine) to be used in a cell constructor
*/
template <Dim_t DimS, Dim_t DimM,
typename FFTEngine=FFTWMPIEngine<DimS, DimM>>
inline
std::unique_ptr<ProjectionBase<DimS, DimM>>
parallel_cell_input(Ccoord_t<DimS> resolutions,
Rcoord_t<DimS> lengths,
Formulation form,
const Communicator & comm) {
- auto fft_ptr{std::make_unique<FFTEngine>(resolutions, lengths, comm)};
+ auto fft_ptr{std::make_unique<FFTEngine>(resolutions, comm)};
switch (form)
{
case Formulation::finite_strain: {
using Projection = ProjectionFiniteStrainFast<DimS, DimM>;
- return std::make_unique<Projection>(std::move(fft_ptr));
+ return std::make_unique<Projection>(std::move(fft_ptr), lengths);
break;
}
case Formulation::small_strain: {
using Projection = ProjectionSmallStrain<DimS, DimM>;
- return std::make_unique<Projection>(std::move(fft_ptr));
+ return std::make_unique<Projection>(std::move(fft_ptr), lengths);
break;
}
default: {
throw std::runtime_error("unknown formulation");
break;
}
}
}
/**
* convenience function to create a parallel cell (avoids having to build
* and move the chain of unique_ptrs
*/
template <size_t DimS, size_t DimM=DimS,
typename Cell=CellBase<DimS, DimM>,
typename FFTEngine=FFTWMPIEngine<DimS, DimM>>
inline
Cell make_parallel_cell(Ccoord_t<DimS> resolutions,
Rcoord_t<DimS> lengths,
Formulation form,
const Communicator & comm) {
auto && input = parallel_cell_input<DimS, DimM, FFTEngine>(resolutions,
lengths,
form, comm);
auto cell{Cell{std::move(input)}};
return cell;
}
#endif /* WITH_MPI */
} // muSpectre
#endif /* CELL_FACTORY_H */
diff --git a/src/fft/fft_engine_base.cc b/src/fft/fft_engine_base.cc
index 6440cfc..1036a88 100644
--- a/src/fft/fft_engine_base.cc
+++ b/src/fft/fft_engine_base.cc
@@ -1,66 +1,65 @@
/**
* @file fft_engine_base.cc
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 03 Dec 2017
*
* @brief implementation for FFT engine base class
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include "fft/fft_engine_base.hh"
namespace muSpectre {
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
- FFTEngineBase<DimS, DimM>::FFTEngineBase(Ccoord resolutions, Rcoord lengths,
+ FFTEngineBase<DimS, DimM>::FFTEngineBase(Ccoord resolutions,
Communicator comm)
- :comm{comm}, resolutions{resolutions}, locations{},
+ :comm{comm}, subdomain_resolutions{resolutions}, subdomain_locations{},
fourier_resolutions{CcoordOps::get_hermitian_sizes(resolutions)},
fourier_locations{}, domain_resolutions{resolutions},
- lengths{lengths},
work{make_field<Workspace_t>("work space", work_space_container)},
norm_factor{1./CcoordOps::get_size(domain_resolutions)}
{}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void FFTEngineBase<DimS, DimM>::initialise(FFT_PlanFlags /*plan_flags*/) {
this->work_space_container.initialise();
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
size_t FFTEngineBase<DimS, DimM>::size() const {
- return CcoordOps::get_size(this->resolutions);
+ return CcoordOps::get_size(this->subdomain_resolutions);
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
size_t FFTEngineBase<DimS, DimM>::workspace_size() const {
return this->work_space_container.size();
}
template class FFTEngineBase<twoD, twoD>;
template class FFTEngineBase<twoD, threeD>;
template class FFTEngineBase<threeD, threeD>;
} // muSpectre
diff --git a/src/fft/fft_engine_base.hh b/src/fft/fft_engine_base.hh
index 59740a7..450b9b7 100644
--- a/src/fft/fft_engine_base.hh
+++ b/src/fft/fft_engine_base.hh
@@ -1,163 +1,159 @@
/**
* @file fft_engine_base.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 01 Dec 2017
*
* @brief Interface for FFT engines
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef FFT_ENGINE_BASE_H
#define FFT_ENGINE_BASE_H
#include "common/common.hh"
#include "common/communicator.hh"
#include "common/field_collection.hh"
namespace muSpectre {
/**
* Virtual base class for FFT engines. To be implemented by all
* FFT_engine implementations.
*/
template <Dim_t DimS, Dim_t DimM>
class FFTEngineBase
{
public:
constexpr static Dim_t sdim{DimS}; //!< spatial dimension of the cell
constexpr static Dim_t mdim{DimM}; //!< material dimension of the cell
//! cell coordinates type
using Ccoord = Ccoord_t<DimS>;
- //! spatial coordinates type
- using Rcoord = std::array<Real, DimS>;
//! global FieldCollection
using GFieldCollection_t = GlobalFieldCollection<DimS>;
//! local FieldCollection (for Fourier-space pixels)
using LFieldCollection_t = LocalFieldCollection<DimS>;
//! Field type on which to apply the projection
using Field_t = TensorField<GFieldCollection_t, Real, 2, DimM>;
/**
* Field type holding a Fourier-space representation of a
* real-valued second-order tensor field
*/
using Workspace_t = TensorField<LFieldCollection_t, Complex, 2, DimM>;
/**
* iterator over Fourier-space discretisation point
*/
using iterator = typename LFieldCollection_t::iterator;
//! Default constructor
FFTEngineBase() = delete;
//! Constructor with cell resolutions
- FFTEngineBase(Ccoord resolutions, Rcoord lengths,
- Communicator comm=Communicator());
+ FFTEngineBase(Ccoord resolutions, Communicator comm=Communicator());
//! Copy constructor
FFTEngineBase(const FFTEngineBase &other) = delete;
//! Move constructor
FFTEngineBase(FFTEngineBase &&other) = default;
//! Destructor
virtual ~FFTEngineBase() = default;
//! Copy assignment operator
FFTEngineBase& operator=(const FFTEngineBase &other) = delete;
//! Move assignment operator
FFTEngineBase& operator=(FFTEngineBase &&other) = default;
//! compute the plan, etc
virtual void initialise(FFT_PlanFlags /*plan_flags*/);
//! forward transform (dummy for interface)
virtual Workspace_t & fft(Field_t & /*field*/) = 0;
//! inverse transform (dummy for interface)
virtual void ifft(Field_t & /*field*/) const = 0;
/**
* iterators over only those pixels that exist in frequency space
* (i.e. about half of all pixels, see rfft)
*/
//! returns an iterator to the first pixel in Fourier space
inline iterator begin() {return this->work_space_container.begin();}
//! returns an iterator past to the last pixel in Fourier space
inline iterator end() {return this->work_space_container.end();}
//! nb of pixels (mostly for debugging)
size_t size() const;
//! nb of pixels in Fourier space
size_t workspace_size() const;
//! return the communicator object
const Communicator & get_communicator() const {return this->comm;}
//! returns the process-local resolutions of the cell
- const Ccoord & get_resolutions() const {return this->resolutions;}
+ const Ccoord & get_subdomain_resolutions() const {
+ return this->subdomain_resolutions;}
//! returns the process-local locations of the cell
- const Ccoord & get_locations() const {return this->locations;}
+ const Ccoord & get_subdomain_locations() const {
+ return this->subdomain_locations;}
//! returns the process-local resolutions of the cell in Fourier space
const Ccoord & get_fourier_resolutions() const {return this->fourier_resolutions;}
//! returns the process-local locations of the cell in Fourier space
const Ccoord & get_fourier_locations() const {return this->fourier_locations;}
//! returns the resolutions of the cell
const Ccoord & get_domain_resolutions() const {return this->domain_resolutions;}
- //! returns the physical sizes of the cell
- const Rcoord & get_lengths() const {return this->lengths;}
//! only required for testing and debugging
LFieldCollection_t & get_field_collection() {
return this->work_space_container;}
//! only required for testing and debugging
Workspace_t& get_work_space() {return this->work;}
//! factor by which to multiply projection before inverse transform (this is
//! typically 1/nb_pixels for so-called unnormalized transforms (see,
//! e.g. http://www.fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html#Multi_002dDimensional-DFTs-of-Real-Data
//! or https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.fft.html
//! . Rather than scaling the inverse transform (which would cost one more
//! loop), FFT engines provide this value so it can be used in the
//! projection operator (where no additional loop is required)
inline Real normalisation() const {return norm_factor;};
protected:
/**
* Field collection in which to store fields associated with
* Fourier-space points
*/
Communicator comm; //!< communicator
LFieldCollection_t work_space_container{};
- Ccoord resolutions; //!< resolutions of the process-local (subdomain) portion of the cell
- Ccoord locations; // !< location of the process-local (subdomain) portion of the cell
+ Ccoord subdomain_resolutions; //!< resolutions of the process-local (subdomain) portion of the cell
+ Ccoord subdomain_locations; // !< location of the process-local (subdomain) portion of the cell
Ccoord fourier_resolutions; //!< resolutions of the process-local (subdomain) portion of the Fourier transformed data
Ccoord fourier_locations; // !< location of the process-local (subdomain) portion of the Fourier transformed data
const Ccoord domain_resolutions; //!< resolutions of the full domain of the cell
- const Rcoord lengths; //!< physical sizes of the cell
Workspace_t & work; //!< field to store the Fourier transform of P
const Real norm_factor; //!< normalisation coefficient of fourier transform
private:
};
} // muSpectre
#endif /* FFT_ENGINE_BASE_H */
diff --git a/src/fft/fftw_engine.cc b/src/fft/fftw_engine.cc
index 0bfabc5..3c98987 100644
--- a/src/fft/fftw_engine.cc
+++ b/src/fft/fftw_engine.cc
@@ -1,153 +1,155 @@
/**
* @file fftw_engine.cc
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 03 Dec 2017
*
* @brief implements the fftw engine
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include "common/ccoord_operations.hh"
#include "fft/fftw_engine.hh"
namespace muSpectre {
template <Dim_t DimS, Dim_t DimM>
- FFTWEngine<DimS, DimM>::FFTWEngine(Ccoord resolutions, Rcoord lengths,
- Communicator comm)
- :Parent{resolutions, lengths, comm}
+ FFTWEngine<DimS, DimM>::FFTWEngine(Ccoord resolutions, Communicator comm)
+ :Parent{resolutions, comm}
{
for (auto && pixel: CcoordOps::Pixels<DimS>(this->fourier_resolutions)) {
this->work_space_container.add_pixel(pixel);
}
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void FFTWEngine<DimS, DimM>::initialise(FFT_PlanFlags plan_flags) {
if (this->initialised) {
throw std::runtime_error("double initialisation, will leak memory");
}
Parent::initialise(plan_flags);
const int & rank = DimS;
std::array<int, DimS> narr;
const int * const n = &narr[0];
- std::copy(this->resolutions.begin(), this->resolutions.end(), narr.begin());
+ std::copy(this->subdomain_resolutions.begin(),
+ this->subdomain_resolutions.end(),
+ narr.begin());
int howmany = Field_t::nb_components;
//temporary buffer for plan
- size_t alloc_size = CcoordOps::get_size(this->resolutions) *howmany;
+ size_t alloc_size = (CcoordOps::get_size(this->subdomain_resolutions)*
+ howmany);
Real * r_work_space = fftw_alloc_real(alloc_size);
Real * in = r_work_space;
const int * const inembed = nullptr;//nembed are tricky: they refer to physical layout
int istride = howmany;
int idist = 1;
fftw_complex * out = reinterpret_cast<fftw_complex*>(this->work.data());
const int * const onembed = nullptr;
int ostride = howmany;
int odist = idist;
unsigned int flags;
switch (plan_flags) {
case FFT_PlanFlags::estimate: {
flags = FFTW_ESTIMATE;
break;
}
case FFT_PlanFlags::measure: {
flags = FFTW_MEASURE;
break;
}
case FFT_PlanFlags::patient: {
flags = FFTW_PATIENT;
break;
}
default:
throw std::runtime_error("unknown planner flag type");
break;
}
this->plan_fft = fftw_plan_many_dft_r2c(rank, n, howmany, in, inembed, istride,
idist, out, onembed, ostride, odist,
FFTW_PRESERVE_INPUT | flags);
if (this->plan_fft == nullptr) {
throw std::runtime_error("Plan failed");
}
fftw_complex * i_in = reinterpret_cast<fftw_complex*>(this->work.data());
Real * i_out = r_work_space;
this->plan_ifft = fftw_plan_many_dft_c2r(rank, n, howmany, i_in, inembed,
istride, idist, i_out, onembed,
ostride, odist, flags);
if (this->plan_ifft == nullptr) {
throw std::runtime_error("Plan failed");
}
fftw_free(r_work_space);
this->initialised = true;
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
FFTWEngine<DimS, DimM>::~FFTWEngine<DimS, DimM>() noexcept {
fftw_destroy_plan(this->plan_fft);
fftw_destroy_plan(this->plan_ifft);
// TODO: We cannot run fftw_cleanup since subsequent FFTW calls will fail
// but multiple FFT engines can be active at the same time.
//fftw_cleanup();
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
typename FFTWEngine<DimS, DimM>::Workspace_t &
FFTWEngine<DimS, DimM>::fft (Field_t & field) {
if (this->plan_fft == nullptr) {
throw std::runtime_error("fft plan not initialised");
}
- if (field.size() != CcoordOps::get_size(this->resolutions)) {
+ if (field.size() != CcoordOps::get_size(this->subdomain_resolutions)) {
throw std::runtime_error("size mismatch");
}
fftw_execute_dft_r2c(this->plan_fft,
field.data(),
reinterpret_cast<fftw_complex*>(this->work.data()));
return this->work;
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void
FFTWEngine<DimS, DimM>::ifft (Field_t & field) const {
if (this->plan_ifft == nullptr) {
throw std::runtime_error("ifft plan not initialised");
}
- if (field.size() != CcoordOps::get_size(this->resolutions)) {
+ if (field.size() != CcoordOps::get_size(this->subdomain_resolutions)) {
throw std::runtime_error("size mismatch");
}
fftw_execute_dft_c2r(this->plan_ifft,
reinterpret_cast<fftw_complex*>(this->work.data()),
field.data());
}
template class FFTWEngine<twoD, twoD>;
template class FFTWEngine<twoD, threeD>;
template class FFTWEngine<threeD, threeD>;
} // muSpectre
diff --git a/src/fft/fftw_engine.hh b/src/fft/fftw_engine.hh
index b453687..13907f6 100644
--- a/src/fft/fftw_engine.hh
+++ b/src/fft/fftw_engine.hh
@@ -1,92 +1,90 @@
/**
* @file fftw_engine.hh
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 03 Dec 2017
*
* @brief FFT engine using FFTW
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef FFTW_ENGINE_H
#define FFTW_ENGINE_H
#include "fft/fft_engine_base.hh"
#include <fftw3.h>
namespace muSpectre {
/**
* implements the `muSpectre::FftEngine_Base` interface using the
* FFTW library
*/
template <Dim_t DimS, Dim_t DimM>
class FFTWEngine: public FFTEngineBase<DimS, DimM>
{
public:
using Parent = FFTEngineBase<DimS, DimM>; //!< base class
using Ccoord = typename Parent::Ccoord; //!< cell coordinates type
- using Rcoord = typename Parent::Rcoord; //!< spatial coordinates type
//! field for Fourier transform of second-order tensor
using Workspace_t = typename Parent::Workspace_t;
//! real-valued second-order tensor
using Field_t = typename Parent::Field_t;
//! Default constructor
FFTWEngine() = delete;
//! Constructor with cell resolutions
- FFTWEngine(Ccoord resolutions, Rcoord lengths,
- Communicator comm=Communicator());
+ FFTWEngine(Ccoord resolutions, Communicator comm=Communicator());
//! Copy constructor
FFTWEngine(const FFTWEngine &other) = delete;
//! Move constructor
FFTWEngine(FFTWEngine &&other) = default;
//! Destructor
virtual ~FFTWEngine() noexcept;
//! Copy assignment operator
FFTWEngine& operator=(const FFTWEngine &other) = delete;
//! Move assignment operator
FFTWEngine& operator=(FFTWEngine &&other) = default;
// compute the plan, etc
virtual void initialise(FFT_PlanFlags plan_flags) override;
//! forward transform
virtual Workspace_t & fft(Field_t & field) override;
//! inverse transform
virtual void ifft(Field_t & field) const override;
protected:
fftw_plan plan_fft{}; //!< holds the plan for forward fourier transform
fftw_plan plan_ifft{}; //!< holds the plan for inverse fourier transform
bool initialised{false}; //!< to prevent double initialisation
private:
};
} // muSpectre
#endif /* FFTW_ENGINE_H */
diff --git a/src/fft/fftwmpi_engine.cc b/src/fft/fftwmpi_engine.cc
index 1bca7db..f5d6bf2 100644
--- a/src/fft/fftwmpi_engine.cc
+++ b/src/fft/fftwmpi_engine.cc
@@ -1,229 +1,235 @@
/**
* @file fftwmpi_engine.cc
*
* @author Lars Pastewka <lars.pastewka@imtek.uni-freiburg.de>
*
* @date 06 Mar 2017
*
* @brief implements the MPI-parallel fftw engine
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include "common/ccoord_operations.hh"
#include "fft/fftwmpi_engine.hh"
namespace muSpectre {
template <Dim_t DimsS, Dim_t DimM>
int FFTWMPIEngine<DimsS, DimM>::nb_engines{0};
template <Dim_t DimS, Dim_t DimM>
- FFTWMPIEngine<DimS, DimM>::FFTWMPIEngine(Ccoord resolutions, Rcoord lengths,
+ FFTWMPIEngine<DimS, DimM>::FFTWMPIEngine(Ccoord resolutions,
Communicator comm)
- :Parent{resolutions, lengths, comm}
+ :Parent{resolutions, comm}
{
if (!this->nb_engines) fftw_mpi_init();
this->nb_engines++;
std::array<ptrdiff_t, DimS> narr;
std::copy(this->domain_resolutions.begin(), this->domain_resolutions.end(),
narr.begin());
narr[DimS-1] = this->domain_resolutions[DimS-1]/2+1;
ptrdiff_t res_x, loc_x, res_y, loc_y;
this->workspace_size =
fftw_mpi_local_size_many_transposed(DimS, narr.data(),
Field_t::nb_components,
FFTW_MPI_DEFAULT_BLOCK,
FFTW_MPI_DEFAULT_BLOCK,
this->comm.get_mpi_comm(),
&res_x, &loc_x, &res_y, &loc_y);
this->fourier_resolutions[1] = this->fourier_resolutions[0];
this->fourier_locations[1] = this->fourier_locations[0];
- this->resolutions[0] = res_x;
- this->locations[0] = loc_x;
+ this->subdomain_resolutions[0] = res_x;
+ this->subdomain_locations[0] = loc_x;
this->fourier_resolutions[0] = res_y;
this->fourier_locations[0] = loc_y;
- for (auto & n: this->resolutions) {
+ for (auto & n: this->subdomain_resolutions) {
if (n == 0) {
throw std::runtime_error("FFTW MPI planning returned zero resolution. "
"You may need to run on fewer processes.");
}
}
for (auto & n: this->fourier_resolutions) {
if (n == 0) {
throw std::runtime_error("FFTW MPI planning returned zero Fourier "
"resolution. You may need to run on fewer "
"processes.");
}
}
for (auto && pixel:
std::conditional_t<
DimS==2,
CcoordOps::Pixels<DimS, 1, 0>,
CcoordOps::Pixels<DimS, 1, 0, 2>
>(this->fourier_resolutions, this->fourier_locations)) {
this->work_space_container.add_pixel(pixel);
}
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void FFTWMPIEngine<DimS, DimM>::initialise(FFT_PlanFlags plan_flags) {
if (this->initialised) {
throw std::runtime_error("double initialisation, will leak memory");
}
// Initialize parent after local resolutions have been determined and
// work space has been initialized
Parent::initialise(plan_flags);
this->real_workspace = fftw_alloc_real(2*this->workspace_size);
// We need to check whether the workspace provided by our field is large
// enough. MPI parallel FFTW may request a workspace size larger than the
// nominal size of the complex buffer.
if (long(this->work.size()*Field_t::nb_components) < this->workspace_size) {
this->work.set_pad_size(this->workspace_size -
Field_t::nb_components*this->work.size());
}
unsigned int flags;
switch (plan_flags) {
case FFT_PlanFlags::estimate: {
flags = FFTW_ESTIMATE;
break;
}
case FFT_PlanFlags::measure: {
flags = FFTW_MEASURE;
break;
}
case FFT_PlanFlags::patient: {
flags = FFTW_PATIENT;
break;
}
default:
throw std::runtime_error("unknown planner flag type");
break;
}
std::array<ptrdiff_t, DimS> narr;
std::copy(this->domain_resolutions.begin(), this->domain_resolutions.end(),
narr.begin());
Real * in{this->real_workspace};
fftw_complex * out{reinterpret_cast<fftw_complex*>(this->work.data())};
this->plan_fft = fftw_mpi_plan_many_dft_r2c(
DimS, narr.data(), Field_t::nb_components, FFTW_MPI_DEFAULT_BLOCK,
FFTW_MPI_DEFAULT_BLOCK, in, out, this->comm.get_mpi_comm(),
FFTW_MPI_TRANSPOSED_OUT | flags);
if (this->plan_fft == nullptr) {
throw std::runtime_error("r2c plan failed");
}
fftw_complex * i_in = reinterpret_cast<fftw_complex*>(this->work.data());
Real * i_out = this->real_workspace;
this->plan_ifft = fftw_mpi_plan_many_dft_c2r(
DimS, narr.data(), Field_t::nb_components, FFTW_MPI_DEFAULT_BLOCK,
FFTW_MPI_DEFAULT_BLOCK, i_in, i_out, this->comm.get_mpi_comm(),
FFTW_MPI_TRANSPOSED_IN | flags);
if (this->plan_ifft == nullptr) {
throw std::runtime_error("c2r plan failed");
}
this->initialised = true;
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
FFTWMPIEngine<DimS, DimM>::~FFTWMPIEngine<DimS, DimM>() noexcept {
if (this->real_workspace != nullptr) fftw_free(this->real_workspace);
if (this->plan_fft != nullptr) fftw_destroy_plan(this->plan_fft);
if (this->plan_ifft != nullptr) fftw_destroy_plan(this->plan_ifft);
// TODO: We cannot run fftw_mpi_cleanup since also calls fftw_cleanup
// and any running FFTWEngine will fail afterwards.
//this->nb_engines--;
//if (!this->nb_engines) fftw_mpi_cleanup();
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
typename FFTWMPIEngine<DimS, DimM>::Workspace_t &
FFTWMPIEngine<DimS, DimM>::fft (Field_t & field) {
if (this->plan_fft == nullptr) {
throw std::runtime_error("fft plan not initialised");
}
- if (field.size() != CcoordOps::get_size(this->resolutions)) {
+ if (field.size() != CcoordOps::get_size(this->subdomain_resolutions)) {
throw std::runtime_error("size mismatch");
}
// Copy non-padded field to padded real_workspace.
// Transposed output of M x N x L transform for >= 3 dimensions is padded
// M x N x 2*(L/2+1).
- ptrdiff_t fstride = Field_t::nb_components*this->resolutions[DimS-1];
- ptrdiff_t wstride = Field_t::nb_components*2*(this->resolutions[DimS-1]/2+1);
- ptrdiff_t n = field.size()/this->resolutions[DimS-1];
+ ptrdiff_t fstride = (Field_t::nb_components*
+ this->subdomain_resolutions[DimS-1]);
+ ptrdiff_t wstride = (Field_t::nb_components*2*
+ (this->subdomain_resolutions[DimS-1]/2+1));
+ ptrdiff_t n = field.size()/this->subdomain_resolutions[DimS-1];
auto fdata = field.data();
auto wdata = this->real_workspace;
for (int i = 0; i < n; i++) {
std::copy(fdata, fdata+fstride, wdata);
fdata += fstride;
wdata += wstride;
}
// Compute FFT
fftw_mpi_execute_dft_r2c(
this->plan_fft, this->real_workspace,
reinterpret_cast<fftw_complex*>(this->work.data()));
return this->work;
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void
FFTWMPIEngine<DimS, DimM>::ifft (Field_t & field) const {
if (this->plan_ifft == nullptr) {
throw std::runtime_error("ifft plan not initialised");
}
- if (field.size() != CcoordOps::get_size(this->resolutions)) {
+ if (field.size() != CcoordOps::get_size(this->subdomain_resolutions)) {
throw std::runtime_error("size mismatch");
}
// Compute inverse FFT
fftw_mpi_execute_dft_c2r(
this->plan_ifft, reinterpret_cast<fftw_complex*>(this->work.data()),
this->real_workspace);
// Copy non-padded field to padded real_workspace.
// Transposed output of M x N x L transform for >= 3 dimensions is padded
// M x N x 2*(L/2+1).
- ptrdiff_t fstride{Field_t::nb_components*this->resolutions[DimS-1]};
- ptrdiff_t wstride{Field_t::nb_components*2*(this->resolutions[DimS-1]/2+1)};
- ptrdiff_t n(field.size()/this->resolutions[DimS-1]);
+ ptrdiff_t fstride{
+ Field_t::nb_components*this->subdomain_resolutions[DimS-1]
+ };
+ ptrdiff_t wstride{
+ Field_t::nb_components*2*(this->subdomain_resolutions[DimS-1]/2+1)
+ };
+ ptrdiff_t n(field.size()/this->subdomain_resolutions[DimS-1]);
auto fdata{field.data()};
auto wdata{this->real_workspace};
for (int i = 0; i < n; i++) {
std::copy(wdata, wdata+fstride, fdata);
fdata += fstride;
wdata += wstride;
}
}
template class FFTWMPIEngine<twoD, twoD>;
template class FFTWMPIEngine<twoD, threeD>;
template class FFTWMPIEngine<threeD, threeD>;
} // muSpectre
diff --git a/src/fft/fftwmpi_engine.hh b/src/fft/fftwmpi_engine.hh
index 5f16b5a..7ce4fe7 100644
--- a/src/fft/fftwmpi_engine.hh
+++ b/src/fft/fftwmpi_engine.hh
@@ -1,95 +1,93 @@
/**
* @file fftwmpi_engine.hh
*
* @author Lars Pastewka <lars.pastewka@imtek.uni-freiburg.de>
*
* @date 06 Mar 2017
*
* @brief FFT engine using MPI-parallel FFTW
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef FFTWMPI_ENGINE_H
#define FFTWMPI_ENGINE_H
#include "fft/fft_engine_base.hh"
#include <fftw3-mpi.h>
namespace muSpectre {
/**
* implements the `muSpectre::FFTEngineBase` interface using the
* FFTW library
*/
template <Dim_t DimS, Dim_t DimM>
class FFTWMPIEngine: public FFTEngineBase<DimS, DimM>
{
public:
using Parent = FFTEngineBase<DimS, DimM>; //!< base class
using Ccoord = typename Parent::Ccoord; //!< cell coordinates type
- using Rcoord = typename Parent::Rcoord; //!< spatial coordinates type
//! field for Fourier transform of second-order tensor
using Workspace_t = typename Parent::Workspace_t;
//! real-valued second-order tensor
using Field_t = typename Parent::Field_t;
//! Default constructor
FFTWMPIEngine() = delete;
//! Constructor with system resolutions
- FFTWMPIEngine(Ccoord resolutions, Rcoord lengths,
- Communicator comm=Communicator());
+ FFTWMPIEngine(Ccoord resolutions, Communicator comm=Communicator());
//! Copy constructor
FFTWMPIEngine(const FFTWMPIEngine &other) = delete;
//! Move constructor
FFTWMPIEngine(FFTWMPIEngine &&other) = default;
//! Destructor
virtual ~FFTWMPIEngine() noexcept;
//! Copy assignment operator
FFTWMPIEngine& operator=(const FFTWMPIEngine &other) = delete;
//! Move assignment operator
FFTWMPIEngine& operator=(FFTWMPIEngine &&other) = default;
// compute the plan, etc
virtual void initialise(FFT_PlanFlags plan_flags) override;
//! forward transform
virtual Workspace_t & fft(Field_t & field) override;
//! inverse transform
virtual void ifft(Field_t & field) const override;
protected:
static int nb_engines; //!< number of times this engine has been instatiated
fftw_plan plan_fft{}; //!< holds the plan for forward fourier transform
fftw_plan plan_ifft{}; //!< holds the plan for inverse fourier transform
ptrdiff_t workspace_size{}; //!< size of workspace buffer returned by planner
Real *real_workspace{}; //!< temporary real workspace that is correctly padded
bool initialised{false}; //!< to prevent double initialisation
private:
};
} // muSpectre
#endif /* FFTWMPI_ENGINE_H */
diff --git a/src/fft/pfft_engine.cc b/src/fft/pfft_engine.cc
index 9afbe4f..ffe56f7 100644
--- a/src/fft/pfft_engine.cc
+++ b/src/fft/pfft_engine.cc
@@ -1,248 +1,247 @@
/**
* @file pfft_engine.cc
*
* @author Lars Pastewka <lars.pastewka@imtek.uni-freiburg.de>
*
* @date 06 Mar 2017
*
* @brief implements the MPI-parallel pfft engine
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include "common/ccoord_operations.hh"
#include "fft/pfft_engine.hh"
namespace muSpectre {
template <Dim_t DimsS, Dim_t DimM>
int PFFTEngine<DimsS, DimM>::nb_engines{0};
template <Dim_t DimS, Dim_t DimM>
- PFFTEngine<DimS, DimM>::PFFTEngine(Ccoord resolutions, Rcoord lengths,
- Communicator comm)
- :Parent{resolutions, lengths, comm}, mpi_comm{comm.get_mpi_comm()}
+ PFFTEngine<DimS, DimM>::PFFTEngine(Ccoord resolutions, Communicator comm)
+ :Parent{resolutions, comm}, mpi_comm{comm.get_mpi_comm()}
{
if (!this->nb_engines) pfft_init();
this->nb_engines++;
int size{comm.size()};
int dim_x{size};
int dim_y{1};
// Note: All TODOs below don't affect the function of the PFFT engine. It
// presently uses slab decompositions, the TODOs are what needs to be done
// to get stripe decomposition to work - but it does not work yet. Slab
// vs stripe decomposition may have an impact on how the code scales.
// TODO: Enable this to enable 2d process mesh. This does not pass tests.
//if (DimS > 2) {
if (false) {
dim_y = int(sqrt(size));
while ((size/dim_y)*dim_y != size) dim_y--;
dim_x = size/dim_y;
}
// TODO: Enable this to enable 2d process mesh. This does not pass tests.
//if (DimS > 2) {
if (false) {
if (pfft_create_procmesh_2d(this->comm.get_mpi_comm(), dim_x, dim_y,
&this->mpi_comm)) {
throw std::runtime_error("Failed to create 2d PFFT process mesh.");
}
}
std::array<ptrdiff_t, DimS> narr;
std::copy(this->domain_resolutions.begin(), this->domain_resolutions.end(),
narr.begin());
ptrdiff_t res[DimS], loc[DimS], fres[DimS], floc[DimS];
this->workspace_size =
pfft_local_size_many_dft_r2c(DimS, narr.data(), narr.data(), narr.data(),
Field_t::nb_components,
PFFT_DEFAULT_BLOCKS, PFFT_DEFAULT_BLOCKS,
this->mpi_comm,
PFFT_TRANSPOSED_OUT,
res, loc, fres, floc);
- std::copy(res, res+DimS, this->resolutions.begin());
- std::copy(loc, loc+DimS, this->locations.begin());
+ std::copy(res, res+DimS, this->subdomain_resolutions.begin());
+ std::copy(loc, loc+DimS, this->subdomain_locations.begin());
std::copy(fres, fres+DimS, this->fourier_resolutions.begin());
std::copy(floc, floc+DimS, this->fourier_locations.begin());
// TODO: Enable this to enable 2d process mesh. This does not pass tests.
//for (int i = 0; i < DimS-1; ++i) {
for (int i = 0; i < 1; ++i) {
std::swap(this->fourier_resolutions[i], this->fourier_resolutions[i+1]);
std::swap(this->fourier_locations[i], this->fourier_locations[i+1]);
}
- for (auto & n: this->resolutions) {
+ for (auto & n: this->subdomain_resolutions) {
if (n == 0) {
throw std::runtime_error("PFFT planning returned zero resolution. "
"You may need to run on fewer processes.");
}
}
for (auto & n: this->fourier_resolutions) {
if (n == 0) {
throw std::runtime_error("PFFT planning returned zero Fourier "
"resolution. You may need to run on fewer "
"processes.");
}
}
for (auto && pixel:
std::conditional_t<
DimS==2,
CcoordOps::Pixels<DimS, 1, 0>,
// TODO: This should be the correct order of dimension for a 2d
// process mesh, but tests don't pass.
//CcoordOps::Pixels<DimS, 1, 2, 0>
CcoordOps::Pixels<DimS, 1, 0, 2>
>(this->fourier_resolutions, this->fourier_locations)) {
this->work_space_container.add_pixel(pixel);
}
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void PFFTEngine<DimS, DimM>::initialise(FFT_PlanFlags plan_flags) {
if (this->initialised) {
throw std::runtime_error("double initialisation, will leak memory");
}
// Initialize parent after local resolutions have been determined and
// work space has been initialized
Parent::initialise(plan_flags);
this->real_workspace = pfft_alloc_real(2*this->workspace_size);
// We need to check whether the workspace provided by our field is large
// enough. PFFT may request a workspace size larger than the nominal size
// of the complex buffer.
if (long(this->work.size()*Field_t::nb_components) < this->workspace_size) {
this->work.set_pad_size(this->workspace_size -
Field_t::nb_components*this->work.size());
}
unsigned int flags;
switch (plan_flags) {
case FFT_PlanFlags::estimate: {
flags = PFFT_ESTIMATE;
break;
}
case FFT_PlanFlags::measure: {
flags = PFFT_MEASURE;
break;
}
case FFT_PlanFlags::patient: {
flags = PFFT_PATIENT;
break;
}
default:
throw std::runtime_error("unknown planner flag type");
break;
}
std::array<ptrdiff_t, DimS> narr;
std::copy(this->domain_resolutions.begin(), this->domain_resolutions.end(),
narr.begin());
Real * in{this->real_workspace};
pfft_complex * out{reinterpret_cast<pfft_complex*>(this->work.data())};
this->plan_fft = pfft_plan_many_dft_r2c(DimS, narr.data(), narr.data(),
narr.data(), Field_t::nb_components,
PFFT_DEFAULT_BLOCKS,
PFFT_DEFAULT_BLOCKS,
in, out, this->mpi_comm,
PFFT_FORWARD,
PFFT_TRANSPOSED_OUT | flags);
if (this->plan_fft == nullptr) {
throw std::runtime_error("r2c plan failed");
}
pfft_complex * i_in{reinterpret_cast<pfft_complex*>(this->work.data())};
Real * i_out{this->real_workspace};
this->plan_ifft = pfft_plan_many_dft_c2r(DimS, narr.data(), narr.data(),
narr.data(),
Field_t::nb_components,
PFFT_DEFAULT_BLOCKS,
PFFT_DEFAULT_BLOCKS,
i_in, i_out,
this->mpi_comm,
PFFT_BACKWARD,
PFFT_TRANSPOSED_IN | flags);
if (this->plan_ifft == nullptr) {
throw std::runtime_error("c2r plan failed");
}
this->initialised = true;
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
PFFTEngine<DimS, DimM>::~PFFTEngine<DimS, DimM>() noexcept {
if (this->real_workspace != nullptr) pfft_free(this->real_workspace);
if (this->plan_fft != nullptr) pfft_destroy_plan(this->plan_fft);
if (this->plan_ifft != nullptr) pfft_destroy_plan(this->plan_ifft);
if (this->mpi_comm != this->comm.get_mpi_comm()) {
MPI_Comm_free(&this->mpi_comm);
}
// TODO: We cannot run fftw_mpi_cleanup since also calls fftw_cleanup
// and any running FFTWEngine will fail afterwards.
//this->nb_engines--;
//if (!this->nb_engines) pfft_cleanup();
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
typename PFFTEngine<DimS, DimM>::Workspace_t &
PFFTEngine<DimS, DimM>::fft (Field_t & field) {
if (!this->plan_fft) {
throw std::runtime_error("fft plan not allocated");
}
- if (field.size() != CcoordOps::get_size(this->resolutions)) {
+ if (field.size() != CcoordOps::get_size(this->subdomain_resolutions)) {
throw std::runtime_error("size mismatch");
}
// Copy field data to workspace buffer. This is necessary because workspace
// buffer is larger than field buffer.
std::copy(field.data(), field.data()+Field_t::nb_components*field.size(),
this->real_workspace);
pfft_execute_dft_r2c(
this->plan_fft, this->real_workspace,
reinterpret_cast<pfft_complex*>(this->work.data()));
return this->work;
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void
PFFTEngine<DimS, DimM>::ifft (Field_t & field) const {
if (!this->plan_ifft) {
throw std::runtime_error("ifft plan not allocated");
}
- if (field.size() != CcoordOps::get_size(this->resolutions)) {
+ if (field.size() != CcoordOps::get_size(this->subdomain_resolutions)) {
throw std::runtime_error("size mismatch");
}
pfft_execute_dft_c2r(
this->plan_ifft, reinterpret_cast<pfft_complex*>(this->work.data()),
this->real_workspace);
std::copy(this->real_workspace,
this->real_workspace+Field_t::nb_components*field.size(),
field.data());
}
template class PFFTEngine<twoD, twoD>;
template class PFFTEngine<twoD, threeD>;
template class PFFTEngine<threeD, threeD>;
} // muSpectre
diff --git a/src/fft/pfft_engine.hh b/src/fft/pfft_engine.hh
index ac2c1f6..b7d5f13 100644
--- a/src/fft/pfft_engine.hh
+++ b/src/fft/pfft_engine.hh
@@ -1,98 +1,96 @@
/**
* @file pfft_engine.hh
*
* @author Lars Pastewka <lars.pastewka@imtek.uni-freiburg.de>
*
* @date 06 Mar 2017
*
* @brief FFT engine using MPI-parallel PFFT
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef PFFT_ENGINE_H
#define PFFT_ENGINE_H
#include "common/communicator.hh"
#include "fft/fft_engine_base.hh"
#include <pfft.h>
namespace muSpectre {
/**
* implements the `muSpectre::FFTEngineBase` interface using the
* FFTW library
*/
template <Dim_t DimS, Dim_t DimM>
class PFFTEngine: public FFTEngineBase<DimS, DimM>
{
public:
using Parent = FFTEngineBase<DimS, DimM>; //!< base class
using Ccoord = typename Parent::Ccoord; //!< cell coordinates type
- using Rcoord = typename Parent::Rcoord; //!< spatial coordinates type
//! field for Fourier transform of second-order tensor
using Workspace_t = typename Parent::Workspace_t;
//! real-valued second-order tensor
using Field_t = typename Parent::Field_t;
//! Default constructor
PFFTEngine() = delete;
//! Constructor with system resolutions
- PFFTEngine(Ccoord resolutions, Rcoord lengths,
- Communicator comm=Communicator());
+ PFFTEngine(Ccoord resolutions, Communicator comm=Communicator());
//! Copy constructor
PFFTEngine(const PFFTEngine &other) = delete;
//! Move constructor
PFFTEngine(PFFTEngine &&other) = default;
//! Destructor
virtual ~PFFTEngine() noexcept;
//! Copy assignment operator
PFFTEngine& operator=(const PFFTEngine &other) = delete;
//! Move assignment operator
PFFTEngine& operator=(PFFTEngine &&other) = default;
// compute the plan, etc
virtual void initialise(FFT_PlanFlags plan_flags) override;
//! forward transform
virtual Workspace_t & fft(Field_t & field) override;
//! inverse transform
virtual void ifft(Field_t & field) const override;
protected:
MPI_Comm mpi_comm; //! < MPI communicator
static int nb_engines; //!< number of times this engine has been instatiated
pfft_plan plan_fft{}; //!< holds the plan for forward fourier transform
pfft_plan plan_ifft{}; //!< holds the plan for inverse fourier transform
ptrdiff_t workspace_size{}; //!< size of workspace buffer returned by planner
Real *real_workspace{}; //!< temporary real workspace that is correctly padded
bool initialised{false}; //!< to prevent double initialisation
private:
};
} // muSpectre
#endif /* PFFT_ENGINE_H */
diff --git a/src/fft/projection_base.cc b/src/fft/projection_base.cc
index e063fb7..c291fc6 100644
--- a/src/fft/projection_base.cc
+++ b/src/fft/projection_base.cc
@@ -1,56 +1,57 @@
/**
* @file projection_base.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 06 Dec 2017
*
* @brief implementation of base class for projections
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include "fft/projection_base.hh"
namespace muSpectre {
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
ProjectionBase<DimS, DimM>::ProjectionBase(FFTEngine_ptr engine,
+ Rcoord domain_lengths,
Formulation form)
- : fft_engine{std::move(engine)}, form{form},
+ : fft_engine{std::move(engine)}, domain_lengths{domain_lengths}, form{form},
projection_container{this->fft_engine->get_field_collection()}
{
static_assert((DimS == FFTEngine::sdim),
"spatial dimensions are incompatible");
static_assert((DimM == FFTEngine::mdim),
"material dimensions are incompatible");
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void ProjectionBase<DimS, DimM>::
initialise(FFT_PlanFlags flags) {
fft_engine->initialise(flags);
}
template class ProjectionBase<twoD, twoD>;
template class ProjectionBase<twoD, threeD>;
template class ProjectionBase<threeD, threeD>;
} // muSpectre
diff --git a/src/fft/projection_base.hh b/src/fft/projection_base.hh
index 57a7768..96bad7b 100644
--- a/src/fft/projection_base.hh
+++ b/src/fft/projection_base.hh
@@ -1,169 +1,169 @@
/**
* @file projection_base.hh
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 03 Dec 2017
*
* @brief Base class for Projection operators
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef PROJECTION_BASE_H
#define PROJECTION_BASE_H
#include "common/common.hh"
#include "common/field_collection.hh"
#include "common/field.hh"
#include "fft/fft_engine_base.hh"
#include <memory>
namespace muSpectre {
/**
* base class for projection related exceptions
*/
class ProjectionError: public std::runtime_error {
public:
//! constructor
explicit ProjectionError(const std::string& what)
:std::runtime_error(what){}
//! constructor
explicit ProjectionError(const char * what)
:std::runtime_error(what){}
};
template<class Projection>
struct Projection_traits {
};
/**
* defines the interface which must be implemented by projection operators
*/
template <Dim_t DimS, Dim_t DimM>
class ProjectionBase
{
public:
//! type of fft_engine used
using FFTEngine = FFTEngineBase<DimS, DimM>;
//! reference to fft engine is safely managed through a `std::unique_ptr`
using FFTEngine_ptr = std::unique_ptr<FFTEngine>;
//! cell coordinates type
using Ccoord = typename FFTEngine::Ccoord;
//! spatial coordinates type
- using Rcoord = typename FFTEngine::Rcoord;
+ using Rcoord = Rcoord_t<DimS>;
//! global FieldCollection
using GFieldCollection_t = typename FFTEngine::GFieldCollection_t;
//! local FieldCollection (for Fourier-space pixels)
using LFieldCollection_t = typename FFTEngine::LFieldCollection_t;
//! Field type on which to apply the projection
using Field_t = typename FFTEngine::Field_t;
/**
* iterator over all pixels. This is taken from the FFT engine,
* because depending on the real-to-complex FFT employed, only
* roughly half of the pixels are present in Fourier space
* (because of the hermitian nature of the transform)
*/
using iterator = typename FFTEngine::iterator;
//! Default constructor
ProjectionBase() = delete;
//! Constructor with cell sizes
- ProjectionBase(FFTEngine_ptr engine, Formulation form);
+ ProjectionBase(FFTEngine_ptr engine, Rcoord domain_lengths, Formulation form);
//! Copy constructor
ProjectionBase(const ProjectionBase &other) = delete;
//! Move constructor
ProjectionBase(ProjectionBase &&other) = default;
//! Destructor
virtual ~ProjectionBase() = default;
//! Copy assignment operator
ProjectionBase& operator=(const ProjectionBase &other) = delete;
//! Move assignment operator
ProjectionBase& operator=(ProjectionBase &&other) = default;
//! initialises the fft engine (plan the transform)
virtual void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate);
//! apply the projection operator to a field
virtual void apply_projection(Field_t & field) = 0;
//! returns the process-local resolutions of the cell
- const Ccoord & get_resolutions() const {
- return this->fft_engine->get_resolutions();}
+ const Ccoord & get_subdomain_resolutions() const {
+ return this->fft_engine->get_subdomain_resolutions();}
//! returns the process-local locations of the cell
- const Ccoord & get_locations() const {
- return this->fft_engine->get_locations();}
+ const Ccoord & get_subdomain_locations() const {
+ return this->fft_engine->get_subdomain_locations();}
//! returns the resolutions of the cell
const Ccoord & get_domain_resolutions() const {
return this->fft_engine->get_domain_resolutions();}
//! returns the physical sizes of the cell
- const Rcoord & get_lengths() const {
- return this->fft_engine->get_lengths();}
+ const Rcoord & get_domain_lengths() const {return this->domain_lengths;}
/**
* return the `muSpectre::Formulation` that is used in solving
* this cell. This allows tho check whether a projection is
* compatible with the chosen formulation
*/
const Formulation & get_formulation() const {return this->form;}
//! return the raw projection operator. This is mainly intended
//! for maintenance and debugging and should never be required in
//! regular use
virtual Eigen::Map<Eigen::ArrayXXd> get_operator() = 0;
//! return the communicator object
const Communicator & get_communicator() const {
return this->fft_engine->get_communicator();
}
protected:
//! handle on the fft_engine used
FFTEngine_ptr fft_engine;
+ const Rcoord domain_lengths; //!< physical sizes of the cell
/**
* formulation this projection can be applied to (determines
* whether the projection enforces gradients, small strain tensor
* or symmetric smal strain tensor
*/
const Formulation form;
/**
* A local `muSpectre::FieldCollection` to store the projection
* operator per k-space point. This is a local rather than a
* global collection, since the pixels considered depend on the
* FFT implementation. See
* http://www.fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html#Multi_002dDimensional-DFTs-of-Real-Data
* for an example
*/
LFieldCollection_t & projection_container{};
private:
};
} // muSpectre
#endif /* PROJECTION_BASE_H */
diff --git a/src/fft/projection_default.cc b/src/fft/projection_default.cc
index 62a9b39..426ddf1 100644
--- a/src/fft/projection_default.cc
+++ b/src/fft/projection_default.cc
@@ -1,65 +1,66 @@
/**
* @file projection_default.cc
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 14 Jan 2018
*
* @brief Implementation default projection implementation
*
* Copyright © 2018 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include "fft/projection_default.hh"
#include "fft/fft_engine_base.hh"
namespace muSpectre {
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
ProjectionDefault<DimS, DimM>::ProjectionDefault(FFTEngine_ptr engine,
+ Rcoord lengths,
Formulation form)
- :Parent{std::move(engine), form},
+ :Parent{std::move(engine), lengths, form},
Gfield{make_field<Proj_t>("Projection Operator",
this->projection_container)},
Ghat{Gfield}
{}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void ProjectionDefault<DimS, DimM>::apply_projection(Field_t & field) {
Vector_map field_map{this->fft_engine->fft(field)};
Real factor = this->fft_engine->normalisation();
for (auto && tup: akantu::zip(this->Ghat, field_map)) {
auto & G{std::get<0>(tup)};
auto & f{std::get<1>(tup)};
f = factor * (G*f).eval();
}
this->fft_engine->ifft(field);
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
Eigen::Map<Eigen::ArrayXXd> ProjectionDefault<DimS, DimM>::get_operator() {
return this->Gfield.dyn_eigen();
}
template class ProjectionDefault<twoD, twoD>;
template class ProjectionDefault<threeD, threeD>;
} // muSpectre
diff --git a/src/fft/projection_default.hh b/src/fft/projection_default.hh
index d8edb63..9c2cbab 100644
--- a/src/fft/projection_default.hh
+++ b/src/fft/projection_default.hh
@@ -1,98 +1,99 @@
/**
* @file projection_default.hh
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 14 Jan 2018
*
* @brief virtual base class for default projection implementation, where the
* projection operator is stored as a full fourth-order tensor per
* k-space point (as opposed to 'smart' faster implementations, such as
* ProjectionFiniteStrainFast
*
* Copyright (C) 2018 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef PROJECTION_DEFAULT_H
#define PROJECTION_DEFAULT_H
#include "fft/projection_base.hh"
namespace muSpectre {
/**
* base class to inherit from if one implements a projection
* operator that is stored in form of a fourth-order tensor of real
* values per k-grid point
*/
template <Dim_t DimS, Dim_t DimM>
class ProjectionDefault: public ProjectionBase<DimS, DimM>
{
public:
using Parent = ProjectionBase<DimS, DimM>; //!< base class
//! polymorphic FFT pointer type
using FFTEngine_ptr = typename Parent::FFTEngine_ptr;
using Ccoord = typename Parent::Ccoord; //!< cell coordinates type
+ using Rcoord = typename Parent::Rcoord; //!< spatial coordinates type
//! global field collection
using GFieldCollection_t = GlobalFieldCollection<DimS>;
//! local field collection for Fourier-space fields
using LFieldCollection_t = LocalFieldCollection<DimS>;
//! Real space second order tensor fields (to be projected)
using Field_t = TensorField<GFieldCollection_t, Real, secondOrder, DimM>;
//! Fourier-space field containing the projection operator itself
using Proj_t = TensorField<LFieldCollection_t, Real, fourthOrder, DimM>;
//! iterable form of the operator
using Proj_map = T4MatrixFieldMap<LFieldCollection_t, Real, DimM>;
//! vectorized version of the Fourier-space second-order tensor field
using Vector_map = MatrixFieldMap<LFieldCollection_t, Complex, DimM*DimM, 1>;
//! Default constructor
ProjectionDefault() = delete;
//! Constructor with cell sizes and formulation
- ProjectionDefault(FFTEngine_ptr engine, Formulation form);
+ ProjectionDefault(FFTEngine_ptr engine, Rcoord lengths, Formulation form);
//! Copy constructor
ProjectionDefault(const ProjectionDefault &other) = delete;
//! Move constructor
ProjectionDefault(ProjectionDefault &&other) = default;
//! Destructor
virtual ~ProjectionDefault() = default;
//! Copy assignment operator
ProjectionDefault& operator=(const ProjectionDefault &other) = delete;
//! Move assignment operator
ProjectionDefault& operator=(ProjectionDefault &&other) = delete;
//! apply the projection operator to a field
void apply_projection(Field_t & field) override final;
Eigen::Map<Eigen::ArrayXXd> get_operator() override final;
protected:
Proj_t & Gfield; //!< field holding the operator
Proj_map Ghat; //!< iterable version of operator
private:
};
} // muSpectre
#endif /* PROJECTION_DEFAULT_H */
diff --git a/src/fft/projection_finite_strain.cc b/src/fft/projection_finite_strain.cc
index 292ab27..55695ef 100644
--- a/src/fft/projection_finite_strain.cc
+++ b/src/fft/projection_finite_strain.cc
@@ -1,92 +1,92 @@
/**
* @file projection_finite_strain.cc
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 05 Dec 2017
*
* @brief implementation of standard finite strain projection operator
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include "fft/projection_finite_strain.hh"
#include "fft/fftw_engine.hh"
#include "fft/fft_utils.hh"
#include "common/field_map.hh"
#include "common/tensor_algebra.hh"
#include "common/iterators.hh"
#include "Eigen/Dense"
namespace muSpectre {
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
ProjectionFiniteStrain<DimS, DimM>::
- ProjectionFiniteStrain(FFTEngine_ptr engine)
- :Parent{std::move(engine), Formulation::finite_strain}
+ ProjectionFiniteStrain(FFTEngine_ptr engine, Rcoord lengths)
+ :Parent{std::move(engine), lengths, Formulation::finite_strain}
{
for (auto res: this->fft_engine->get_domain_resolutions()) {
if (res % 2 == 0) {
throw ProjectionError
("Only an odd number of gridpoints in each direction is supported");
}
}
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void ProjectionFiniteStrain<DimS, DimM>::
initialise(FFT_PlanFlags flags) {
Parent::initialise(flags);
FFT_freqs<DimS> fft_freqs(this->fft_engine->get_domain_resolutions(),
- this->fft_engine->get_lengths());
+ this->domain_lengths);
for (auto && tup: akantu::zip(*this->fft_engine, this->Ghat)) {
const auto & ccoord = std::get<0> (tup);
auto & G = std::get<1>(tup);
auto xi = fft_freqs.get_unit_xi(ccoord);
//! this is simplifiable using Curnier's Méthodes numériques, 6.69(c)
G = Matrices::outer_under(Matrices::I2<DimM>(), xi*xi.transpose());
// The commented block below corresponds to the original
// definition of the operator in de Geus et
// al. (https://doi.org/10.1016/j.cma.2016.12.032). However,
// they use a bizarre definition of the double contraction
// between fourth-order and second-order tensors that has a
// built-in transpose operation (i.e., C = A:B <-> AᵢⱼₖₗBₗₖ =
// Cᵢⱼ , note the inverted ₗₖ instead of ₖₗ), here, we define
// the double contraction without the transposition. As a
// result, the Projection operator produces the transpose of de
// Geus's
// for (Dim_t im = 0; im < DimS; ++im) {
// for (Dim_t j = 0; j < DimS; ++j) {
// for (Dim_t l = 0; l < DimS; ++l) {
// get(G, im, j, l, im) = xi(j)*xi(l);
// }
// }
// }
}
- if (this->get_locations() == Ccoord{}) {
+ if (this->get_subdomain_locations() == Ccoord{}) {
this->Ghat[0].setZero();
}
}
template class ProjectionFiniteStrain<twoD, twoD>;
template class ProjectionFiniteStrain<threeD, threeD>;
} // muSpectre
diff --git a/src/fft/projection_finite_strain.hh b/src/fft/projection_finite_strain.hh
index 529037e..fcefe46 100644
--- a/src/fft/projection_finite_strain.hh
+++ b/src/fft/projection_finite_strain.hh
@@ -1,90 +1,91 @@
/**
* @file projection_finite_strain.hh
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 05 Dec 2017
*
* @brief Class for standard finite-strain gradient projections see de Geus et
* al. (https://doi.org/10.1016/j.cma.2016.12.032) for derivation
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef PROJECTION_FINITE_STRAIN_H
#define PROJECTION_FINITE_STRAIN_H
#include "fft/projection_default.hh"
#include "common/common.hh"
#include "common/field_collection.hh"
#include "common/field_map.hh"
namespace muSpectre {
/**
* Implements the finite strain gradient projection operator as
* defined in de Geus et
* al. (https://doi.org/10.1016/j.cma.2016.12.032) for derivation
*/
template <Dim_t DimS, Dim_t DimM>
class ProjectionFiniteStrain: public ProjectionDefault<DimS, DimM>
{
public:
using Parent = ProjectionDefault<DimS, DimM>; //!< base class
//! polymorphic pointer to FFT engines
using FFTEngine_ptr = typename Parent::FFTEngine_ptr;
using Ccoord = typename Parent::Ccoord; //!< cell coordinates type
+ using Rcoord = typename Parent::Rcoord; //!< spatial coordinates type
//! local field collection (for Fourier-space representations)
using LFieldCollection_t = LocalFieldCollection<DimS>;
//! iterable operator
using Proj_map = T4MatrixFieldMap<LFieldCollection_t, Real, DimM>;
//! iterable vectorised version of the Fourier-space tensor field
using Vector_map = MatrixFieldMap<LFieldCollection_t, Complex, DimM*DimM, 1>;
//! Default constructor
ProjectionFiniteStrain() = delete;
//! Constructor with fft_engine
- ProjectionFiniteStrain(FFTEngine_ptr engine);
+ ProjectionFiniteStrain(FFTEngine_ptr engine, Rcoord lengths);
//! Copy constructor
ProjectionFiniteStrain(const ProjectionFiniteStrain &other) = delete;
//! Move constructor
ProjectionFiniteStrain(ProjectionFiniteStrain &&other) = default;
//! Destructor
virtual ~ProjectionFiniteStrain() = default;
//! Copy assignment operator
ProjectionFiniteStrain& operator=(const ProjectionFiniteStrain &other) = delete;
//! Move assignment operator
ProjectionFiniteStrain& operator=(ProjectionFiniteStrain &&other) = default;
//! initialises the fft engine (plan the transform)
virtual void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate) override final;
protected:
private:
};
} // muSpectre
#endif /* PROJECTION_FINITE_STRAIN_H */
diff --git a/src/fft/projection_finite_strain_fast.cc b/src/fft/projection_finite_strain_fast.cc
index 08337cf..dd3300c 100644
--- a/src/fft/projection_finite_strain_fast.cc
+++ b/src/fft/projection_finite_strain_fast.cc
@@ -1,92 +1,92 @@
/**
* @file projection_finite_strain_fast.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 12 Dec 2017
*
* @brief implementation for fast projection in finite strain
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include "fft/projection_finite_strain_fast.hh"
#include "fft/fft_utils.hh"
#include "common/tensor_algebra.hh"
#include "common/iterators.hh"
namespace muSpectre {
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
ProjectionFiniteStrainFast<DimS, DimM>::
- ProjectionFiniteStrainFast(FFTEngine_ptr engine)
- :Parent{std::move(engine), Formulation::finite_strain},
+ ProjectionFiniteStrainFast(FFTEngine_ptr engine, Rcoord lengths)
+ :Parent{std::move(engine), lengths, Formulation::finite_strain},
xiField{make_field<Proj_t>("Projection Operator",
this->projection_container)},
xis(xiField)
{
for (auto res: this->fft_engine->get_domain_resolutions()) {
if (res % 2 == 0) {
throw ProjectionError
("Only an odd number of gridpoints in each direction is supported");
}
}
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void ProjectionFiniteStrainFast<DimS, DimM>::
initialise(FFT_PlanFlags flags) {
Parent::initialise(flags);
FFT_freqs<DimS> fft_freqs(this->fft_engine->get_domain_resolutions(),
- this->fft_engine->get_lengths());
+ this->domain_lengths);
for (auto && tup: akantu::zip(*this->fft_engine, this->xis)) {
const auto & ccoord = std::get<0> (tup);
auto & xi = std::get<1>(tup);
xi = fft_freqs.get_unit_xi(ccoord);
}
- if (this->get_locations() == Ccoord{}) {
+ if (this->get_subdomain_locations() == Ccoord{}) {
this->xis[0].setZero();
}
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void ProjectionFiniteStrainFast<DimS, DimM>::apply_projection(Field_t & field) {
Grad_map field_map{this->fft_engine->fft(field)};
Real factor = this->fft_engine->normalisation();
for (auto && tup: akantu::zip(this->xis, field_map)) {
auto & xi{std::get<0>(tup)};
auto & f{std::get<1>(tup)};
f = factor * ((f*xi).eval()*xi.transpose());
}
this->fft_engine->ifft(field);
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
Eigen::Map<Eigen::ArrayXXd> ProjectionFiniteStrainFast<DimS, DimM>::
get_operator() {
return this->xiField.dyn_eigen();
}
template class ProjectionFiniteStrainFast<twoD, twoD>;
template class ProjectionFiniteStrainFast<threeD, threeD>;
} // muSpectre
diff --git a/src/fft/projection_finite_strain_fast.hh b/src/fft/projection_finite_strain_fast.hh
index 1b99bbe..8a09407 100644
--- a/src/fft/projection_finite_strain_fast.hh
+++ b/src/fft/projection_finite_strain_fast.hh
@@ -1,105 +1,106 @@
/**
* @file projection_finite_strain_fast.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 12 Dec 2017
*
* @brief Faster alternative to ProjectionFinitestrain
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef PROJECTION_FINITE_STRAIN_FAST_H
#define PROJECTION_FINITE_STRAIN_FAST_H
#include "fft/projection_base.hh"
#include "common/common.hh"
#include "common/field_collection.hh"
#include "common/field_map.hh"
namespace muSpectre {
/**
* replaces `muSpectre::ProjectionFiniteStrain` with a faster and
* less memory-hungry alternative formulation. Use this if you don't
* have a very good reason not to (and tell me (author) about it,
* I'd be interested to hear it).
*/
template <Dim_t DimS, Dim_t DimM>
class ProjectionFiniteStrainFast: public ProjectionBase<DimS, DimM>
{
public:
using Parent = ProjectionBase<DimS, DimM>; //!< base class
//! polymorphic pointer to FFT engines
using FFTEngine_ptr = typename Parent::FFTEngine_ptr;
using Ccoord = typename Parent::Ccoord; //!< cell coordinates type
+ using Rcoord = typename Parent::Rcoord; //!< spatial coordinates type
//! global field collection (for real-space representations)
using GFieldCollection_t = GlobalFieldCollection<DimS>;
//! local field collection (for Fourier-space representations)
using LFieldCollection_t = LocalFieldCollection<DimS>;
//! Real space second order tensor fields (to be projected)
using Field_t = TensorField<GFieldCollection_t, Real, secondOrder, DimM>;
//! Fourier-space field containing the projection operator itself
using Proj_t = TensorField<LFieldCollection_t, Real, firstOrder, DimM>;
//! iterable form of the operator
using Proj_map = MatrixFieldMap<LFieldCollection_t, Real, DimM, 1>;
//! iterable Fourier-space second-order tensor field
using Grad_map = MatrixFieldMap<LFieldCollection_t, Complex, DimM, DimM>;
//! Default constructor
ProjectionFiniteStrainFast() = delete;
//! Constructor with fft_engine
- ProjectionFiniteStrainFast(FFTEngine_ptr engine);
+ ProjectionFiniteStrainFast(FFTEngine_ptr engine, Rcoord lengths);
//! Copy constructor
ProjectionFiniteStrainFast(const ProjectionFiniteStrainFast &other) = delete;
//! Move constructor
ProjectionFiniteStrainFast(ProjectionFiniteStrainFast &&other) = default;
//! Destructor
virtual ~ProjectionFiniteStrainFast() = default;
//! Copy assignment operator
ProjectionFiniteStrainFast& operator=(const ProjectionFiniteStrainFast &other) = delete;
//! Move assignment operator
ProjectionFiniteStrainFast& operator=(ProjectionFiniteStrainFast &&other) = default;
//! initialises the fft engine (plan the transform)
virtual void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate) override final;
//! apply the projection operator to a field
void apply_projection(Field_t & field) override final;
Eigen::Map<Eigen::ArrayXXd> get_operator() override final;
protected:
Proj_t & xiField; //!< field of normalised wave vectors
Proj_map xis; //!< iterable normalised wave vectors
private:
};
} // muSpectre
#endif /* PROJECTION_FINITE_STRAIN_FAST_H */
diff --git a/src/fft/projection_small_strain.cc b/src/fft/projection_small_strain.cc
index 7279fc6..9d86625 100644
--- a/src/fft/projection_small_strain.cc
+++ b/src/fft/projection_small_strain.cc
@@ -1,83 +1,83 @@
/**
* @file projection_small_strain.cc
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 14 Jan 2018
*
* @brief Implementation for ProjectionSmallStrain
*
* Copyright © 2018 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include "fft/projection_small_strain.hh"
#include "fft/fft_utils.hh"
namespace muSpectre {
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
ProjectionSmallStrain<DimS, DimM>::
- ProjectionSmallStrain(FFTEngine_ptr engine)
- : Parent{std::move(engine), Formulation::small_strain}
+ ProjectionSmallStrain(FFTEngine_ptr engine, Rcoord lengths)
+ : Parent{std::move(engine), lengths, Formulation::small_strain}
{
for (auto res: this->fft_engine->get_domain_resolutions()) {
if (res % 2 == 0) {
throw ProjectionError
("Only an odd number of gridpoints in each direction is supported");
}
}
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
void ProjectionSmallStrain<DimS, DimM>::initialise(FFT_PlanFlags flags) {
Parent::initialise(flags);
FFT_freqs<DimS> fft_freqs(this->fft_engine->get_domain_resolutions(),
- this->fft_engine->get_lengths());
+ this->domain_lengths);
for (auto && tup: akantu::zip(*this->fft_engine, this->Ghat)) {
const auto & ccoord = std::get<0> (tup);
auto & G = std::get<1>(tup);
auto xi = fft_freqs.get_unit_xi(ccoord);
auto kron = [](const Dim_t i, const Dim_t j) -> Real{
return (i==j) ? 1. : 0.;
};
for (Dim_t i{0}; i < DimS; ++i) {
for (Dim_t j{0}; j < DimS; ++j) {
for (Dim_t l{0}; l < DimS; ++l) {
for (Dim_t m{0}; m < DimS; ++m ) {
Real & g = get(G, i, j, l, m);
g = 0.5* (xi(i) * kron(j, l) * xi(m) +
xi(i) * kron(j, m) * xi(l) +
xi(j) * kron(i, l) * xi(m) +
xi(j) * kron(i, m) * xi(l)) -
xi(i)*xi(j)*xi(l)*xi(m);
}
}
}
}
}
- if (this->get_locations() == Ccoord{}) {
+ if (this->get_subdomain_locations() == Ccoord{}) {
this->Ghat[0].setZero();
}
}
template class ProjectionSmallStrain<twoD, twoD>;
template class ProjectionSmallStrain<threeD, threeD>;
} // muSpectre
diff --git a/src/fft/projection_small_strain.hh b/src/fft/projection_small_strain.hh
index 2ccb835..ba4e5f5 100644
--- a/src/fft/projection_small_strain.hh
+++ b/src/fft/projection_small_strain.hh
@@ -1,93 +1,94 @@
/**
* @file projection_small_strain.cc
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 14 Jan 2018
*
* @brief Small strain projection operator as defined in Appendix A1 of
* DOI: 10.1002/nme.5481 ("A finite element perspective on nonlinear
* FFT-based micromechanical simulations", Int. J. Numer. Meth. Engng
* 2017; 111 :903–926)
*
* Copyright © 2018 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef PROJECTION_SMALL_STRAIN_H
#define PROJECTION_SMALL_STRAIN_H
#include "fft/projection_default.hh"
namespace muSpectre {
/**
* Implements the small strain projection operator as defined in
* Appendix A1 of DOI: 10.1002/nme.5481 ("A finite element
* perspective on nonlinear FFT-based micromechanical
* simulations", Int. J. Numer. Meth. Engng 2017; 111
* :903–926)
*/
template <Dim_t DimS, Dim_t DimM>
class ProjectionSmallStrain: public ProjectionDefault<DimS, DimM>
{
public:
using Parent = ProjectionDefault<DimS, DimM>; //!< base class
//! polymorphic pointer to FFT engines
using FFTEngine_ptr = typename Parent::FFTEngine_ptr;
using Ccoord = typename Parent::Ccoord; //!< cell coordinates type
+ using Rcoord = typename Parent::Rcoord; //!< spatial coordinates type
//! local field collection (for Fourier-space representations)
using LFieldCollection_t = LocalFieldCollection<DimS>;
//! Fourier-space field containing the projection operator itself
using Proj_t = TensorField<LFieldCollection_t, Real, fourthOrder, DimM>;
//! iterable operator
using Proj_map = T4MatrixFieldMap<LFieldCollection_t, Real, DimM>;
//! iterable vectorised version of the Fourier-space tensor field
using Vector_map = MatrixFieldMap<LFieldCollection_t, Complex, DimM*DimM, 1>;
//! Default constructor
ProjectionSmallStrain() = delete;
//! Constructor with fft_engine
- ProjectionSmallStrain(FFTEngine_ptr engine);
+ ProjectionSmallStrain(FFTEngine_ptr engine, Rcoord lengths);
//! Copy constructor
ProjectionSmallStrain(const ProjectionSmallStrain &other) = delete;
//! Move constructor
ProjectionSmallStrain(ProjectionSmallStrain &&other) = default;
//! Destructor
virtual ~ProjectionSmallStrain() = default;
//! Copy assignment operator
ProjectionSmallStrain& operator=(const ProjectionSmallStrain &other) = delete;
//! Move assignment operator
ProjectionSmallStrain& operator=(ProjectionSmallStrain &&other) = delete;
//! initialises the fft engine (plan the transform)
virtual void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate) override final;
protected:
private:
};
} // muSpectre
#endif /* PROJECTION_SMALL_STRAIN_H */
diff --git a/src/solver/solver_base.hh b/src/solver/solver_base.hh
index 8a51f56..1d82b33 100644
--- a/src/solver/solver_base.hh
+++ b/src/solver/solver_base.hh
@@ -1,148 +1,148 @@
/**
* @file solver_base.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 18 Dec 2017
*
* @brief Base class for solvers
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef SOLVER_BASE_H
#define SOLVER_BASE_H
#include "solver/solver_error.hh"
#include "common/common.hh"
#include "cell/cell_base.hh"
#include "common/tensor_algebra.hh"
#include <Eigen/Dense>
#include <vector>
namespace muSpectre {
/* ---------------------------------------------------------------------- */
/**
* Virtual base class for solvers. Any implementation of this interface can be used with the solver functions prototyped in solvers.hh
*/
template <Dim_t DimS, Dim_t DimM=DimS>
class SolverBase
{
public:
/**
* Enum to describe in what kind the solver relies tangent stiffnesses
*/
enum class TangentRequirement{NoNeed, NeedEffect, NeedTangents};
using Cell_t = CellBase<DimS, DimM>; //!< Cell type
using Ccoord = Ccoord_t<DimS>; //!< cell coordinates type
//! Field collection to store temporary fields in
using Collection_t = GlobalFieldCollection<DimS>;
//! Input vector for solvers
using SolvVectorIn = Eigen::Ref<Eigen::VectorXd>;
//! Input vector for solvers
using SolvVectorInC = Eigen::Ref<const Eigen::VectorXd>;
//! Output vector for solvers
using SolvVectorOut = Eigen::VectorXd;
//! Default constructor
SolverBase() = delete;
//! Constructor with domain resolutions
SolverBase(Cell_t & cell, Real tol, Uint maxiter=0, bool verbose =false);
//! Copy constructor
SolverBase(const SolverBase &other) = delete;
//! Move constructor
SolverBase(SolverBase &&other) = default;
//! Destructor
virtual ~SolverBase() = default;
//! Copy assignment operator
SolverBase& operator=(const SolverBase &other) = delete;
//! Move assignment operator
SolverBase& operator=(SolverBase &&other) = default;
//! Allocate fields used during the solution
virtual void initialise() {
- this->collection.initialise(this->cell.get_resolutions(),
- this->cell.get_locations());
+ this->collection.initialise(this->cell.get_subdomain_resolutions(),
+ this->cell.get_subdomain_locations());
}
//! determine whether this solver requires full tangent stiffnesses
bool need_tangents() const {
return (this->get_tangent_req() == TangentRequirement::NeedTangents);}
//! determine whether this solver requires evaluation of directional tangent
bool need_effect() const {
return (this->get_tangent_req() == TangentRequirement::NeedEffect);}
//! determine whether this solver has no need for tangents
bool no_need_tangent() const {
return (this->get_tangent_req() == TangentRequirement::NoNeed);}
//! returns whether the solver has converged
virtual bool has_converged() const = 0;
//! reset the iteration counter to zero
void reset_counter();
//! get the count of how many solve steps have been executed since
//! construction of most recent counter reset
Uint get_counter() const;
//! executes the solver
virtual SolvVectorOut solve(const SolvVectorInC rhs, SolvVectorIn x_0) = 0;
//! return a reference to the cell
Cell_t & get_cell() {return cell;}
//! read the current maximum number of iterations setting
Uint get_maxiter() const {return this->maxiter;}
//! set the maximum number of iterations
void set_maxiter(Uint val) {this->maxiter = val;}
//! read the current tolerance setting
Real get_tol() const {return this->tol;}
//! set the torelance setting
void set_tol(Real val) {this->tol = val;}
//! returns the name of the solver
virtual std::string name() const = 0;
protected:
//! returns the tangent requirements of this solver
virtual TangentRequirement get_tangent_req() const = 0;
Cell_t & cell; //!< reference to the cell
Real tol; //!< convergence tolerance
Uint maxiter;//!< maximum number of iterations
bool verbose;//!< whether or not to write information to the std output
Uint counter{0}; //!< iteration counter
//! storage for internal fields to avoid reallocations between calls
Collection_t collection{};
private:
};
} // muSpectre
#endif /* SOLVER_BASE_H */
diff --git a/src/solver/solvers.cc b/src/solver/solvers.cc
index 954c3ba..94e3b60 100644
--- a/src/solver/solvers.cc
+++ b/src/solver/solvers.cc
@@ -1,430 +1,432 @@
/**
* @file solvers.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 20 Dec 2017
*
* @brief implementation of solver functions
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include "solver/solvers.hh"
#include "solver/solver_cg.hh"
#include "common/iterators.hh"
#include <Eigen/IterativeLinearSolvers>
#include <iomanip>
#include <cmath>
namespace muSpectre {
template <Dim_t DimS, Dim_t DimM>
std::vector<OptimizeResult>
de_geus (CellBase<DimS, DimM> & cell, const GradIncrements<DimM> & delFs,
SolverBase<DimS, DimM> & solver, Real newton_tol,
Real equil_tol,
Dim_t verbose) {
using Field_t = typename MaterialBase<DimS, DimM>::StrainField_t;
const Communicator & comm = cell.get_communicator();
auto solver_fields{std::make_unique<GlobalFieldCollection<DimS>>()};
- solver_fields->initialise(cell.get_resolutions(), cell.get_locations());
+ solver_fields->initialise(cell.get_subdomain_resolutions(),
+ cell.get_subdomain_locations());
// Corresponds to symbol δF or δε
auto & incrF{make_field<Field_t>("δF", *solver_fields)};
// Corresponds to symbol ΔF or Δε
auto & DeltaF{make_field<Field_t>("ΔF", *solver_fields)};
// field to store the rhs for cg calculations
auto & rhs{make_field<Field_t>("rhs", *solver_fields)};
solver.initialise();
if (solver.get_maxiter() == 0) {
solver.set_maxiter(cell.size()*DimM*DimM*10);
}
size_t count_width{};
const auto form{cell.get_formulation()};
std::string strain_symb{};
if (verbose > 0 && comm.rank() == 0) {
//setup of algorithm 5.2 in Nocedal, Numerical Optimization (p. 111)
std::cout << "de Geus-" << solver.name() << " for ";
switch (form) {
case Formulation::small_strain: {
strain_symb = "ε";
std::cout << "small";
break;
}
case Formulation::finite_strain: {
strain_symb = "F";
std::cout << "finite";
break;
}
default:
throw SolverError("unknown formulation");
break;
}
std::cout << " strain with" << std::endl
<< "newton_tol = " << newton_tol << ", cg_tol = "
<< solver.get_tol() << " maxiter = " << solver.get_maxiter() << " and Δ"
<< strain_symb << " =" <<std::endl;
for (auto&& tup: akantu::enumerate(delFs)) {
auto && counter{std::get<0>(tup)};
auto && grad{std::get<1>(tup)};
std::cout << "Step " << counter + 1 << ":" << std::endl
<< grad << std::endl;
}
count_width = size_t(std::log10(solver.get_maxiter()))+1;
}
// initialise F = I or ε = 0
auto & F{cell.get_strain()};
switch (form) {
case Formulation::finite_strain: {
F.get_map() = Matrices::I2<DimM>();
break;
}
case Formulation::small_strain: {
F.get_map() = Matrices::I2<DimM>().Zero();
for (const auto & delF: delFs) {
if (!check_symmetry(delF)) {
throw SolverError("all Δε must be symmetric!");
}
}
break;
}
default:
throw SolverError("Unknown formulation");
break;
}
// initialise return value
std::vector<OptimizeResult> ret_val{};
// initialise materials
constexpr bool need_tangent{true};
cell.initialise_materials(need_tangent);
Grad_t<DimM> previous_grad{Grad_t<DimM>::Zero()};
for (const auto & delF: delFs) { //incremental loop
std::string message{"Has not converged"};
Real incrNorm{2*newton_tol}, gradNorm{1};
Real stressNorm{2*equil_tol};
bool has_converged{false};
auto convergence_test = [&incrNorm, &gradNorm, &newton_tol,
&stressNorm, &equil_tol, &message, &has_converged] () {
bool incr_test = incrNorm/gradNorm <= newton_tol;
bool stress_test = stressNorm < equil_tol;
if (incr_test) {
message = "Residual tolerance reached";
} else if (stress_test) {
message = "Reached stress divergence tolerance";
}
has_converged = incr_test || stress_test;
return has_converged;
};
Uint newt_iter{0};
for (;
(newt_iter < solver.get_maxiter()) && (!has_converged ||
(newt_iter==1));
++newt_iter) {
// obtain material response
auto res_tup{cell.evaluate_stress_tangent(F)};
auto & P{std::get<0>(res_tup)};
auto & K{std::get<1>(res_tup)};
auto tangent_effect = [&cell, &K] (const Field_t & dF, Field_t & dP) {
cell.directional_stiffness(K, dF, dP);
};
if (newt_iter == 0) {
DeltaF.get_map() = -(delF-previous_grad); // neg sign because rhs
tangent_effect(DeltaF, rhs);
stressNorm = std::sqrt(comm.sum(rhs.eigen().matrix().squaredNorm()));
if (convergence_test()) {
break;
}
incrF.eigenvec() = solver.solve(rhs.eigenvec(), incrF.eigenvec());
F.eigen() -= DeltaF.eigen();
} else {
rhs.eigen() = -P.eigen();
cell.project(rhs);
stressNorm = std::sqrt(comm.sum(rhs.eigen().matrix().squaredNorm()));
if (convergence_test()) {
break;
}
incrF.eigen() = 0;
incrF.eigenvec() = solver.solve(rhs.eigenvec(), incrF.eigenvec());
}
F.eigen() += incrF.eigen();
incrNorm = std::sqrt(comm.sum(incrF.eigen().matrix().squaredNorm()));
gradNorm = std::sqrt(comm.sum(F.eigen().matrix().squaredNorm()));
if (verbose > 0 && comm.rank() == 0) {
std::cout << "at Newton step " << std::setw(count_width) << newt_iter
<< ", |δ" << strain_symb << "|/|Δ" << strain_symb
<< "| = " << std::setw(17) << incrNorm/gradNorm
<< ", tol = " << newton_tol << std::endl;
if (verbose-1>1) {
std::cout << "<" << strain_symb << "> =" << std::endl
<< F.get_map().mean() << std::endl;
}
}
convergence_test();
}
// update previous gradient
previous_grad = delF;
ret_val.push_back(OptimizeResult{F.eigen(), cell.get_stress().eigen(),
has_converged, Int(has_converged),
message,
newt_iter, solver.get_counter()});
// store history variables here
cell.save_history_variables();
}
return ret_val;
}
//! instantiation for two-dimensional cells
template std::vector<OptimizeResult>
de_geus (CellBase<twoD, twoD> & cell, const GradIncrements<twoD>& delF0,
SolverBase<twoD, twoD> & solver, Real newton_tol,
Real equil_tol,
Dim_t verbose);
// template typename CellBase<twoD, threeD>::StrainField_t &
// de_geus (CellBase<twoD, threeD> & cell, const GradIncrements<threeD>& delF0,
// const Real cg_tol, const Real newton_tol, Uint maxiter,
// Dim_t verbose);
//! instantiation for three-dimensional cells
template std::vector<OptimizeResult>
de_geus (CellBase<threeD, threeD> & cell, const GradIncrements<threeD>& delF0,
SolverBase<threeD, threeD> & solver, Real newton_tol,
Real equil_tol,
Dim_t verbose);
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
std::vector<OptimizeResult>
newton_cg (CellBase<DimS, DimM> & cell, const GradIncrements<DimM> & delFs,
SolverBase<DimS, DimM> & solver, Real newton_tol,
Real equil_tol,
Dim_t verbose) {
using Field_t = typename MaterialBase<DimS, DimM>::StrainField_t;
const Communicator & comm = cell.get_communicator();
auto solver_fields{std::make_unique<GlobalFieldCollection<DimS>>()};
- solver_fields->initialise(cell.get_resolutions(), cell.get_locations());
+ solver_fields->initialise(cell.get_subdomain_resolutions(),
+ cell.get_subdomain_locations());
// Corresponds to symbol δF or δε
auto & incrF{make_field<Field_t>("δF", *solver_fields)};
// field to store the rhs for cg calculations
auto & rhs{make_field<Field_t>("rhs", *solver_fields)};
solver.initialise();
if (solver.get_maxiter() == 0) {
solver.set_maxiter(cell.size()*DimM*DimM*10);
}
size_t count_width{};
const auto form{cell.get_formulation()};
std::string strain_symb{};
if (verbose > 0 && comm.rank() == 0) {
//setup of algorithm 5.2 in Nocedal, Numerical Optimization (p. 111)
std::cout << "Newton-" << solver.name() << " for ";
switch (form) {
case Formulation::small_strain: {
strain_symb = "ε";
std::cout << "small";
break;
}
case Formulation::finite_strain: {
strain_symb = "Fy";
std::cout << "finite";
break;
}
default:
throw SolverError("unknown formulation");
break;
}
std::cout << " strain with" << std::endl
<< "newton_tol = " << newton_tol << ", cg_tol = "
<< solver.get_tol() << " maxiter = " << solver.get_maxiter() << " and Δ"
<< strain_symb << " =" <<std::endl;
for (auto&& tup: akantu::enumerate(delFs)) {
auto && counter{std::get<0>(tup)};
auto && grad{std::get<1>(tup)};
std::cout << "Step " << counter + 1 << ":" << std::endl
<< grad << std::endl;
}
count_width = size_t(std::log10(solver.get_maxiter()))+1;
}
// initialise F = I or ε = 0
auto & F{cell.get_strain()};
switch (cell.get_formulation()) {
case Formulation::finite_strain: {
F.get_map() = Matrices::I2<DimM>();
break;
}
case Formulation::small_strain: {
F.get_map() = Matrices::I2<DimM>().Zero();
for (const auto & delF: delFs) {
if (!check_symmetry(delF)) {
throw SolverError("all Δε must be symmetric!");
}
}
break;
}
default:
throw SolverError("Unknown formulation");
break;
}
// initialise return value
std::vector<OptimizeResult> ret_val{};
// initialise materials
constexpr bool need_tangent{true};
cell.initialise_materials(need_tangent);
Grad_t<DimM> previous_grad{Grad_t<DimM>::Zero()};
for (const auto & delF: delFs) { //incremental loop
// apply macroscopic strain increment
for (auto && grad: F.get_map()) {
grad += delF - previous_grad;
}
std::string message{"Has not converged"};
Real incrNorm{2*newton_tol}, gradNorm{1};
Real stressNorm{2*equil_tol};
bool has_converged{false};
auto convergence_test = [&incrNorm, &gradNorm, &newton_tol,
&stressNorm, &equil_tol, &message, &has_converged] () {
bool incr_test = incrNorm/gradNorm <= newton_tol;
bool stress_test = stressNorm < equil_tol;
if (incr_test) {
message = "Residual tolerance reached";
} else if (stress_test) {
message = "Reached stress divergence tolerance";
}
has_converged = incr_test || stress_test;
return has_converged;
};
Uint newt_iter{0};
for (;
newt_iter < solver.get_maxiter() && !has_converged;
++newt_iter) {
// obtain material response
auto res_tup{cell.evaluate_stress_tangent(F)};
auto & P{std::get<0>(res_tup)};
rhs.eigen() = -P.eigen();
cell.project(rhs);
stressNorm = std::sqrt(comm.sum(rhs.eigen().matrix().squaredNorm()));
if (convergence_test()) {
break;
}
incrF.eigen() = 0;
incrF.eigenvec() = solver.solve(rhs.eigenvec(), incrF.eigenvec());
F.eigen() += incrF.eigen();
incrNorm = std::sqrt(comm.sum(incrF.eigen().matrix().squaredNorm()));
gradNorm = std::sqrt(comm.sum(F.eigen().matrix().squaredNorm()));
if (verbose > 0 && comm.rank() == 0) {
std::cout << "at Newton step " << std::setw(count_width) << newt_iter
<< ", |δ" << strain_symb << "|/|Δ" << strain_symb
<< "| = " << std::setw(17) << incrNorm/gradNorm
<< ", tol = " << newton_tol << std::endl;
if (verbose-1>1) {
std::cout << "<" << strain_symb << "> =" << std::endl
<< F.get_map().mean() << std::endl;
}
}
convergence_test();
}
// update previous gradient
previous_grad = delF;
ret_val.push_back(OptimizeResult{F.eigen(), cell.get_stress().eigen(),
convergence_test(), Int(convergence_test()),
message,
newt_iter, solver.get_counter()});
//store history variables for next load increment
cell.save_history_variables();
}
return ret_val;
}
//! instantiation for two-dimensional cells
template std::vector<OptimizeResult>
newton_cg (CellBase<twoD, twoD> & cell, const GradIncrements<twoD>& delF0,
SolverBase<twoD, twoD> & solver, Real newton_tol,
Real equil_tol,
Dim_t verbose);
// template typename CellBase<twoD, threeD>::StrainField_t &
// newton_cg (CellBase<twoD, threeD> & cell, const GradIncrements<threeD>& delF0,
// const Real cg_tol, const Real newton_tol, Uint maxiter,
// Dim_t verbose);
//! instantiation for three-dimensional cells
template std::vector<OptimizeResult>
newton_cg (CellBase<threeD, threeD> & cell, const GradIncrements<threeD>& delF0,
SolverBase<threeD, threeD> & solver, Real newton_tol,
Real equil_tol,
Dim_t verbose);
/* ---------------------------------------------------------------------- */
bool check_symmetry(const Eigen::Ref<const Eigen::ArrayXXd>& eps,
Real rel_tol){
return (rel_tol >= (eps-eps.transpose()).matrix().norm()/eps.matrix().norm() ||
rel_tol >= eps.matrix().norm());
}
} // muSpectre
diff --git a/tests/mpi_test_fft_engine.cc b/tests/mpi_test_fft_engine.cc
index 21d1f56..962a11e 100644
--- a/tests/mpi_test_fft_engine.cc
+++ b/tests/mpi_test_fft_engine.cc
@@ -1,174 +1,175 @@
/**
* @file mpi_test_fft_engine.cc
*
* @author Lars Pastewka <lars.pastewka@imtek.uni-freiburg.de>
*
* @date 06 Mar 2017
*
* @brief tests for MPI-parallel fft engine implementations
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#define BOOST_MPL_CFG_NO_PREPROCESSED_HEADERS
#define BOOST_MPL_LIMIT_LIST_SIZE 50
#include <boost/mpl/list.hpp>
#include "tests.hh"
#include "mpi_context.hh"
#include "fft/fftw_engine.hh"
#ifdef WITH_FFTWMPI
#include "fft/fftwmpi_engine.hh"
#endif
#ifdef WITH_PFFT
#include "fft/pfft_engine.hh"
#endif
#include "common/ccoord_operations.hh"
#include "common/field_collection.hh"
#include "common/field_map.hh"
#include "common/iterators.hh"
namespace muSpectre {
BOOST_AUTO_TEST_SUITE(mpi_fft_engine);
/* ---------------------------------------------------------------------- */
template <typename Engine, Dim_t resolution, bool serial=false>
struct FFTW_fixture {
constexpr static Dim_t box_resolution{resolution};
constexpr static Dim_t serial_engine{serial};
constexpr static Real box_length{4.5};
constexpr static Dim_t sdim{Engine::sdim};
constexpr static Dim_t mdim{Engine::mdim};
constexpr static Ccoord_t<sdim> res() {
return CcoordOps::get_cube<sdim>(box_resolution);
}
- FFTW_fixture(): engine(res(), CcoordOps::get_cube<sdim>(box_length),
- MPIContext::get_context().comm) {}
+ FFTW_fixture(): engine(res(), MPIContext::get_context().comm) {}
Engine engine;
};
template <typename Engine>
struct FFTW_fixture_python_segfault{
constexpr static Dim_t serial_engine{false};
constexpr static Dim_t dim{twoD};
constexpr static Dim_t sdim{twoD};
constexpr static Dim_t mdim{twoD};
constexpr static Ccoord_t<sdim> res() {return {6, 4};}
FFTW_fixture_python_segfault():
- engine{res(), {3., 3}, MPIContext::get_context().comm} {}
+ engine{res(), MPIContext::get_context().comm} {}
Engine engine;
};
using fixlist = boost::mpl::list<
#ifdef WITH_FFTWMPI
FFTW_fixture<FFTWMPIEngine< twoD, twoD>, 3>,
FFTW_fixture<FFTWMPIEngine< twoD, threeD>, 3>,
FFTW_fixture<FFTWMPIEngine<threeD, threeD>, 3>,
FFTW_fixture<FFTWMPIEngine< twoD, twoD>, 4>,
FFTW_fixture<FFTWMPIEngine< twoD, threeD>, 4>,
FFTW_fixture<FFTWMPIEngine<threeD, threeD>, 4>,
FFTW_fixture_python_segfault<FFTWMPIEngine<twoD, twoD>>,
#endif
#ifdef WITH_PFFT
FFTW_fixture<PFFTEngine< twoD, twoD>, 3>,
FFTW_fixture<PFFTEngine< twoD, threeD>, 3>,
FFTW_fixture<PFFTEngine<threeD, threeD>, 3>,
FFTW_fixture<PFFTEngine< twoD, twoD>, 4>,
FFTW_fixture<PFFTEngine< twoD, threeD>, 4>,
FFTW_fixture<PFFTEngine<threeD, threeD>, 4>,
FFTW_fixture_python_segfault<PFFTEngine<twoD, twoD>>,
#endif
FFTW_fixture<FFTWEngine< twoD, twoD>, 3, true>>;
/* ---------------------------------------------------------------------- */
BOOST_FIXTURE_TEST_CASE_TEMPLATE(Constructor_test, Fix, fixlist, Fix) {
Communicator &comm = MPIContext::get_context().comm;
if (Fix::serial_engine && comm.size() > 1) {
return;
}
else {
BOOST_CHECK_NO_THROW(Fix::engine.initialise(FFT_PlanFlags::estimate));
}
BOOST_CHECK_EQUAL(comm.sum(Fix::engine.size()),
CcoordOps::get_size(Fix::res()));
}
/* ---------------------------------------------------------------------- */
BOOST_FIXTURE_TEST_CASE_TEMPLATE(fft_test, Fix, fixlist, Fix) {
if (Fix::serial_engine && Fix::engine.get_communicator().size() > 1) {
// dont test serial engies in parallel
return;
}
else {
Fix::engine.initialise(FFT_PlanFlags::estimate);
}
constexpr Dim_t order{2};
using FC_t = GlobalFieldCollection<Fix::sdim>;
FC_t fc;
auto & input{make_field<TensorField<FC_t, Real, order, Fix::mdim>>("input", fc)};
auto & ref {make_field<TensorField<FC_t, Real, order, Fix::mdim>>("reference", fc)};
auto & result{make_field<TensorField<FC_t, Real, order, Fix::mdim>>("result", fc)};
- fc.initialise(Fix::engine.get_resolutions(), Fix::engine.get_locations());
+ fc.initialise(Fix::engine.get_subdomain_resolutions(),
+ Fix::engine.get_subdomain_locations());
using map_t = MatrixFieldMap<FC_t, Real, Fix::mdim, Fix::mdim>;
map_t inmap{input};
auto refmap{map_t{ref}};
auto resultmap{map_t{result}};
size_t cntr{0};
for (auto tup: akantu::zip(inmap, refmap)) {
cntr++;
auto & in_{std::get<0>(tup)};
auto & ref_{std::get<1>(tup)};
in_.setRandom();
ref_ = in_;
}
auto & complex_field = Fix::engine.fft(input);
using cmap_t = MatrixFieldMap<LocalFieldCollection<Fix::sdim>, Complex, Fix::mdim, Fix::mdim>;
cmap_t complex_map(complex_field);
- if (Fix::engine.get_locations() == CcoordOps::get_cube<Fix::sdim>(0)) {
+ if (Fix::engine.get_subdomain_locations() ==
+ CcoordOps::get_cube<Fix::sdim>(0)) {
// Check that 0,0 location has no imaginary part.
Real error = complex_map[0].imag().norm();
BOOST_CHECK_LT(error, tol);
}
/* make sure, the engine has not modified input (which is
unfortunately const-casted internally, hence this test) */
for (auto && tup: akantu::zip(inmap, refmap)) {
Real error{(std::get<0>(tup) - std::get<1>(tup)).norm()};
BOOST_CHECK_LT(error, tol);
}
/* make sure that the ifft of fft returns the original*/
Fix::engine.ifft(result);
for (auto && tup: akantu::zip(resultmap, refmap)) {
Real error{(std::get<0>(tup)*Fix::engine.normalisation() - std::get<1>(tup)).norm()};
BOOST_CHECK_LT(error, tol);
if (error > tol) {
std::cout << std::get<0>(tup).array()/std::get<1>(tup).array() << std::endl << std::endl;
}
}
}
BOOST_AUTO_TEST_SUITE_END();
} // muSpectre
diff --git a/tests/mpi_test_projection.hh b/tests/mpi_test_projection.hh
index 438d39e..dae7c39 100644
--- a/tests/mpi_test_projection.hh
+++ b/tests/mpi_test_projection.hh
@@ -1,88 +1,88 @@
/**
* @file test_projection.hh
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 16 Jan 2018
*
* @brief common declarations for testing both the small and finite strain
* projection operators
*
* Copyright © 2018 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include "tests.hh"
#include "mpi_context.hh"
#include <boost/mpl/list.hpp>
#include <Eigen/Dense>
namespace muSpectre {
/* ---------------------------------------------------------------------- */
template <Dim_t DimS>
struct Sizes {
};
template<>
struct Sizes<twoD> {
constexpr static Ccoord_t<twoD> get_resolution() {
return Ccoord_t<twoD>{3, 5};}
constexpr static Rcoord_t<twoD> get_lengths() {
return Rcoord_t<twoD>{3.4, 5.8};}
};
template<>
struct Sizes<threeD> {
constexpr static Ccoord_t<threeD> get_resolution() {
return Ccoord_t<threeD>{3, 5, 7};}
constexpr static Rcoord_t<threeD> get_lengths() {
return Rcoord_t<threeD>{3.4, 5.8, 6.7};}
};
template <Dim_t DimS>
struct Squares {
};
template<>
struct Squares<twoD> {
constexpr static Ccoord_t<twoD> get_resolution() {
return Ccoord_t<twoD>{5, 5};}
constexpr static Rcoord_t<twoD> get_lengths() {
return Rcoord_t<twoD>{5, 5};}
};
template<>
struct Squares<threeD> {
constexpr static Ccoord_t<threeD> get_resolution() {
return Ccoord_t<threeD>{7, 7, 7};}
constexpr static Rcoord_t<threeD> get_lengths() {
return Rcoord_t<threeD>{7, 7, 7};}
};
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM, class SizeGiver, class Proj, class Engine,
bool parallel=true>
struct ProjectionFixture {
using Parent = Proj;
constexpr static Dim_t sdim{DimS};
constexpr static Dim_t mdim{DimM};
constexpr static bool is_parallel{parallel};
ProjectionFixture()
:projector(std::make_unique<Engine>(SizeGiver::get_resolution(),
- SizeGiver::get_lengths(),
- MPIContext::get_context().comm)){}
+ MPIContext::get_context().comm),
+ SizeGiver::get_lengths()){}
Parent projector;
};
} // muSpectre
diff --git a/tests/mpi_test_projection_finite.cc b/tests/mpi_test_projection_finite.cc
index 5811186..6a12061 100644
--- a/tests/mpi_test_projection_finite.cc
+++ b/tests/mpi_test_projection_finite.cc
@@ -1,185 +1,185 @@
/**
* @file test_projection_finite.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 07 Dec 2017
*
* @brief tests for standard finite strain projection operator
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#define BOOST_MPL_CFG_NO_PREPROCESSED_HEADERS
#define BOOST_MPL_LIMIT_LIST_SIZE 50
#include "fft/projection_finite_strain.hh"
#include "fft/projection_finite_strain_fast.hh"
#include "fft/fft_utils.hh"
#include "mpi_test_projection.hh"
#include "fft/fftw_engine.hh"
#ifdef WITH_FFTWMPI
#include "fft/fftwmpi_engine.hh"
#endif
#ifdef WITH_PFFT
#include "fft/pfft_engine.hh"
#endif
#include <Eigen/Dense>
namespace muSpectre {
BOOST_AUTO_TEST_SUITE(mpi_projection_finite_strain);
/* ---------------------------------------------------------------------- */
using fixlist = boost::mpl::list<
#ifdef WITH_FFTWMPI
ProjectionFixture<twoD, twoD, Squares<twoD>,
ProjectionFiniteStrain<twoD, twoD>,
FFTWMPIEngine<twoD, twoD>>,
ProjectionFixture<threeD, threeD, Squares<threeD>,
ProjectionFiniteStrain<threeD, threeD>,
FFTWMPIEngine<threeD, threeD>>,
ProjectionFixture<twoD, twoD, Sizes<twoD>,
ProjectionFiniteStrain<twoD, twoD>,
FFTWMPIEngine<twoD, twoD>>,
ProjectionFixture<threeD, threeD, Sizes<threeD>,
ProjectionFiniteStrain<threeD, threeD>,
FFTWMPIEngine<threeD, threeD>>,
ProjectionFixture<twoD, twoD, Squares<twoD>,
ProjectionFiniteStrainFast<twoD, twoD>,
FFTWMPIEngine<twoD, twoD>>,
ProjectionFixture<threeD, threeD, Squares<threeD>,
ProjectionFiniteStrainFast<threeD, threeD>,
FFTWMPIEngine<threeD, threeD>>,
ProjectionFixture<twoD, twoD, Sizes<twoD>,
ProjectionFiniteStrainFast<twoD, twoD>,
FFTWMPIEngine<twoD, twoD>>,
ProjectionFixture<threeD, threeD, Sizes<threeD>,
ProjectionFiniteStrainFast<threeD, threeD>,
FFTWMPIEngine<threeD, threeD>>,
#endif
#ifdef WITH_PFFT
ProjectionFixture<twoD, twoD, Squares<twoD>,
ProjectionFiniteStrain<twoD, twoD>,
PFFTEngine<twoD, twoD>>,
ProjectionFixture<threeD, threeD, Squares<threeD>,
ProjectionFiniteStrain<threeD, threeD>,
PFFTEngine<threeD, threeD>>,
ProjectionFixture<twoD, twoD, Sizes<twoD>,
ProjectionFiniteStrain<twoD, twoD>,
PFFTEngine<twoD, twoD>>,
ProjectionFixture<threeD, threeD, Sizes<threeD>,
ProjectionFiniteStrain<threeD, threeD>,
PFFTEngine<threeD, threeD>>,
ProjectionFixture<twoD, twoD, Squares<twoD>,
ProjectionFiniteStrainFast<twoD, twoD>,
PFFTEngine<twoD, twoD>>,
ProjectionFixture<threeD, threeD, Squares<threeD>,
ProjectionFiniteStrainFast<threeD, threeD>,
PFFTEngine<threeD, threeD>>,
ProjectionFixture<twoD, twoD, Sizes<twoD>,
ProjectionFiniteStrainFast<twoD, twoD>,
PFFTEngine<twoD, twoD>>,
ProjectionFixture<threeD, threeD, Sizes<threeD>,
ProjectionFiniteStrainFast<threeD, threeD>,
PFFTEngine<threeD, threeD>>,
#endif
ProjectionFixture<twoD, twoD, Squares<twoD>,
ProjectionFiniteStrain<twoD, twoD>,
FFTWEngine<twoD, twoD>,
false>
>;
/* ---------------------------------------------------------------------- */
BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, fix, fixlist, fix) {
if (fix::is_parallel || fix::projector.get_communicator().size() == 1) {
BOOST_CHECK_NO_THROW(fix::projector.initialise(FFT_PlanFlags::estimate));
}
}
/* ---------------------------------------------------------------------- */
BOOST_FIXTURE_TEST_CASE_TEMPLATE(Gradient_preservation_test,
fix, fixlist, fix) {
if (!fix::is_parallel || fix::projector.get_communicator().size() > 1) {
return;
}
// create a gradient field with a zero mean gradient and verify
// that the projection preserves it
constexpr Dim_t dim{fix::sdim}, sdim{fix::sdim}, mdim{fix::mdim};
static_assert(dim == fix::mdim,
"These tests assume that the material and spatial dimension are "
"identical");
using Fields = GlobalFieldCollection<sdim>;
using FieldT = TensorField<Fields, Real, secondOrder, mdim>;
using FieldMap = MatrixFieldMap<Fields, Real, mdim, mdim>;
using Vector = Eigen::Matrix<Real, dim, 1>;
Fields fields{};
FieldT & f_grad{make_field<FieldT>("gradient", fields)};
FieldT & f_var{make_field<FieldT>("working field", fields)};
FieldMap grad(f_grad);
FieldMap var(f_var);
- fields.initialise(fix::projector.get_resolutions(),
- fix::projector.get_locations());
+ fields.initialise(fix::projector.get_subdomain_resolutions(),
+ fix::projector.get_subdomain_locations());
Vector k; for (Dim_t i = 0; i < dim; ++i) {
// the wave vector has to be such that it leads to an integer
// number of periods in each length of the domain
- k(i) = (i+1)*2*pi/fix::projector.get_lengths()[i];
+ k(i) = (i+1)*2*pi/fix::projector.get_domain_lengths()[i];
}
for (auto && tup: akantu::zip(fields, grad, var)) {
auto & ccoord = std::get<0>(tup);
auto & g = std::get<1>(tup);
auto & v = std::get<2>(tup);
Vector vec = CcoordOps::get_vector(ccoord,
- fix::projector.get_lengths()/
+ fix::projector.get_domain_lengths()/
fix::projector.get_domain_resolutions());
g.row(0) = k.transpose() * cos(k.dot(vec));
v.row(0) = g.row(0);
}
fix::projector.initialise(FFT_PlanFlags::estimate);
fix::projector.apply_projection(f_var);
for (auto && tup: akantu::zip(fields, grad, var)) {
auto & ccoord = std::get<0>(tup);
auto & g = std::get<1>(tup);
auto & v = std::get<2>(tup);
Real error = (g-v).norm();
BOOST_CHECK_LT(error, tol);
if (error >=tol) {
Vector vec = CcoordOps::get_vector(ccoord,
- fix::projector.get_lengths()/
+ fix::projector.get_domain_lengths()/
fix::projector.get_domain_resolutions());
std::cout << std::endl << "grad_ref :" << std::endl << g << std::endl;
std::cout << std::endl << "grad_proj :" << std::endl << v << std::endl;
std::cout << std::endl << "ccoord :" << std::endl << ccoord << std::endl;
std::cout << std::endl << "vector :" << std::endl << vec.transpose() << std::endl;
}
}
}
BOOST_AUTO_TEST_SUITE_END();
} // muSpectre
diff --git a/tests/mpi_test_projection_small.cc b/tests/mpi_test_projection_small.cc
index b745f5e..a87ea52 100644
--- a/tests/mpi_test_projection_small.cc
+++ b/tests/mpi_test_projection_small.cc
@@ -1,170 +1,170 @@
/**
* @file test_projection_small.cc
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 16 Jan 2018
*
* @brief tests for standard small strain projection operator
*
* Copyright © 2018 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#define BOOST_MPL_CFG_NO_PREPROCESSED_HEADERS
#define BOOST_MPL_LIMIT_LIST_SIZE 50
#include "fft/projection_small_strain.hh"
#include "mpi_test_projection.hh"
#include "fft/fft_utils.hh"
#include "fft/fftw_engine.hh"
#ifdef WITH_FFTWMPI
#include "fft/fftwmpi_engine.hh"
#endif
#ifdef WITH_PFFT
#include "fft/pfft_engine.hh"
#endif
#include <Eigen/Dense>
namespace muSpectre {
BOOST_AUTO_TEST_SUITE(mpi_projection_small_strain);
using fixlist = boost::mpl::list<
#ifdef WITH_FFTWMPI
ProjectionFixture<twoD, twoD, Squares<twoD>,
ProjectionSmallStrain<twoD, twoD>,
FFTWMPIEngine<twoD, twoD>>,
ProjectionFixture<threeD, threeD, Squares<threeD>,
ProjectionSmallStrain<threeD, threeD>,
FFTWMPIEngine<threeD, threeD>>,
ProjectionFixture<twoD, twoD, Sizes<twoD>,
ProjectionSmallStrain<twoD, twoD>,
FFTWMPIEngine<twoD, twoD>>,
ProjectionFixture<threeD, threeD, Sizes<threeD>,
ProjectionSmallStrain<threeD, threeD>,
FFTWMPIEngine<threeD, threeD>>,
#endif
#ifdef WITH_PFFT
ProjectionFixture<twoD, twoD, Squares<twoD>,
ProjectionSmallStrain<twoD, twoD>,
PFFTEngine<twoD, twoD>>,
ProjectionFixture<threeD, threeD, Squares<threeD>,
ProjectionSmallStrain<threeD, threeD>,
PFFTEngine<threeD, threeD>>,
ProjectionFixture<twoD, twoD, Sizes<twoD>,
ProjectionSmallStrain<twoD, twoD>,
PFFTEngine<twoD, twoD>>,
ProjectionFixture<threeD, threeD, Sizes<threeD>,
ProjectionSmallStrain<threeD, threeD>,
PFFTEngine<threeD, threeD>>,
#endif
ProjectionFixture<twoD, twoD, Squares<twoD>,
ProjectionSmallStrain<twoD, twoD>,
FFTWEngine<twoD, twoD>,
false>
>;
/* ---------------------------------------------------------------------- */
BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, fix, fixlist, fix) {
if (fix::is_parallel || fix::projector.get_communicator().size() == 1) {
BOOST_CHECK_NO_THROW(fix::projector.initialise(FFT_PlanFlags::estimate));
}
}
/* ---------------------------------------------------------------------- */
BOOST_FIXTURE_TEST_CASE_TEMPLATE(Gradient_preservation_test,
fix, fixlist, fix) {
if (!fix::is_parallel || fix::projector.get_communicator().size() > 1) {
return;
}
// create a gradient field with a zero mean gradient and verify
// that the projection preserves it
constexpr Dim_t dim{fix::sdim}, sdim{fix::sdim}, mdim{fix::mdim};
static_assert(dim == fix::mdim,
"These tests assume that the material and spatial dimension are "
"identical");
using Fields = GlobalFieldCollection<sdim>;
using FieldT = TensorField<Fields, Real, secondOrder, mdim>;
using FieldMap = MatrixFieldMap<Fields, Real, mdim, mdim>;
using Vector = Eigen::Matrix<Real, dim, 1>;
Fields fields{};
FieldT & f_grad{make_field<FieldT>("strain", fields)};
FieldT & f_var{make_field<FieldT>("working field", fields)};
FieldMap grad(f_grad);
FieldMap var(f_var);
- fields.initialise(fix::projector.get_resolutions(),
- fix::projector.get_locations());
+ fields.initialise(fix::projector.get_subdomain_resolutions(),
+ fix::projector.get_subdomain_locations());
Vector k;
for (Dim_t i = 0; i < dim; ++i) {
// the wave vector has to be such that it leads to an integer
// number of periods in each length of the domain
- k(i) = (i+1)*2*pi/fix::projector.get_lengths()[i];
+ k(i) = (i+1)*2*pi/fix::projector.get_domain_lengths()[i];
}
for (auto && tup: akantu::zip(fields, grad, var)) {
auto & ccoord = std::get<0>(tup);
auto & g = std::get<1>(tup);
auto & v = std::get<2>(tup);
Vector vec = CcoordOps::get_vector(ccoord,
- fix::projector.get_lengths()/
- fix::projector.get_resolutions());
+ fix::projector.get_domain_lengths()/
+ fix::projector.get_domain_resolutions());
g.row(0) << k.transpose() * cos(k.dot(vec));
// We need to add I to the term, because this field has a net
// zero gradient, which leads to a net -I strain
g = 0.5*((g-g.Identity()).transpose() + (g-g.Identity())).eval()+g.Identity();
v = g;
}
fix::projector.initialise(FFT_PlanFlags::estimate);
fix::projector.apply_projection(f_var);
constexpr bool verbose{false};
for (auto && tup: akantu::zip(fields, grad, var)) {
auto & ccoord = std::get<0>(tup);
auto & g = std::get<1>(tup);
auto & v = std::get<2>(tup);
Vector vec = CcoordOps::get_vector(ccoord,
- fix::projector.get_lengths()/
- fix::projector.get_resolutions());
+ fix::projector.get_domain_lengths()/
+ fix::projector.get_domain_resolutions());
Real error = (g-v).norm();
BOOST_CHECK_LT(error, tol);
if ((error >=tol) || verbose) {
std::cout << std::endl << "grad_ref :" << std::endl << g << std::endl;
std::cout << std::endl << "grad_proj :" << std::endl << v << std::endl;
std::cout << std::endl << "ccoord :" << std::endl << ccoord << std::endl;
std::cout << std::endl << "vector :" << std::endl << vec.transpose() << std::endl;
std::cout << "means:" << std::endl
<< "<strain>:" << std::endl
<< grad.mean() << std::endl
<< "<proj>:" << std::endl
<< var.mean();
}
}
}
BOOST_AUTO_TEST_SUITE_END();
} // muSpectre
diff --git a/tests/mpi_test_solver_newton_cg.cc b/tests/mpi_test_solver_newton_cg.cc
index 0c3f8d4..e7ea9c9 100644
--- a/tests/mpi_test_solver_newton_cg.cc
+++ b/tests/mpi_test_solver_newton_cg.cc
@@ -1,201 +1,200 @@
/**
* @file test_solver_newton_cg.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 20 Dec 2017
*
* @brief Tests for the standard Newton-Raphson + Conjugate Gradient solver
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include "tests.hh"
#include "mpi_context.hh"
#include "solver/solvers.hh"
#include "solver/solver_cg.hh"
#include "solver/solver_cg_eigen.hh"
#include "fft/fftwmpi_engine.hh"
#include "fft/projection_finite_strain_fast.hh"
#include "materials/material_linear_elastic1.hh"
#include "common/iterators.hh"
#include "common/ccoord_operations.hh"
#include "cell/cell_factory.hh"
namespace muSpectre {
BOOST_AUTO_TEST_SUITE(newton_cg_tests);
BOOST_AUTO_TEST_CASE(manual_construction_test) {
const Communicator & comm = MPIContext::get_context().comm;
// constexpr Dim_t dim{twoD};
constexpr Dim_t dim{threeD};
// constexpr Ccoord_t<dim> resolutions{3, 3};
// constexpr Rcoord_t<dim> lengths{2.3, 2.7};
constexpr Ccoord_t<dim> resolutions{5, 5, 5};
constexpr Rcoord_t<dim> lengths{5, 5, 5};
- auto fft_ptr{std::make_unique<FFTWMPIEngine<dim, dim>>(resolutions, lengths,
- comm)};
- auto proj_ptr{std::make_unique<ProjectionFiniteStrainFast<dim, dim>>(std::move(fft_ptr))};
+ auto fft_ptr{std::make_unique<FFTWMPIEngine<dim, dim>>(resolutions, comm)};
+ auto proj_ptr{std::make_unique<ProjectionFiniteStrainFast<dim, dim>>(std::move(fft_ptr), lengths)};
CellBase<dim, dim> sys(std::move(proj_ptr));
using Mat_t = MaterialLinearElastic1<dim, dim>;
//const Real Young{210e9}, Poisson{.33};
const Real Young{1.0030648180242636}, Poisson{0.29930675909878679};
// const Real lambda{Young*Poisson/((1+Poisson)*(1-2*Poisson))};
// const Real mu{Young/(2*(1+Poisson))};
auto& Material_hard = Mat_t::make(sys, "hard", 10*Young, Poisson);
auto& Material_soft = Mat_t::make(sys, "soft", Young, Poisson);
- auto& loc = sys.get_locations();
+ auto& loc = sys.get_subdomain_locations();
for (auto && tup: akantu::enumerate(sys)) {
auto && pixel = std::get<1>(tup);
if (loc == Ccoord_t<threeD>{0, 0} && std::get<0>(tup) == 0) {
Material_hard.add_pixel(pixel);
} else {
Material_soft.add_pixel(pixel);
}
}
sys.initialise();
Grad_t<dim> delF0;
delF0 << 0, 1., 0, 0, 0, 0, 0, 0, 0;
constexpr Real cg_tol{1e-8}, newton_tol{1e-5};
constexpr Uint maxiter{CcoordOps::get_size(resolutions)*ipow(dim, secondOrder)*10};
constexpr bool verbose{false};
GradIncrements<dim> grads; grads.push_back(delF0);
SolverCG<dim> cg{sys, cg_tol, maxiter, bool(verbose)};
Eigen::ArrayXXd res1{de_geus(sys, grads, cg, newton_tol, verbose)[0].grad};
SolverCG<dim> cg2{sys, cg_tol, maxiter, bool(verbose)};
Eigen::ArrayXXd res2{newton_cg(sys, grads, cg2, newton_tol, verbose)[0].grad};
BOOST_CHECK_LE(abs(res1-res2).mean(), cg_tol);
}
BOOST_AUTO_TEST_CASE(small_strain_patch_test) {
const Communicator & comm = MPIContext::get_context().comm;
constexpr Dim_t dim{twoD};
using Ccoord = Ccoord_t<dim>;
using Rcoord = Rcoord_t<dim>;
constexpr Ccoord resolutions{CcoordOps::get_cube<dim>(3)};
constexpr Rcoord lengths{CcoordOps::get_cube<dim>(1.)};
constexpr Formulation form{Formulation::small_strain};
// number of layers in the hard material
constexpr Uint nb_lays{1};
constexpr Real contrast{2};
static_assert(nb_lays < resolutions[0],
"the number or layers in the hard material must be smaller "
"than the total number of layers in dimension 0");
auto sys{make_parallel_cell(resolutions, lengths, form, comm)};
using Mat_t = MaterialLinearElastic1<dim, dim>;
constexpr Real Young{2.}, Poisson{.33};
auto material_hard{std::make_unique<Mat_t>("hard", contrast*Young, Poisson)};
auto material_soft{std::make_unique<Mat_t>("soft", Young, Poisson)};
for (const auto & pixel: sys) {
if (pixel[0] < Dim_t(nb_lays)) {
material_hard->add_pixel(pixel);
} else {
material_soft->add_pixel(pixel);
}
}
sys.add_material(std::move(material_hard));
sys.add_material(std::move(material_soft));
sys.initialise();
Grad_t<dim> delEps0{Grad_t<dim>::Zero()};
constexpr Real eps0 = 1.;
//delEps0(0, 1) = delEps0(1, 0) = eps0;
delEps0(0, 0) = eps0;
constexpr Real cg_tol{1e-8}, newton_tol{1e-5}, equil_tol{1e-10};
constexpr Uint maxiter{dim*10};
constexpr Dim_t verbose{0};
SolverCGEigen<dim> cg{sys, cg_tol, maxiter, bool(verbose)};
auto result = de_geus(sys, delEps0, cg, newton_tol,
equil_tol, verbose);
if (verbose) {
std::cout << "result:" << std::endl << result.grad << std::endl;
std::cout << "mean strain = " << std::endl
<< sys.get_strain().get_map().mean() << std::endl;
}
/**
* verification of resultant strains: subscript ₕ for hard and ₛ
* for soft, Nₕ is nb_lays and Nₜₒₜ is resolutions, k is contrast
*
* Δl = εl = Δlₕ + Δlₛ = εₕlₕ+εₛlₛ
* => ε = εₕ Nₕ/Nₜₒₜ + εₛ (Nₜₒₜ-Nₕ)/Nₜₒₜ
*
* σ is constant across all layers
* σₕ = σₛ
* => Eₕ εₕ = Eₛ εₛ
* => εₕ = 1/k εₛ
* => ε / (1/k Nₕ/Nₜₒₜ + (Nₜₒₜ-Nₕ)/Nₜₒₜ) = εₛ
*/
constexpr Real factor{1/contrast * Real(nb_lays)/resolutions[0]
+ 1.-nb_lays/Real(resolutions[0])};
constexpr Real eps_soft{eps0/factor};
constexpr Real eps_hard{eps_soft/contrast};
if (verbose) {
std::cout << "εₕ = " << eps_hard << ", εₛ = " << eps_soft << std::endl;
std::cout << "ε = εₕ Nₕ/Nₜₒₜ + εₛ (Nₜₒₜ-Nₕ)/Nₜₒₜ" << std::endl;
}
Grad_t<dim> Eps_hard; Eps_hard << eps_hard, 0, 0, 0;
Grad_t<dim> Eps_soft; Eps_soft << eps_soft, 0, 0, 0;
// verify uniaxial tension patch test
for (const auto & pixel: sys) {
if (pixel[0] < Dim_t(nb_lays)) {
BOOST_CHECK_LE((Eps_hard-sys.get_strain().get_map()[pixel]).norm(), tol);
} else {
BOOST_CHECK_LE((Eps_soft-sys.get_strain().get_map()[pixel]).norm(), tol);
}
}
delEps0 = Grad_t<dim>::Zero();
delEps0(0, 1) = delEps0(1, 0) = eps0;
SolverCG<dim> cg2{sys, cg_tol, maxiter, bool(verbose)};
result = newton_cg(sys, delEps0, cg2, newton_tol,
equil_tol, verbose);
Eps_hard << 0, eps_hard, eps_hard, 0;
Eps_soft << 0, eps_soft, eps_soft, 0;
// verify pure shear patch test
for (const auto & pixel: sys) {
if (pixel[0] < Dim_t(nb_lays)) {
BOOST_CHECK_LE((Eps_hard-sys.get_strain().get_map()[pixel]).norm(), tol);
} else {
BOOST_CHECK_LE((Eps_soft-sys.get_strain().get_map()[pixel]).norm(), tol);
}
}
}
BOOST_AUTO_TEST_SUITE_END();
} // muSpectre
diff --git a/tests/python_fft_tests.py b/tests/python_fft_tests.py
index 924bb5b..de93da2 100644
--- a/tests/python_fft_tests.py
+++ b/tests/python_fft_tests.py
@@ -1,100 +1,97 @@
#!/usr/bin/env python3
# -*- coding:utf-8 -*-
"""
file python_fft_tests.py
@author Till Junge <till.junge@altermail.ch>
@date 17 Jan 2018
@brief Compare µSpectre's fft implementations to numpy reference
@section LICENSE
Copyright © 2018 Till Junge
µSpectre 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, or (at
your option) any later version.
µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
Free Software Foundation, Inc., 59 Temple Place - Suite 330,
Boston, MA 02111-1307, USA.
"""
import unittest
import numpy as np
from python_test_imports import µ
try:
from mpi4py import MPI
comm = MPI.COMM_WORLD
except ImportError:
comm = None
class FFT_Check(unittest.TestCase):
def setUp(self):
self.resolution = [6, 4]
self.dim = len(self.resolution)
- self.lengths = [3., 3.]
self.engines = [('fftw', False),
('fftwmpi', True),
('pfft', True)]
self.tol = 1e-14 * np.prod(self.resolution)
def test_forward_transform(self):
for engine_str, transposed in self.engines:
try:
- engine = µ.fft.FFT(self.resolution, self.lengths,
- fft=engine_str)
+ engine = µ.fft.FFT(self.resolution, fft=engine_str)
except KeyError:
# This FFT engine has not been compiled into the code. Skip
# test.
continue
engine.initialise()
in_arr = np.random.random([*self.resolution, self.dim, self.dim])
out_ref = np.fft.rfftn(in_arr, axes=(0, 1))
if transposed:
out_ref = out_ref.swapaxes(0, 1)
out_msp = engine.fft(in_arr.reshape(-1, self.dim**2).T).T
err = np.linalg.norm(out_ref -
out_msp.reshape(out_ref.shape))
self.assertTrue(err < self.tol)
def test_reverse_transform(self):
for engine_str, transposed in self.engines:
try:
- engine = µ.fft.FFT(self.resolution, self.lengths,
- fft=engine_str)
+ engine = µ.fft.FFT(self.resolution, fft=engine_str)
except KeyError:
# This FFT engine has not been compiled into the code. Skip
# test.
continue
engine.initialise()
complex_res = µ.get_hermitian_sizes(self.resolution)
in_arr = np.zeros([*complex_res, self.dim, self.dim],
dtype=complex)
in_arr.real = np.random.random(in_arr.shape)
in_arr.imag = np.random.random(in_arr.shape)
out_ref = np.fft.irfftn(in_arr, axes=(0, 1))
if transposed:
in_arr = in_arr.swapaxes(0, 1)
out_msp = engine.ifft(in_arr.reshape(-1, self.dim**2).T).T
out_msp *= engine.normalisation()
err = np.linalg.norm(out_ref -
out_msp.reshape(out_ref.shape))
self.assertTrue(err < self.tol)
if __name__ == '__main__':
unittest.main()
diff --git a/tests/test_cell_base.cc b/tests/test_cell_base.cc
index ae204b2..0dcc84d 100644
--- a/tests/test_cell_base.cc
+++ b/tests/test_cell_base.cc
@@ -1,194 +1,194 @@
/**
* @file test_cell_base.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 14 Dec 2017
*
* @brief Tests for the basic cell class
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include <boost/mpl/list.hpp>
#include <Eigen/Dense>
#include "tests.hh"
#include "common/common.hh"
#include "common/iterators.hh"
#include "common/field_map.hh"
#include "tests/test_goodies.hh"
#include "cell/cell_factory.hh"
#include "materials/material_linear_elastic1.hh"
namespace muSpectre {
BOOST_AUTO_TEST_SUITE(cell_base);
template <Dim_t DimS>
struct Sizes {
};
template<>
struct Sizes<twoD> {
constexpr static Dim_t sdim{twoD};
constexpr static Ccoord_t<sdim> get_resolution() {
return Ccoord_t<sdim>{3, 5};}
constexpr static Rcoord_t<sdim> get_lengths() {
return Rcoord_t<sdim>{3.4, 5.8};}
};
template<>
struct Sizes<threeD> {
constexpr static Dim_t sdim{threeD};
constexpr static Ccoord_t<sdim> get_resolution() {
return Ccoord_t<sdim>{3, 5, 7};}
constexpr static Rcoord_t<sdim> get_lengths() {
return Rcoord_t<sdim>{3.4, 5.8, 6.7};}
};
template <Dim_t DimS, Dim_t DimM, Formulation form>
struct CellBaseFixture: CellBase<DimS, DimM> {
constexpr static Dim_t sdim{DimS};
constexpr static Dim_t mdim{DimM};
constexpr static Formulation formulation{form};
CellBaseFixture()
:CellBase<DimS, DimM>{
std::move(cell_input<DimS, DimM>(Sizes<DimS>::get_resolution(),
Sizes<DimS>::get_lengths(),
form))} {}
};
using fixlist = boost::mpl::list<CellBaseFixture<twoD, twoD,
Formulation::finite_strain>,
CellBaseFixture<threeD, threeD,
Formulation::finite_strain>,
CellBaseFixture<twoD, twoD,
Formulation::small_strain>,
CellBaseFixture<threeD, threeD,
Formulation::small_strain>>;
BOOST_AUTO_TEST_CASE(manual_construction) {
constexpr Dim_t dim{twoD};
Ccoord_t<dim> resolutions{3, 3};
Rcoord_t<dim> lengths{2.3, 2.7};
Formulation form{Formulation::finite_strain};
- auto fft_ptr{std::make_unique<FFTWEngine<dim, dim>>(resolutions, lengths)};
- auto proj_ptr{std::make_unique<ProjectionFiniteStrainFast<dim, dim>>(std::move(fft_ptr))};
+ auto fft_ptr{std::make_unique<FFTWEngine<dim, dim>>(resolutions)};
+ auto proj_ptr{std::make_unique<ProjectionFiniteStrainFast<dim, dim>>(std::move(fft_ptr), lengths)};
CellBase<dim, dim> sys{std::move(proj_ptr)};
auto sys2{make_cell<dim, dim>(resolutions, lengths, form)};
auto sys2b{std::move(sys2)};
BOOST_CHECK_EQUAL(sys2b.size(), sys.size());
}
BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, fix, fixlist, fix) {
BOOST_CHECK_THROW(fix::check_material_coverage(), std::runtime_error);
BOOST_CHECK_THROW(fix::initialise(), std::runtime_error);
}
BOOST_FIXTURE_TEST_CASE_TEMPLATE(add_material_test, fix, fixlist, fix) {
constexpr Dim_t dim{fix::sdim};
using Material_t = MaterialLinearElastic1<dim, dim>;
auto Material_hard = std::make_unique<Material_t>("hard", 210e9, .33);
BOOST_CHECK_NO_THROW(fix::add_material(std::move(Material_hard)));
}
BOOST_FIXTURE_TEST_CASE_TEMPLATE(simple_evaluation_test, fix, fixlist, fix) {
constexpr Dim_t dim{fix::sdim};
constexpr Formulation form{fix::formulation};
using Mat_t = MaterialLinearElastic1<dim, dim>;
const Real Young{210e9}, Poisson{.33};
const Real lambda{Young*Poisson/((1+Poisson)*(1-2*Poisson))};
const Real mu{Young/(2*(1+Poisson))};
auto Material_hard = std::make_unique<Mat_t>("hard", Young, Poisson);
for (auto && pixel: *this) {
Material_hard->add_pixel(pixel);
}
fix::add_material(std::move(Material_hard));
auto & F = fix::get_strain();
auto F_map = F.get_map();
// finite strain formulation expects the deformation gradient F,
// while small strain expects infinitesimal strain ε
for (auto grad: F_map) {
switch (form) {
case Formulation::finite_strain: {
grad = grad.Identity();
break;
}
case Formulation::small_strain: {
grad = grad.Zero();
break;
}
default:
BOOST_CHECK(false);
break;
}
}
fix::initialise_materials();
auto res_tup{fix::evaluate_stress_tangent(F)};
auto stress{std::get<0>(res_tup).get_map()};
auto tangent{std::get<1>(res_tup).get_map()};
auto tup = testGoodies::objective_hooke_explicit
(lambda, mu, Matrices::I2<dim>());
auto P_ref = std::get<0>(tup);
for (auto mat: stress) {
Real norm = (mat - P_ref).norm();
BOOST_CHECK_EQUAL(norm, 0.);
}
auto tan_ref = std::get<1>(tup);
for (const auto tan: tangent) {
Real norm = (tan - tan_ref).norm();
BOOST_CHECK_EQUAL(norm, 0.);
}
}
BOOST_FIXTURE_TEST_CASE_TEMPLATE(evaluation_test, fix, fixlist, fix) {
constexpr Dim_t dim{fix::sdim};
using Mat_t = MaterialLinearElastic1<dim, dim>;
auto Material_hard = std::make_unique<Mat_t>("hard", 210e9, .33);
auto Material_soft = std::make_unique<Mat_t>("soft", 70e9, .3);
for (auto && cnt_pixel: akantu::enumerate(*this)) {
auto counter = std::get<0>(cnt_pixel);
auto && pixel = std::get<1>(cnt_pixel);
if (counter < 5) {
Material_hard->add_pixel(pixel);
} else {
Material_soft->add_pixel(pixel);
}
}
fix::add_material(std::move(Material_hard));
fix::add_material(std::move(Material_soft));
auto & F = fix::get_strain();
fix::initialise_materials();
fix::evaluate_stress_tangent(F);
fix::evaluate_stress_tangent(F);
}
BOOST_AUTO_TEST_SUITE_END();
} // muSpectre
diff --git a/tests/test_fftw_engine.cc b/tests/test_fftw_engine.cc
index cafa4a4..3fdb161 100644
--- a/tests/test_fftw_engine.cc
+++ b/tests/test_fftw_engine.cc
@@ -1,133 +1,132 @@
/**
* @file test_fftw_engine.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 05 Dec 2017
*
* @brief tests for the fftw fft engine implementation
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include <boost/mpl/list.hpp>
#include "tests.hh"
#include "fft/fftw_engine.hh"
#include "common/ccoord_operations.hh"
#include "common/field_collection.hh"
#include "common/field_map.hh"
#include "common/iterators.hh"
namespace muSpectre {
BOOST_AUTO_TEST_SUITE(fftw_engine);
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM, Dim_t resolution>
struct FFTW_fixture {
constexpr static Dim_t box_resolution{resolution};
constexpr static Real box_length{4.5};
constexpr static Dim_t sdim{DimS};
constexpr static Dim_t mdim{DimM};
constexpr static Ccoord_t<sdim> res() {
return CcoordOps::get_cube<DimS>(box_resolution);
}
constexpr static Ccoord_t<sdim> loc() {
return CcoordOps::get_cube<DimS>(0);
}
- FFTW_fixture() :engine(res(),
- CcoordOps::get_cube<DimS>(box_length)){}
+ FFTW_fixture() : engine(res()) {}
FFTWEngine<DimS, DimM> engine;
};
struct FFTW_fixture_python_segfault{
constexpr static Dim_t dim{twoD};
constexpr static Dim_t sdim{twoD};
constexpr static Dim_t mdim{twoD};
constexpr static Ccoord_t<sdim> res() {return {6, 4};}
constexpr static Ccoord_t<sdim> loc() {return {0, 0};}
- FFTW_fixture_python_segfault():engine{res(), {3., 3}} {}
+ FFTW_fixture_python_segfault() : engine{res()} {}
FFTWEngine<sdim, mdim> engine;
};
using fixlist = boost::mpl::list<FFTW_fixture< twoD, twoD, 3>,
FFTW_fixture< twoD, threeD, 3>,
FFTW_fixture<threeD, threeD, 3>,
FFTW_fixture< twoD, twoD, 4>,
FFTW_fixture< twoD, threeD, 4>,
FFTW_fixture<threeD, threeD, 4>,
FFTW_fixture_python_segfault>;
/* ---------------------------------------------------------------------- */
BOOST_FIXTURE_TEST_CASE_TEMPLATE(Constructor_test, Fix, fixlist, Fix) {
BOOST_CHECK_NO_THROW(Fix::engine.initialise(FFT_PlanFlags::estimate));
BOOST_CHECK_EQUAL(Fix::engine.size(), CcoordOps::get_size(Fix::res()));
}
/* ---------------------------------------------------------------------- */
BOOST_FIXTURE_TEST_CASE_TEMPLATE(fft_test, Fix, fixlist, Fix) {
Fix::engine.initialise(FFT_PlanFlags::estimate);
constexpr Dim_t order{2};
using FC_t = GlobalFieldCollection<Fix::sdim>;
FC_t fc;
auto & input{make_field<TensorField<FC_t, Real, order, Fix::mdim>>("input", fc)};
auto & ref {make_field<TensorField<FC_t, Real, order, Fix::mdim>>("reference", fc)};
auto & result{make_field<TensorField<FC_t, Real, order, Fix::mdim>>("result", fc)};
fc.initialise(Fix::res(), Fix::loc());
using map_t = MatrixFieldMap<FC_t, Real, Fix::mdim, Fix::mdim>;
map_t inmap{input};
auto refmap{map_t{ref}};
auto resultmap{map_t{result}};
size_t cntr{0};
for (auto tup: akantu::zip(inmap, refmap)) {
cntr++;
auto & in_{std::get<0>(tup)};
auto & ref_{std::get<1>(tup)};
in_.setRandom();
ref_ = in_;
}
auto & complex_field = Fix::engine.fft(input);
using cmap_t = MatrixFieldMap<LocalFieldCollection<Fix::sdim>, Complex, Fix::mdim, Fix::mdim>;
cmap_t complex_map(complex_field);
Real error = complex_map[0].imag().norm();
BOOST_CHECK_LT(error, tol);
/* make sure, the engine has not modified input (which is
unfortunately const-casted internally, hence this test) */
for (auto && tup: akantu::zip(inmap, refmap)) {
Real error{(std::get<0>(tup) - std::get<1>(tup)).norm()};
BOOST_CHECK_LT(error, tol);
}
/* make sure that the ifft of fft returns the original*/
Fix::engine.ifft(result);
for (auto && tup: akantu::zip(resultmap, refmap)) {
Real error{(std::get<0>(tup)*Fix::engine.normalisation() - std::get<1>(tup)).norm()};
BOOST_CHECK_LT(error, tol);
if (error > tol) {
std::cout << std::get<0>(tup).array()/std::get<1>(tup).array() << std::endl << std::endl;
}
}
}
BOOST_AUTO_TEST_SUITE_END();
} // muSpectre
diff --git a/tests/test_projection.hh b/tests/test_projection.hh
index 0b2004a..a6371e2 100644
--- a/tests/test_projection.hh
+++ b/tests/test_projection.hh
@@ -1,86 +1,86 @@
/**
* @file test_projection.hh
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 16 Jan 2018
*
* @brief common declarations for testing both the small and finite strain
* projection operators
*
* Copyright © 2018 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include "tests.hh"
#include "fft/fftw_engine.hh"
#include <boost/mpl/list.hpp>
#include <Eigen/Dense>
namespace muSpectre {
/* ---------------------------------------------------------------------- */
template <Dim_t DimS>
struct Sizes {
};
template<>
struct Sizes<twoD> {
constexpr static Ccoord_t<twoD> get_resolution() {
return Ccoord_t<twoD>{3, 5};}
constexpr static Rcoord_t<twoD> get_lengths() {
return Rcoord_t<twoD>{3.4, 5.8};}
};
template<>
struct Sizes<threeD> {
constexpr static Ccoord_t<threeD> get_resolution() {
return Ccoord_t<threeD>{3, 5, 7};}
constexpr static Rcoord_t<threeD> get_lengths() {
return Rcoord_t<threeD>{3.4, 5.8, 6.7};}
};
template <Dim_t DimS>
struct Squares {
};
template<>
struct Squares<twoD> {
constexpr static Ccoord_t<twoD> get_resolution() {
return Ccoord_t<twoD>{5, 5};}
constexpr static Rcoord_t<twoD> get_lengths() {
return Rcoord_t<twoD>{5, 5};}
};
template<>
struct Squares<threeD> {
constexpr static Ccoord_t<threeD> get_resolution() {
return Ccoord_t<threeD>{7, 7, 7};}
constexpr static Rcoord_t<threeD> get_lengths() {
return Rcoord_t<threeD>{7, 7, 7};}
};
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM, class SizeGiver, class Proj>
struct ProjectionFixture {
using Engine = FFTWEngine<DimS, DimM>;
using Parent = Proj;
constexpr static Dim_t sdim{DimS};
constexpr static Dim_t mdim{DimM};
ProjectionFixture()
- :projector(std::make_unique<Engine>(SizeGiver::get_resolution(),
- SizeGiver::get_lengths())){}
+ : projector(std::make_unique<Engine>(SizeGiver::get_resolution()),
+ SizeGiver::get_lengths()) {}
Parent projector;
};
} // muSpectre
diff --git a/tests/test_projection_finite.cc b/tests/test_projection_finite.cc
index 60700af..916f0b3 100644
--- a/tests/test_projection_finite.cc
+++ b/tests/test_projection_finite.cc
@@ -1,138 +1,137 @@
/**
* @file test_projection_finite.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 07 Dec 2017
*
* @brief tests for standard finite strain projection operator
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include "fft/projection_finite_strain.hh"
#include "fft/projection_finite_strain_fast.hh"
#include "fft/fft_utils.hh"
#include "test_projection.hh"
#include <Eigen/Dense>
namespace muSpectre {
BOOST_AUTO_TEST_SUITE(projection_finite_strain);
/* ---------------------------------------------------------------------- */
using fixlist = boost::mpl::list<
ProjectionFixture<twoD, twoD, Squares<twoD>,
ProjectionFiniteStrain<twoD, twoD>>,
ProjectionFixture<threeD, threeD, Squares<threeD>,
ProjectionFiniteStrain<threeD, threeD>>,
ProjectionFixture<twoD, twoD, Sizes<twoD>,
ProjectionFiniteStrain<twoD, twoD>>,
ProjectionFixture<threeD, threeD, Sizes<threeD>,
ProjectionFiniteStrain<threeD, threeD>>,
ProjectionFixture<twoD, twoD, Squares<twoD>,
ProjectionFiniteStrainFast<twoD, twoD>>,
ProjectionFixture<threeD, threeD, Squares<threeD>,
ProjectionFiniteStrainFast<threeD, threeD>>,
ProjectionFixture<twoD, twoD, Sizes<twoD>,
ProjectionFiniteStrainFast<twoD, twoD>>,
ProjectionFixture<threeD, threeD, Sizes<threeD>,
ProjectionFiniteStrainFast<threeD, threeD>>>;
/* ---------------------------------------------------------------------- */
BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, fix, fixlist, fix) {
BOOST_CHECK_NO_THROW(fix::projector.initialise(FFT_PlanFlags::estimate));
}
/* ---------------------------------------------------------------------- */
BOOST_AUTO_TEST_CASE(even_grid_test) {
using Engine = FFTWEngine<twoD, twoD>;
using proj = ProjectionFiniteStrainFast<twoD, twoD>;
- auto engine = std::make_unique<Engine>(Ccoord_t<twoD>{2, 2},
- Rcoord_t<twoD>{4.3, 4.3});
- BOOST_CHECK_THROW(proj (std::move(engine)),
- std::runtime_error);
+ auto engine = std::make_unique<Engine>(Ccoord_t<twoD>{2, 2});
+ BOOST_CHECK_THROW(proj(std::move(engine), Rcoord_t<twoD>{4.3, 4.3}),
+ std::runtime_error);
}
/* ---------------------------------------------------------------------- */
BOOST_FIXTURE_TEST_CASE_TEMPLATE(Gradient_preservation_test,
fix, fixlist, fix) {
// create a gradient field with a zero mean gradient and verify
// that the projection preserves it
constexpr Dim_t dim{fix::sdim}, sdim{fix::sdim}, mdim{fix::mdim};
static_assert(dim == fix::mdim,
"These tests assume that the material and spatial dimension are "
"identical");
using Fields = GlobalFieldCollection<sdim>;
using FieldT = TensorField<Fields, Real, secondOrder, mdim>;
using FieldMap = MatrixFieldMap<Fields, Real, mdim, mdim>;
using Vector = Eigen::Matrix<Real, dim, 1>;
Fields fields{};
FieldT & f_grad{make_field<FieldT>("gradient", fields)};
FieldT & f_var{make_field<FieldT>("working field", fields)};
FieldMap grad(f_grad);
FieldMap var(f_var);
- fields.initialise(fix::projector.get_resolutions(),
- fix::projector.get_locations());
+ fields.initialise(fix::projector.get_subdomain_resolutions(),
+ fix::projector.get_subdomain_locations());
FFT_freqs<dim> freqs{fix::projector.get_domain_resolutions(),
- fix::projector.get_lengths()};
+ fix::projector.get_domain_lengths()};
Vector k; for (Dim_t i = 0; i < dim; ++i) {
// the wave vector has to be such that it leads to an integer
// number of periods in each length of the domain
- k(i) = (i+1)*2*pi/fix::projector.get_lengths()[i]; ;
+ k(i) = (i+1)*2*pi/fix::projector.get_domain_lengths()[i]; ;
}
for (auto && tup: akantu::zip(fields, grad, var)) {
auto & ccoord = std::get<0>(tup);
auto & g = std::get<1>(tup);
auto & v = std::get<2>(tup);
Vector vec = CcoordOps::get_vector(ccoord,
- fix::projector.get_lengths()/
- fix::projector.get_resolutions());
+ fix::projector.get_domain_lengths()/
+ fix::projector.get_domain_resolutions());
g.row(0) = k.transpose() * cos(k.dot(vec));
v.row(0) = g.row(0);
}
fix::projector.initialise(FFT_PlanFlags::estimate);
fix::projector.apply_projection(f_var);
for (auto && tup: akantu::zip(fields, grad, var)) {
auto & ccoord = std::get<0>(tup);
auto & g = std::get<1>(tup);
auto & v = std::get<2>(tup);
Vector vec = CcoordOps::get_vector(ccoord,
- fix::projector.get_lengths()/
- fix::projector.get_resolutions());
+ fix::projector.get_domain_lengths()/
+ fix::projector.get_domain_resolutions());
Real error = (g-v).norm();
BOOST_CHECK_LT(error, tol);
if (error >=tol) {
std::cout << std::endl << "grad_ref :" << std::endl << g << std::endl;
std::cout << std::endl << "grad_proj :" << std::endl << v << std::endl;
std::cout << std::endl << "ccoord :" << std::endl << ccoord << std::endl;
std::cout << std::endl << "vector :" << std::endl << vec.transpose() << std::endl;
}
}
}
BOOST_AUTO_TEST_SUITE_END();
} // muSpectre
diff --git a/tests/test_projection_small.cc b/tests/test_projection_small.cc
index 60e4e27..f8ae988 100644
--- a/tests/test_projection_small.cc
+++ b/tests/test_projection_small.cc
@@ -1,130 +1,130 @@
/**
* @file test_projection_small.cc
*
* @author Till Junge <till.junge@altermail.ch>
*
* @date 16 Jan 2018
*
* @brief tests for standard small strain projection operator
*
* Copyright © 2018 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include "fft/projection_small_strain.hh"
#include "test_projection.hh"
#include "fft/fft_utils.hh"
#include <Eigen/Dense>
namespace muSpectre {
BOOST_AUTO_TEST_SUITE(projection_small_strain);
using fixlist = boost::mpl::list<
ProjectionFixture<twoD, twoD, Squares<twoD>,
ProjectionSmallStrain<twoD, twoD>>,
ProjectionFixture<threeD, threeD, Squares<threeD>,
ProjectionSmallStrain<threeD, threeD>>,
ProjectionFixture<twoD, twoD, Sizes<twoD>,
ProjectionSmallStrain<twoD, twoD>>,
ProjectionFixture<threeD, threeD, Sizes<threeD>,
ProjectionSmallStrain<threeD, threeD>>>;
/* ---------------------------------------------------------------------- */
BOOST_FIXTURE_TEST_CASE_TEMPLATE(constructor_test, fix, fixlist, fix) {
BOOST_CHECK_NO_THROW(fix::projector.initialise(FFT_PlanFlags::estimate));
}
/* ---------------------------------------------------------------------- */
BOOST_FIXTURE_TEST_CASE_TEMPLATE(Gradient_preservation_test,
fix, fixlist, fix) {
// create a gradient field with a zero mean gradient and verify
// that the projection preserves it
constexpr Dim_t dim{fix::sdim}, sdim{fix::sdim}, mdim{fix::mdim};
static_assert(dim == fix::mdim,
"These tests assume that the material and spatial dimension are "
"identical");
using Fields = GlobalFieldCollection<sdim>;
using FieldT = TensorField<Fields, Real, secondOrder, mdim>;
using FieldMap = MatrixFieldMap<Fields, Real, mdim, mdim>;
using Vector = Eigen::Matrix<Real, dim, 1>;
Fields fields{};
FieldT & f_grad{make_field<FieldT>("strain", fields)};
FieldT & f_var{make_field<FieldT>("working field", fields)};
FieldMap grad(f_grad);
FieldMap var(f_var);
- fields.initialise(fix::projector.get_resolutions(),
- fix::projector.get_locations());
- FFT_freqs<dim> freqs{fix::projector.get_resolutions(),
- fix::projector.get_lengths()};
+ fields.initialise(fix::projector.get_subdomain_resolutions(),
+ fix::projector.get_subdomain_locations());
+ FFT_freqs<dim> freqs{fix::projector.get_domain_resolutions(),
+ fix::projector.get_domain_lengths()};
Vector k;
for (Dim_t i = 0; i < dim; ++i) {
// the wave vector has to be such that it leads to an integer
// number of periods in each length of the domain
- k(i) = (i+1)*2*pi/fix::projector.get_lengths()[i];
+ k(i) = (i+1)*2*pi/fix::projector.get_domain_lengths()[i];
}
for (auto && tup: akantu::zip(fields, grad, var)) {
auto & ccoord = std::get<0>(tup);
auto & g = std::get<1>(tup);
auto & v = std::get<2>(tup);
Vector vec = CcoordOps::get_vector(ccoord,
- fix::projector.get_lengths()/
- fix::projector.get_resolutions());
+ fix::projector.get_domain_lengths()/
+ fix::projector.get_domain_resolutions());
g.row(0) << k.transpose() * cos(k.dot(vec));
// We need to add I to the term, because this field has a net
// zero gradient, which leads to a net -I strain
g = 0.5*((g-g.Identity()).transpose() + (g-g.Identity())).eval()+g.Identity();
v = g;
}
fix::projector.initialise(FFT_PlanFlags::estimate);
fix::projector.apply_projection(f_var);
constexpr bool verbose{false};
for (auto && tup: akantu::zip(fields, grad, var)) {
auto & ccoord = std::get<0>(tup);
auto & g = std::get<1>(tup);
auto & v = std::get<2>(tup);
Vector vec = CcoordOps::get_vector(ccoord,
- fix::projector.get_lengths()/
- fix::projector.get_resolutions());
+ fix::projector.get_domain_lengths()/
+ fix::projector.get_domain_resolutions());
Real error = (g-v).norm();
BOOST_CHECK_LT(error, tol);
if ((error >=tol) || verbose) {
std::cout << std::endl << "grad_ref :" << std::endl << g << std::endl;
std::cout << std::endl << "grad_proj :" << std::endl << v << std::endl;
std::cout << std::endl << "ccoord :" << std::endl << ccoord << std::endl;
std::cout << std::endl << "vector :" << std::endl << vec.transpose() << std::endl;
std::cout << "means:" << std::endl
<< "<strain>:" << std::endl
<< grad.mean() << std::endl
<< "<proj>:" << std::endl
<< var.mean();
}
}
}
BOOST_AUTO_TEST_SUITE_END();
} // muSpectre
diff --git a/tests/test_solver_newton_cg.cc b/tests/test_solver_newton_cg.cc
index 92e1749..76c228f 100644
--- a/tests/test_solver_newton_cg.cc
+++ b/tests/test_solver_newton_cg.cc
@@ -1,195 +1,195 @@
/**
* @file test_solver_newton_cg.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 20 Dec 2017
*
* @brief Tests for the standard Newton-Raphson + Conjugate Gradient solver
*
* Copyright © 2017 Till Junge
*
* µSpectre 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, or (at
* your option) any later version.
*
* µSpectre 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include "tests.hh"
#include "solver/solvers.hh"
#include "solver/solver_cg.hh"
#include "solver/solver_cg_eigen.hh"
#include "fft/fftw_engine.hh"
#include "fft/projection_finite_strain_fast.hh"
#include "materials/material_linear_elastic1.hh"
#include "common/iterators.hh"
#include "common/ccoord_operations.hh"
#include "cell/cell_factory.hh"
namespace muSpectre {
BOOST_AUTO_TEST_SUITE(newton_cg_tests);
BOOST_AUTO_TEST_CASE(manual_construction_test) {
// constexpr Dim_t dim{twoD};
constexpr Dim_t dim{threeD};
// constexpr Ccoord_t<dim> resolutions{3, 3};
// constexpr Rcoord_t<dim> lengths{2.3, 2.7};
constexpr Ccoord_t<dim> resolutions{5, 5, 5};
constexpr Rcoord_t<dim> lengths{5, 5, 5};
- auto fft_ptr{std::make_unique<FFTWEngine<dim, dim>>(resolutions, lengths)};
- auto proj_ptr{std::make_unique<ProjectionFiniteStrainFast<dim, dim>>(std::move(fft_ptr))};
+ auto fft_ptr{std::make_unique<FFTWEngine<dim, dim>>(resolutions)};
+ auto proj_ptr{std::make_unique<ProjectionFiniteStrainFast<dim, dim>>(std::move(fft_ptr), lengths)};
CellBase<dim, dim> sys(std::move(proj_ptr));
using Mat_t = MaterialLinearElastic1<dim, dim>;
//const Real Young{210e9}, Poisson{.33};
const Real Young{1.0030648180242636}, Poisson{0.29930675909878679};
// const Real lambda{Young*Poisson/((1+Poisson)*(1-2*Poisson))};
// const Real mu{Young/(2*(1+Poisson))};
auto& Material_hard = Mat_t::make(sys, "hard", 10*Young, Poisson);
auto& Material_soft = Mat_t::make(sys, "soft", Young, Poisson);
for (auto && tup: akantu::enumerate(sys)) {
auto && pixel = std::get<1>(tup);
if (std::get<0>(tup) == 0) {
Material_hard.add_pixel(pixel);
} else {
Material_soft.add_pixel(pixel);
}
}
sys.initialise();
Grad_t<dim> delF0;
delF0 << 0, 1., 0, 0, 0, 0, 0, 0, 0;
constexpr Real cg_tol{1e-8}, newton_tol{1e-5};
constexpr Uint maxiter{CcoordOps::get_size(resolutions)*ipow(dim, secondOrder)*10};
constexpr bool verbose{false};
GradIncrements<dim> grads; grads.push_back(delF0);
SolverCG<dim> cg{sys, cg_tol, maxiter, bool(verbose)};
Eigen::ArrayXXd res1{de_geus(sys, grads, cg, newton_tol, verbose)[0].grad};
SolverCG<dim> cg2{sys, cg_tol, maxiter, bool(verbose)};
Eigen::ArrayXXd res2{newton_cg(sys, grads, cg2, newton_tol, verbose)[0].grad};
BOOST_CHECK_LE(abs(res1-res2).mean(), cg_tol);
}
BOOST_AUTO_TEST_CASE(small_strain_patch_test) {
constexpr Dim_t dim{twoD};
using Ccoord = Ccoord_t<dim>;
using Rcoord = Rcoord_t<dim>;
constexpr Ccoord resolutions{CcoordOps::get_cube<dim>(3)};
constexpr Rcoord lengths{CcoordOps::get_cube<dim>(1.)};
constexpr Formulation form{Formulation::small_strain};
// number of layers in the hard material
constexpr Uint nb_lays{1};
constexpr Real contrast{2};
static_assert(nb_lays < resolutions[0],
"the number or layers in the hard material must be smaller "
"than the total number of layers in dimension 0");
auto sys{make_cell(resolutions, lengths, form)};
using Mat_t = MaterialLinearElastic1<dim, dim>;
constexpr Real Young{2.}, Poisson{.33};
auto material_hard{std::make_unique<Mat_t>("hard", contrast*Young, Poisson)};
auto material_soft{std::make_unique<Mat_t>("soft", Young, Poisson)};
for (const auto & pixel: sys) {
if (pixel[0] < Dim_t(nb_lays)) {
material_hard->add_pixel(pixel);
} else {
material_soft->add_pixel(pixel);
}
}
sys.add_material(std::move(material_hard));
sys.add_material(std::move(material_soft));
sys.initialise();
Grad_t<dim> delEps0{Grad_t<dim>::Zero()};
constexpr Real eps0 = 1.;
//delEps0(0, 1) = delEps0(1, 0) = eps0;
delEps0(0, 0) = eps0;
constexpr Real cg_tol{1e-8}, newton_tol{1e-5}, equil_tol{1e-10};
constexpr Uint maxiter{dim*10};
constexpr Dim_t verbose{0};
SolverCGEigen<dim> cg{sys, cg_tol, maxiter, bool(verbose)};
auto result = de_geus(sys, delEps0, cg, newton_tol,
equil_tol, verbose);
if (verbose) {
std::cout << "result:" << std::endl << result.grad << std::endl;
std::cout << "mean strain = " << std::endl
<< sys.get_strain().get_map().mean() << std::endl;
}
/**
* verification of resultant strains: subscript ₕ for hard and ₛ
* for soft, Nₕ is nb_lays and Nₜₒₜ is resolutions, k is contrast
*
* Δl = εl = Δlₕ + Δlₛ = εₕlₕ+εₛlₛ
* => ε = εₕ Nₕ/Nₜₒₜ + εₛ (Nₜₒₜ-Nₕ)/Nₜₒₜ
*
* σ is constant across all layers
* σₕ = σₛ
* => Eₕ εₕ = Eₛ εₛ
* => εₕ = 1/k εₛ
* => ε / (1/k Nₕ/Nₜₒₜ + (Nₜₒₜ-Nₕ)/Nₜₒₜ) = εₛ
*/
constexpr Real factor{1/contrast * Real(nb_lays)/resolutions[0]
+ 1.-nb_lays/Real(resolutions[0])};
constexpr Real eps_soft{eps0/factor};
constexpr Real eps_hard{eps_soft/contrast};
if (verbose) {
std::cout << "εₕ = " << eps_hard << ", εₛ = " << eps_soft << std::endl;
std::cout << "ε = εₕ Nₕ/Nₜₒₜ + εₛ (Nₜₒₜ-Nₕ)/Nₜₒₜ" << std::endl;
}
Grad_t<dim> Eps_hard; Eps_hard << eps_hard, 0, 0, 0;
Grad_t<dim> Eps_soft; Eps_soft << eps_soft, 0, 0, 0;
// verify uniaxial tension patch test
for (const auto & pixel: sys) {
if (pixel[0] < Dim_t(nb_lays)) {
BOOST_CHECK_LE((Eps_hard-sys.get_strain().get_map()[pixel]).norm(), tol);
} else {
BOOST_CHECK_LE((Eps_soft-sys.get_strain().get_map()[pixel]).norm(), tol);
}
}
delEps0 = Grad_t<dim>::Zero();
delEps0(0, 1) = delEps0(1, 0) = eps0;
SolverCG<dim> cg2{sys, cg_tol, maxiter, bool(verbose)};
result = newton_cg(sys, delEps0, cg2, newton_tol,
equil_tol, verbose);
Eps_hard << 0, eps_hard, eps_hard, 0;
Eps_soft << 0, eps_soft, eps_soft, 0;
// verify pure shear patch test
for (const auto & pixel: sys) {
if (pixel[0] < Dim_t(nb_lays)) {
BOOST_CHECK_LE((Eps_hard-sys.get_strain().get_map()[pixel]).norm(), tol);
} else {
BOOST_CHECK_LE((Eps_soft-sys.get_strain().get_map()[pixel]).norm(), tol);
}
}
}
BOOST_AUTO_TEST_SUITE_END();
} // muSpectre