ReciprocalMatDiffusion

Calculates the divergence of flux for a variable mobility in reciprocal space.

ReciprocalMatDiffusion computes the divergence of the flux $var element = document.getElementById("moose-equation-6b5766e9-3c80-4be3-8c8a-6fa0f0fc73c1");katex.render("\\mathbf{J} = - M \\, \\nabla \\mu", element, {displayMode:false,throwOnError:false});$ with variable mobility $var element = document.getElementById("moose-equation-fef30af8-d795-4a19-95c5-9965e7fbf915");katex.render("M", element, {displayMode:false,throwOnError:false});$ using FFTs, and adds a no-flux correction on masked boundaries psi>0.

Overview

Inputs: - chemical_potential: field $var element = document.getElementById("moose-equation-ee00f37d-6b09-402d-9a92-0b7d58a807b3");katex.render("\\mu", element, {displayMode:false,throwOnError:false});$ . - mobility: mobility $var element = document.getElementById("moose-equation-87354c5a-7368-4440-b75b-9469bd333172");katex.render("M(\\cdot)", element, {displayMode:false,throwOnError:false});$ . - psi: mask for imposing Neumann no-flux boundaries. Set always_update_psi = true if the mask changes.

The object computes $var element = document.getElementById("moose-equation-16b91d53-07bd-4e49-b970-6835730df828");katex.render("\\nabla \\cdot \\mathbf{J}", element, {displayMode:false,throwOnError:false});$ in reciprocal space and adds $var element = document.getElementById("moose-equation-a44c03ac-92f9-4072-9f3b-94f124a2315d");katex.render("\\nabla psi / psi \\cdot \\mathbf{J}", element, {displayMode:false,throwOnError:false});$ inside the mask to enforce $var element = document.getElementById("moose-equation-d825eefc-1578-4e6e-b1c1-4181dea1e46b");katex.render("\\mathbf{n} \\cdot \\mathbf{J} = 0", element, {displayMode:false,throwOnError:false});$ .

Example Input File Syntax

[TensorComputes<<<{"href": "../../syntax/TensorComputes/index.html"}>>>]
  [Solve<<<{"href": "../../syntax/TensorComputes/Solve/index.html"}>>>]
    [kappa_grad_eta]
      type = ReciprocalMatDiffusion<<<{"description": "Calculates the divergence of flux for a variable mobility in reciprocal space.", "href": "ReciprocalMatDiffusion.html"}>>>
      buffer<<<{"description": "The buffer this compute is writing to"}>>> = 'kappa_grad_eta'
      chemical_potential<<<{"description": "Chemical potential buffer name"}>>> = 'eta'
      mobility<<<{"description": "Mobility buffer name"}>>> = 'L_kappa'
      psi<<<{"description": "Variable to impose Neuamnn BC."}>>> = psi
    []
  []
[]

Input Parameters

bufferThe buffer this compute is writing to
C++ Type:std::string
Controllable:No
Description:The buffer this compute is writing to
chemical_potentialChemical potential buffer name
C++ Type:std::string
Controllable:No
Description:Chemical potential buffer name
mobilityMobility buffer name
C++ Type:std::string
Controllable:No
Description:Mobility buffer name

Required Parameters

always_update_psiFalseSet to true if the BC changes .
Default:False
C++ Type:bool
Controllable:No
Description:Set to true if the BC changes .
psiVariable to impose Neuamnn BC.
C++ Type:std::string
Controllable:No
Description:Variable to impose Neuamnn BC.

Optional Parameters

control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Controllable:No
Description:Set the enabled status of the MooseObject.

Advanced Parameters

(test/tests/kks/KKS_no_flux_bc.i)

#
# Kim-Kim-Suzuki with no-flux BC imposed using the smooth boundary method (SBM), solved on a 2D grid.
# Mask tensor 'psi' supplies the mask for the solve region to the system.
# Note: c is not directly conserved here - the masked value (psi > 0.0)*c will however be conserved.
#


# Constants for Initial Conditions
r = 30
l = 4.2

# Initial condition function for order parameter
eta_IC = '0.5*(1-tanh(2*(sqrt(x^2+y^2)-${r})/${l}))'

# Phase-field model parameters
kappa_eta = 5
rho_sq = 2
w = 1
M = 5
L = 5
c0_a = 0.3
c0_b = 0.7

# Expressions for switching function and bulk Gibbs energy
h_eta = 'eta^3*(6*eta^2-15*eta+10)'
F = '${h_eta}*(${rho_sq}*((c - (1-${h_eta})*(${c0_b} - ${c0_a}))-${c0_a})^2) + (1-${h_eta})*(${rho_sq}*((c + (${h_eta})*(${c0_b} - ${c0_a}))-${c0_b})^2 ) + ${w}*(eta^2)*(1-eta)^2'


[Domain]
  dim = 2
  nx = 20
  ny = 20

  xmin = -50
  xmax = 50
  ymin = -50
  ymax = 50

  # run on a CUDA device (adjust this to `cpu` if not available)
  device_names = 'cpu'

  # automatically create a matching mesh
  mesh_mode = DUMMY
[]

[Functions]
  [psi_func]
    type = ParsedFunction
    expression = 'if(x<x_min-${l},0,if(x>x_min+${l},1,0.5-0.5*cos(pi*(x-(x_min-${l}))/2/${l}) )) * if(x<x_max-${l},1,if(x>x_max+${l},0,0.5+0.5*cos(pi*(x-(x_max-${l}))/2/${l}) ))'
    symbol_names = 'x_min x_max y_min y_max'
    symbol_values = '30 70 0 100'
  []
[]

[TensorComputes]
  [Initialize]
    [c_IC]
      type = ParsedCompute
      buffer = c
      expression = '0.6 + (${c0_a}-0.6)*${eta_IC}'
      extra_symbols = 'true'
    []
    [eta_IC]
      type = ParsedCompute
      buffer = eta
      expression = '${eta_IC}'
      extra_symbols = 'true'
    []
    [psi_init]
       type = MooseFunctionTensor
      function = psi_func
      buffer = psi
    []
    [zero]
      type = ConstantReciprocalTensor
      buffer = zero
    []
    [M]
      type = ConstantTensor
      buffer = M
      real = ${M}
    []
    [L]
      type = ConstantTensor
      buffer = L
      real = ${L}
    []
    [L_kappa]
      type = ConstantTensor
      buffer = L_kappa
      real = ${fparse  ${L} * ${kappa_eta} }
    []
  []
  [Solve]
    [cbar]
      type = ForwardFFT
      buffer = cbar
      input = c
    []
    [etabar]
      type = ForwardFFT
      buffer = etabar
      input = eta
    []
    [mu]
      type = ParsedCompute
      buffer = 'mu'
      expression = '${F}'
      inputs = 'c eta'
      derivatives = 'c'
    []
    [div_J]
      type = ReciprocalMatDiffusion
      buffer = 'div_J'
      chemical_potential = mu
      mobility = M
      psi = psi
    []
    [domega_chem_deta]
      type = ParsedCompute
      buffer = 'domega_chem_deta'
      expression = '${F} - mu*c'
      inputs = 'mu c eta'
      derivatives = 'eta'
    []
    [AC_bulk]
      type = ReciprocalAllenCahn
      buffer = AC_bulk
      dF_chem_deta = domega_chem_deta
      L = L
      psi = psi
    []
    [kappa_grad_eta]
      type = ReciprocalMatDiffusion
      buffer = 'kappa_grad_eta'
      chemical_potential = 'eta'
      mobility = 'L_kappa'
      psi = psi
    []
    [AC_bar]
      type = ParsedCompute
      buffer = AC_bar
      expression = 'kappa_grad_eta + AC_bulk'
      inputs = 'AC_bulk kappa_grad_eta'

    []
  []
[]

[TensorSolver]
  type = AdamsBashforthMoulton
  buffer = 'c eta'
  reciprocal_buffer = 'cbar etabar'
  linear_reciprocal = 'zero zero'
  nonlinear_reciprocal = 'div_J AC_bar'
  substeps = 1e3
  predictor_order = 3
[]

[Postprocessors]
  [total_C]
    type = TensorIntegralPostprocessor
    buffer = c
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [total_eta]
    type = TensorIntegralPostprocessor
    buffer = eta
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

[TensorOutputs]
  [xdmf]
    type = XDMFTensorOutput
    buffer = 'eta c mu psi'
    enable_hdf5 = true
    transpose = false
  []
[]

[Executioner]
  type = Transient
  dt = 0.1
  num_steps = 10
[]

[Outputs]
  csv = true
  perf_graph = true
  execute_on = 'INITIAL TIMESTEP_END'
[]

(examples/libtorch_kks/KKS_libtorch.i)

#
# Kim-Kim-Suzuki with Gibbs energy supplied by a torch model, solved on a 2D grid.
#


# Constants for Initial Conditions
r = 30
l = 4.2

# Initial condition function for order parameter
eta_IC = '0.5*(1-tanh(2*(sqrt(x^2+y^2)-${r})/${l}))'

# Phase-field model parameters
kappa_eta = 5
w = 1
M = 5
L = 5

# Expressions for switching function and bulk Gibbs energy
h_eta = 'eta^3*(6*eta^2-15*eta+10)'


[Domain]
    dim = 2
    nx = 100
    ny = 100

    xmin = -50
    xmax = 50
    ymin = -50
    ymax = 50

    # automatically create a matching mesh
    mesh_mode = DUMMY
[]

[TensorComputes]
    [Initialize]
        [c_IC]
            type = ParsedCompute
            buffer = c
            expression = '0.6 + (0.3-0.6)*${eta_IC}'
            extra_symbols = 'true'
        []
        [eta_IC]
            type = ParsedCompute
            buffer = eta
            expression = '${eta_IC}'
            extra_symbols = 'true'
        []
        [psi_init]
            type = ConstantTensor
            buffer = psi
            real = 1
        []

        [M]
            type = ConstantTensor
            buffer = M
            real = ${M}
        []
        [L]
            type = ConstantTensor
            buffer = L
            real = ${L}
        []
        [L_kappa]
            type = ReciprocalLaplacianFactor
            buffer = L_kappa
            factor = ${fparse  ${L} * ${kappa_eta} }
        []
        [h_eta_IC]
            type = ParsedCompute
            buffer = h_eta
            expression = '${h_eta}'
            inputs = eta
        []
        [G_func_IC]
            type = LibtorchGibbsEnergy
            buffer = 'G'
            phase_fractions = 'h_eta'
            concentrations = 'c'
            domega_detas = 'dG_dh'
            chem_pots = 'mu'
            libtorch_model_file = 'torch_NN_gibbs_model.pt'
        []
        [smooth]
            type = DeAliasingTensor
            method = HOULI
            buffer = smooth
        []
    []
    [Solve]
        [h_eta]
            type = ParsedCompute
            buffer = h_eta
            expression = '${h_eta}'
            inputs = eta
        []
        [G_func]
            type = LibtorchGibbsEnergy
            buffer = 'G'
            phase_fractions = 'h_eta'
            concentrations = 'c'
            domega_detas = 'dG_dh'
            chem_pots = 'mu'
            libtorch_model_file = 'torch_NN_gibbs_model.pt'
        []
        [dG_deta]
            type = ParsedCompute
            buffer = 'dG_deta'
            inputs = 'eta dG_dh'
            expression = 'dG_dh * ${h_eta} + ${w} * eta^2 * (1-eta^2)^2'
            derivatives = 'eta'
        []

        [etabar]
            type = ForwardFFT
            buffer = etabar
            input = eta
        []
        [AC_bulk]
            type = ReciprocalAllenCahn
            L = L
            buffer = AC_bulk
            dF_chem_deta = dG_deta
            psi = psi
        []
        [NL_eta]
            type = ParsedCompute
            buffer = NL_eta
            expression = 'AC_bulk '
            inputs = 'AC_bulk'
        []
        [cbar]
            type = ForwardFFT
            buffer = cbar
            input = c
        []
        [div_J]
            type = ReciprocalMatDiffusion
            buffer = 'div_J'
            chemical_potential = mu
            mobility = M
            psi = psi
        []
        [NL_c]
            type = ParsedCompute
            buffer = 'NL_c'
            inputs = 'div_J smooth'
            expression = 'smooth * div_J'
        []
    []
[]

[TensorSolver]
    type = AdamsBashforthMoulton
    buffer = 'c eta'
    reciprocal_buffer = 'cbar etabar'
    linear_reciprocal = '0 L_kappa'
    nonlinear_reciprocal = 'NL_c NL_eta'
    substeps = 1e3
    predictor_order = 3
    corrector_order = 1
    corrector_steps = 1
[]

[Postprocessors]
    [total_c]
        type = TensorIntegralPostprocessor
        buffer = c
    []
[]

[TensorOutputs]
    [xdmf]
        type = XDMFTensorOutput
        buffer = 'eta c mu psi dG_deta dG_dh G'
        enable_hdf5 = true
        transpose = false
    []
[]

[Executioner]
    type = Transient
    num_steps = 100
  [TimeStepper]
    type = IterationAdaptiveDT
    growth_factor = 1.25
    dt = 0.1
  []
  dtmax = 10
[]

[Outputs]
    csv = true
    perf_graph = true
    execute_on = 'INITIAL TIMESTEP_END'
[]

(test/tests/kks/KKS_libtorch.i)

#
# Kim-Kim-Suzuki with Gibbs energy supplied by a torch model, solved on a 2D grid.
#


# Constants for Initial Conditions
r = 30
l = 4.2

# Initial condition function for order parameter
eta_IC = '0.5*(1-tanh(2*(sqrt(x^2+y^2)-${r})/${l}))'

# Phase-field model parameters
kappa_eta = 5
w = 1
M = 5
L = 5

# Expressions for switching function and bulk Gibbs energy
h_eta = 'eta^3*(6*eta^2-15*eta+10)'


[Domain]
    dim = 2
    nx = 50
    ny = 50

    xmin = -50
    xmax = 50
    ymin = -50
    ymax = 50

    # automatically create a matching mesh
    mesh_mode = DUMMY
[]

[TensorComputes]
    [Initialize]
        [c_IC]
            type = ParsedCompute
            buffer = c
            expression = '0.7 + (0.3-0.6)*${eta_IC}'
            extra_symbols = 'true'
        []
        [eta_IC]
            type = ParsedCompute
            buffer = eta
            expression = '${eta_IC}'
            extra_symbols = 'true'
        []
        [psi_init]
            type = ConstantTensor
            buffer = psi
            real = 1
        []

        [M]
            type = ConstantTensor
            buffer = M
            real = ${M}
        []
        [L]
            type = ConstantTensor
            buffer = L
            real = ${L}
        []
        [L_kappa]
            type = ReciprocalLaplacianFactor
            buffer = L_kappa
            factor = ${fparse  ${L} * ${kappa_eta} }
        []
        [h_eta_IC]
            type = ParsedCompute
            buffer = h_eta
            expression = '${h_eta}'
            inputs = eta
        []
        [G_func_IC]
            type = LibtorchGibbsEnergy
            buffer = 'G'
            phase_fractions = 'h_eta'
            concentrations = 'c'
            domega_detas = 'dG_dh'
            chem_pots = 'mu'
            libtorch_model_file = 'marlin:libtorch_gibbs_energy/torch_NN_gibbs_model.pt'
        []
    []
    [Solve]
        [h_eta]
            type = ParsedCompute
            buffer = h_eta
            expression = '${h_eta}'
            inputs = eta
        []
        [G_func]
            type = LibtorchGibbsEnergy
            buffer = 'G'
            phase_fractions = 'h_eta'
            concentrations = 'c'
            domega_detas = 'dG_dh'
            chem_pots = 'mu'
            libtorch_model_file = 'marlin:libtorch_gibbs_energy/torch_NN_gibbs_model.pt'
        []
        [dG_deta]
            type = ParsedCompute
            buffer = 'dG_deta'
            inputs = 'eta dG_dh'
            expression = 'dG_dh * ${h_eta} + ${w} * eta^2 * (1-eta^2)^2'
            derivatives = 'eta'
        []

        [etabar]
            type = ForwardFFT
            buffer = etabar
            input = eta
        []
        [AC_bulk]
            type = ReciprocalAllenCahn
            L = L
            buffer = AC_bulk
            dF_chem_deta = dG_deta
            psi = psi
        []
        [NL_eta]
            type = ParsedCompute
            buffer = NL_eta
            expression = 'AC_bulk '
            inputs = 'AC_bulk'
        []
        [cbar]
            type = ForwardFFT
            buffer = cbar
            input = c
        []
        [div_J]
            type = ReciprocalMatDiffusion
            buffer = 'div_J'
            chemical_potential = mu
            mobility = M
            psi = psi
        []
    []
[]

[TensorSolver]
    type = AdamsBashforthMoulton
    buffer = 'c eta'
    reciprocal_buffer = 'cbar etabar'
    linear_reciprocal = '0 L_kappa'
    nonlinear_reciprocal = 'div_J NL_eta'
    substeps = 1e3
    predictor_order = 3
[]

[TensorOutputs]
    [xdmf]
        type = XDMFTensorOutput
        buffer = 'eta c mu psi dG_deta dG_dh G'
        enable_hdf5 = true
        transpose = false
    []
[]

[Executioner]
    type = Transient
    dt = 0.1
    num_steps = 10
[]

[Outputs]
    csv = true
    perf_graph = true
    execute_on = 'INITIAL TIMESTEP_END'
[]

(test/tests/kks/KKS_no_flux_bc.i)

#
# Kim-Kim-Suzuki with no-flux BC imposed using the smooth boundary method (SBM), solved on a 2D grid.
# Mask tensor 'psi' supplies the mask for the solve region to the system.
# Note: c is not directly conserved here - the masked value (psi > 0.0)*c will however be conserved.
#


# Constants for Initial Conditions
r = 30
l = 4.2

# Initial condition function for order parameter
eta_IC = '0.5*(1-tanh(2*(sqrt(x^2+y^2)-${r})/${l}))'

# Phase-field model parameters
kappa_eta = 5
rho_sq = 2
w = 1
M = 5
L = 5
c0_a = 0.3
c0_b = 0.7

# Expressions for switching function and bulk Gibbs energy
h_eta = 'eta^3*(6*eta^2-15*eta+10)'
F = '${h_eta}*(${rho_sq}*((c - (1-${h_eta})*(${c0_b} - ${c0_a}))-${c0_a})^2) + (1-${h_eta})*(${rho_sq}*((c + (${h_eta})*(${c0_b} - ${c0_a}))-${c0_b})^2 ) + ${w}*(eta^2)*(1-eta)^2'


[Domain]
  dim = 2
  nx = 20
  ny = 20

  xmin = -50
  xmax = 50
  ymin = -50
  ymax = 50

  # run on a CUDA device (adjust this to `cpu` if not available)
  device_names = 'cpu'

  # automatically create a matching mesh
  mesh_mode = DUMMY
[]

[Functions]
  [psi_func]
    type = ParsedFunction
    expression = 'if(x<x_min-${l},0,if(x>x_min+${l},1,0.5-0.5*cos(pi*(x-(x_min-${l}))/2/${l}) )) * if(x<x_max-${l},1,if(x>x_max+${l},0,0.5+0.5*cos(pi*(x-(x_max-${l}))/2/${l}) ))'
    symbol_names = 'x_min x_max y_min y_max'
    symbol_values = '30 70 0 100'
  []
[]

[TensorComputes]
  [Initialize]
    [c_IC]
      type = ParsedCompute
      buffer = c
      expression = '0.6 + (${c0_a}-0.6)*${eta_IC}'
      extra_symbols = 'true'
    []
    [eta_IC]
      type = ParsedCompute
      buffer = eta
      expression = '${eta_IC}'
      extra_symbols = 'true'
    []
    [psi_init]
       type = MooseFunctionTensor
      function = psi_func
      buffer = psi
    []
    [zero]
      type = ConstantReciprocalTensor
      buffer = zero
    []
    [M]
      type = ConstantTensor
      buffer = M
      real = ${M}
    []
    [L]
      type = ConstantTensor
      buffer = L
      real = ${L}
    []
    [L_kappa]
      type = ConstantTensor
      buffer = L_kappa
      real = ${fparse  ${L} * ${kappa_eta} }
    []
  []
  [Solve]
    [cbar]
      type = ForwardFFT
      buffer = cbar
      input = c
    []
    [etabar]
      type = ForwardFFT
      buffer = etabar
      input = eta
    []
    [mu]
      type = ParsedCompute
      buffer = 'mu'
      expression = '${F}'
      inputs = 'c eta'
      derivatives = 'c'
    []
    [div_J]
      type = ReciprocalMatDiffusion
      buffer = 'div_J'
      chemical_potential = mu
      mobility = M
      psi = psi
    []
    [domega_chem_deta]
      type = ParsedCompute
      buffer = 'domega_chem_deta'
      expression = '${F} - mu*c'
      inputs = 'mu c eta'
      derivatives = 'eta'
    []
    [AC_bulk]
      type = ReciprocalAllenCahn
      buffer = AC_bulk
      dF_chem_deta = domega_chem_deta
      L = L
      psi = psi
    []
    [kappa_grad_eta]
      type = ReciprocalMatDiffusion
      buffer = 'kappa_grad_eta'
      chemical_potential = 'eta'
      mobility = 'L_kappa'
      psi = psi
    []
    [AC_bar]
      type = ParsedCompute
      buffer = AC_bar
      expression = 'kappa_grad_eta + AC_bulk'
      inputs = 'AC_bulk kappa_grad_eta'

    []
  []
[]

[TensorSolver]
  type = AdamsBashforthMoulton
  buffer = 'c eta'
  reciprocal_buffer = 'cbar etabar'
  linear_reciprocal = 'zero zero'
  nonlinear_reciprocal = 'div_J AC_bar'
  substeps = 1e3
  predictor_order = 3
[]

[Postprocessors]
  [total_C]
    type = TensorIntegralPostprocessor
    buffer = c
    execute_on = 'INITIAL TIMESTEP_END'
  []
  [total_eta]
    type = TensorIntegralPostprocessor
    buffer = eta
    execute_on = 'INITIAL TIMESTEP_END'
  []
[]

[TensorOutputs]
  [xdmf]
    type = XDMFTensorOutput
    buffer = 'eta c mu psi'
    enable_hdf5 = true
    transpose = false
  []
[]

[Executioner]
  type = Transient
  dt = 0.1
  num_steps = 10
[]

[Outputs]
  csv = true
  perf_graph = true
  execute_on = 'INITIAL TIMESTEP_END'
[]

Overview
Example Input File Syntax
Input Parameters
Input Files