Homec4science
Diffusion µSpectre prototype implementation (add_eigenstrain) Browse

muSpectre_proto/bin/demonstrator2.cc
8353160ae1bdadd_eigenstrain

Last Change
Raw File

demonstrator2.cc
View Options

/**
* file demonstator2.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 15 May 2017
*
* @brief for comparison with akantu demonstrator
*
* @section LICENCE
*
* Copyright (C) 2017 Till Junge
*
* µSpectr 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.
*
* µSpectr 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 <iostream>
#include <vector>
#include <string>
#include "common/common.hh"
#include "external/cxxopts.hpp"
#include <memory>
#include "system/system_base.hh"
#include "system/fftw_engine_r2c.hh"
#include "materials/material_hyper_elastic.hh"
#include "chrono"
using opt_ptr = std::unique_ptr<cxxopts::Options>;
opt_ptr parse_args(int argc, char **argv) {
opt_ptr options =
std::make_unique<cxxopts::Options>(argv[0], "Tests MPI fft scalability");
try {
options->add_options()
("0,N0", "number of rows", cxxopts::value<int>(), "N0")
("h,help", "print help")
("positional",
"Positional arguments: these are the arguments that are entered "
"without an option", cxxopts::value<std::vector<std::string>>());
options->parse_positional(std::vector<std::string>{"N0", "positional"});
options->parse(argc, argv);
if (options->count("help")) {
std::cout << options->help({"", "Group"}) << std::endl;
exit(0);
}
if (options->count("N0") != 1 ) {
throw cxxopts::OptionException("Parameter N0 missing");
} else if ((*options)["N0"].as<int>()%2 != 1) {
throw cxxopts::OptionException("N0 must be odd");
} else if (options->count("positional") > 0) {
throw cxxopts::OptionException("There are too many positional arguments");
}
} catch (const cxxopts::OptionException & e) {
std::cout << "Error parsing options: " << e.what() << std::endl;
exit(1);
}
return options;
}
using namespace muSpectre;
int main(int argc, char *argv[])
{
auto options = parse_args(argc, argv);
auto & opt = *options;
const Dim_t size = opt["N0"].as<int>();
const Real fsize = 1;
const Dim_t dim = 3;
std::cout << "Number of dofs: " << size*size*size*dim*dim << std::endl;
const std::array<Real, dim> sizes({fsize, fsize, fsize});
const std::array<Dim_t, dim> nb_pixels{size, size, size};
using Sys_t = SystemBase<dim, dim>;
SystemBase<dim, dim> system(sizes, nb_pixels);
auto start = std::chrono::high_resolution_clock::now();
typename Sys_t::FFT_Ptr fft_eng =
std::make_shared<FFTW_EngineR2C<dim, dim>>
(nb_pixels, FFT_PlanFlags::measure);
std::chrono::duration<Real> dur_prep = std::chrono::high_resolution_clock::now() - start;
const Real E = 1.0030648180242636;
const Real nu = 0.29930675909878679;
typename Sys_t::mat_Ptr hyperelastic_soft =
std::make_shared<MaterialHyperElastic<dim, dim>>
("soft", E, nu);
typename Sys_t::mat_Ptr hyperelastic_hard =
std::make_shared<MaterialHyperElastic<dim, dim>>
("hard", 10*E, nu);
system.set_fft_engine(fft_eng);
system.add_material(hyperelastic_soft);
system.add_material(hyperelastic_hard);
Real newton_tol = 1e-4;
Real cg_tol = 1e-7;
size_t maxiter = 0;
int counter= 0;
for (const auto && pixel:system) {
int sum = 0;
for (Dim_t i = 0; i < dim; ++i) {
sum += pixel[i]*2 / nb_pixels[i];
}
if (sum == 0) {
hyperelastic_hard->add_pixel(pixel);
counter ++;
} else {
hyperelastic_soft->add_pixel(pixel);
}
}
typename Sys_t::T2 DeltaF;
DeltaF.setZero();
DeltaF(0, 1) = .1;
bool verbose = false;
start = std::chrono::high_resolution_clock::now();
system.solve(DeltaF, newton_tol, cg_tol, maxiter, verbose);
std::chrono::duration<Real> dur = std::chrono::high_resolution_clock::now() - start;
std::cout << "FFT prep time = " << dur_prep.count() << "s" << std::endl;
std::cout << "Resolution time per iter= " << dur.count()/system.get_nb_cg_calls() << "s" << std::endl;
std::cout << system.get_nb_linsolv()/dur.count() << " linear solves per second" << std::endl;
return 0;
}
c4science · Help