diff --git a/tests/test_model.cpp b/tests/test_model.cpp
index 00f6cfc..ff91260 100644
--- a/tests/test_model.cpp
+++ b/tests/test_model.cpp
@@ -1,181 +1,181 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2022 Lucas Frérot
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "hooke.hh"
 #include "isotropic_hardening.hh"
 #include "logger.hh"
 #include "model_factory.hh"
 #include "mpi_interface.hh"
 #include "test.hh"
 #include <random>
 /* -------------------------------------------------------------------------- */
 
 using namespace tamaas;
 
 TEST(TestModel, applyElasticity) {
   auto model =
       ModelFactory::createModel(model_type::volume_2d, {1., 1., 1.}, {1, 1, 1});
 
   Grid<Real, 3> gradient({1, 1, 1}, 9), stress({1, 1, 1}, 9);
 
   // Random data objects
   std::random_device rnd;
   std::mt19937 mt(rnd());
   std::normal_distribution<> dis(0, 1);
 
   // Filling gradient with random data (sequential)
   for (auto& du : gradient)
     du = dis(mt);
 
   std::uniform_real_distribution<> unif(0, 0.5);
   model->setElasticity(std::abs(dis(mt)) * unif(mt), unif(mt));
 
   // Computing stresses with seperate array
   model->applyElasticity(stress, gradient);
 
   // Checking correct isotropic elasticity
   auto mu = model->getShearModulus(), nu = model->getPoissonRatio();
   auto lambda = 2 * mu * nu / (1 - 2 * nu);
 
   MatrixProxy<Real, 3, 3> grad(gradient(0));
   auto trace = grad.trace();
 
   Matrix<Real, 3, 3> sigma;
 
   for (UInt i = 0; i < 3; ++i)
     for (UInt j = 0; j < 3; ++j)
       sigma(i, j) = (i == j) * lambda * trace + mu * (grad(i, j) + grad(j, i));
 
   EXPECT_TRUE(compare(stress, sigma, AreFloatEqual())) << "Elasticity fail";
 
   // Computing stress with same array
   model->applyElasticity(gradient, gradient);
 
   EXPECT_TRUE(compare(gradient, stress, AreFloatEqual()))
       << "Applying elasticity in-place fail";
 }
 
 TEST(TestModel, hookeField) {
+  mpi::sequential_guard guard;
   auto model =
       ModelFactory::createModel(model_type::volume_2d, {1., 1., 1.}, {2, 1, 1});
 
   Grid<Real, 3> gradient({2, 1, 1}, 9), stress({2, 1, 1}, 9);
 
   auto hooke_field = HookeField<model_type::volume_2d>(model.get());
   auto& mu = hooke_field.template field<Real>("mu");
   mu(0) *= .5;
 
-  Loop::loop([](MatrixProxy<Real, 3, 3> g) {
-    g(0, 0) = g(1, 1) = g(2, 2) = 1;
-  }, range<MatrixProxy<Real, 3, 3>>(gradient));
+  Loop::loop([](MatrixProxy<Real, 3, 3> g) { g(0, 0) = g(1, 1) = g(2, 2) = 1; },
+             range<MatrixProxy<Real, 3, 3>>(gradient));
 
   hooke_field.apply(gradient, stress);
 
   ASSERT_NEAR(stress(0, 0, 0, 0), stress(1, 0, 0, 0) * 0.5, 1e-12);
 }
 
 TEST(TestModel, applyElasticitySym) {
   auto model =
       ModelFactory::createModel(model_type::volume_2d, {1., 1., 1.}, {1, 1, 1});
 
   Grid<Real, 3> gradient({1, 1, 1}, 6), stress({1, 1, 1}, 6);
 
   // Random data objects
   std::random_device rnd;
   std::mt19937 mt(rnd());
   std::normal_distribution<> dis(0, 1);
 
   // Filling gradient with random data (sequential)
   for (auto& du : gradient)
     du = dis(mt);
 
   std::uniform_real_distribution<> unif(0, 0.5);
   model->setElasticity(std::abs(dis(mt)) * unif(mt), unif(mt));
 
   // Computing stresses with seperate array
   model->applyElasticity(stress, gradient);
 
   // Checking correct isotropic elasticity
   auto mu = model->getShearModulus(), nu = model->getPoissonRatio();
   auto lambda = 2 * mu * nu / (1 - 2 * nu);
 
   SymMatrixProxy<Real, 3> sym_grad(gradient(0));
   auto grad = dense(sym_grad);
   auto trace = grad.trace();
 
   Matrix<Real, 3, 3> sigma;
 
   for (UInt i = 0; i < 3; ++i)
     for (UInt j = 0; j < 3; ++j)
       sigma(i, j) = (i == j) * lambda * trace + mu * (grad(i, j) + grad(j, i));
 
   EXPECT_TRUE(compare(stress, symmetrize(sigma), AreFloatEqual()))
       << "Elasticity fail";
 
   // Computing stress with same array
   model->applyElasticity(gradient, gradient);
 
   EXPECT_TRUE(compare(gradient, stress, AreFloatEqual()))
       << "Applying elasticity in-place fail";
 }
 
 TEST(TestIsotropicHardening, computePlasticIncrement) {
   mpi::sequential_guard guard;
 
   auto model =
       ModelFactory::createModel(model_type::volume_2d, {1., 1., 1.}, {1, 1, 1});
   auto sigma_0 = 1., h = 0.1;
   IsotropicHardening<model_type::volume_2d> hardening(model.get(), sigma_0, h);
   Grid<Real, 3> strain({1, 1, 1}, 6), strain_increment({1, 1, 1}, 6),
       plastic_strain_increment({1, 1, 1}, 6), solution({1, 1, 1}, 6);
   auto E = model->getYoungModulus();
   auto nu = model->getPoissonRatio();
   auto mu = model->getShearModulus();
 
   // uniform tension state
   Real sigma = 0.9;
   strain_increment(0, 0, 0, 0) = sigma / E;
   strain_increment(0, 0, 0, 1) = -nu * sigma / E;
   strain_increment(0, 0, 0, 2) = -nu * sigma / E;
 
   hardening.computePlasticIncrement<false>(plastic_strain_increment, strain,
                                            strain_increment);
   // just checking everything is zero
   EXPECT_TRUE(compare(solution, plastic_strain_increment, AreFloatEqual()))
       << "Elastic radial return fail (plastic increment)";
 
   // uniform tension state
   sigma = 1.1;
   strain_increment(0, 0, 0, 0) = sigma / E;
   strain_increment(0, 0, 0, 1) = -nu * sigma / E;
   strain_increment(0, 0, 0, 2) = -nu * sigma / E;
 
   hardening.computePlasticIncrement<false>(plastic_strain_increment, strain,
                                            strain_increment);
   Real a = (sigma - sigma_0) / (3 * mu + h);
   solution(0, 0, 0, 0) = a;
   solution(0, 0, 0, 1) = -a / 2;
   solution(0, 0, 0, 2) = -a / 2;
 
   EXPECT_TRUE(compare(solution, plastic_strain_increment, AreFloatEqual()))
       << "Plastic radial return fail (plastic increment)";
 }