DiffusionCG.C

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
 
#include "DiffusionCG.h"
#include "MooseVariableBase.h"
 
// Register the actions for the objects actually used
registerDiffusionPhysicsBaseTasks("MooseApp", DiffusionCG);
registerMooseAction("MooseApp", DiffusionCG, "add_kernel");
registerMooseAction("MooseApp", DiffusionCG, "add_bc");
registerMooseAction("MooseApp", DiffusionCG, "add_variables_physics");
 
InputParameters
DiffusionCG::validParams()
{
  InputParameters params = DiffusionPhysicsBase::validParams();
  params.addClassDescription(
      "Discretizes a diffusion equation with the continuous Galerkin finite element method");
  params.addParam<bool>(
      "use_automatic_differentiation",
      true,
      "Whether to use automatic differentiation for all the terms in the equation");
  params.transferParam<MooseEnum>(MooseVariableBase::validParams(), "order", "variable_order");
  params.addParam<MaterialPropertyName>("source_matprop", "Source term in the diffusion problem");
 
  return params;
}
 
DiffusionCG::DiffusionCG(const InputParameters & parameters)
  : PhysicsBase(parameters),
    PhysicsComponentInterface(parameters),
    DiffusionPhysicsBase(parameters),
    _use_ad(getParam<bool>("use_automatic_differentiation"))
{
  checkAtMostOneParamSetByUser({"source_functor", "source_matprop"});
}
 
void
DiffusionCG::addFEKernels()
{
  // Diffusion term(s)
  unsigned int num_diffusion_terms = 1;
  if (isParamValid("diffusivity_matprop"))
    num_diffusion_terms = getParam<std::vector<MaterialPropertyName>>("diffusivity_matprop").size();
  else if (isParamValid("diffusivity_functor"))
    num_diffusion_terms = getParam<std::vector<MooseFunctorName>>("diffusivity_functor").size();
 
  for (const auto i_diff_block_group : make_range(num_diffusion_terms))
  {
    // Select the kernel type based on the user parameters
    std::string kernel_type;
    if (isParamValid("diffusivity_matprop"))
      kernel_type = _use_ad ? "ADMatDiffusion" : "MatDiffusion";
    else if (isParamValid("diffusivity_functor"))
    {
      const auto & d =
          getParam<std::vector<MooseFunctorName>>("diffusivity_functor")[i_diff_block_group];
      if (getProblem().hasFunction(d) || MooseUtils::parsesToReal(d))
        kernel_type = "FunctionDiffusion";
      else
        paramError("diffusivity_functor",
                   "No diffusion kernel implemented for the diffusivity type of ",
                   d);
    }
    else
      kernel_type = _use_ad ? "ADDiffusion" : "Diffusion";
    InputParameters params = getFactory().getValidParams(kernel_type);
    params.set<NonlinearVariableName>("variable") = _var_name;
    if (isParamValid("diffusivity_blocks"))
      assignBlocks(params,
                   getParam<std::vector<std::vector<SubdomainName>>>(
                       "diffusivity_blocks")[i_diff_block_group]);
    else
      assignBlocks(params, _blocks);
 
    // Transfer the diffusivity parameter from the Physics to the kernel
    // From parameters
    if (isParamValid("diffusivity_matprop"))
      params.set<MaterialPropertyName>("diffusivity") =
          getParam<std::vector<MaterialPropertyName>>("diffusivity_matprop")[i_diff_block_group];
    else if (isParamValid("diffusivity_functor"))
      params.set<FunctionName>("function") =
          getParam<std::vector<MooseFunctorName>>("diffusivity_functor")[i_diff_block_group];
 
    const auto kernel_name =
        prefix() + _var_name + "_diffusion" +
        ((num_diffusion_terms > 1) ? "_" + std::to_string(i_diff_block_group) : "");
    getProblem().addKernel(kernel_type, kernel_name, params);
  }
 
  // Source term
  if (isParamValid("source_functor") || isParamValid("source_matprop"))
  {
    // Select the kernel type based on the user parameters
    std::string kernel_type;
    std::string source =
        isParamValid("source_functor") ? getParam<MooseFunctorName>("source_functor") : "";
    if (isParamValid("source_matprop"))
      source = getParam<MaterialPropertyName>("source_matprop");
 
    if (MooseUtils::parsesToReal(source) || getProblem().hasFunction(source) ||
        getProblem().hasPostprocessorValueByName(source))
      kernel_type = _use_ad ? "ADBodyForce" : "BodyForce";
    else if (getProblem().hasVariable(source))
      kernel_type = _use_ad ? "ADCoupledForce" : "CoupledForce";
    else if (getProblem().hasFunctor(source, 0))
      kernel_type = "FunctorKernel";
    else if (getProblem().getMaterialPropertyRegistry().hasProperty(source))
      kernel_type = _use_ad ? "ADMatBodyForce" : "MatBodyForce";
    else
      paramError(isParamValid("source_functor") ? "source_functor" : "source_matprop",
                 "No kernel defined for a source term in CG for the type of '",
                 source,
                 "'");
 
    InputParameters params = getFactory().getValidParams(kernel_type);
    params.set<NonlinearVariableName>("variable") = _var_name;
    assignBlocks(params, _blocks);
 
    // Transfer the source and coefficient parameter from the Physics to the kernel
    const auto coef = getParam<Real>("source_coef");
    if (MooseUtils::parsesToReal(source))
      params.set<Real>("value") = MooseUtils::convert<Real>(source) * coef;
    else if (getProblem().hasFunction(source))
    {
      params.set<Real>("value") = coef;
      params.set<FunctionName>("function") = source;
    }
    else if (getProblem().hasPostprocessorValueByName(source))
    {
      params.set<Real>("value") = coef;
      params.set<PostprocessorName>("postprocessor") = source;
    }
    else if (getProblem().hasFunctor(source, 0))
    {
      params.set<MooseFunctorName>("functor") = source;
      params.set<bool>("functor_on_rhs") = true;
      if (isParamSetByUser("source_coef"))
        paramError("source_coef",
                   "Setting a coefficient is not implemented with a source functor. Use an "
                   "intermediate (parsed) functor material to define a new functors which "
                   "multiplies the functor by the coefficient");
    }
    else if (getProblem().getMaterialPropertyRegistry().hasProperty(source))
    {
      params.set<MaterialPropertyName>("material_property") = source;
      if (isParamSetByUser("source_coef"))
        paramError(
            "source_coef",
            "Setting a coefficient is not implemented with a material property source. Use an "
            "intermediate (parsed) material to define a new functors which "
            "multiplies the property by the coefficient");
    }
    else
    {
      params.set<Real>("coef") = coef;
      params.set<std::vector<VariableName>>("v") = {source};
    }
 
    getProblem().addKernel(kernel_type, prefix() + _var_name + "_source", params);
  }
 
  // Time derivative term
  if (shouldCreateTimeDerivative(_var_name, _blocks, false))
  {
    const std::string kernel_type = _use_ad ? "ADTimeDerivative" : "TimeDerivative";
    InputParameters params = getFactory().getValidParams(kernel_type);
    params.set<NonlinearVariableName>("variable") = _var_name;
    assignBlocks(params, _blocks);
    getProblem().addKernel(kernel_type, prefix() + _var_name + "_time", params);
  }
}
 
void
DiffusionCG::addFEBCs()
{
  if (isParamSetByUser("neumann_boundaries"))
  {
    const auto & boundary_fluxes = getParam<std::vector<MooseFunctorName>>("boundary_fluxes");
    for (const auto i : index_range(_neumann_boundaries))
    {
      const auto & bc_flux = boundary_fluxes[i];
      // Select the boundary type based on the user parameters and what we know to be most efficient
      // We could actually just use the very last option for everything but the performance is
      // better if one uses the specialized objects
      std::string bc_type = "";
      if (MooseUtils::parsesToReal(bc_flux))
        bc_type = _use_ad ? "ADNeumannBC" : "NeumannBC";
      else if (getProblem().hasVariable(bc_flux))
        bc_type = "CoupledVarNeumannBC"; // not AD, but still perfect Jacobian
      else if (getProblem().hasFunction(bc_flux))
        bc_type = _use_ad ? "ADFunctionNeumannBC" : "FunctionNeumannBC";
      else if (getProblem().hasPostprocessorValueByName(bc_flux))
        bc_type = "PostprocessorNeumannBC";
      else // this is AD, but we can mix AD and non-AD
        bc_type = "FunctorNeumannBC";
 
      // Get the parameters for the object type chosen and set the common parameters
      InputParameters params = getFactory().getValidParams(bc_type);
      params.set<NonlinearVariableName>("variable") = _var_name;
      params.set<std::vector<BoundaryName>>("boundary") = {_neumann_boundaries[i]};
 
      // Set the flux parameter for the specific type of NeumannBC used
      if (MooseUtils::parsesToReal(bc_flux))
        params.set<Real>("value") = MooseUtils::convert<Real>(bc_flux);
      else if (getProblem().hasVariable(bc_flux))
        params.set<std::vector<VariableName>>("v") = {bc_flux};
      else if (getProblem().hasFunction(bc_flux))
        params.set<FunctionName>("function") = bc_flux;
      else if (getProblem().hasPostprocessorValueByName(bc_flux))
        params.set<PostprocessorName>("postprocessor") = bc_flux;
      else
        params.set<MooseFunctorName>("functor") = bc_flux;
 
      getProblem().addBoundaryCondition(
          bc_type, prefix() + _var_name + "_neumann_bc_" + _neumann_boundaries[i], params);
    }
  }
  if (isParamSetByUser("dirichlet_boundaries"))
  {
    const auto & boundary_values = getParam<std::vector<MooseFunctorName>>("boundary_values");
    for (const auto i : index_range(_dirichlet_boundaries))
    {
      const auto & bc_value = boundary_values[i];
      // Select the boundary type based on the user parameters and what we know to be most efficient
      std::string bc_type = "";
      if (MooseUtils::parsesToReal(bc_value))
        bc_type = _use_ad ? "ADDirichletBC" : "DirichletBC";
      else if (getProblem().hasVariable(bc_value))
        bc_type = _use_ad ? "ADMatchedValueBC" : "MatchedValueBC";
      else if (getProblem().hasFunction(bc_value))
        bc_type = _use_ad ? "ADFunctionDirichletBC" : "FunctionDirichletBC";
      else if (getProblem().hasPostprocessorValueByName(bc_value))
        bc_type = "PostprocessorDirichletBC";
      else // this is AD, but we can mix AD and non-AD
        bc_type = "FunctorDirichletBC";
 
      InputParameters params = getFactory().getValidParams(bc_type);
      params.set<NonlinearVariableName>("variable") = _var_name;
      params.set<std::vector<BoundaryName>>("boundary") = {_dirichlet_boundaries[i]};
 
      // Set the flux parameter for the specific type of DirichletBC used
      if (MooseUtils::parsesToReal(bc_value))
        params.set<Real>("value") = MooseUtils::convert<Real>(bc_value);
      else if (getProblem().hasVariable(bc_value))
        params.set<std::vector<VariableName>>("v") = {bc_value};
      else if (getProblem().hasFunction(bc_value))
        params.set<FunctionName>("function") = bc_value;
      else if (getProblem().hasPostprocessorValueByName(bc_value))
        params.set<PostprocessorName>("postprocessor") = bc_value;
      else
        params.set<MooseFunctorName>("functor") = bc_value;
 
      getProblem().addBoundaryCondition(
          bc_type, prefix() + _var_name + "_dirichlet_bc_" + _dirichlet_boundaries[i], params);
    }
  }
}
 
void
DiffusionCG::addSolverVariables()
{
  // If the variable was added outside the Physics
  if (!shouldCreateVariable(_var_name, _blocks, true))
  {
    reportPotentiallyMissedParameters({"variable_order", "system_names"}, "MooseVariable");
    return;
  }
 
  const std::string variable_type = "MooseVariable";
  InputParameters params = getFactory().getValidParams(variable_type);
  params.set<MooseEnum>("order") = getParam<MooseEnum>("variable_order");
  assignBlocks(params, _blocks);
  params.set<SolverSystemName>("solver_sys") = getSolverSystem(_var_name);
 
  getProblem().addVariable(variable_type, _var_name, params);
}
 
void
DiffusionCG::addBoundaryConditionsFromComponents()
{
  const auto bc_type_fixed = "FunctorDirichletBC";
  InputParameters params_fixed = getFactory().getValidParams(bc_type_fixed);
  const auto bc_type_flux = "FunctorNeumannBC";
  InputParameters params_flux = getFactory().getValidParams(bc_type_flux);
 
  for (const auto & comp_pair : _components_boundary_conditions)
  {
    for (const auto & bc_pair : comp_pair.second)
    {
      params_fixed.set<NonlinearVariableName>("variable") = _var_name;
      params_flux.set<NonlinearVariableName>("variable") = _var_name;
 
      const auto bc_type = bc_pair.second.second;
      const auto functor_name = bc_pair.second.first;
      params_fixed.set<MooseFunctorName>("functor") = functor_name;
      params_flux.set<MooseFunctorName>("functor") = functor_name;
 
      const auto surface_name = bc_pair.first.second;
      params_fixed.set<std::vector<BoundaryName>>("boundary") = {surface_name};
      params_flux.set<std::vector<BoundaryName>>("boundary") = {surface_name};
 
      if (bc_type == ComponentBoundaryConditionInterface::BoundaryConditionType::FIXED_VALUE)
        getProblem().addBoundaryCondition(bc_type_fixed,
                                          prefix() + _var_name + "_comp_" + comp_pair.first +
                                              "_fixed_value_" + surface_name,
                                          params_fixed);
      else if (bc_type == ComponentBoundaryConditionInterface::BoundaryConditionType::FLUX)
        getProblem().addBoundaryCondition(bc_type_flux,
                                          prefix() + _var_name + "_comp_" + comp_pair.first +
                                              "_flux_" + surface_name,
                                          params_flux);
      else
        mooseError("Boundary condition type not implemented");
    }
  }
}