diff --git a/tests/mpi_test_solver_newton_cg.cc b/tests/mpi_test_solver_newton_cg.cc index 819f9e5..0c3f8d4 100644 --- a/tests/mpi_test_solver_newton_cg.cc +++ b/tests/mpi_test_solver_newton_cg.cc @@ -1,202 +1,201 @@ /** * @file test_solver_newton_cg.cc * * @author Till Junge * * @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 resolutions{3, 3}; // constexpr Rcoord_t lengths{2.3, 2.7}; constexpr Ccoord_t resolutions{5, 5, 5}; constexpr Rcoord_t lengths{5, 5, 5}; auto fft_ptr{std::make_unique>(resolutions, lengths, - MPIContext::get_context().comm)}; + comm)}; auto proj_ptr{std::make_unique>(std::move(fft_ptr))}; CellBase sys(std::move(proj_ptr)); using Mat_t = MaterialLinearElastic1; //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(); for (auto && tup: akantu::enumerate(sys)) { auto && pixel = std::get<1>(tup); if (loc == Ccoord_t{0, 0} && std::get<0>(tup) == 0) { Material_hard.add_pixel(pixel); } else { Material_soft.add_pixel(pixel); } } sys.initialise(); Grad_t 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{true}; + constexpr bool verbose{false}; GradIncrements grads; grads.push_back(delF0); SolverCG cg{sys, cg_tol, maxiter, bool(verbose)}; Eigen::ArrayXXd res1{de_geus(sys, grads, cg, newton_tol, verbose)[0].grad}; -#if 0 SolverCG 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); -#endif } -#if 0 BOOST_AUTO_TEST_CASE(small_strain_patch_test) { + const Communicator & comm = MPIContext::get_context().comm; constexpr Dim_t dim{twoD}; using Ccoord = Ccoord_t; using Rcoord = Rcoord_t; constexpr Ccoord resolutions{CcoordOps::get_cube(3)}; constexpr Rcoord lengths{CcoordOps::get_cube(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)}; + auto sys{make_parallel_cell(resolutions, lengths, form, comm)}; using Mat_t = MaterialLinearElastic1; constexpr Real Young{2.}, Poisson{.33}; auto material_hard{std::make_unique("hard", contrast*Young, Poisson)}; auto material_soft{std::make_unique("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 delEps0{Grad_t::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 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 Eps_hard; Eps_hard << eps_hard, 0, 0, 0; Grad_t 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::Zero(); delEps0(0, 1) = delEps0(1, 0) = eps0; SolverCG 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); } } } -#endif BOOST_AUTO_TEST_SUITE_END(); } // muSpectre