SwitchingFunctionPenalty

Penalty kernel to constrain the sum of all switching functions in a multiphase system.

The penalty is added in the form of a free energy term

var element = document.getElementById("moose-equation-fee15889-f3be-4c7f-b2bf-5c64582d439f");katex.render("F_p=\\gamma k(\\vec \\eta)^2,", element, {displayMode:true,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});

where $var element = document.getElementById("moose-equation-404c849c-4e3c-48b0-b584-cdef08087af7");katex.render("\\gamma", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is the penalty prefactor (penalty) and $var element = document.getElementById("moose-equation-038be530-d438-44b2-a773-cd7ff45fc189");katex.render("k = 1-\\sum_i\\eta_i", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ is the constraint. The constraint is enforced approximately to a tolerance of $var element = document.getElementById("moose-equation-f29ca895-d9eb-4cdf-8cdb-a1466c30e167");katex.render("\\frac1\\gamma", element, {displayMode:false,throwOnError:false,macros:{"\\pf":"\\frac{\\partial #1}{\\partial #2}"}});$ (depending on the shape and units of the free energy).

Also see: Multiphase models

Input Parameters

etaseta_i order parameters, one for each h
C++ Type:std::vector<VariableName>
Options:
Description:eta_i order parameters, one for each h
variableThe name of the variable that this residual object operates on
C++ Type:NonlinearVariableName
Options:
Description:The name of the variable that this residual object operates on

Required Parameters

blockThe list of block ids (SubdomainID) that this object will be applied
C++ Type:std::vector<SubdomainName>
Options:
Description:The list of block ids (SubdomainID) that this object will be applied
displacementsThe displacements
C++ Type:std::vector<VariableName>
Options:
Description:The displacements
h_namesSwitching Function Materials that provide h(eta_i)
C++ Type:std::vector<MaterialPropertyName>
Options:
Description:Switching Function Materials that provide h(eta_i)
penalty1Penalty scaling factor
Default:1
C++ Type:double
Options:
Description:Penalty scaling factor

Optional Parameters

control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Options:
Description:Adds user-defined labels for accessing object parameters via control logic.
diag_save_inThe name of auxiliary variables to save this Kernel's diagonal Jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Options:
Description:The name of auxiliary variables to save this Kernel's diagonal Jacobian contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Options:
Description:Set the enabled status of the MooseObject.
implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Options:
Description:Determines whether this object is calculated using an implicit or explicit form
save_inThe name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
C++ Type:std::vector<AuxVariableName>
Options:
Description:The name of auxiliary variables to save this Kernel's residual contributions to. Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)
seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Options:
Description:The seed for the master random number generator
use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Options:
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.

Advanced Parameters

extra_matrix_tagsThe extra tags for the matrices this Kernel should fill
C++ Type:std::vector<TagName>
Options:
Description:The extra tags for the matrices this Kernel should fill
extra_vector_tagsThe extra tags for the vectors this Kernel should fill
C++ Type:std::vector<TagName>
Options:
Description:The extra tags for the vectors this Kernel should fill
matrix_tagssystemThe tag for the matrices this Kernel should fill
Default:system
C++ Type:MultiMooseEnum
Options:nontime, system
Description:The tag for the matrices this Kernel should fill
vector_tagsnontimeThe tag for the vectors this Kernel should fill
Default:nontime
C++ Type:MultiMooseEnum
Options:nontime, time
Description:The tag for the vectors this Kernel should fill

Tagging Parameters

Input Files

(modules/phase_field/examples/multiphase/DerivativeMultiPhaseMaterial.i)
(modules/phase_field/test/tests/MultiPhase/penalty.i)

(modules/phase_field/examples/multiphase/DerivativeMultiPhaseMaterial.i)

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 40
  ny = 40
  nz = 0
  xmin = -12
  xmax = 12
  ymin = -12
  ymax = 12
  elem_type = QUAD4
[]

[GlobalParams]
  # let's output all material properties for demonstration purposes
  outputs = exodus

  # prefactor on the penalty function kernels. The higher this value is, the
  # more rigorously the constraint is enforced
  penalty = 1e3
[]

#
# These AuxVariables hold the directly calculated free energy density in the
# simulation cell. They are provided for visualization purposes.
#
[AuxVariables]
  [./local_energy]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./cross_energy]
    order = CONSTANT
    family = MONOMIAL
  [../]
[]

[AuxKernels]
  [./local_free_energy]
    type = TotalFreeEnergy
    variable = local_energy
    interfacial_vars = 'c'
    kappa_names = 'kappa_c'
    additional_free_energy = cross_energy
  [../]

  #
  # Helper kernel to cpompute the gradient contribution from interfaces of order
  # parameters evolved using the ACMultiInterface kernel
  #
  [./cross_terms]
    type = CrossTermGradientFreeEnergy
    variable = cross_energy
    interfacial_vars = 'eta1 eta2 eta3'
    #
    # The interface coefficient matrix. This should be symmetrical!
    #
    kappa_names = 'kappa11 kappa12 kappa13
                   kappa21 kappa22 kappa23
                   kappa31 kappa32 kappa33'
  [../]
[]

[Variables]
  [./c]
    order = FIRST
    family = LAGRANGE
    #
    # We set up a smooth cradial concentrtaion gradient
    # The concentration will quickly change to adapt to the preset order
    # parameters eta1, eta2, and eta3
    #
    [./InitialCondition]
      type = SmoothCircleIC
      x1 = 0.0
      y1 = 0.0
      radius = 5.0
      invalue = 1.0
      outvalue = 0.01
      int_width = 10.0
    [../]
  [../]

  [./eta1]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = FunctionIC
      #
      # Note: this initial conditions sets up a _sharp_ interface. Ideally
      # we should start with a smooth interface with a width consistent
      # with the kappa parameter supplied for the given interface.
      #
      function = 'r:=sqrt(x^2+y^2);if(r<=4,1,0)'
    [../]
  [../]
  [./eta2]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = FunctionIC
      function = 'r:=sqrt(x^2+y^2);if(r>4&r<=7,1,0)'
    [../]
  [../]
  [./eta3]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = FunctionIC
      function = 'r:=sqrt(x^2+y^2);if(r>7,1,0)'
    [../]
  [../]
[]

[Kernels]
  #
  # Cahn-Hilliard kernel for the concentration variable.
  # Note that we are not using an interfcae kernel on this variable, but rather
  # rely on the interface width enforced on the order parameters. This allows us
  # to use a direct solve using the CahnHilliard kernel _despite_ only using first
  # order elements.
  #
  [./c_res]
    type = CahnHilliard
    variable = c
    f_name = F
    args = 'eta1 eta2 eta3'
  [../]
  [./time]
    type = TimeDerivative
    variable = c
  [../]

  #
  # Order parameter eta1
  # Each order parameter is acted on by 4 kernels:
  #  1. The stock time derivative deta_i/dt kernel
  #  2. The Allen-Cahn kernel that takes a Dervative Material for the free energy
  #  3. A gradient interface kernel that includes cross terms
  #     see http://mooseframework.org/wiki/PhysicsModules/PhaseField/DevelopingModels/MultiPhaseModels/ACMultiInterface/
  #  4. A penalty contribution that forces the interface contributions h(eta)
  #     to sum up to unity
  #
  [./deta1dt]
    type = TimeDerivative
    variable = eta1
  [../]
  [./ACBulk1]
    type = AllenCahn
    variable = eta1
    args = 'eta2 eta3 c'
    mob_name = L1
    f_name = F
  [../]
  [./ACInterface1]
    type = ACMultiInterface
    variable = eta1
    etas = 'eta1 eta2 eta3'
    mob_name = L1
    kappa_names = 'kappa11 kappa12 kappa13'
  [../]
  [./penalty1]
    type = SwitchingFunctionPenalty
    variable = eta1
    etas    = 'eta1 eta2 eta3'
    h_names = 'h1   h2   h3'
  [../]

  #
  # Order parameter eta2
  #
  [./deta2dt]
    type = TimeDerivative
    variable = eta2
  [../]
  [./ACBulk2]
    type = AllenCahn
    variable = eta2
    args = 'eta1 eta3 c'
    mob_name = L2
    f_name = F
  [../]
  [./ACInterface2]
    type = ACMultiInterface
    variable = eta2
    etas = 'eta1 eta2 eta3'
    mob_name = L2
    kappa_names = 'kappa21 kappa22 kappa23'
  [../]
  [./penalty2]
    type = SwitchingFunctionPenalty
    variable = eta2
    etas    = 'eta1 eta2 eta3'
    h_names = 'h1   h2   h3'
  [../]

  #
  # Order parameter eta3
  #
  [./deta3dt]
    type = TimeDerivative
    variable = eta3
  [../]
  [./ACBulk3]
    type = AllenCahn
    variable = eta3
    args = 'eta1 eta2 c'
    mob_name = L3
    f_name = F
  [../]
  [./ACInterface3]
    type = ACMultiInterface
    variable = eta3
    etas = 'eta1 eta2 eta3'
    mob_name = L3
    kappa_names = 'kappa31 kappa32 kappa33'
  [../]
  [./penalty3]
    type = SwitchingFunctionPenalty
    variable = eta3
    etas    = 'eta1 eta2 eta3'
    h_names = 'h1   h2   h3'
  [../]
[]

[BCs]
  [./Periodic]
    [./All]
      auto_direction = 'x y'
    [../]
  [../]
[]

[Materials]
  # here we declare some of the model parameters: the mobilities and interface
  # gradient prefactors. For this example we use arbitrary numbers. In an actual simulation
  # physical mobilities would be used, and the interface gradient prefactors would
  # be readjusted to the free energy magnitudes.
  [./consts]
    type = GenericConstantMaterial
    prop_names  = 'M   kappa_c L1 L2 L3  kappa11 kappa12 kappa13 kappa21 kappa22 kappa23 kappa31 kappa32 kappa33'
    prop_values = '0.2 0.75    1  1  1   0.75    0.75    0.75    0.75    0.75    0.75    0.75    0.75    0.75   '
  [../]

  # This material sums up the individual phase contributions. It is written to the output file
  # (see GlobalParams section above) and can be used to check the constraint enforcement.
  [./etasummat]
    type = ParsedMaterial
    f_name = etasum
    args = 'eta1 eta2 eta3'
    material_property_names = 'h1 h2 h3'
    function = 'h1+h2+h3'
  [../]

  # The phase contribution factors for each material point are computed using the
  # SwitchingFunctionMaterials. Each phase with an order parameter eta contributes h(eta)
  # to the global free energy density. h is a function that switches smoothly from 0 to 1
  [./switching1]
    type = SwitchingFunctionMaterial
    function_name = h1
    eta = eta1
    h_order = SIMPLE
  [../]
  [./switching2]
    type = SwitchingFunctionMaterial
    function_name = h2
    eta = eta2
    h_order = SIMPLE
  [../]
  [./switching3]
    type = SwitchingFunctionMaterial
    function_name = h3
    eta = eta3
    h_order = SIMPLE
  [../]

  # The barrier function adds a phase transformation energy barrier. It also
  # Drives order parameters toward the [0:1] interval to avoid negative or larger than 1
  # order parameters (these are set to 0 and 1 contribution by the switching functions
  # above)
  [./barrier]
    type = MultiBarrierFunctionMaterial
    etas = 'eta1 eta2 eta3'
  [../]

  # We use DerivativeParsedMaterials to specify three (very) simple free energy
  # expressions for the three phases. All necessary derivatives are built automatically.
  # In a real problem these expressions can be arbitrarily complex (or even provided
  # by custom kernels).
  [./phase_free_energy_1]
    type = DerivativeParsedMaterial
    f_name = F1
    function = '(c-1)^2'
    args = 'c'
  [../]
  [./phase_free_energy_2]
    type = DerivativeParsedMaterial
    f_name = F2
    function = '(c-0.5)^2'
    args = 'c'
  [../]
  [./phase_free_energy_3]
    type = DerivativeParsedMaterial
    f_name = F3
    function = 'c^2'
    args = 'c'
  [../]

  # The DerivativeMultiPhaseMaterial ties the phase free energies together into a global free energy.
  # http://mooseframework.org/wiki/PhysicsModules/PhaseField/DevelopingModels/MultiPhaseModels/
  [./free_energy]
    type = DerivativeMultiPhaseMaterial
    f_name = F
    # we use a constant free energy (GeneriConstantmaterial property Fx)
    fi_names = 'F1  F2  F3'
    hi_names = 'h1  h2  h3'
    etas     = 'eta1 eta2 eta3'
    args = 'c'
    W = 1
  [../]
[]

[Postprocessors]
  # The total free energy of the simulation cell to observe the energy reduction.
  [./total_free_energy]
    type = ElementIntegralVariablePostprocessor
    variable = local_energy
  [../]

  # for testing we also monitor the total solute amount, which should be conserved.
  [./total_solute]
    type = ElementIntegralVariablePostprocessor
    variable = c
  [../]
[]

[Preconditioning]
  # This preconditioner makes sure the Jacobian Matrix is fully populated. Our
  # kernels compute all Jacobian matrix entries.
  # This allows us to use the Newton solver below.
  [./SMP]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient
  scheme = 'bdf2'

  # Automatic differentiation provedes a _full_ Jacobian in this example
  # so we can safely use NEWTON for a fast solve
  solve_type = 'NEWTON'

  l_max_its = 15
  l_tol = 1.0e-6

  nl_max_its = 50
  nl_rel_tol = 1.0e-6
  nl_abs_tol = 1.0e-6

  start_time = 0.0
  end_time   = 150.0

  [./TimeStepper]
    type = SolutionTimeAdaptiveDT
    dt = 0.1
  [../]
[]

[Debug]
  # show_var_residual_norms = true
[]

[Outputs]
  execute_on = 'timestep_end'
  exodus = true
  [./table]
    type = CSV
    delimiter = ' '
  [../]
[]

(modules/phase_field/test/tests/MultiPhase/penalty.i)

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 14
  ny = 10
  nz = 0
  xmin = 10
  xmax = 40
  ymin = 15
  ymax = 35
  elem_type = QUAD4
[]

[GlobalParams]
  penalty = 5
[]

[Variables]
  [./c]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = SmoothCircleIC
      x1 = 25.0
      y1 = 25.0
      radius = 6.0
      invalue = 0.9
      outvalue = 0.1
      int_width = 3.0
    [../]
  [../]
  [./w]
    order = FIRST
    family = LAGRANGE
  [../]

  [./eta1]
    order = FIRST
    family = LAGRANGE
    [./InitialCondition]
      type = SmoothCircleIC
      x1 = 30.0
      y1 = 25.0
      radius = 4.0
      invalue = 0.9
      outvalue = 0.1
      int_width = 2.0
    [../]
  [../]
  [./eta2]
    order = FIRST
    family = LAGRANGE
    initial_condition = 0.5
  [../]
[]

[Kernels]
  [./deta1dt]
    type = TimeDerivative
    variable = eta1
  [../]
  [./ACBulk1]
    type = AllenCahn
    variable = eta1
    args = 'c eta2'
    f_name = F
  [../]
  [./ACInterface1]
    type = ACInterface
    variable = eta1
    kappa_name = kappa_eta
  [../]
  [./penalty1]
    type = SwitchingFunctionPenalty
    variable = eta1
    etas    = 'eta1 eta2'
    h_names = 'h1   h2'
  [../]

  [./deta2dt]
    type = TimeDerivative
    variable = eta2
  [../]
  [./ACBulk2]
    type = AllenCahn
    variable = eta2
    args = 'c eta1'
    f_name = F
  [../]
  [./ACInterface2]
    type = ACInterface
    variable = eta2
    kappa_name = kappa_eta
  [../]
  [./penalty2]
    type = SwitchingFunctionPenalty
    variable = eta2
    etas    = 'eta1 eta2'
    h_names = 'h1   h2'
  [../]

  [./c_res]
    type = SplitCHParsed
    variable = c
    f_name = F
    kappa_name = kappa_c
    w = w
    args = 'eta1 eta2'
  [../]
  [./w_res]
    type = SplitCHWRes
    variable = w
    mob_name = M
  [../]
  [./time1]
    type = CoupledTimeDerivative
    variable = w
    v = c
  [../]
[]

[BCs]
  [./Periodic]
    [./All]
      auto_direction = 'x y'
    [../]
  [../]
[]

[Materials]
  [./consts]
    type = GenericConstantMaterial
    prop_names  = 'L kappa_eta'
    prop_values = '1 1        '
  [../]
  [./consts2]
    type = GenericConstantMaterial
    prop_names  = 'M kappa_c'
    prop_values = '1 1'
  [../]

  [./hsum]
    type = ParsedMaterial
    function = h1+h2
    f_name = hsum
    material_property_names = 'h1 h2'
    args = 'c'
    outputs = exodus
  [../]

  [./switching1]
    type = SwitchingFunctionMaterial
    function_name = h1
    eta = eta1
    h_order = SIMPLE
  [../]
  [./switching2]
    type = SwitchingFunctionMaterial
    function_name = h2
    eta = eta2
    h_order = SIMPLE
  [../]

  [./barrier]
    type = MultiBarrierFunctionMaterial
    etas = 'eta1 eta2'
  [../]

  [./free_energy_A]
    type = DerivativeParsedMaterial
    f_name = Fa
    args = 'c'
    function = '(c-0.1)^2'
    derivative_order = 2
  [../]
  [./free_energy_B]
    type = DerivativeParsedMaterial
    f_name = Fb
    args = 'c'
    function = '(c-0.9)^2'
    derivative_order = 2
  [../]

  [./free_energy]
    type = DerivativeMultiPhaseMaterial
    f_name = F
    fi_names = 'Fa   Fb'
    hi_names = 'h1   h2'
    etas     = 'eta1 eta2'
    args = 'c'
    derivative_order = 2
  [../]
[]

[Preconditioning]
  [./SMP]
    type = SMP
    full = true
  [../]
[]

[Executioner]
  type = Transient
  scheme = 'bdf2'

  solve_type = 'PJFNK'
  petsc_options_iname = '-pc_type -sub_pc_type'
  petsc_options_value = 'asm       lu'

  l_max_its = 15
  l_tol = 1.0e-6

  nl_max_its = 50
  nl_rel_tol = 1.0e-7
  nl_abs_tol = 1.0e-9

  start_time = 0.0
  num_steps = 2
  dt = 0.05
  dtmin = 0.01
[]

[Debug]
  # show_var_residual_norms = true
[]

[Outputs]
  execute_on = 'timestep_end'
  exodus = true
[]

Input Parameters
Input Files

Install MOOSE

New Users

Examples and Tutorials

Application Usage

Physics and Syntax

Application Development

Framework Development

MOOSEDocs

Infrastructure

Questions

Information and Tools

INL Applications and Remote Access