runConvergence.cpp

Go to the documentation of this file.

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   This file is part of CEPS.
 
   CEPS 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 of the License, or
   (at your option) any later version.
 
   CEPS 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 CEPS (see file LICENSE at root of project).
   If not, see <https://www.gnu.org/licenses/>.
 
 
   Copyright 2019-2024 Inria, Universite de Bordeaux
 
   Authors, in alphabetical order:
 
     Pierre-Elliott BECUE, Florian CARO, Yves COUDIERE(*), Andjela DAVIDOVIC, 
     Charlie DOUANLA-LONTSI, Marc FUENTES, Mehdi JUHOOR, Michael LEGUEBE(*), 
     Pauline MIGERDITICHAN, Valentin PANNETIER(*), Nejib ZEMZEMI.
     * : currently active authors
 
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
#include "applications/runConvergence.hpp"
#include "applications/AllProblems.hpp"
 
template<typename PbType>
void
copyActivationMap(ConvergenceStudy<PbType>* cs, CepsReal dt, CepsReal hmax)
{
  InputParameters* p    = cs->getParameters();
  CepsString       base = ceps::getBaseName(p->getString("output file name"));
  CepsString       dir  = ceps::getDir     (p->getString("output file name"));
  if (p->hasKey("activation time data"))
  {
    CepsString dir2 = dir + "/snapshots/";
    CepsString amap = dir2 + base + "_APD-0";
    CepsString ext  = (ceps::fileExists(amap+".pvtu")) ? ".pvtu" : ".vtk";
    CepsString command = "cp " + amap + ext + " " + dir + "/" + base + "_APD_" 
                       + ceps::toString(dt) + "_" + ceps::toString(hmax) + ext;
    ceps::execute(command);
  }
  if (p->hasKey("probe points"))
  {
    CepsString command = "cp " + dir + "/" + base + "_probes.dat "
                       + dir + "/" + base + "_probes_" + ceps::toString(dt) + "_" + ceps::toString(hmax) + ".dat";
    ceps::execute(command);
  }
  return;
}
 
template<typename PbType>
void
runConvergenceT(InputParameters* p, CepsBool refSol, CepsBool copyMap)
{
  ConvergenceStudy<PbType>* cs = ceps::getNew<ConvergenceStudy<PbType>>(p,refSol);
 
  // Set callback if required
  if (copyMap and ceps::isMaster()) {
    // // Callback method for cardiac convergence tests that copies the activation map in the top directory
    // // When this function is called, the output directory is already set to dir/cvNum/
    std::function<void(ConvergenceStudy<PbType>*,CepsReal,CepsReal)> copyFunc = copyActivationMap<PbType>; 
    cs->setCallBackMethod(copyFunc);
  }
 
  // Run the convergence analysis
 
  cs->run();
 
  ceps::destroyObject(cs);
}
 
void
doConvergence(InputParameters* p)
{
  CepsString s = ceps::toKey(p->getString("problem type"));
  if (s == ceps::toKey("dirichlet anode cathode"))
    runConvergenceT<DirichletAnodeCathodeProblem>(p);
  else if (s == ceps::toKey("flux anode cathode"))
    runConvergenceT<FluxAnodeCathodeProblem>(p);
  else if (s == ceps::toKey("monodomain"))
    runConvergenceT<MonodomainProblem>(p,true,true);
  else if (s == ceps::toKey("cl monodomain"))
    runConvergenceT<CLMonodomainProblem>(p,true,true);
  else if (s == ceps::toKey("bidomain"))
    runConvergenceT<BidomainProblem>(p,true,true);
  else if (s == ceps::toKey("extended bidomain"))
    runConvergenceT<ExtendedBidomainProblem>(p, true, true);
  else if (s == ceps::toKey("bilayer monodomain"))
    runConvergenceT<BilayerMonodomainProblem>(p,true,true);
  else if (s == ceps::toKey("pacemaker bidomain"))
    runConvergenceT<PacemakerBidomainProblem>(p, true, true);
  else if (s == ceps::toKey("pacemaker poisson"))
    runConvergenceT<PacemakerPoissonProblem>(p, true, true);
  else 
    CEPS_ABORT(
      "Unknown type of problem " << s << "\n" <<
      "   Valid types are:\n" << validProblems
    );
 
  return; 
}