TensorRobinBC

c_function

C++ Type:FunctionName

Unit:(no unit assumed)

Controllable:No

Description:Function giving c along the face, evaluated at the face position at the solve time.

c_buffer

C++ Type:std::string

Controllable:No

Description:Buffer whose first owned cell layer next to the face holds c along the face. The compute that writes it must evaluate c at the face position, not at the cell centres.

(test/tests/real_space/robin_steady_state.i)

# Verification against an independent analytic solution: the steady 1D elliptic
# Robin-BC benchmark from OpenGeoSys's documentation,
# https://www.opengeosys.org/6.5.7/docs/benchmarks/elliptic/elliptic-robin/
# ("second benchmark": Dirichlet at x=0, Robin at x=1). The governing equation is
# k*Laplacian(h) = 0 on [0,1] with k=1, boundary conditions
#   h(0) = g_D,                 dh/dn|_{x=1} = alpha*(h0 - h(1)),
# and exact linear solution
#   h(x) = [alpha*(h0 - g_D)/(1 + alpha)] * x + g_D.
# TensorRobinBC's convention is a*u + b*du/dn = c; the OGS Robin condition
# du/dn = alpha*(h0 - u) rearranges to alpha*u + du/dn = alpha*h0, i.e.
# a = alpha, b = 1, c = alpha*h0.
#
# Uses alpha = +2 rather than OGS's alpha = -2: the ghost coefficient on u0,
# -(a/2 - b/dx)/(a/2 + b/dx), has magnitude > 1 for a < 0, so explicit relaxation to
# the steady state diverges (OGS solves the steady system directly). With alpha = 2,
# h0 = 1.5 and g_D = 1 the exact solution is h(x) = x/3 + 1.
#
# Y is periodic, and the IC matches neither the steady profile nor the BCs.
[Domain]
  dim = 2
  nx = 20
  ny = 8
  parallel_mode = REAL_SPACE
  periodic_directions = 'Y'
  xmin = 0
  xmax = 1
  ymin = -1
  ymax = 1
[]

[Functions]
  [u_ic_func]
    type = ParsedFunction
    expression = '0.5'
  []
  [u_exact_func]
    type = ParsedFunction
    expression = '1 + x / 3'
  []
[]

[TensorComputes]
  [Initialize]
    [u_ic]
      type = MooseFunctionTensor
      buffer = u
      function = u_ic_func
    []
  []
  [BCs]
    [left]
      type = TensorDirichletBC
      variable = u
      boundary = left
      value = 1
    []
    [right]
      # alpha*u + du/dn = alpha*h0, with alpha = 2, h0 = 1.5.
      type = TensorRobinBC
      variable = u
      boundary = right
      a = 2
      b = 1
      c = 3
    []
  []
  [Solve]
    [dt]
      type = FiniteDifferenceLaplacian
      input = u
      buffer = dudt
    []
    [u_exact]
      type = MooseFunctionTensor
      buffer = u_exact
      function = u_exact_func
    []
    [diff]
      type = ParsedCompute
      buffer = diff
      expression = 'abs(u - u_exact)'
      inputs = 'u u_exact'
    []
  []
[]

[TensorSolver]
  type = RealSpaceForwardEuler
  buffer = u
  time_derivative = dudt
  substeps = 1000
[]

[Postprocessors]
  [max_error]
    type = TensorExtremeValuePostprocessor
    buffer = diff
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
[]

[Executioner]
  type = Transient
  # same grid and CFL margin as dirichlet_linear.i
  dt = 1
  num_steps = 4
[]

[Outputs]
  [out]
    type = CSV
  []
[]

(test/tests/real_space/robin_steady_state.i)

# Verification against an independent analytic solution: the steady 1D elliptic
# Robin-BC benchmark from OpenGeoSys's documentation,
# https://www.opengeosys.org/6.5.7/docs/benchmarks/elliptic/elliptic-robin/
# ("second benchmark": Dirichlet at x=0, Robin at x=1). The governing equation is
# k*Laplacian(h) = 0 on [0,1] with k=1, boundary conditions
#   h(0) = g_D,                 dh/dn|_{x=1} = alpha*(h0 - h(1)),
# and exact linear solution
#   h(x) = [alpha*(h0 - g_D)/(1 + alpha)] * x + g_D.
# TensorRobinBC's convention is a*u + b*du/dn = c; the OGS Robin condition
# du/dn = alpha*(h0 - u) rearranges to alpha*u + du/dn = alpha*h0, i.e.
# a = alpha, b = 1, c = alpha*h0.
#
# Uses alpha = +2 rather than OGS's alpha = -2: the ghost coefficient on u0,
# -(a/2 - b/dx)/(a/2 + b/dx), has magnitude > 1 for a < 0, so explicit relaxation to
# the steady state diverges (OGS solves the steady system directly). With alpha = 2,
# h0 = 1.5 and g_D = 1 the exact solution is h(x) = x/3 + 1.
#
# Y is periodic, and the IC matches neither the steady profile nor the BCs.
[Domain]
  dim = 2
  nx = 20
  ny = 8
  parallel_mode = REAL_SPACE
  periodic_directions = 'Y'
  xmin = 0
  xmax = 1
  ymin = -1
  ymax = 1
[]

[Functions]
  [u_ic_func]
    type = ParsedFunction
    expression = '0.5'
  []
  [u_exact_func]
    type = ParsedFunction
    expression = '1 + x / 3'
  []
[]

[TensorComputes]
  [Initialize]
    [u_ic]
      type = MooseFunctionTensor
      buffer = u
      function = u_ic_func
    []
  []
  [BCs]
    [left]
      type = TensorDirichletBC
      variable = u
      boundary = left
      value = 1
    []
    [right]
      # alpha*u + du/dn = alpha*h0, with alpha = 2, h0 = 1.5.
      type = TensorRobinBC
      variable = u
      boundary = right
      a = 2
      b = 1
      c = 3
    []
  []
  [Solve]
    [dt]
      type = FiniteDifferenceLaplacian
      input = u
      buffer = dudt
    []
    [u_exact]
      type = MooseFunctionTensor
      buffer = u_exact
      function = u_exact_func
    []
    [diff]
      type = ParsedCompute
      buffer = diff
      expression = 'abs(u - u_exact)'
      inputs = 'u u_exact'
    []
  []
[]

[TensorSolver]
  type = RealSpaceForwardEuler
  buffer = u
  time_derivative = dudt
  substeps = 1000
[]

[Postprocessors]
  [max_error]
    type = TensorExtremeValuePostprocessor
    buffer = diff
    value_type = MAX
    execute_on = 'TIMESTEP_END'
  []
[]

[Executioner]
  type = Transient
  # same grid and CFL margin as dirichlet_linear.i
  dt = 1
  num_steps = 4
[]

[Outputs]
  [out]
    type = CSV
  []
[]

(test/tests/implicit_fd/bc_mms.i)

# Implicit solve of u_t - D Laplacian(u) = f on the unit square with boundary faces, against the
# manufactured solution
#
#   u = g + tf t q,  g = cos(3 pi x / 2) cos(pi y) + x y,  q = x^2 - y^2,
#
# with tf = 0 for the steady problem -D Laplacian(u) = f and tf = 1 for backward Euler.
# q is harmonic, so f = tf q - D Laplacian(g) is independent of time and the backward Euler
# step makes no time error: the error is the spatial one, provided the boundary values are
# taken at the time of the unknown, t + dt (at t it would be of order dt).
#
# Each face has a different condition and value source:
#   left   (x = 0)  Dirichlet   u = cos(pi y) - tf t y^2                     function
#   right  (x = 1)  Robin       ra u + rb du/dx = c                         function
#                                                      or buffer (right_buffer, steady)
#   bottom (y = 0)  Neumann     du/dn = -du/dy = -x                          buffer
#   top    (y = 1)  Neumann     du/dn = du/dy = x - 2 tf t                   function
#                                                      or buffer (top_buffer, values at t)

D = 0.5
n = 32
stencil = 3
tf = 0
ops = 'lap'
num_steps = 1
dt = 0.05
# Robin coefficients on the right face (the 5-point stencil supports ra = 1, rb = 0 only)
ra = 1
rb = 1

[Domain]
  dim = 2
  nx = ${n}
  ny = ${n}
  xmax = 1
  ymax = 1
  parallel_mode = REAL_SPACE
[]

[Functions]
  [u_exact]
    type = ParsedFunction
    expression = 'cos(1.5*pi*x) * cos(pi*y) + x*y + ${tf} * t * (x^2 - y^2)'
  []
  [f_func]
    type = ParsedFunction
    # Laplacian(g) = -13/4 pi^2 cos(3 pi x / 2) cos(pi y)
    expression = '${tf} * (x^2 - y^2) + ${D} * 3.25 * pi^2 * cos(1.5*pi*x) * cos(pi*y)'
  []
  [right_c]
    type = ParsedFunction
    # on x = 1: u = y + tf t (1 - y^2), du/dx = 3/2 pi cos(pi y) + y + 2 tf t
    expression = '${ra} * (y + ${tf} * t * (1 - y^2)) + ${rb} * (1.5*pi*cos(pi*y) + y + 2 * ${tf} * t)'
  []
  [bottom_flux]
    type = ParsedFunction
    expression = '-x'
  []
  [top_flux]
    type = ParsedFunction
    expression = 'x - 2 * ${tf} * t'
  []
[]

[TensorBuffers]
  [u]
  []
  [dudt_u]
  []
  [lap_u]
  []
  [f]
  []
  [bottom_flux]
  []
  [right_c]
  []
  [top_flux]
  []
  [exact]
  []
  [diff]
  []
[]

[TensorComputes]
  [Initialize]
    [u_ic]
      type = MooseFunctionTensor
      buffer = u
      function = u_exact
    []
    [f_ic]
      type = MooseFunctionTensor
      buffer = f
      function = f_func
    []
    # a face value: depends on the tangential coordinate x only
    [bottom_flux]
      type = MooseFunctionTensor
      buffer = bottom_flux
      function = bottom_flux
    []
    [right_c]
      type = MooseFunctionTensor
      buffer = right_c
      function = right_c
    []
  []
  [BCs]
    active = 'left right bottom top'
    [left]
      type = TensorDirichletBC
      variable = u
      boundary = left
      value_function = u_exact
    []
    [right]
      type = TensorRobinBC
      variable = u
      boundary = right
      a = ${ra}
      b = ${rb}
      c_function = right_c
    []
    # the Robin c from a buffer instead (active = 'left right_buffer bottom top'), steady only
    [right_buffer]
      type = TensorRobinBC
      variable = u
      boundary = right
      a = ${ra}
      b = ${rb}
      c_buffer = right_c
    []
    [bottom]
      type = TensorNeumannBC
      variable = u
      boundary = bottom
      flux_buffer = bottom_flux
    []
    [top]
      type = TensorNeumannBC
      variable = u
      boundary = top
      flux_function = top_flux
    []
    # the top flux from a buffer instead (active = 'left right bottom top_buffer'); a buffer
    # value is the one its compute wrote at the start of the substep, at t
    [top_buffer]
      type = TensorNeumannBC
      variable = u
      boundary = top
      flux_buffer = top_flux
    []
  []
  [Solve]
    [top_flux]
      type = MooseFunctionTensor
      buffer = top_flux
      function = top_flux
    []
    [dudt]
      type = TimeDerivativeOperator
      buffer = dudt_u
      input = u
    []
    [lap]
      type = FiniteDifferenceLaplacian
      buffer = lap_u
      input = u
      factor = '${fparse -D}'
      stencil_width = ${stencil}
    []
  []
  [Postprocess]
    [exact_compute]
      type = MooseFunctionTensor
      buffer = exact
      function = u_exact
    []
    [diff_compute]
      type = ParsedCompute
      buffer = diff
      expression = 'abs(u - exact)'
      inputs = 'u exact'
    []
  []
[]

[TensorSolver]
  type = ImplicitFDSolver
  buffer = u
  operators = '${ops}'
  source = f
  linear_solver = TORCH_PCG
  l_tol = 1e-10
[]

[Postprocessors]
  [max_error]
    type = TensorExtremeValuePostprocessor
    buffer = diff
    value_type = MAX
  []
[]

[Problem]
  type = TensorProblem
[]

[Executioner]
  type = Transient
  num_steps = ${num_steps}
  dt = ${dt}
[]

[Outputs]
  csv = true
[]

(include/tensor_computes/TensorNeumannBC.h)

/**********************************************************************/
/*                     DO NOT MODIFY THIS HEADER                      */
/*            Marlin, a Fourier spectral solver for MOOSE             */
/*                                                                    */
/*            Copyright 2024 Battelle Energy Alliance, LLC            */
/*                        ALL RIGHTS RESERVED                         */
/**********************************************************************/

#pragma once

#include "TensorRobinBC.h"

/// Neumann preset of TensorRobinBC: a=0, b=1, c=flux (du/dn = flux on the face).
class TensorNeumannBC : public TensorRobinBC
{
public:
  static InputParameters validParams();

  TensorNeumannBC(const InputParameters & parameters);
};

(include/tensor_computes/TensorDirichletBC.h)

/**********************************************************************/
/*                     DO NOT MODIFY THIS HEADER                      */
/*            Marlin, a Fourier spectral solver for MOOSE             */
/*                                                                    */
/*            Copyright 2024 Battelle Energy Alliance, LLC            */
/*                        ALL RIGHTS RESERVED                         */
/**********************************************************************/

#pragma once

#include "TensorRobinBC.h"

/// Dirichlet preset of TensorRobinBC: a=1, b=0, c=value (u = value on the face).
class TensorDirichletBC : public TensorRobinBC
{
public:
  static InputParameters validParams();

  TensorDirichletBC(const InputParameters & parameters);
};