AlephIndexTest.cc

// -*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
//-----------------------------------------------------------------------------
// Copyright 2000-2026 CEA (www.cea.fr) IFPEN (www.ifpenergiesnouvelles.com)
// See the top-level COPYRIGHT file for details.
// SPDX-License-Identifier: Apache-2.0
//-----------------------------------------------------------------------------
/*---------------------------------------------------------------------------*/
/* AlephIndexTest.cc                                           (C) 2000-2015 */
/*                                                                           */
/* Test service for the Aleph+Index service.                                 */
/*---------------------------------------------------------------------------*/
 
#include "arcane/utils/ScopedPtr.h"
 
#include "arcane/aleph/tests/AlephTest.h"
#include "arcane/aleph/tests/AlephIndexTest.h"
#include "arcane/aleph/tests/AlephIndexTest_axl.h"
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
namespace ArcaneTest
{
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
using namespace Arcane;
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
class AlephIndexTest
: public ArcaneAlephIndexTestObject
{
 public:
 
  AlephIndexTest(const ModuleBuildInfo&);
  ~AlephIndexTest(void);
  void init(void);
  void compute(void);
 
 private:
 
  void setValues(const Real, AlephMatrix*);
  void postSolver(const Integer, Real, Array<Real>&, Array<Integer>&);
  static Real geoFaceSurface(Face, VariableItemReal3&);
 
 public:
 
  Integer m_total_nb_cell;
  Integer m_local_nb_cell;
  UniqueArray<Real> m_vector_zeroes;
  AlephKernel* m_aleph_kernel;
  UniqueArray<AlephMatrix*> m_aleph_mat;
  UniqueArray<AlephVector*> m_aleph_rhs;
  UniqueArray<AlephVector*> m_aleph_sol;
  UniqueArray<AlephParams*> m_aleph_params;
  Integer m_get_solution_idx;
  Integer m_nb_solver;
};
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
AlephIndexTest::
AlephIndexTest(const ModuleBuildInfo& mbi)
: ArcaneAlephIndexTestObject(mbi)
, m_total_nb_cell(0)
, m_local_nb_cell(0)
, m_vector_zeroes(0)
, m_aleph_kernel(0)
, m_get_solution_idx(0)
, m_nb_solver(0)
{}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
AlephIndexTest::
~AlephIndexTest(void)
{
  debug() << "[AlephIndexTest::AlephIndexTest] Delete & Free";
}
 
/***************************************************************************
 * AlephIndexTest::init                                                   *
 ***************************************************************************/
void AlephIndexTest::
init(void)
{
  m_nb_solver = options()->alephNumberOfSolvers;
 
  m_aleph_mat.resize(m_nb_solver);
  m_aleph_rhs.resize(m_nb_solver);
  m_aleph_sol.resize(m_nb_solver);
  m_aleph_params.resize(m_nb_solver);
 
  // Set our time step
  m_global_deltat = 1.0;
 
  // initialization of temperature on all cells and faces
  m_cell_temperature.fill(options()->initTemperature());
  m_face_temperature.fill(options()->initTemperature());
 
  // initialization of boundary temperature, loop over boundary conditions
  for (int i = options()->boundaryCondition.size() - 1; i >= 0; --i) {
    Real temperature = options()->boundaryCondition[i]->value();
    FaceGroup face_group = options()->boundaryCondition[i]->surface();
    // loop over faces of the surface
    ENUMERATE_FACE (iFace, face_group)
      m_face_temperature[iFace] = temperature;
  }
  mesh()->checkValidMeshFull();
 
  // We instantiate a minimalistic kernel that will handle the initialization for us
  // AlephKernel(subDomain(),0,0); seems to cause issues, but 2,0, or 0,1 do not:
  // pre-kernel topology initialization problem?
  // The value 2 corresponds to the Kernel for Hypre.
  m_aleph_kernel = new AlephKernel(subDomain(), 2, 1);
 
  // Another initialization of vector indices
  ENUMERATE_CELL (cell, ownCells()) {
    m_vector_zeroes.add(0.0);
  }
}
 
/***************************************************************************
 * Entry points leads us here
 ***************************************************************************/
void AlephIndexTest::
compute(void)
{
  UniqueArray<SharedArray<Integer>> indexs;
  UniqueArray<SharedArray<Real>> values;
  Integer nb_iteration;
  Real residual_norm[4];
 
  if (m_nb_solver == 0) {
    subDomain()->timeLoopMng()->stopComputeLoop(true);
    return;
  }
 
  indexs.resize(m_nb_solver);
  values.resize(m_nb_solver);
  for (int i = 0; i < m_nb_solver; ++i)
    postSolver(i, options()->deltaT, values[i], indexs[i]);
 
  // And should be able to get the solutions after the last resolution was fired
  for (int i = m_nb_solver - 1; i >= 0; --i) {
    debug() << "[AlephIndexTest::compute] Getting solution #" << i;
    AlephVector* solution = m_aleph_kernel->syncSolver(i, nb_iteration, &residual_norm[0]);
    if (i != m_get_solution_idx)
      continue;
    debug() << "Solved in \33[7m" << nb_iteration << "\33[m iterations,"
            << "residuals=[\33[1m" << residual_norm[0] << "\33[m," << residual_norm[3] << "]";
    debug() << "[AlephIndexTest::compute] Applying solution #" << m_get_solution_idx;
    solution->getLocalComponents(values[i]);
    m_get_solution_idx += 1;
  }
  m_get_solution_idx %= m_nb_solver;
 
  // Otherwise, we copy the results
  debug() << "[AlephIndexTest::compute] Now get our results";
  ENUMERATE_CELL (cell, ownCells())
    m_cell_temperature[cell] = values[m_get_solution_idx][m_aleph_kernel->indexing()->get(m_cell_temperature, *cell)];
 
  // If we have reached our maximum number of iterations, we exit
  if (subDomain()->commonVariables().globalIteration() >= options()->iterations)
    subDomain()->timeLoopMng()->stopComputeLoop(true);
 
  debug() << "[AlephIndexTest::compute] done";
}
 
/***************************************************************************
 * job
 ***************************************************************************/
void AlephIndexTest::
postSolver(const Integer i, Real optionDeltaT,
           Array<Real>& values,
           Array<Integer>& indexs)
{
  // We force the deltaTs to be different to have calculation times that can be scheduled
  Real deltaT = (1.0 + (Real)i) * optionDeltaT;
  Integer fake_nb_iteration = 0;
  Real fake_residual_norm[4];
 
  m_aleph_params[i] = new AlephParams(traceMng(),
                                      1.0e-10, // m_param_epsilon convergence epsilon
                                      2000, // m_param_max_iteration max iterations
                                      TypesSolver::AMG, // m_param_preconditioner_method preconditioning: DIAGONAL, AMG, IC
                                      TypesSolver::PCG, // m_param_solver_method solving method
                                      -1, // m_param_gamma
                                      -1.0, // m_param_alpha
                                      false, // m_param_xo_user by default Xo is not equal to 0
                                      false, // m_param_check_real_residue
                                      false, // m_param_print_real_residue
                                      false, // m_param_debug_info
                                      1.e-40, // m_param_min_rhs_norm
                                      false, // m_param_convergence_analyse
                                      true, // m_param_stop_error_strategy
                                      false, // m_param_write_matrix_to_file_error_strategy
                                      "SolveErrorAlephMatrix.dbg", // m_param_write_matrix_name_error_strategy
                                      false, // m_param_listing_output
                                      0., // m_param_threshold
                                      false, // m_param_print_cpu_time_resolution
                                      0, // m_param_amg_coarsening_method: by default Sloop's,
                                      0, // m_param_output_level
                                      1, // m_param_amg_cycle: 1-cycle amg in V, 2= amg cycle in W, 3=cycle in Full Multigrid V
                                      1, // m_param_amg_solver_iterations
                                      1, // m_param_amg_smoother_iterations
                                      TypesSolver::SymHybGSJ_smoother, // m_param_amg_smootherOption
                                      TypesSolver::ParallelRugeStuben, // m_param_amg_coarseningOption
                                      TypesSolver::CG_coarse_solver, // m_param_amg_coarseSolverOption
                                      false, // m_param_keep_solver_structure
                                      false, // m_param_sequential_solver
                                      TypesSolver::RB); // m_param_criteria_stop
 
  // Filling the right-hand side: boundary conditions + right-hand side
  debug() << "[AlephIndexTest::compute] Resize (" << m_vector_zeroes.size() << ") + filling the right-hand side";
  values.resize(m_vector_zeroes.size());
  indexs.resize(m_vector_zeroes.size());
  ENUMERATE_CELL (cell, ownCells()) {
    Integer row = m_aleph_kernel->indexing()->get(m_cell_temperature, cell);
    indexs[row] = row;
    values[row] = m_cell_temperature[cell];
  }
 
  debug() << "[AlephIndexTest::compute] ENUMERATE_FACE";
  ENUMERATE_FACE (iFace, allCells().outerFaceGroup()) {
    if (!iFace->cell(0).isOwn())
      continue;
    values[m_aleph_kernel->indexing()->get(m_cell_temperature, iFace->cell(0))] +=
    deltaT * (m_face_temperature[iFace]) / geoFaceSurface(*iFace, nodesCoordinates());
  }
 
  // Creation of the MatVec matrix and Aleph requirements
  m_aleph_mat[i] = m_aleph_kernel->createSolverMatrix();
  m_aleph_rhs[i] = m_aleph_kernel->createSolverVector(); // First vector returned IS the rhs
  m_aleph_sol[i] = m_aleph_kernel->createSolverVector(); // Next one IS the solution
 
  m_aleph_mat[i]->create();
  m_aleph_rhs[i]->create();
  m_aleph_sol[i]->create();
  m_aleph_mat[i]->reset();
 
  // Filling the matrix and assembly
  setValues(deltaT, m_aleph_mat[i]);
  m_aleph_mat[i]->assemble();
 
  debug() << "[AlephIndexTest::job] setLocalComponents";
  //m_aleph_rhs[i]->setLocalComponents(indexs.size(), indexs.view(), values.view());
  m_aleph_rhs[i]->setLocalComponents(values.view());
  m_aleph_rhs[i]->assemble();
 
  //m_aleph_sol[i]->setLocalComponents(indexs.size(), indexs.view(), m_vector_zeroes.view());
  m_aleph_sol[i]->setLocalComponents(m_vector_zeroes.view());
  m_aleph_sol[i]->assemble();
  debug() << "[AlephIndexTest::job] done setLocalComponents";
  traceMng()->flush();
 
  // Now solve with Aleph
  debug() << "[AlephIndexTest::job] Now solve with Aleph";
  m_aleph_mat[i]->solve(m_aleph_sol[i],
                        m_aleph_rhs[i],
                        fake_nb_iteration,
                        &fake_residual_norm[0],
                        m_aleph_params[i],
                        true); // We wish to post asynchronously
  debug() << "[AlephIndexTest::job] done with Aleph solve";
  traceMng()->flush();
}
 
/***************************************************************************
 * AlephTestModule::_FaceComputeSetValues                                  *
 ***************************************************************************/
void AlephIndexTest::setValues(const Real deltaT, AlephMatrix* aleph_mat)
{
  // Flush the coefficients
  ENUMERATE_CELL (iCell, ownCells())
    m_cell_coefs[iCell] = 0.;
  // 'inner' faces
  debug() << "[AlephIndexTest::setValues] inner-faces";
  ENUMERATE_FACE (iFace, allCells().innerFaceGroup()) {
    if (iFace->backCell().isOwn()) {
      const Real surface = geoFaceSurface(*iFace, nodesCoordinates());
      aleph_mat->setValue(m_cell_temperature, iFace->backCell(),
                          m_cell_temperature, iFace->frontCell(),
                          -deltaT / surface);
      m_cell_coefs[iFace->backCell()] += 1.0 / surface;
    }
    if (iFace->frontCell().isOwn()) {
      const Real surface = geoFaceSurface(*iFace, nodesCoordinates());
      aleph_mat->setValue(m_cell_temperature, iFace->frontCell(),
                          m_cell_temperature, iFace->backCell(),
                          -deltaT / surface);
      m_cell_coefs[iFace->frontCell()] += 1.0 / surface;
    }
  }
  debug() << "[AlephIndexTest::setValues] outer-faces";
  ENUMERATE_FACE (iFace, allCells().outerFaceGroup()) {
    if (!iFace->cell(0).isOwn())
      continue;
    m_cell_coefs[iFace->cell(0)] += 1.0 / geoFaceSurface(*iFace, nodesCoordinates());
  }
  debug() << "[AlephIndexTest::setValues] diagonal";
  ENUMERATE_CELL (cell, ownCells()) {
    aleph_mat->setValue(m_cell_temperature, cell,
                        m_cell_temperature, cell,
                        1.0 + deltaT * m_cell_coefs[cell]);
  }
  debug() << "[AlephIndexTest::setValues] done";
}
 
/***************************************************************************
 * geoFaceSurface
 ***************************************************************************/
Real AlephIndexTest::geoFaceSurface(Face face, VariableItemReal3& nodes_coords)
{
  if (face.nbNode() == 4) {
    const Real3 xyz0 = nodes_coords[face.node(0)];
    const Real3 xyz1 = nodes_coords[face.node(1)];
    const Real3 xyz3 = nodes_coords[face.node(3)];
    return (xyz0 - xyz1).normL2() * (xyz0 - xyz3).normL2();
  }
  if (face.nbNode() == 2) {
    const Real3 xyz0 = nodes_coords[face.node(0)];
    const Real3 xyz1 = nodes_coords[face.node(1)];
    return (xyz0 - xyz1).normL2();
  }
  throw FatalErrorException("geoFace", "Nb nodes != 4 !=2");
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
ARCANE_DEFINE_STANDARD_MODULE(AlephIndexTest, AlephIndexTest);
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
} // namespace ArcaneTest
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/