List of tutorials
All tutorials are written assuming that you are reasonably familiar with MOOSE. If you find most of the tutorials difficult to follow, please refer to the official MOOSE website for learning resources.
Tutorial 1: Multidomain state of a FE block
This tutorial (and others) covers the basic usage of the ferroelectric phase field method implemented in FERRET.
This simple problem considers solving the time-dependent LGD equation with just bulk and gradient contributions to the free energy in the presence of a depolarization phenomena (via solving the Poisson equation iteratively).
Consider a computational domain with a geometry which we define with the Mesh block.
[Mesh]
type = GeneratedMesh
dim = 3
nx = 32
ny = 32
nz = 6
xmin = -15
xmax = 15
ymin = -15
ymax = 15
zmin = -2
zmax = 2
elem_type = HEX8
[]
In general, the geometry defined in the 'Mesh' block never carries units. The length scale is introduced through Materials, Kernels, or other MOOSE objects. For this problem, the length scale is introduced through the units in the Materials objects that connect to the Kernels. Ferret uses a special base units system for a number of problems. This reduces the load quite extensively on the PETSc solvers to iterate the problem.
In the total free energy density, we include the bulk free energy density,
(1)
to the sixth order and the wall free energy density, ,
(2)
where and are suitable coefficients for lead-titanate. To solve the problem, we need the time-dependent Landau-Ginzburg equation,
(3)
which has variational derivatives of the total free energy,
(4)
The computations of the two relevant energy terms are introduced through a number of Kernels in the input files as below.
[Kernels]
#---------------------------------------#
# #
# Bulk (homogeneous) free energy #
# #
#---------------------------------------#
[./bed_x]
type = BulkEnergyDerivativeEighth
variable = polar_x
component = 0
[../]
[./bed_y]
type = BulkEnergyDerivativeEighth
variable = polar_y
component = 1
[../]
[./bed_z]
type = BulkEnergyDerivativeEighth
variable = polar_z
component = 2
[../]
#---------------------------------------#
# #
# Gradient free energy #
# #
#---------------------------------------#
[./walled_x]
type = WallEnergyDerivative
variable = polar_x
component = 0
[../]
[./walled_y]
type = WallEnergyDerivative
variable = polar_y
component = 1
[../]
[./walled_z]
type = WallEnergyDerivative
variable = polar_z
component = 2
[../]
[./walled2_x]
type = Wall2EnergyDerivative
variable = polar_x
component = 0
[../]
[./walled2_y]
type = Wall2EnergyDerivative
variable = polar_y
component = 1
[../]
[./walled2_z]
type = Wall2EnergyDerivative
variable = polar_z
component = 2
[../]
#---------------------------------------#
# #
# Poissons equation and P*E #
# #
#---------------------------------------#
[./polar_x_electric_E]
type = PolarElectricEStrong
variable = potential_int
[../]
[./FE_E_int]
type = Electrostatics
variable = potential_int
[../]
[./polar_electric_px]
type = PolarElectricPStrong
variable = polar_x
component = 0
[../]
[./polar_electric_py]
type = PolarElectricPStrong
variable = polar_y
component = 1
[../]
[./polar_electric_pz]
type = PolarElectricPStrong
variable = polar_z
component = 2
[../]
#---------------------------------------#
# #
# Time dependence of Pj #
# #
#---------------------------------------#
[./polar_x_time]
type = TimeDerivativeScaled
variable = polar_x
time_scale = 1.0
[../]
[./polar_y_time]
type = TimeDerivativeScaled
variable = polar_y
time_scale = 1.0
[../]
[./polar_z_time]
type = TimeDerivativeScaled
variable = polar_z
time_scale = 1.0
[../]
[]
We need one Kernel for each variable (, and ). The computation of Eq. \ref{varderiv} is split into three Kernel contributions for each of the three components of . The FERRET Syntax page details the algebra needed to turn these terms (BulkEnergyDerivativeEighth, WallEnergyDerivative, and Wall2EnergyDerivative) into weak forms.
The reader will notice that we also have included terms related to the free energy due to an applied field in PolarElectricPStrong and also terms related to the coupling to the electric field. The Poisson equation for our system reads,
(5)
which is satisfied at every step of the time evolution. The left hand side is computed by Electrostatics and the right hand side computed by PolarElectricPStrong. Note that corresponds to a background dielectric strength which is assigned to high frequency (core) electrons. Next, we have the time derivatives of , and with TimeDerivativeScaled. We implement some reduction of input file length by using the GlobalParams block.
[GlobalParams]
len_scale = 1.0
polar_x = polar_x
polar_y = polar_y
polar_z = polar_z
potential_E_int = potential_int
[]
which adds these lines automatically in the relevant Kernel and Materials objects. To seed the relevant materials coefficients to the problem, we use the Materials system in MOOSE,
[Materials]
## Use coefficients for lead-titanate in this convienient unit system (aC, kg, sec, nm)
[./Landau_P_FE]
type = GenericConstantMaterial
prop_names = 'alpha1 alpha11 alpha12 alpha111 alpha112 alpha123 alpha1111 alpha1112 alpha1122 alpha1123'
prop_values = '-0.1722883 -0.073 0.75 0.26 0.61 -3.67 0.0 0.0 0.0 0.0'
block = '0'
[../]
[./Landau_G_FE]
type = GenericConstantMaterial
prop_names = 'G110 G11_G110 G12_G110 G44_G110 G44P_G110'
prop_values = '0.173 0.6 0.0 0.3 0.3'
block = '0'
[../]
[./permitivitty_1]
###############################################
##
## so-called background dielectric constant
## (it encapsulates the motion of core electrons
## at high frequency) = e_b*e_0 (here we use
## e_b = 10), see PRB. 74, 104014, (2006)
## This does not influence results much.
##
###############################################
type = GenericConstantMaterial
prop_names = 'permittivity'
prop_values = '0.08854187'
[../]
[]
In this problem, the boundary conditions are left empty. This is equivalently known as a natural boundary condition in MOOSe. As such, large depolarizing fields are generated from the surface charges () where is some surface normal. The action of the depolarization drives the polar topology into a classic flux-closure domain pattern as shown in the below figure.
Figure 1: Classic flux-closure domain pattern solution to the first tutorial.
A few other FERRET and MOOSE objects exist in the input file. We use the Postprocessors system to track various aspects of the simulation.
[Postprocessors]
[./avePz]
type = ElementAverageValue
variable = polar_z
execute_on = 'timestep_end'
[../]
[./FbP]
type = BulkEnergyEighth
execute_on = 'initial timestep_end'
[../]
[./FgP]
type = WallEnergy
execute_on = 'initial timestep_end'
[../]
[./Ftot]
type = LinearCombinationPostprocessor
pp_names = 'FbP FgP'
pp_coefs = ' 1 1'
execute_on = 'timestep_end'
##########################################
#
# NOTE: Ferret output is in attojoules
#
##########################################
[../]
[./perc_change]
type = PercentChangePostprocessor
postprocessor = Ftot
[../]
[]
The bulk and gradient energy contributions are listed (computed volume integral over the simulation box). Their sum is computed with the LinearCombinationPostprocessor to be equal to . An astute user should ask why the electric field contribution is also not included here. We leave that as an exercise for the user to verify that the total energy with this contribution still is decreasing with every time step as gauranteed by the mathematical form of Eq. \refLGD.
Then, the relative change of this energy between adjacent time steps is calculated with PercentChangePostprocessor as,
(6)
The Terminator located in the UserObjects block
[UserObjects]
[./kill]
type = Terminator
expression = 'perc_change <= 1.0e-4'
[../]
[]
which tracks when the postprocessed value is less than . This is a useful tool to end the simulation when the suspected ground state has been reached. In general, this value seems to be sufficient for these types of problems.
Finally, we share our PETSc and Executioner options that seem to be the most efficient for polar domain prediction problems.
[Preconditioning]
[./smp]
type = SMP
full = true
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-ksp_gmres_restart -snes_atol -snes_rtol -ksp_rtol -pc_type'
petsc_options_value = ' 120 1e-10 1e-8 1e-4 bjacobi'
[../]
[]
and
[Executioner]
type = Transient
solve_type = 'NEWTON'
scheme = 'implicit-euler'
dtmax = 0.7
[]
These flags can be optimized to lead to less wall clock times but this is a subject of future work.