first_order_reaction_base.cc

Go to the documentation of this file.

/*!
 *
 * Copyright (C) 2015 Technical University of Liberec.  All rights reserved.
 * 
 * This program is free software; you can redistribute it and/or modify it under
 * the terms of the GNU General Public License version 3 as published by the
 * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
 * 
 * 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 General Public License for more details.
 *
 * 
 * @file    first_order_reaction_base.cc
 * @brief   
 */

#include "reaction/first_order_reaction_base.hh"
#include "reaction/reaction_term.hh"

#include "reaction/linear_ode_solver.hh"


#include "system/global_defs.h"
#include "system/sys_profiler.hh"

#include "mesh/mesh.h"
#include "la/distribution.hh"
#include "input/accessors.hh"



using namespace Input::Type;


FirstOrderReactionBase::FirstOrderReactionBase(Mesh &init_mesh, Input::Record in_rec)
    : ReactionTerm(init_mesh, in_rec)
{
    linear_ode_solver_ = std::make_shared<LinearODESolver>();
}

FirstOrderReactionBase::~FirstOrderReactionBase()
{
}

void FirstOrderReactionBase::initialize()
{
    ASSERT(time_ != nullptr).error("Time governor has not been set yet.\n");
    ASSERT_LT(0, substances_.size()).error("No substances for rection term.\n");
    
    n_substances_ = substances_.size();
    initialize_from_input();

    // allocation
    prev_conc_.resize(n_substances_);
    reaction_matrix_.resize(n_substances_, n_substances_);
    molar_matrix_.resize(n_substances_, n_substances_);
    molar_mat_inverse_.resize(n_substances_, n_substances_);

    // initialize diagonal matrices with molar masses
    molar_matrix_.zeros();
    molar_mat_inverse_.zeros();
    for (unsigned int i=0; i<n_substances_; ++i)
    {
        molar_matrix_(i,i) = substances_[i].molar_mass();
        molar_mat_inverse_(i,i) = 1./substances_[i].molar_mass();
    }
}


void FirstOrderReactionBase::zero_time_step()
{
    ASSERT(time_ != nullptr).error("Time governor has not been set yet.\n");
    ASSERT_LT(0, substances_.size()).error("No substances for rection term.\n");

    assemble_ode_matrix();
    // make scaling that takes into account different molar masses of substances
    reaction_matrix_ = molar_matrix_ * reaction_matrix_ * molar_mat_inverse_;
    
    linear_ode_solver_->set_system_matrix(reaction_matrix_);
}


void FirstOrderReactionBase::compute_reaction(const DHCellAccessor& dh_cell)
{      
    unsigned int sbi;  // row in the concentration matrix, regards the substance index
    arma::vec new_conc;
    
    IntIdx dof_p0 = dh_cell.get_loc_dof_indices()[0];

    // save previous concentrations to column vector
    for(sbi = 0; sbi < n_substances_; sbi++)
        prev_conc_(sbi) = conc_mobile_fe[sbi]->vec()[dof_p0];
    
    // compute new concetrations R*c
    linear_ode_solver_->update_solution(prev_conc_, new_conc);
    
    // save new concentrations to the concentration matrix
    for(sbi = 0; sbi < n_substances_; sbi++)
        conc_mobile_fe[sbi]->vec()[dof_p0] = new_conc(sbi);
}

void FirstOrderReactionBase::update_solution(void)
{
    //DebugOut() << "FirstOrderReactionBases - update solution\n";
    if(time_->is_changed_dt())
    {
        linear_ode_solver_->set_step(time_->dt());
    }

    START_TIMER("linear reaction step");

    for ( DHCellAccessor dh_cell : dof_handler_->own_range() )
    {
        compute_reaction(dh_cell);
    }
    END_TIMER("linear reaction step");
}


unsigned int FirstOrderReactionBase::find_subst_name(const string &name)
{
    unsigned int k=0;
        for(; k < n_substances_; k++)
                if (name == substances_[k].name()) return k;

        return k;
}


bool FirstOrderReactionBase::evaluate_time_constraint(double &time_constraint)
{
    if (!linear_ode_solver_->evaluate_time_constraint(time_constraint)) return false;
    
    DebugOut().fmt("CFL constraint(first order reaction): {}.\n", time_constraint);
    
    return true;
}