FVAdvectedMinmodWeightBased.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 "FVAdvectedMinmodWeightBased.h"

#include "MathFVUtils.h"

#include <algorithm>
#include <limits>

registerMooseObject("MooseApp", FVAdvectedMinmodWeightBased);

InputParameters
FVAdvectedMinmodWeightBased::validParams()
{
  InputParameters params = FVInterpolationMethod::validParams();
  params.addClassDescription(
      "Minmod interpolation for advected quantities implemented as limited blending weights "
      "(no MUSCL reconstruction, no deferred correction).");
  params.addParam<bool>(
      "limit_to_linear",
      true,
      "Whether to limit the scheme to be no more downwind-biased than linear interpolation.");
  params.addRangeCheckedParam<Real>(
      "blending_factor",
      1.0,
      "blending_factor>=0 & blending_factor<=1",
      "Scales the high-order blending strength; 0 gives pure upwind, 1 gives the full limited "
      "blending. Values < 1 can improve linear solver robustness for fully implicit assembly.");
  return params;
}

FVAdvectedMinmodWeightBased::FVAdvectedMinmodWeightBased(const InputParameters & params)
  : FVInterpolationMethod(params),
    _limit_to_linear(getParam<bool>("limit_to_linear")),
    _blending_factor(getParam<Real>("blending_factor"))
{
}

FVAdvectedInterpolationMethod::AdvectedSystemContribution
FVAdvectedMinmodWeightBased::advectedInterpolate(const FaceInfo & face,
                                                 Real elem_value,
                                                 Real neighbor_value,
                                                 const VectorValue<Real> * elem_grad,
                                                 const VectorValue<Real> * neighbor_grad,
                                                 Real mass_flux) const
{
  mooseAssert(elem_grad && neighbor_grad,
              "Minmod advected interpolation requires both element and neighbor gradients.");

  const bool upwind_is_elem = mass_flux >= 0.0;

  const Real phi_upwind = upwind_is_elem ? elem_value : neighbor_value;
  const Real phi_downwind = upwind_is_elem ? neighbor_value : elem_value;

  const VectorValue<Real> grad_upwind = upwind_is_elem ? *elem_grad : *neighbor_grad;
  const Point upwind_to_downwind = upwind_is_elem ? face.dCN() : Point(-face.dCN());

  const auto r_f = Moose::FV::rF(phi_upwind, phi_downwind, grad_upwind, upwind_to_downwind);
  const Real beta = std::max(Real(0.0), std::min(Real(1.0), r_f));

  // Geometric weight associated with the upwind cell for this face.
  const Real w_f = upwind_is_elem ? face.gC() : (1.0 - face.gC());

  // Use a Greenshields-style blending:
  //   phi_f = (1-g)*phi_upwind + g*phi_downwind, where g = beta*(1-w_f)
  const Real g_unclamped = _blending_factor * beta * (1.0 - w_f);
  const Real g_clamped = std::min(std::max(g_unclamped, 0.0), 1.0 - w_f);
  // Clamp to [0, 1-w_f] so the weights do not become more downwind-biased than linear.
  const Real g = _limit_to_linear ? g_clamped : g_unclamped;

  const Real w_upwind = 1.0 - g;
  const Real w_downwind = g;

  // Map (upwind, downwind) weights back to (elem, neighbor) ordering.
  const Real w_elem = upwind_is_elem ? w_upwind : w_downwind;
  const Real w_neighbor = upwind_is_elem ? w_downwind : w_upwind;

  AdvectedSystemContribution result;
  result.weights_matrix = std::make_pair(w_elem, w_neighbor);

  return result;
}

Real
FVAdvectedMinmodWeightBased::advectedInterpolateValue(const FaceInfo & face,
                                                      Real elem_value,
                                                      Real neighbor_value,
                                                      const VectorValue<Real> * elem_grad,
                                                      const VectorValue<Real> * neighbor_grad,
                                                      Real mass_flux) const
{
  const auto result =
      advectedInterpolate(face, elem_value, neighbor_value, elem_grad, neighbor_grad, mass_flux);
  return result.weights_matrix.first * elem_value + result.weights_matrix.second * neighbor_value;
}