/* SPDX-FileCopyrightText: Copyright (c) 2016 - 2025, the adamantine authors.
 * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
 */

#include "MaterialStates.hh"
#include <Geometry.hh>
#include <MaterialProperty.hh>
#include <ThermalPhysics.hh>
#include <Timer.hh>

#include <deal.II/base/quadrature_lib.h>

#include <fstream>

namespace tt = boost::test_tools;

boost::property_tree::ptree basic_geometry_database()
{
  boost::property_tree::ptree geometry_database;
  geometry_database.put("import_mesh", false);
  geometry_database.put("length", 12e-3);
  geometry_database.put("length_divisions", 1);
  geometry_database.put("height", 6e-3);
  geometry_database.put("height_divisions", 1);

  return geometry_database;
}

boost::property_tree::ptree basic_material_properies_database()
{
  // MaterialProperty database
  boost::property_tree::ptree material_property_database;
  material_property_database.put("property_format", "polynomial");
  material_property_database.put("n_materials", 1);
  material_property_database.put("material_0.solid.density", 1.);
  material_property_database.put("material_0.powder.density", 10.);
  material_property_database.put("material_0.liquid.density", 1.);
  material_property_database.put("material_0.solid.specific_heat", 1.);
  material_property_database.put("material_0.powder.specific_heat", 2.);
  material_property_database.put("material_0.liquid.specific_heat", 1.);
  material_property_database.put("material_0.solid.thermal_conductivity_x", 1.);
  material_property_database.put("material_0.solid.thermal_conductivity_z", 1.);
  material_property_database.put("material_0.powder.thermal_conductivity_x",
                                 1.);
  material_property_database.put("material_0.powder.thermal_conductivity_z",
                                 1.);
  material_property_database.put("material_0.liquid.thermal_conductivity_x",
                                 1.);
  material_property_database.put("material_0.liquid.thermal_conductivity_z",
                                 1.);

  return material_property_database;
}

boost::property_tree::ptree basic_input_database()
{
  boost::property_tree::ptree database;
  // Source database
  database.put("sources.n_beams", 1);
  database.put("sources.beam_0.depth", 1e100);
  database.put("sources.beam_0.diameter", 1e100);
  database.put("sources.beam_0.max_power", 1e300);
  database.put("sources.beam_0.scan_path_file", "scan_path.txt");
  database.put("sources.beam_0.absorption_efficiency", 0.3);
  database.put("sources.beam_0.type", "electron_beam");
  database.put("sources.beam_0.scan_path_file",
               "scan_path_test_thermal_physics.txt");
  database.put("sources.beam_0.scan_path_file_format", "segment");
  // Boundary database
  database.put("boundary.type", "adiabatic");

  // Time-stepping database
  database.put("time_stepping.method", "rk_fourth_order");
  return database;
}

template <typename MemorySpaceType>
void thermal_2d(boost::property_tree::ptree &database, double time_step)
{
  MPI_Comm communicator = MPI_COMM_WORLD;

  // Geometry database
  boost::property_tree::ptree geometry_database;
  geometry_database.put("import_mesh", false);
  geometry_database.put("length", 12e-3);
  geometry_database.put("length_divisions", 4);
  geometry_database.put("height", 6e-3);
  geometry_database.put("height_divisions", 5);
  boost::optional<boost::property_tree::ptree const &> units_optional_database;
  // Build Geometry
  adamantine::Geometry<2> geometry(communicator, geometry_database,
                                   units_optional_database);
  // Create the Boundary
  boost::property_tree::ptree boundary_database;
  boundary_database.put("type", "adiabatic");
  adamantine::Boundary boundary(
      boundary_database, geometry.get_triangulation().get_boundary_ids());
  // MaterialProperty database
  boost::property_tree::ptree material_property_database;
  material_property_database.put("property_format", "polynomial");
  material_property_database.put("n_materials", 1);
  material_property_database.put("material_0.solid.density", 1.);
  material_property_database.put("material_0.powder.density", 1.);
  material_property_database.put("material_0.liquid.density", 1.);
  material_property_database.put("material_0.solid.specific_heat", 1.);
  material_property_database.put("material_0.powder.specific_heat", 1.);
  material_property_database.put("material_0.liquid.specific_heat", 1.);
  material_property_database.put("material_0.solid.thermal_conductivity_x", 1.);
  material_property_database.put("material_0.solid.thermal_conductivity_z", 1.);
  material_property_database.put("material_0.powder.thermal_conductivity_x",
                                 1.);
  material_property_database.put("material_0.powder.thermal_conductivity_z",
                                 1.);
  material_property_database.put("material_0.liquid.thermal_conductivity_x",
                                 1.);
  material_property_database.put("material_0.liquid.thermal_conductivity_z",
                                 1.);
  // Build MaterialProperty
  adamantine::MaterialProperty<2, 1, 2, adamantine::SolidLiquidPowder,
                               MemorySpaceType>
      material_properties(communicator, geometry.get_triangulation(),
                          material_property_database);
  // Source database
  database.put("sources.n_beams", 1);
  database.put("sources.beam_0.depth", 1e100);
  database.put("sources.beam_0.diameter", 1e100);
  database.put("sources.beam_0.max_power", 1e300);
  database.put("sources.beam_0.absorption_efficiency", 0.1);
  database.put("sources.beam_0.type", "electron_beam");
  database.put("sources.beam_0.scan_path_file",
               "scan_path_test_thermal_physics.txt");
  database.put("sources.beam_0.scan_path_file_format", "segment");

  // Build ThermalPhysics
  adamantine::ThermalPhysics<2, 1, 2, 2, adamantine::SolidLiquidPowder,
                             MemorySpaceType, dealii::QGauss<1>>
      physics(communicator, database, geometry, boundary, material_properties);
  physics.setup();
  dealii::LA::distributed::Vector<double, MemorySpaceType> solution;
  physics.initialize_dof_vector(0., solution);

  std::vector<adamantine::Timer> timers(adamantine::Timing::n_timers);
  double time = 0;
  while (time < 0.1)
  {
    time = physics.evolve_one_time_step(time, time_step, solution, timers);
  }

  double const tolerance = 1e-3;
  BOOST_TEST(time == 0.1, tt::tolerance(tolerance));
  BOOST_TEST(solution.l2_norm() == 0.291705, tt::tolerance(tolerance));

  physics.initialize_dof_vector(1000., solution);
  BOOST_TEST(solution.l1_norm() == 1000. * solution.size());
}

template <typename MemorySpaceType>
void thermal_2d_manufactured_solution()
{
  MPI_Comm communicator = MPI_COMM_WORLD;

  // Geometry database
  boost::property_tree::ptree geometry_database;
  geometry_database.put("import_mesh", false);
  geometry_database.put("length", 1e3);
  geometry_database.put("length_divisions", 4);
  geometry_database.put("height", 6e3);
  geometry_database.put("height_divisions", 5);
  boost::optional<boost::property_tree::ptree const &> units_optional_database;
  // Build Geometry
  adamantine::Geometry<2> geometry(communicator, geometry_database,
                                   units_optional_database);
  // Create the Boundary
  boost::property_tree::ptree boundary_database;
  boundary_database.put("type", "adiabatic");
  adamantine::Boundary boundary(
      boundary_database, geometry.get_triangulation().get_boundary_ids());

  // MaterialProperty database
  boost::property_tree::ptree material_property_database;
  material_property_database.put("property_format", "polynomial");
  material_property_database.put("n_materials", 1);
  material_property_database.put("material_0.solid.density", 1.);
  material_property_database.put("material_0.powder.density", 1.);
  material_property_database.put("material_0.liquid.density", 1.);
  material_property_database.put("material_0.solid.specific_heat", 1.);
  material_property_database.put("material_0.powder.specific_heat", 1.);
  material_property_database.put("material_0.liquid.specific_heat", 1.);
  material_property_database.put("material_0.solid.thermal_conductivity_x", 1.);
  material_property_database.put("material_0.solid.thermal_conductivity_z", 1.);
  material_property_database.put("material_0.powder.thermal_conductivity_x",
                                 1.);
  material_property_database.put("material_0.powder.thermal_conductivity_z",
                                 1.);
  material_property_database.put("material_0.liquid.thermal_conductivity_x",
                                 1.);
  material_property_database.put("material_0.liquid.thermal_conductivity_z",
                                 1.);
  // Build MaterialProperty
  adamantine::MaterialProperty<2, 1, 1, adamantine::SolidLiquidPowder,
                               MemorySpaceType>
      material_properties(communicator, geometry.get_triangulation(),
                          material_property_database);

  boost::property_tree::ptree database;
  // Source database
  database.put("sources.n_beams", 1);
  database.put("sources.beam_0.depth", 1e100);
  database.put("sources.beam_0.diameter", 1e100);
  database.put("sources.beam_0.max_power", 1e300);
  database.put("sources.beam_0.absorption_efficiency",
               0.1 / 0.29317423955177113);
  database.put("sources.beam_0.scan_path_file",
               "scan_path_test_thermal_physics.txt");
  database.put("sources.beam_0.type", "electron_beam");
  database.put("sources.beam_0.scan_path_file_format", "segment");

  // Time-stepping database
  database.put("time_stepping.method", "rk_fourth_order");
  // Build ThermalPhysics
  adamantine::ThermalPhysics<2, 1, 1, 2, adamantine::SolidLiquidPowder,
                             MemorySpaceType, dealii::QGauss<1>>
      physics(communicator, database, geometry, boundary, material_properties);
  physics.setup();
  dealii::LA::distributed::Vector<double, MemorySpaceType> solution;
  physics.initialize_dof_vector(0., solution);

  std::vector<adamantine::Timer> timers(adamantine::Timing::n_timers);
  double time = physics.evolve_one_time_step(0., 0.1, solution, timers);

  double const tolerance = 1e-5;

  BOOST_TEST(time == 0.1);

  if (std::is_same<MemorySpaceType, dealii::MemorySpace::Host>::value)
  {
    for (unsigned int i = 0; i < solution.locally_owned_size(); ++i)
      BOOST_TEST(solution.local_element(i) == 0.1, tt::tolerance(tolerance));
  }
  else
  {
    dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host>
        solution_host(solution.get_partitioner());
    solution_host.import_elements(solution, dealii::VectorOperation::insert);
    for (unsigned int i = 0; i < solution.size(); ++i)
      BOOST_TEST(solution_host[i] == 0.1, tt::tolerance(tolerance));
  }
}

template <typename MemorySpaceType>
void initial_temperature()
{
  MPI_Comm communicator = MPI_COMM_WORLD;

  // Build Geometry
  auto geometry_database = basic_geometry_database();
  boost::optional<boost::property_tree::ptree const &> units_optional_database;
  adamantine::Geometry<2> geometry(communicator, geometry_database,
                                   units_optional_database);

  // Build MaterialProperty
  auto material_property_database = basic_material_properies_database();
  adamantine::MaterialProperty<2, -1, 4, adamantine::SolidLiquidPowder,
                               MemorySpaceType>
      material_properties(communicator, geometry.get_triangulation(),
                          material_property_database);

  // Other generic input parameters
  auto database = basic_input_database();

  // Build Boundary
  adamantine::Boundary boundary(
      database.get_child("boundary"),
      geometry.get_triangulation().get_boundary_ids());

  // Build ThermalPhysics
  adamantine::ThermalPhysics<2, -1, 4, 2, adamantine::SolidLiquidPowder,
                             MemorySpaceType, dealii::QGauss<1>>
      physics(communicator, database, geometry, boundary, material_properties);
  physics.setup();
  dealii::LA::distributed::Vector<double, MemorySpaceType> solution;
  physics.initialize_dof_vector(1000., solution);

  BOOST_TEST(solution.l1_norm() == 1000. * solution.size());
}

template <typename MemorySpaceType>
void energy_conservation()
{
  MPI_Comm communicator = MPI_COMM_WORLD;

  // Geometry database
  boost::property_tree::ptree geometry_database;
  geometry_database.put("import_mesh", false);
  geometry_database.put("length", 10);
  geometry_database.put("length_divisions", 10);
  geometry_database.put("height", 10);
  geometry_database.put("height_divisions", 10);
  boost::optional<boost::property_tree::ptree const &> units_optional_database;
  // Build Geometry
  adamantine::Geometry<2> geometry(communicator, geometry_database,
                                   units_optional_database);
  boost::property_tree::ptree material_property_database;
  // Create the Boundary
  boost::property_tree::ptree boundary_database;
  boundary_database.put("type", "adiabatic");
  adamantine::Boundary boundary(
      boundary_database, geometry.get_triangulation().get_boundary_ids());
  // MaterialProperty database
  material_property_database.put("property_format", "polynomial");
  material_property_database.put("n_materials", 1);
  material_property_database.put("material_0.solid.density", 0.5);
  material_property_database.put("material_0.powder.density", 0.5);
  material_property_database.put("material_0.liquid.density", 0.5);
  material_property_database.put("material_0.solid.specific_heat", 4.);
  material_property_database.put("material_0.powder.specific_heat", 4.);
  material_property_database.put("material_0.liquid.specific_heat", 4.);
  material_property_database.put("material_0.solid.thermal_conductivity_x", 2.);
  material_property_database.put("material_0.solid.thermal_conductivity_z", 2.);
  material_property_database.put("material_0.powder.thermal_conductivity_x",
                                 2.);
  material_property_database.put("material_0.powder.thermal_conductivity_z",
                                 2.);
  material_property_database.put("material_0.liquid.thermal_conductivity_x",
                                 2.);
  material_property_database.put("material_0.liquid.thermal_conductivity_z",
                                 2.);
  // Build MaterialProperty
  adamantine::MaterialProperty<2, 1, 0, adamantine::SolidLiquidPowder,
                               MemorySpaceType>
      material_properties(communicator, geometry.get_triangulation(),
                          material_property_database);
  boost::property_tree::ptree database;
  // Source database
  database.put("sources.n_beams", 1);
  database.put("sources.beam_0.type", "cube");
  database.put("sources.beam_0.start_time", 0);
  database.put("sources.beam_0.end_time", 5);
  database.put("sources.beam_0.value", 5);
  database.put("sources.beam_0.min_x", 4);
  database.put("sources.beam_0.min_y", 4);
  database.put("sources.beam_0.max_x", 6);
  database.put("sources.beam_0.max_y", 6);
  // Time-stepping database
  database.put("time_stepping.method", "forward_euler");
  // Build ThermalPhysics
  adamantine::ThermalPhysics<2, 1, 0, 2, adamantine::SolidLiquidPowder,
                             MemorySpaceType, dealii::QGauss<1>>
      physics(communicator, database, geometry, boundary, material_properties);
  physics.setup();
  dealii::LA::distributed::Vector<double, MemorySpaceType> solution;
  double constexpr initial_temperature = 10;
  double constexpr final_temperature = 10.5;
  physics.initialize_dof_vector(initial_temperature, solution);

  std::vector<adamantine::Timer> timers(adamantine::Timing::n_timers);
  double time = 0;
  while (time < 100)
  {
    time = physics.evolve_one_time_step(time, 0.05, solution, timers);
  }

  double max = -1;
  double min = 1e4;
  if (std::is_same<MemorySpaceType, dealii::MemorySpace::Host>::value)
  {
    for (auto v : solution)
    {
      if (max < v)
        max = v;
      if (min > v)
        min = v;
    }
  }
  else
  {
    dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host>
        solution_host(solution.get_partitioner());
    solution_host.import_elements(solution, dealii::VectorOperation::insert);
    for (auto v : solution_host)
    {
      if (max < v)
        max = v;
      if (min > v)
        min = v;
    }
  }

  double constexpr tolerance = 1e-9;
  BOOST_TEST(solution.mean_value() == final_temperature,
             tt::tolerance(tolerance));
  BOOST_TEST(min == max, tt::tolerance(tolerance));
}

template <typename MemorySpaceType>
void radiation_bcs()
{
  MPI_Comm communicator = MPI_COMM_WORLD;

  // Geometry database
  boost::property_tree::ptree geometry_database;
  geometry_database.put("import_mesh", false);
  geometry_database.put("length", 5);
  geometry_database.put("length_divisions", 5);
  geometry_database.put("height", 5);
  geometry_database.put("height_divisions", 5);
  boost::optional<boost::property_tree::ptree const &> units_optional_database;
  // Build Geometry
  adamantine::Geometry<2> geometry(communicator, geometry_database,
                                   units_optional_database);
  // Create the Boundary
  boost::property_tree::ptree boundary_database;
  boundary_database.put("type", "radiative");
  adamantine::Boundary boundary(
      boundary_database, geometry.get_triangulation().get_boundary_ids());
  // MaterialProperty database
  boost::property_tree::ptree material_property_database;
  material_property_database.put("property_format", "polynomial");
  material_property_database.put("n_materials", 1);
  material_property_database.put("material_0.solid.density", 1.);
  material_property_database.put("material_0.powder.density", 1.);
  material_property_database.put("material_0.liquid.density", 1.);
  material_property_database.put("material_0.solid.specific_heat", 1.);
  material_property_database.put("material_0.powder.specific_heat", 1.);
  material_property_database.put("material_0.liquid.specific_heat", 1.);
  material_property_database.put("material_0.solid.thermal_conductivity_x", 1.);
  material_property_database.put("material_0.solid.thermal_conductivity_z", 1.);
  material_property_database.put("material_0.powder.thermal_conductivity_x",
                                 1.);
  material_property_database.put("material_0.powder.thermal_conductivity_z",
                                 1.);
  material_property_database.put("material_0.liquid.thermal_conductivity_x",
                                 1.);
  material_property_database.put("material_0.liquid.thermal_conductivity_z",
                                 1.);
  material_property_database.put("material_0.solid.emissivity", 1.);
  material_property_database.put("material_0.powder.emissivity", 1.);
  material_property_database.put("material_0.liquid.emissivity", 1.);
  material_property_database.put(
      "material_0.solid.radiation_heat_transfer_coef", 1.);
  material_property_database.put(
      "material_0.powder.radiation_heat_transfer_coef", 1.);
  material_property_database.put(
      "material_0.liquid.radiation_heat_transfer_coef", 1.);
  material_property_database.put(
      "material_0.solid.convection_heat_transfer_coef", 1.);
  material_property_database.put(
      "material_0.powder.convection_heat_transfer_coef", 1.);
  material_property_database.put(
      "material_0.liquid.convection_heat_transfer_coef", 1.);
  material_property_database.put("material_0.radiation_temperature_infty",
                                 20.0);
  material_property_database.put("material_0.convection_temperature_infty",
                                 0.0);
  // Build MaterialProperty
  adamantine::MaterialProperty<2, 1, 1, adamantine::SolidLiquidPowder,
                               MemorySpaceType>
      material_properties(communicator, geometry.get_triangulation(),
                          material_property_database);
  boost::property_tree::ptree database;
  // Source database
  database.put("sources.n_beams", 0);
  // Time-stepping database
  database.put("time_stepping.method", "forward_euler");
  // Build ThermalPhysics
  adamantine::ThermalPhysics<2, 1, 1, 2, adamantine::SolidLiquidPowder,
                             dealii::MemorySpace::Host, dealii::QGauss<1>>
      physics(communicator, database, geometry, boundary, material_properties);
  physics.setup();
  dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host> solution;
  double constexpr initial_temperature = 10;
  physics.initialize_dof_vector(initial_temperature, solution);
  std::vector<adamantine::Timer> timers(adamantine::Timing::n_timers);
  double time = 0;
  while (time < 100)
  {
    time = physics.evolve_one_time_step(time, 0.05, solution, timers);
  }

  double max = -1;
  double min = 1e4;
  if (std::is_same<MemorySpaceType, dealii::MemorySpace::Host>::value)
  {
    for (auto v : solution)
    {
      if (max < v)
        max = v;
      if (min > v)
        min = v;
    }
  }
  else
  {
    dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host>
        solution_host(solution.get_partitioner());
    solution_host.import_elements(solution, dealii::VectorOperation::insert);
    for (auto v : solution_host)
    {
      if (max < v)
        max = v;
      if (min > v)
        min = v;
    }
  }

  BOOST_TEST(min >= 10.);
  BOOST_TEST(min <= 20.);
  BOOST_TEST(max > 10.);
  BOOST_TEST(max <= 20.);
}

template <typename MemorySpaceType>
void convection_bcs()
{
  MPI_Comm communicator = MPI_COMM_WORLD;

  // Geometry database
  boost::property_tree::ptree geometry_database;
  geometry_database.put("import_mesh", false);
  geometry_database.put("length", 0.1);
  geometry_database.put("length_divisions", 100);
  geometry_database.put("height", 0.1);
  geometry_database.put("height_divisions", 1);
  geometry_database.put("width", 0.1);
  geometry_database.put("width_divisions", 1);
  boost::optional<boost::property_tree::ptree const &> units_optional_database;
  // Build Geometry
  adamantine::Geometry<3> geometry(communicator, geometry_database,
                                   units_optional_database);
  boost::property_tree::ptree material_property_database;
  // Create the Boundary
  boost::property_tree::ptree boundary_database;
  boundary_database.put("type", "adiabatic");
  boundary_database.put("boundary_1.type", "convective");
  adamantine::Boundary boundary(
      boundary_database, geometry.get_triangulation().get_boundary_ids());
  // MaterialProperty database
  material_property_database.put("property_format", "polynomial");
  material_property_database.put("n_materials", 1);
  material_property_database.put("material_0.solid.density", 100.);
  material_property_database.put("material_0.powder.density", 100.);
  material_property_database.put("material_0.liquid.density", 100.);
  material_property_database.put("material_0.solid.specific_heat", 100.);
  material_property_database.put("material_0.powder.specific_heat", 100.);
  material_property_database.put("material_0.liquid.specific_heat", 100.);
  material_property_database.put("material_0.solid.thermal_conductivity_x", 1.);
  material_property_database.put("material_0.solid.thermal_conductivity_z", 1.);
  material_property_database.put("material_0.powder.thermal_conductivity_x",
                                 1.);
  material_property_database.put("material_0.powder.thermal_conductivity_z",
                                 1.);
  material_property_database.put("material_0.liquid.thermal_conductivity_x",
                                 1.);
  material_property_database.put("material_0.liquid.thermal_conductivity_z",
                                 1.);
  material_property_database.put("material_0.solid.emissivity", 1.);
  material_property_database.put("material_0.powder.emissivity", 1.);
  material_property_database.put("material_0.liquid.emissivity", 1.);
  material_property_database.put(
      "material_0.solid.convection_heat_transfer_coef", 50.);
  material_property_database.put(
      "material_0.powder.convection_heat_transfer_coef", 50.);
  material_property_database.put(
      "material_0.liquid.convection_heat_transfer_coef", 50.);
  material_property_database.put("material_0.radiation_temperature_infty", 0.0);
  material_property_database.put("material_0.convection_temperature_infty",
                                 300.0);
  // Build MaterialProperty
  adamantine::MaterialProperty<3, -1, 0, adamantine::SolidLiquidPowder,
                               MemorySpaceType>
      material_properties(communicator, geometry.get_triangulation(),
                          material_property_database);
  boost::property_tree::ptree database;
  // Source database
  database.put("sources.n_beams", 0);
  // Time-stepping database
  database.put("time_stepping.method", "forward_euler");
  // Build ThermalPhysics
  adamantine::ThermalPhysics<3, -1, 0, 2, adamantine::SolidLiquidPowder,
                             MemorySpaceType, dealii::QGauss<1>>
      physics(communicator, database, geometry, boundary, material_properties);
  physics.setup();
  dealii::LA::distributed::Vector<double, MemorySpaceType> solution;
  double constexpr initial_temperature = 400;
  physics.initialize_dof_vector(initial_temperature, solution);
  std::vector<adamantine::Timer> timers(adamantine::Timing::n_timers);
  double time = 0;
  while (time < 1.0)
  {
    time = physics.evolve_one_time_step(time, 2.5e-4, solution, timers);
  }

  if (dealii::Utilities::MPI::n_mpi_processes(communicator) == 1)
  {
    // Compare the solution to a gold solution. Comparison to an analytical
    // solution, use analytical_1D_transient.py
    std::ifstream gold_file("convective_bc_gold.txt");
    if (std::is_same<MemorySpaceType, dealii::MemorySpace::Host>::value)
    {
      for (unsigned int i = 0; i < solution.locally_owned_size(); ++i)
      {
        double gold_value = -1.;
        gold_file >> gold_value;
        BOOST_TEST(solution.local_element(i) == gold_value);
      }
    }
    else
    {
      dealii::LA::distributed::Vector<double, dealii::MemorySpace::Host>
          solution_host(solution.get_partitioner());
      solution_host.import_elements(solution, dealii::VectorOperation::insert);
      for (unsigned int i = 0; i < solution.locally_owned_size(); ++i)
      {
        double gold_value = -1.;
        gold_file >> gold_value;
        BOOST_TEST(solution_host.local_element(i) == gold_value);
      }
    }
  }
  else
  {
    BOOST_TEST(solution.l2_norm() == 16860.163757);
  }
}

template <typename MemorySpaceType>
void reference_temperature()
{
  MPI_Comm communicator = MPI_COMM_WORLD;

  // Build Geometry
  auto geometry_database = basic_geometry_database();
  boost::optional<boost::property_tree::ptree const &> units_optional_database;
  adamantine::Geometry<2> geometry(communicator, geometry_database,
                                   units_optional_database);

  // Build MaterialProperty
  auto material_property_database = basic_material_properies_database();
  adamantine::MaterialProperty<2, -1, 4, adamantine::SolidLiquidPowder,
                               MemorySpaceType>
      material_properties(communicator, geometry.get_triangulation(),
                          material_property_database);

  // Other generic input parameters
  auto database = basic_input_database();

  // Build Boundary
  adamantine::Boundary boundary(
      database.get_child("boundary"),
      geometry.get_triangulation().get_boundary_ids());

  // Build ThermalPhysics
  adamantine::ThermalPhysics<2, -1, 4, 2, adamantine::SolidLiquidPowder,
                             MemorySpaceType, dealii::QGauss<1>>
      physics(communicator, database, geometry, boundary, material_properties);
  physics.setup();
  dealii::LA::distributed::Vector<double, MemorySpaceType> solution;
  physics.initialize_dof_vector(1000., solution);

  // Now check that the melting indicator works as expected
  std::vector<double> reference_temperatures({1500.0, 300.0});

  auto has_melted = physics.get_has_melted_vector();

  // Check that has_melted is the correct size
  BOOST_TEST(has_melted.size() ==
             geometry.get_triangulation().n_locally_owned_active_cells());

  // Check that nothing is marked as melted
  for (auto indicator : has_melted)
    BOOST_CHECK(indicator == false);

  // Mark cells above the melting temperature, expect none to be marked
  physics.mark_has_melted(reference_temperatures[0], solution);
  has_melted = physics.get_has_melted_vector();

  for (auto indicator : has_melted)
    BOOST_CHECK(indicator == false);

  // Increase the temperature above the reference temperature, expect all to be
  // marked now
  for (unsigned int i = 0; i < solution.locally_owned_size(); ++i)
    solution.local_element(i) = 1600.0;

  physics.mark_has_melted(reference_temperatures[0], solution);
  has_melted = physics.get_has_melted_vector();

  for (auto indicator : has_melted)
    BOOST_CHECK(indicator == true);

  // Decrease the temperature back below the reference temperature, expect all
  // to still be marked
  for (unsigned int i = 0; i < solution.locally_owned_size(); ++i)
    solution.local_element(i) = 1100.0;

  physics.mark_has_melted(reference_temperatures[0], solution);
  has_melted = physics.get_has_melted_vector();

  for (auto indicator : has_melted)
    BOOST_CHECK(indicator == true);
}