PressureJumpInterfaceTest.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
 
#include "gtest/gtest.h"
 
#include "PressureJumpInterface.h"
 
#include <cmath>
#include <limits>
 
TEST(PressureJumpInterface, OrthogonalIsotropic)
{
  const RealVectorValue normal(1.0, 0.0, 0.0);
  const Point half_cell(0.5, 0.0, 0.0);
  const RealVectorValue diffusion(2.0, 2.0, 2.0);
  const RealVectorValue gradient(3.0, -4.0, 0.0);
 
  const auto data = NS::FV::pressureJumpInterfaceData(normal,
                                                      half_cell,
                                                      half_cell,
                                                      diffusion,
                                                      diffusion,
                                                      gradient,
                                                      gradient,
                                                      {},
                                                      {},
                                                      3.0,
                                                      true,
                                                      0.0);
 
  ASSERT_TRUE(data.valid);
  EXPECT_TRUE(data.is_baffle);
  EXPECT_DOUBLE_EQ(data.elem_conductance, 12.0);
  EXPECT_DOUBLE_EQ(data.neighbor_conductance, 12.0);
  EXPECT_DOUBLE_EQ(data.transmissibility, 6.0);
  EXPECT_DOUBLE_EQ(data.elem_correction, 0.0);
  EXPECT_DOUBLE_EQ(data.neighbor_correction, 0.0);
  EXPECT_DOUBLE_EQ(data.correction, 0.0);
}
 
TEST(PressureJumpInterface, ConservesFluxAndPreservesJump)
{
  const RealVectorValue normal(0.8, 0.6, 0.0);
  const Point elem_to_face(0.5, 0.1, 0.0);
  const Point face_to_neighbor(0.35, 0.3, 0.0);
  const RealVectorValue elem_diffusion(2.0, 5.0, 1.0);
  const RealVectorValue neighbor_diffusion(7.0, 3.0, 1.0);
  const RealVectorValue elem_gradient(1.2, -0.4, 0.0);
  const RealVectorValue neighbor_gradient(-0.3, 2.1, 0.0);
  const Real elem_pressure = 11.0;
  const Real neighbor_pressure = 4.0;
  const Real jump = 2.5;
 
  const auto data = NS::FV::pressureJumpInterfaceData(normal,
                                                      elem_to_face,
                                                      face_to_neighbor,
                                                      elem_diffusion,
                                                      neighbor_diffusion,
                                                      elem_gradient,
                                                      neighbor_gradient,
                                                      RealVectorValue(),
                                                      RealVectorValue(),
                                                      1.7,
                                                      true,
                                                      jump);
  ASSERT_TRUE(data.valid);
  EXPECT_DOUBLE_EQ(data.elem_jump, jump);
 
  const Real flux = NS::FV::pressureJumpFlux(data, elem_pressure, neighbor_pressure);
  const Real assembled_flux = data.transmissibility * elem_pressure -
                              data.transmissibility * neighbor_pressure -
                              (data.correction + data.transmissibility * data.elem_jump);
  const Real elem_face_pressure =
      elem_pressure - (flux + data.elem_correction) / data.elem_conductance;
  const Real neighbor_face_pressure =
      neighbor_pressure + (flux + data.neighbor_correction) / data.neighbor_conductance;
  const Real elem_half_flux =
      data.elem_conductance * (elem_pressure - elem_face_pressure) - data.elem_correction;
  const Real neighbor_half_flux =
      data.neighbor_conductance * (neighbor_face_pressure - neighbor_pressure) -
      data.neighbor_correction;
  const Real recovered_elem_face_pressure =
      NS::FV::pressureJumpOneSidedFaceValue(data, elem_pressure, flux, /*elem_side=*/true);
  const Real recovered_neighbor_face_pressure =
      NS::FV::pressureJumpOneSidedFaceValue(data, neighbor_pressure, flux, /*elem_side=*/false);
 
  EXPECT_NEAR(elem_face_pressure - neighbor_face_pressure, jump, 1e-13);
  EXPECT_NEAR(assembled_flux, flux, 1e-13);
  EXPECT_NEAR(elem_half_flux, flux, 1e-13);
  EXPECT_NEAR(neighbor_half_flux, flux, 1e-13);
  EXPECT_NEAR(recovered_elem_face_pressure, elem_face_pressure, 1e-13);
  EXPECT_NEAR(recovered_neighbor_face_pressure, neighbor_face_pressure, 1e-13);
}
 
TEST(PressureJumpInterface, ManufacturedPiecewiseLinearPressure)
{
  const RealVectorValue normal(0.8, 0.6, 0.0);
  const Point elem_to_face(0.5, 0.1, 0.0);
  const Point face_to_neighbor(0.35, 0.3, 0.0);
  const RealVectorValue elem_diffusion(2.0, 5.0, 1.0);
  const RealVectorValue neighbor_diffusion(7.0, 3.0, 1.0);
  const RealVectorValue elem_gradient(1.2, -0.4, 0.0);
  const RealVectorValue neighbor_gradient(-0.3, 4.0 / 3.0, 0.0);
  const Real face_area = 1.7;
  const Real elem_face_pressure = 11.0;
  const Real jump = 2.5;
  const Real neighbor_face_pressure = elem_face_pressure - jump;
  const Real elem_pressure = elem_face_pressure - elem_gradient * elem_to_face;
  const Real neighbor_pressure = neighbor_face_pressure + neighbor_gradient * face_to_neighbor;
  const Real expected_flux = -face_area * normal *
                             RealVectorValue(elem_diffusion(0) * elem_gradient(0),
                                             elem_diffusion(1) * elem_gradient(1),
                                             elem_diffusion(2) * elem_gradient(2));
 
  const auto data = NS::FV::pressureJumpInterfaceData(normal,
                                                      elem_to_face,
                                                      face_to_neighbor,
                                                      elem_diffusion,
                                                      neighbor_diffusion,
                                                      elem_gradient,
                                                      neighbor_gradient,
                                                      RealVectorValue(),
                                                      RealVectorValue(),
                                                      face_area,
                                                      true,
                                                      jump);
  ASSERT_TRUE(data.valid);
 
  const Real flux = NS::FV::pressureJumpFlux(data, elem_pressure, neighbor_pressure);
  const Real recovered_elem_face_pressure =
      NS::FV::pressureJumpOneSidedFaceValue(data, elem_pressure, flux, /*elem_side=*/true);
  const Real recovered_neighbor_face_pressure =
      NS::FV::pressureJumpOneSidedFaceValue(data, neighbor_pressure, flux, /*elem_side=*/false);
 
  EXPECT_NEAR(flux, expected_flux, 1e-13);
  EXPECT_NEAR(recovered_elem_face_pressure, elem_face_pressure, 1e-13);
  EXPECT_NEAR(recovered_neighbor_face_pressure, neighbor_face_pressure, 1e-13);
}
 
TEST(PressureJumpInterface, HydrostaticBalanceAcrossJump)
{
  // At rest each side balances its body force with its own pressure gradient,
  // rho_s (H/A)_s = -(rho_s / A_s) grad(p)_s. With jumps in the diffusion coefficients and the
  // gradients, the one-sided H/A fluxes h_s differ. Each side's own flux is still zero, so the
  // face mass flux must vanish and the recovered face pressures must be the one-sided linear
  // extrapolations.
  const RealVectorValue normal(0.8, 0.6, 0.0);
  const Point elem_to_face(0.5, 0.1, 0.0);
  const Point face_to_neighbor(0.35, 0.3, 0.0);
  const RealVectorValue elem_diffusion(2.0, 5.0, 1.0);
  const RealVectorValue neighbor_diffusion(7.0, 3.0, 1.0);
  const RealVectorValue elem_gradient(1.2, -0.4, 0.0);
  const RealVectorValue neighbor_gradient(-0.3, 2.1, 0.0);
  const Real face_area = 1.7;
  const Real elem_face_pressure = 11.0;
  const Real jump = 2.5;
  const Real neighbor_face_pressure = elem_face_pressure - jump;
  const Real elem_pressure = elem_face_pressure - elem_gradient * elem_to_face;
  const Real neighbor_pressure = neighbor_face_pressure + neighbor_gradient * face_to_neighbor;
 
  RealVectorValue elem_rho_hbya;
  RealVectorValue neighbor_rho_hbya;
  for (const auto i : make_range(Moose::dim))
  {
    elem_rho_hbya(i) = -elem_diffusion(i) * elem_gradient(i);
    neighbor_rho_hbya(i) = -neighbor_diffusion(i) * neighbor_gradient(i);
  }
 
  const auto data = NS::FV::pressureJumpInterfaceData(normal,
                                                      elem_to_face,
                                                      face_to_neighbor,
                                                      elem_diffusion,
                                                      neighbor_diffusion,
                                                      elem_gradient,
                                                      neighbor_gradient,
                                                      elem_rho_hbya,
                                                      neighbor_rho_hbya,
                                                      face_area,
                                                      true,
                                                      jump);
  ASSERT_TRUE(data.valid);
  // The test is meaningful only if the side H/A fluxes differ; with equal fluxes any weighting of
  // them gives the same face value.
  ASSERT_GT(std::abs(data.elem_hbya_flux - data.neighbor_hbya_flux), 1.0);
 
  const Real mass_flux = NS::FV::pressureJumpFlux(data, elem_pressure, neighbor_pressure) -
                         NS::FV::pressureJumpHbyAFlux(data);
 
  EXPECT_NEAR(mass_flux, 0.0, 1e-13);
  EXPECT_NEAR(
      NS::FV::pressureJumpOneSidedFaceValue(data, elem_pressure, mass_flux, /*elem_side=*/true),
      elem_face_pressure,
      1e-13);
  EXPECT_NEAR(NS::FV::pressureJumpOneSidedFaceValue(
                  data, neighbor_pressure, mass_flux, /*elem_side=*/false),
              neighbor_face_pressure,
              1e-13);
}
 
TEST(PressureJumpInterface, OrientationReversal)
{
  const RealVectorValue normal(0.8, 0.6, 0.0);
  const Point elem_to_face(0.5, 0.1, 0.0);
  const Point face_to_neighbor(0.35, 0.3, 0.0);
  const RealVectorValue elem_diffusion(2.0, 5.0, 1.0);
  const RealVectorValue neighbor_diffusion(7.0, 3.0, 1.0);
  const RealVectorValue elem_gradient(1.2, -0.4, 0.0);
  const RealVectorValue neighbor_gradient(-0.3, 2.1, 0.0);
 
  const auto forward = NS::FV::pressureJumpInterfaceData(normal,
                                                         elem_to_face,
                                                         face_to_neighbor,
                                                         elem_diffusion,
                                                         neighbor_diffusion,
                                                         elem_gradient,
                                                         neighbor_gradient,
                                                         RealVectorValue(),
                                                         RealVectorValue(),
                                                         1.7,
                                                         true,
                                                         2.5);
  const auto reverse = NS::FV::pressureJumpInterfaceData(-normal,
                                                         -face_to_neighbor,
                                                         -elem_to_face,
                                                         neighbor_diffusion,
                                                         elem_diffusion,
                                                         neighbor_gradient,
                                                         elem_gradient,
                                                         RealVectorValue(),
                                                         RealVectorValue(),
                                                         1.7,
                                                         true,
                                                         -2.5);
  ASSERT_TRUE(forward.valid);
  ASSERT_TRUE(reverse.valid);
 
  const Real forward_flux = NS::FV::pressureJumpFlux(forward, 11.0, 4.0);
  const Real reverse_flux = NS::FV::pressureJumpFlux(reverse, 4.0, 11.0);
 
  EXPECT_NEAR(reverse.elem_conductance, forward.neighbor_conductance, 1e-13);
  EXPECT_NEAR(reverse.neighbor_conductance, forward.elem_conductance, 1e-13);
  EXPECT_NEAR(reverse.correction, -forward.correction, 1e-13);
  EXPECT_NEAR(reverse_flux, -forward_flux, 1e-13);
}
 
TEST(PressureJumpInterface, RejectsInvalidGeometryAndCoefficients)
{
  const auto build_operator =
      [](const Point & elem_to_face, const RealVectorValue & elem_diffusion, const Real face_area)
  {
    return NS::FV::pressureJumpInterfaceData(RealVectorValue(1.0, 0.0, 0.0),
                                             elem_to_face,
                                             Point(0.5, 0.0, 0.0),
                                             elem_diffusion,
                                             RealVectorValue(1.0, 1.0, 1.0),
                                             RealVectorValue(),
                                             RealVectorValue(),
                                             RealVectorValue(),
                                             RealVectorValue(),
                                             face_area,
                                             true,
                                             0.0);
  };
 
  EXPECT_FALSE(build_operator(Point(0.0, 1.0, 0.0), RealVectorValue(1.0, 1.0, 1.0), 1.0).valid);
  EXPECT_FALSE(build_operator(Point(0.5, 0.0, 0.0), RealVectorValue(0.0, 1.0, 1.0), 1.0).valid);
  EXPECT_FALSE(build_operator(Point(0.5, 0.0, 0.0),
                              RealVectorValue(std::numeric_limits<Real>::infinity(), 1.0, 1.0),
                              1.0)
                   .valid);
  EXPECT_FALSE(build_operator(Point(0.5, 0.0, 0.0), RealVectorValue(1.0, 1.0, 1.0), 0.0).valid);
}