Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "AddCoupledSolidKinSpeciesAction.h"
      11             : #include "MooseUtils.h"
      12             : #include "FEProblem.h"
      13             : #include "Factory.h"
      14             : #include "MooseError.h"
      15             : #include "Parser.h"
      16             : #include <algorithm>
      17             : 
      18             : // Regular expression includes
      19             : #include "pcrecpp.h"
      20             : 
      21             : registerMooseAction("ChemicalReactionsApp", AddCoupledSolidKinSpeciesAction, "add_kernel");
      22             : 
      23             : registerMooseAction("ChemicalReactionsApp", AddCoupledSolidKinSpeciesAction, "add_aux_kernel");
      24             : 
      25             : InputParameters
      26          41 : AddCoupledSolidKinSpeciesAction::validParams()
      27             : {
      28          41 :   InputParameters params = Action::validParams();
      29          82 :   params.addRequiredParam<std::vector<NonlinearVariableName>>("primary_species",
      30             :                                                               "The list of primary species to add");
      31          82 :   params.addParam<std::vector<AuxVariableName>>(
      32             :       "secondary_species", "The list of solid kinetic species to be output as aux variables");
      33          82 :   params.addRequiredParam<std::string>("kin_reactions", "The list of solid kinetic reactions");
      34          82 :   params.addRequiredParam<std::vector<Real>>("log10_keq",
      35             :                                              "The list of equilibrium constants for all reactions");
      36          82 :   params.addRequiredParam<std::vector<Real>>(
      37             :       "specific_reactive_surface_area",
      38             :       "The list of specific reactive surface area for all minerals (m^2/L)");
      39          82 :   params.addRequiredParam<std::vector<Real>>(
      40             :       "kinetic_rate_constant", "The list of kinetic rate constant for all reactions (mol/m^2/s)");
      41          82 :   params.addRequiredParam<std::vector<Real>>(
      42             :       "activation_energy", "The list of activation energy values for all reactions (J/mol)");
      43          82 :   params.addParam<Real>("gas_constant", 8.314, "Gas constant. Default is 8.314 (J/mol/K)");
      44          82 :   params.addRequiredParam<std::vector<Real>>(
      45             :       "reference_temperature", "The list of reference temperatures for all reactions (K)");
      46          82 :   params.addRequiredCoupledVar("system_temperature",
      47             :                                "The system temperature for all reactions (K)");
      48          41 :   params.addClassDescription("Adds solid kinetic Kernels and AuxKernels for primary species");
      49          41 :   return params;
      50           0 : }
      51             : 
      52          41 : AddCoupledSolidKinSpeciesAction::AddCoupledSolidKinSpeciesAction(const InputParameters & params)
      53             :   : Action(params),
      54          41 :     _primary_species(getParam<std::vector<NonlinearVariableName>>("primary_species")),
      55         123 :     _secondary_species(getParam<std::vector<AuxVariableName>>("secondary_species")),
      56          41 :     _kinetic_species_involved(_primary_species.size()),
      57          41 :     _weights(_primary_species.size()),
      58         123 :     _input_reactions(getParam<std::string>("kin_reactions")),
      59          82 :     _logk(getParam<std::vector<Real>>("log10_keq")),
      60          82 :     _r_area(getParam<std::vector<Real>>("specific_reactive_surface_area")),
      61          82 :     _ref_kconst(getParam<std::vector<Real>>("kinetic_rate_constant")),
      62          82 :     _e_act(getParam<std::vector<Real>>("activation_energy")),
      63          82 :     _gas_const(getParam<Real>("gas_constant")),
      64          82 :     _ref_temp(getParam<std::vector<Real>>("reference_temperature")),
      65         164 :     _sys_temp(getParam<std::vector<VariableName>>("system_temperature"))
      66             : {
      67             :   // Note: as the reaction syntax has changed, check to see if the old syntax has
      68             :   // been used and throw an informative error. The number of = signs should be one
      69             :   // more than the number of commas, while the smallest number of spaces possible is 2
      70             :   bool old_syntax = false;
      71          41 :   if (std::count(_input_reactions.begin(), _input_reactions.end(), '=') !=
      72          41 :       std::count(_input_reactions.begin(), _input_reactions.end(), ',') + 1)
      73             :     old_syntax = true;
      74             : 
      75          41 :   if (std::count(_input_reactions.begin(), _input_reactions.end(), ' ') < 2)
      76             :     old_syntax = true;
      77             : 
      78          41 :   if (old_syntax)
      79           0 :     mooseError("Old solid kinetic reaction syntax present.\nReactions should now be comma "
      80             :                "separated, and must have spaces between species and +/-/= operators.\n"
      81             :                "See #9972 for details");
      82             : 
      83             :   // Parse the kinetic reactions
      84             :   pcrecpp::RE re_reactions("(.+?)" // A single reaction (any character until the comma delimiter)
      85             :                            "(?:,\\s*|$)" // comma or end of string
      86             :                            ,
      87          41 :                            pcrecpp::RE_Options().set_extended(true));
      88             : 
      89          41 :   pcrecpp::RE re_terms("(\\S+)");
      90             :   pcrecpp::RE re_coeff_and_species("(?: \\(? (.*?) \\)? )" // match the leading coefficent
      91             :                                    "([A-Za-z].*)"          // match the species
      92             :                                    ,
      93          41 :                                    pcrecpp::RE_Options().set_extended(true));
      94             : 
      95             :   pcrecpp::StringPiece input(_input_reactions);
      96             :   pcrecpp::StringPiece single_reaction, term;
      97             :   std::string single_reaction_str;
      98             : 
      99             :   // Parse reaction network to extract each individual reaction
     100         100 :   while (re_reactions.FindAndConsume(&input, &single_reaction_str))
     101          59 :     _reactions.push_back(single_reaction_str);
     102             : 
     103          41 :   _num_reactions = _reactions.size();
     104             : 
     105          41 :   if (_num_reactions == 0)
     106           0 :     mooseError("No solid kinetic reaction provided!");
     107             : 
     108             :   // Start parsing each reaction
     109         100 :   for (unsigned int i = 0; i < _num_reactions; ++i)
     110             :   {
     111          59 :     single_reaction = _reactions[i];
     112             : 
     113             :     // Capture all of the terms
     114             :     std::string species, coeff_str;
     115             :     Real coeff;
     116             :     int sign = 1;
     117             :     bool secondary = false;
     118             : 
     119             :     std::vector<Real> local_stos;
     120             :     std::vector<VariableName> local_species_list;
     121             : 
     122             :     // Find every single term in this reaction (species and operators)
     123         406 :     while (re_terms.FindAndConsume(&single_reaction, &term))
     124             :     {
     125             :       // Separating the stoichiometric coefficients from species
     126         347 :       if (re_coeff_and_species.PartialMatch(term, &coeff_str, &species))
     127             :       {
     128         203 :         if (coeff_str.length())
     129          45 :           coeff = std::stod(coeff_str);
     130             :         else
     131         158 :           coeff = 1.0;
     132             : 
     133         203 :         coeff *= sign;
     134             : 
     135         203 :         if (secondary)
     136         118 :           _solid_kinetic_species.push_back(species);
     137             :         else
     138             :         {
     139         144 :           local_stos.push_back(coeff);
     140         144 :           local_species_list.push_back(species);
     141             :         }
     142             :       }
     143             :       // Finding the operators and assign value of -1.0 to "-" sign
     144         144 :       else if (term == "+" || term == "=" || term == "-")
     145             :       {
     146         144 :         if (term == "-")
     147             :         {
     148             :           sign = -1;
     149          35 :           term = "+";
     150             :         }
     151             : 
     152         144 :         if (term == "=")
     153             :           secondary = true;
     154             :       }
     155             :       else
     156           0 :         mooseError("Error parsing term: ", term.as_string());
     157             :     }
     158             : 
     159          59 :     _stos.push_back(local_stos);
     160          59 :     _primary_species_involved.push_back(local_species_list);
     161          59 :   }
     162             : 
     163             :   // Start picking out primary species and coupled primary species and assigning
     164             :   // corresponding stoichiomentric coefficients
     165         167 :   for (unsigned int i = 0; i < _primary_species.size(); ++i)
     166         324 :     for (unsigned int j = 0; j < _num_reactions; ++j)
     167             :     {
     168         672 :       for (unsigned int k = 0; k < _primary_species_involved[j].size(); ++k)
     169         474 :         if (_primary_species_involved[j][k] == _primary_species[i])
     170             :         {
     171         144 :           _weights[i].push_back(_stos[j][k]);
     172         144 :           _kinetic_species_involved[i].push_back(_solid_kinetic_species[j]);
     173             :         }
     174             :     }
     175             : 
     176             :   // Print out details of the solid kinetic reactions to the console
     177          41 :   _console << "Solid kinetic reactions:\n";
     178         100 :   for (unsigned int i = 0; i < _num_reactions; ++i)
     179          59 :     _console << "  Reaction " << i + 1 << ": " << _reactions[i] << "\n";
     180          41 :   _console << std::endl;
     181             : 
     182             :   // Check that all secondary species read from the reaction network have been added
     183             :   // as AuxVariables. Note: can't sort the _solid_kinetic_species vector as it throws
     184             :   // out the species and coefficient vectors so use std::is_permutation
     185          41 :   if (!std::is_permutation(
     186             :           _secondary_species.begin(), _secondary_species.end(), _solid_kinetic_species.begin()))
     187           0 :     mooseError("All solid kinetic species must be added as secondary species");
     188             : 
     189             :   // Check that the size of property vectors is equal to the number of reactions
     190          41 :   if (_logk.size() != _num_reactions)
     191           0 :     mooseError("The number of values entered for log10_keq is not equal to the number of solid "
     192             :                "kinetic reactions");
     193          41 :   if (_r_area.size() != _num_reactions)
     194           0 :     mooseError("The number of values entered for specific_reactive_surface_area is not equal to "
     195             :                "the number of solid kinetic reactions");
     196          41 :   if (_ref_kconst.size() != _num_reactions)
     197           0 :     mooseError("The number of values entered for kinetic_rate_constant is not equal to the number "
     198             :                "of solid kinetic reactions");
     199          41 :   if (_e_act.size() != _num_reactions)
     200           0 :     mooseError("The number of values entered for activation_energy is not equal to the number of "
     201             :                "solid kinetic reactions");
     202          41 :   if (_ref_temp.size() != _num_reactions)
     203           0 :     mooseError("The number of values entered for reference_temperature is not equal to the number "
     204             :                "of solid kinetic reactions");
     205          41 : }
     206             : 
     207             : void
     208          82 : AddCoupledSolidKinSpeciesAction::act()
     209             : {
     210          82 :   if (_current_task == "add_kernel")
     211             :   {
     212             :     // Add Kernels for each primary species
     213         167 :     for (unsigned int i = 0; i < _primary_species.size(); ++i)
     214             :     {
     215         252 :       InputParameters params_kin = _factory.getValidParams("CoupledBEKinetic");
     216         252 :       params_kin.set<NonlinearVariableName>("variable") = _primary_species[i];
     217         252 :       params_kin.set<std::vector<Real>>("weight") = _weights[i];
     218         126 :       params_kin.set<std::vector<VariableName>>("v") = _kinetic_species_involved[i];
     219         378 :       _problem->addKernel("CoupledBEKinetic", _primary_species[i] + "_" + "_kin", params_kin);
     220         126 :     }
     221             :   }
     222             : 
     223          82 :   if (_current_task == "add_aux_kernel")
     224             :   {
     225             :     // Add AuxKernels for each solid kinetic species
     226         100 :     for (unsigned int i = 0; i < _num_reactions; ++i)
     227             :     {
     228          59 :       InputParameters params_kin = _factory.getValidParams("KineticDisPreConcAux");
     229         118 :       params_kin.set<AuxVariableName>("variable") = _solid_kinetic_species[i];
     230         118 :       params_kin.defaultCoupledValue("log_k", _logk[i]);
     231          59 :       params_kin.set<Real>("r_area") = _r_area[i];
     232          59 :       params_kin.set<Real>("ref_kconst") = _ref_kconst[i];
     233          59 :       params_kin.set<Real>("e_act") = _e_act[i];
     234          59 :       params_kin.set<Real>("gas_const") = _gas_const;
     235          59 :       params_kin.set<Real>("ref_temp") = _ref_temp[i];
     236         118 :       params_kin.set<std::vector<VariableName>>("sys_temp") = _sys_temp;
     237         118 :       params_kin.set<std::vector<Real>>("sto_v") = _stos[i];
     238          59 :       params_kin.set<std::vector<VariableName>>("v") = _primary_species_involved[i];
     239         118 :       _problem->addAuxKernel(
     240          59 :           "KineticDisPreConcAux", "aux_" + _solid_kinetic_species[i], params_kin);
     241          59 :     }
     242             :   }
     243          82 : }