tstMultiZoneDepleter.cpp

./Core/xf/tests/tstMultiZoneDepleter.cpp

#include <algorithm>
#include <iostream>
#include <string>
#include <vector>

#include "Origen/Core/io/LibraryIO.h"
#include "Origen/Core/xf/MultiZoneDepleter.h"
#include "Origen/Core/xf/Solver_Fake.h"
#include "ScaleUtils/IO/DB.h"

#include "Nemesis/comm/Logger.hh"
#include "Nemesis/gtest/nemesis_gtest.hh"
#include "Nemesis/harness/DBC.hh"
#include "Origen/Core/dc/FakeFactory.h"
#include "Origen/Core/dc/TransitionMatrixP.h"
#include "ScaleSTL/Functions.h"
#include "Standard/Interface/AbstractWriter.h"
#include "Standard/Interface/BasicIOWriter.h"
#include "Standard/Interface/Communicator.h"
#include "Standard/Interface/jdebug.h"

using ScaleUtils::IO::nprintf;
using ScaleSTL::approx_eq;
typedef std::vector<std::string> Vec_Str;

// GLOBAL COMMUNICATOR
Standard::Communicator world;

TEST( MultiZoneDepleter, CommTest )
{
    //    jDebugLine("Host: '" << world.host()
    //            << "', rank: " << world.rank()
    //            << ", size: " << world.size());
    double x = world.rank();
    std::vector<double> vecx;
    world.allgather( x, vecx );

    for( size_t i = 0; i < vecx.size(); ++i )
    {
        // pnout << nprintf("rank: %3d
        // val[%3d]=%8.2f\n",world.rank(),i,vecx[i]);
        EXPECT_EQ( i, vecx[i] );
    }
    world.barrier();

    std::ofstream ofs(
        nprintf( "%s_%04d", world.host().c_str(), world.rank() ) );
    ofs << world.host() << std::endl;
    ofs << world.rank() << std::endl;
    for( size_t i = 0; i < vecx.size(); ++i )
    {
        ofs << nprintf( "val[%3d]=%8.2f\n", i, vecx[i] );
    }
    ofs.close();
}

void get_name_id( int start, std::string& name, int& id )
{
    id = 1000 * ( world.rank() + 1 ) + start;
    name = nprintf( "mix%05d", id );
}

Origen::SP_Material setup_material1( Origen::SP_Library lib )
{
    // Don't be alarmed, this is all just approximate.
    double cm2_barn = 1e-24;              // cm^2/barn
    double at_mol = 6.022e23;             // at/mol
    double dens = 10.4;                   // g/cm^3
    double vol = 12.0;                    // cm^3
    double hm_dens = 0.8815 * dens;       // gHM/cm^3
    double mol_hm = hm_dens / 235.;       // moles/cm^3 HM
    double atom_hm = mol_hm * at_mol;     // atoms/cm^3 HM
    double abcm_hm = atom_hm * cm2_barn;  // atoms/barn*cm HM
    double amt_hm = abcm_hm * vol;        // atoms*cm^2/barn HM

    // library variables
    std::string name0;
    int id0;
    get_name_id( 5, name0, id0 );
    Origen::SP_Material mat(
        new Origen::Material( "GENERAL", lib, name0, id0, vol ) );
    jDebugTaggedLine( "material library assoc finished..." );

    // BOS numden
    std::vector<int> id{92235};
    std::vector<double> num{amt_hm};
    mat->set_amount_bos( num, id );
    jDebugTaggedLine( "set begin-of-step..." );

    // std::cout << mat->to_string() << std::endl;

    jDebugTaggedLine( "leaving setup_material1" );
    return mat;
}

Origen::SP_Material setup_material2( Origen::SP_Library lib )
{
    // Don't be alarmed, this is all just approximate.
    double cm2_barn = 1e-24;              // cm^2/barn
    double at_mol = 6.022e23;             // at/mol
    double dens = 10.4;                   // g/cm^3
    double vol = 7.5;                     // cm^3
    double hm_dens = 0.8815 * dens;       // gHM/cm^3
    double mol_hm = hm_dens / 238.;       // moles/cm^3 HM
    double atom_hm = mol_hm * at_mol;     // atoms/cm^3 HM
    double abcm_hm = atom_hm * cm2_barn;  // atoms/barn*cm HM
    double amt_hm = abcm_hm * vol;        // atoms*cm^2/barn HM

    // library variables
    std::string name0;
    int id0;
    get_name_id( 4, name0, id0 );
    Origen::SP_Material mat(
        new Origen::Material( "GENERAL", lib, name0, id0, vol ) );
    jDebugTaggedLine( "material library assoc finished..." );

    // BOS numden
    std::vector<int> id{92235, 92238};
    std::vector<double> num{0.1 * amt_hm, 0.9 * amt_hm};
    mat->set_amount_bos( num, id );

    jDebugTaggedLine( "set begin-of-step..." );

    // std::cout << mat->to_string() << std::endl;

    jDebugTaggedLine( "leaving setup_material2" );
    return mat;
}

void build_simple_data1( Origen::SP_Library lib,
                         DataContainer& dc,
                         Origen::SP_MaterialMap matmap )
{
    dc.addDouble( "backgroundPower", 0.0 );
    dc.addDouble( "systemPower", 0.0 );
    dc.addBool( "depleteByPower", false );
    dc.addInt( "numSubsteps", 1 );

    // setup material 1
    matmap->add( setup_material1( lib ) );

    jDebugTaggedLine( "leaving build_simple_data1..." );
}

void build_simple_data2( Origen::SP_Library lib,
                         DataContainer& dc,
                         Origen::SP_MaterialMap matmap )
{
    build_simple_data1( lib, dc, matmap );

    // setup material 2
    matmap->add( setup_material2( lib ) );

    jDebugTaggedLine( "leaving build_simple_data2..." );
}

class MultiZoneDepleterTester : public ::testing::Test
{
  protected:
    // data
    Origen::SP_TransitionMatrixP trx;
    Origen::SP_Library lib;
    DataContainer dc;
    Origen::MultiZoneDepleter mzd;
    Origen::SP_MaterialMap matmap;
    Origen::SP_Material mat4, mat5;
    Standard::SP_Communicator comm;

    bool created;

    void SetUp()
    {
        // other stuff
        matmap = Origen::SP_MaterialMap( new Origen::MaterialMap() );

        // get a library
        lib = std::make_shared<Origen::Library>();
        created = Origen::FakeFactory::Library_scale_pwr( *lib );
        trx = lib->newsp_transition_matrix_at( 0 );

        // set solver
        Origen::SP_Solver fake( new Origen::Solver_Fake() );
        mzd.set_solver( fake );

        // create a dc and matmap
        build_simple_data2( lib, dc, matmap );

        // comm = world.split(0); //make a comm with ALL processes
        // mzd.set_comm(comm);
        mzd.set_comm( &world );

        // associate materials
        ASSERT_FALSE( matmap == nullptr );
        mzd.set_materials( matmap );

        // make sure we have mass
        ASSERT_LT( 0.0, matmap->initial_mass() );
        ASSERT_LT( 0.0, matmap->initial_hm_mass() );
        // std::cout << matmap->to_string() << std::endl;

        std::string name0;
        int id0;

        get_name_id( 4, name0, id0 );
        mat4 = matmap->get( id0 );
        ASSERT_EQ( mat4->name(), name0 );
        nemesis::log_local() << name0 << endl;

        get_name_id( 5, name0, id0 );
        mat5 = matmap->get( id0 );
        ASSERT_EQ( mat5->name(), name0 );

        ASSERT_TRUE( mat4 != nullptr );
        ASSERT_TRUE( mat5 != nullptr );

        jDebugTaggedLine( "finished SetUp" );
    }

    void WrapUp()
    {
        dc.deepDelete();
        jDebugTaggedLine( "finished WrapUp" );
    }

    void print( std::ostream& os = std::cout )
    {
        os << matmap->to_string() << std::endl;
    }
};

TEST_F( MultiZoneDepleterTester, DepleteByFlux )
{
    if( !created ) return;

    // set data for this step
    mzd.add_step( 86400.0 * 10 );
    mzd.update( dc );
    mzd.set_substep_dtrel( Origen::Vec_Dbl( 10, 0.1 ) );

    // set data on materials
    mat4->set_transition_matrix( trx );
    mat4->set_flux( 1.0e14 );
    mat5->set_transition_matrix( trx );
    mat5->set_flux( 2.0e14 );

    // execute
    ASSERT_TRUE( mzd.execute() );

    // extract material map
    matmap = mzd.sp_materials();
    ASSERT_FALSE( matmap == nullptr );
    print();
}

TEST_F( MultiZoneDepleterTester, DepleteByPower )
{
    if( !created ) return;

    double stepTime = 86400.0;         // 1 day
    double systemPower = 4400 + 2715;  // material1+material2 assuming ~40 W/g

    // each worker prints
    std::string filename = nprintf( "matmap_rank%02d.txt", world.rank() );
    std::ofstream ofs( filename );

    int nnodes = world.size();
    size_t nstep = 1;
    for( size_t step = 0; step < nstep; ++step )
    {
        // add a new step if beyond the first (added in testing class)
        mzd.add_step( stepTime );

        // Add new step data.
        mat4->set_transition_matrix( trx );
        mat4->set_flux( 3.0e-1 );

        mat5->set_transition_matrix( trx );
        mat5->set_flux( 6.0e-1 );

        // Update with new step data (don't really have to because we aren't
        // changing anything on the dc, but it's how this will be
        // executed in most cases.
        dc.addInt( "numSubsteps", 1 );
        dc.addBool( "depleteByPower", true );
        dc.addDouble( "systemPower", systemPower );  // Watts
        mzd.update( dc );
        mzd.set_constraint_its_max( 2 );

        // check instantaneous powers at bos (trel=0.0)
        EXPECT_NEAR( 5.29545e-12, 0.3 * mat4->power_factor_bos(), 1e-17 );
        EXPECT_NEAR( 1.65698e-10, 0.6 * mat5->power_factor_bos(), 1e-15 );

        // Solve the problem.
        mzd.execute();
        Origen::Vec_Dbl scaling = mzd.scaling_factors();
        ASSERT_EQ( nstep + 1, scaling.size() );
        EXPECT_NEAR( 4.16098e13, scaling[0], 1e8 );

        // Get status.
        const ScaleUtils::IO::DB& stat = mzd.status();

        // Output status to a file.
        ofs << mzd.status() << std::endl;

        // We are replicating materials on all nodes.

        // So with power held constant we can check the following:
        // - local power will scale down proportional to nodes
        EXPECT_FLOAT_EQ( systemPower / nnodes,
                         stat.get( "localMaterialTotalPower", -1.0 ) );
        // - local flux will scale down proportional to nodes
        EXPECT_FLOAT_EQ( 2.4277535e+13 / nnodes,
                         stat.get( "localMaterialAverageFlux", -1.0 ) );
        // - local volume will stay the same
        EXPECT_FLOAT_EQ( 19.5, stat.get( "localMaterialTotalVolume", -1.0 ) );
        // - material 4 will always have half the flux of material 5
        EXPECT_FLOAT_EQ( 0.5, mat4->flux() / mat5->flux() );
        // - global volume will scale up proportionally
        EXPECT_FLOAT_EQ( 19.5 * nnodes,
                         stat.get( "globalMaterialTotalVolume", -1.0 ) );
        // - global average flux will scale down proportionally
        EXPECT_FLOAT_EQ( 2.4277535e+13 / nnodes,
                         stat.get( "globalMaterialAverageFlux", -1.0 ) );
        // - global total power will be constant
        EXPECT_FLOAT_EQ( systemPower,
                         stat.get( "globalMaterialTotalPower", -1.0 ) );
    }

    // print of the whole map
    print( ofs );
}

typedef Origen::Vec_Dbl Vec_Dbl;

class Stepper
{
  public:
    virtual void initialize( Vec_Str ) = 0;  // keys to check for convergence

    // iteration control
    virtual bool done() const = 0;           // check in while loop
    virtual void begin() = 0;                // begin new step
    virtual bool done_resp_its() const = 0;  // uses internal logic to determine
    virtual void begin_resp_its() = 0;  // begin new response iteration loop
    virtual void set_done_resp_its( bool ) = 0;  // overrides any internal logic

    // iteration query
    virtual size_t resp_its_max() const = 0;
    virtual double resp_its_tol() const = 0;

    // time values
    virtual size_t step() const = 0;
    virtual double time0() const = 0;
    virtual double time1() const = 0;
    virtual double dtime() const = 0;

    // update and get what was set
    virtual size_t resp_its_count() const = 0;
    virtual void update_resp( const Vec_Dbl& ) = 0;
    virtual void get_last_resp( Vec_Dbl* ) const = 0;
    virtual void get_resp( Vec_Dbl* ) const = 0;
    virtual void get_delta_resp( Vec_Dbl* ) const = 0;
    virtual void get_rel_delta_resp( Vec_Dbl* ) const = 0;

    // guess value (returns true if interpolating/false if extrapolating)
    virtual bool guess_r( double t, Vec_Dbl* ) const = 0;
    virtual bool guess_dtrel( Vec_Dbl* dtrel ) const = 0;
    virtual size_t guess_dtrel_size() const = 0;

  private:
    friend std::ostream& operator<<( std::ostream& os, const Stepper& s );
};
typedef std::shared_ptr<Stepper> SP_Stepper;

::std::ostream& operator<<(::std::ostream& os, const Stepper& s )
{
    return os << " tol=" << s.resp_its_tol() << " max=" << s.resp_its_max()
              << " dtrel.size()=" << s.guess_dtrel_size();
}

::std::ostream& operator<<(::std::ostream& os, const SP_Stepper& s )
{
    return os << *s;
}

class Stepper_fixed : public virtual Stepper
{
  private:
    Stepper_fixed() {}
  public:
    virtual ~Stepper_fixed() {}
    virtual void initialize( Vec_Str keys )
    {
        Require( keys.size() > 0 );
        b_keys = keys;
        b_resp.resize( b_keys.size(), 0.0 );
        b_resp_last.resize( b_keys.size(), 0.0 );
    }
    Stepper_fixed( Vec_Dbl dtrel, int its_max, double its_tol )
    {
        b_dtrel = dtrel;
        b_resp_its_max = its_max;
        b_resp_its_tol = its_tol;
        b_s = 0;
        b_t0 = 0;
    }

    bool done() const
    {
        Require( b_keys.size() > 0 );
        return !( b_s < b_dtrel.size() );
    }
    void begin()
    {
        ++b_s;
        b_t1 = b_t0 + b_dtrel[b_s - 1];
        b_resp_its_count = 0;
        b_resp_done = false;
    }
    bool done_resp_its() const
    {
        // quick return
        if( b_resp_its_count > b_resp_its_max ) b_resp_done = true;
        if( b_resp_done ) return true;

        Vec_Dbl rel_dr;
        get_delta_resp( &rel_dr );
        for( auto r : rel_dr )
        {
            if( r > b_resp_its_tol ) return false;
        }
        b_resp_done = true;
        return true;
    }
    void begin_resp_its() { ++b_resp_its_count; }
    void set_done_resp_its( bool done ) { b_resp_done = done; }
    size_t resp_its_max() const { return b_resp_its_max; }
    double resp_its_tol() const { return b_resp_its_tol; }
    size_t step() const { return b_s - 1; }
    double time0() const { return b_t0; }
    double time1() const { return b_t1; }
    double dtime() const { return b_t1 - b_t0; }
    size_t resp_its_count() const { return b_resp_its_count; }
    void update_resp( const Vec_Dbl& resp )
    {
        Require( resp.size() == b_keys.size() );
        b_resp_last = b_resp;
        b_resp = resp;
    }
    void get_last_resp( Vec_Dbl* resp ) const { *resp = b_resp_last; }
    void get_resp( Vec_Dbl* resp ) const { *resp = b_resp; }
    void get_delta_resp( Vec_Dbl* resp ) const
    {
        Require( resp->size() == b_keys.size() );
        for( size_t i = 0; i < b_keys.size(); ++i )
        {
            ( *resp )[i] = b_resp[i] - b_resp_last[i];
        }
    }
    void get_rel_delta_resp( Vec_Dbl* resp ) const
    {
        Require( resp->size() == b_keys.size() );
        for( size_t i = 0; i < b_keys.size(); ++i )
        {
            ( *resp )[i] = b_resp[i] - b_resp_last[i];
            if( b_resp_last[i] != 0.0 )
            {
                ( *resp )[i] /= b_resp_last[i];
            }
        }
    }

    // guess value (returns true if interpolating/false if extrapolating)
    bool guess_r( double t, Vec_Dbl* resp ) const
    {
        *resp = b_resp;
        return t == b_t1;
    }
    // true if known/false if guessing
    bool guess_dtrel( Vec_Dbl* dtrel ) const
    {
        *dtrel = b_dtrel;
        return true;
    }
    size_t guess_dtrel_size() const { return b_dtrel.size(); }
  protected:
    Vec_Dbl b_dtrel;
    Vec_Dbl b_resp;
    Vec_Dbl b_resp_last;
    Vec_Str b_keys;
    size_t b_s;
    double b_t0;
    double b_t1;
    int b_resp_its_count;
    int b_resp_its_max;
    double b_resp_its_tol;
    mutable bool b_resp_done;
};

class MultiZoneDepleterTesterP
    : public MultiZoneDepleterTester,
      public virtual ::testing::WithParamInterface<SP_Stepper>
{
  protected:
    void SetUp()
    {
        MultiZoneDepleterTester::SetUp();
        ASSERT_FALSE( matmap == nullptr );
        stepper = GetParam();
    }
    SP_Stepper stepper;
};

// same-size steps, step-constant power constraint
TEST_P( MultiZoneDepleterTesterP, PowerConstraint0 )
{
    if( !created ) return;

    // time steps
    Origen::Vec_Dbl sys_dt( 10, 86400. * 10 );  // 10 day periods

    // power over each step
    Origen::Vec_Dbl sys_p( 10, 7000 );
    {
        size_t c = 0;
        for( auto& p : sys_p )
        {
            p *= ( c + 10. ) / 10.;
            ++c;
        }
    }

    Vec_Dbl dtrel;
    for( size_t n = 0; n < sys_dt.size(); ++n )
    {
        // Allocate memory on all materials for new step.
        // Set global time step (all materials).
        mzd.add_step( sys_dt[n] );

        // Add new step data for each material.
        mat4->set_transition_matrix( trx );
        mat4->set_flux( 0.6 );

        mat5->set_transition_matrix( trx );
        mat5->set_flux( 0.4 );

        // Set global power constraint
        mzd.set_power_constraint( sys_p[n] );

        // Set substep parameters.  Eventually we want to move this into
        // MZD such that it has its own stepper.
        EXPECT_TRUE( stepper->guess_dtrel( &dtrel ) );
        mzd.set_substep_dtrel( dtrel );
        mzd.set_constraint_its_max( stepper->resp_its_max() );
        mzd.set_constraint_its_tol( stepper->resp_its_tol() );

        // Solve the problem.
        mzd.execute();

        // Get & check status.
        const ScaleUtils::IO::DB& stat = mzd.status();

        EXPECT_EQ( stepper->guess_dtrel_size(),
                   stat.get<int>( "numSubsteps" ) );
        EXPECT_EQ( true, stat.get<bool>( "hasConstraint" ) );
        EXPECT_EQ( mzd.dt(), stat.get<double>( "stepTime" ) );
        EXPECT_EQ( "POWER", stat.get<std::string>( "constraintType" ) );
        EXPECT_EQ( sys_p[n], stat.get<double>( "constraintValue" ) );
        EXPECT_EQ( stepper->resp_its_max(),
                   stat.get<int>( "constraintItsMax" ) );
        EXPECT_EQ( stepper->resp_its_tol(),
                   stat.get<double>( "constraintItsTol" ) );

        // check power constraint is met
        EXPECT_FLOAT_EQ( sys_p[n], matmap->power() );
    }
}

const Vec_Dbl dtrel_1( {1.0} );
const Vec_Dbl dtrel_3( {0.1, 0.2, 0.7} );

const std::vector<SP_Stepper> steppers = {
    std::make_shared<Stepper_fixed>(
        dtrel_1, 1, 1e-5 ),  // single step with 1 response iteration
    std::make_shared<Stepper_fixed>(
        dtrel_1, 2, 1e-5 ),  // single step with 2 response iterations
    std::make_shared<Stepper_fixed>(
        dtrel_1,
        2,
        1e-2 ),  // single step with 2 response iterations & reduced tolerance
    std::make_shared<Stepper_fixed>(
        dtrel_3, 1, 1e-5 ),  // multiple steps like above
    std::make_shared<Stepper_fixed>(
        dtrel_3, 2, 1e-5 ),  // single step with 2 response iterations
    std::make_shared<Stepper_fixed>(
        dtrel_3,
        20,
        1e-5 ),  // single step with 2 response iterations & reduced tolerance
};

INSTANTIATE_TEST_CASE_P( StepperTesting,
                         MultiZoneDepleterTesterP,
                         ::testing::ValuesIn( steppers ) );