- componentAn integer corresponding to the direction the variable this kernel acts in. (0 for x, 1 for y, 2 for z)
C++ Type:unsigned int
Controllable:No
Description:An integer corresponding to the direction the variable this kernel acts in. (0 for x, 1 for y, 2 for z)
- polar_xThe x component of the polarization
C++ Type:std::vector<VariableName>
Controllable:No
Description:The x component of the polarization
- polar_yThe y component of the polarization
C++ Type:std::vector<VariableName>
Controllable:No
Description:The y component of the polarization
- volThe name of the volume used for volume normalization
C++ Type:PostprocessorName
Controllable:No
Description:The name of the volume used for volume normalization
DomainVariantPopulation
Calculates the fraction of volume of a given domain population (only works in tetragonal phase at the moment)
Overview
Computes the fraction of the domain variation population in the compuational box of volume ,
This expression is only valid for tetragonal ferroelectrics but we plan to generalize this to orthorhombic and rhombohedral symmetries in the future.
Example Input File Syntax
Input Parameters
- blockThe list of blocks (ids or names) that this object will be applied
C++ Type:std::vector<SubdomainName>
Controllable:No
Description:The list of blocks (ids or names) that this object will be applied
- execute_onTIMESTEP_ENDThe list of flag(s) indicating when this object should be executed, the available options include NONE, INITIAL, LINEAR, NONLINEAR, TIMESTEP_END, TIMESTEP_BEGIN, MULTIAPP_FIXED_POINT_END, MULTIAPP_FIXED_POINT_BEGIN, FINAL, CUSTOM, ALWAYS.
Default:TIMESTEP_END
C++ Type:ExecFlagEnum
Options:NONE, INITIAL, LINEAR, NONLINEAR, TIMESTEP_END, TIMESTEP_BEGIN, MULTIAPP_FIXED_POINT_END, MULTIAPP_FIXED_POINT_BEGIN, FINAL, CUSTOM, ALWAYS, TRANSFER
Controllable:No
Description:The list of flag(s) indicating when this object should be executed, the available options include NONE, INITIAL, LINEAR, NONLINEAR, TIMESTEP_END, TIMESTEP_BEGIN, MULTIAPP_FIXED_POINT_END, MULTIAPP_FIXED_POINT_BEGIN, FINAL, CUSTOM, ALWAYS.
- polar_zThe z component of the polarization
C++ Type:std::vector<VariableName>
Controllable:No
Description:The z component of the polarization
- prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
Optional Parameters
- allow_duplicate_execution_on_initialFalseIn the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable).
Default:False
C++ Type:bool
Controllable:No
Description:In the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable).
- 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:Yes
Description:Set the enabled status of the MooseObject.
- execution_order_group0Execution order groups are executed in increasing order (e.g., the lowest number is executed first). Note that negative group numbers may be used to execute groups before the default (0) group. Please refer to the user object documentation for ordering of user object execution within a group.
Default:0
C++ Type:int
Controllable:No
Description:Execution order groups are executed in increasing order (e.g., the lowest number is executed first). Note that negative group numbers may be used to execute groups before the default (0) group. Please refer to the user object documentation for ordering of user object execution within a group.
- force_postauxFalseForces the UserObject to be executed in POSTAUX
Default:False
C++ Type:bool
Controllable:No
Description:Forces the UserObject to be executed in POSTAUX
- force_preauxFalseForces the UserObject to be executed in PREAUX
Default:False
C++ Type:bool
Controllable:No
Description:Forces the UserObject to be executed in PREAUX
- force_preicFalseForces the UserObject to be executed in PREIC during initial setup
Default:False
C++ Type:bool
Controllable:No
Description:Forces the UserObject to be executed in PREIC during initial setup
- implicitTrueDetermines whether this object is calculated using an implicit or explicit form
Default:True
C++ Type:bool
Controllable:No
Description:Determines whether this object is calculated using an implicit or explicit form
- outputsVector of output names where you would like to restrict the output of variables(s) associated with this object
C++ Type:std::vector<OutputName>
Controllable:No
Description:Vector of output names where you would like to restrict the output of variables(s) associated with this object
- seed0The seed for the master random number generator
Default:0
C++ Type:unsigned int
Controllable:No
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
Controllable:No
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
Input Files
- (test/tests/film/PTO_film3_T298K.i)
- (examples/domain_walls/BTO_wall_T298K.i)
- (examples/monodomain/BTO_monodomain_Tdef.i)
- (tutorial/film.i)
- (examples/films/PZTfilm_e12_T298K_E0.i)
- (examples/monodomain/PZT_monodomain_Tdef.i)
- (examples/films/PTOfilm_e12_T298K_E0.i)
- (examples/monodomain/PTO_monodomain_Tdef.i)
- (test/tests/film/PTO_film3_T298K_noP.i)
(test/tests/film/PTO_film3_T298K.i)
[Mesh]
[gen]
############################################
##=
## Type and dimension of the mesh
##=
############################################
type = GeneratedMeshGenerator
dim = 3
#############################################
##
## Grid definition. Note that it should be
## nJ = 2*(Jmax-Jmin) for J = x, y, z
##=
#############################################
nx = 4
ny = 4
nz = 18 #need approximately 90 processors for this since dof is 953133 and we want 10k dof/proc
#############################################
##
## Actual spatial coordinates of mesh.
## Jmax - Jmin = nJ/2 for J = x, y, z
## Units are in nanometers
##
#############################################
xmin = -2.0
xmax = 2.0
ymin = -2.0
ymax = 2.0
zmin = -10.0
zmax = 8.0
#############################################
##
## FE type/order (hexahedral, tetrahedral
##
#############################################
elem_type = HEX8
[]
[./cnode]
input = gen
############################################
##
## additional boundary sideset (one node)
## to zero one of the elastic displacement vectors
## vectors and eliminates rigid body translations
## from the degrees of freedom
##
## NOTE: This must conform with the about
## [Mesh] block settings
##=
############################################
type = ExtraNodesetGenerator
coord = '-2.0 -2.0 -10.0'
new_boundary = 100
[../]
[subdomains]
type = SubdomainBoundingBoxGenerator
input = cnode
bottom_left = '-2.0 -2.0 -10.0'
block_id = 1
top_right = '2.0 2.0 0'
location = INSIDE
[]
[film_interface]
type = SideSetsBetweenSubdomainsGenerator
input = subdomains
primary_block = 0
paired_block = 1
new_boundary = 52
[]
[]
[GlobalParams]
len_scale = 1.0
polar_x = polar_x
polar_y = polar_y
polar_z = polar_z
potential_E_int = potential_E_int
displacements = 'u_x u_y u_z'
vol = vol
##############################################
##=
## IMPORTANT(!): Units in Ferret are nm, kg,
## seconds, and attocoulombs
##
##############################################
u_x = u_x
u_y = u_y
u_z = u_z
[]
[Variables]
#################################
##
## Variable definitions
## P, u, phi, e^global_ij
## and their initial conditions
##
#################################
[./global_strain]
order = SIXTH
family = SCALAR
[../]
[./polar_x]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.1e-5
max = 0.1e-5
[../]
block = '0 1'
[../]
[./polar_y]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.1e-5
max = 0.1e-5
[../]
block = '0 1'
[../]
[./polar_z]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.1e-5
max = 0.1e-5
[../]
block = '0 1'
[../]
[./potential_E_int]
order = FIRST
family = LAGRANGE
block = '0 1'
[../]
[./u_x]
order = FIRST
family = LAGRANGE
block = '0 1'
[../]
[./u_y]
order = FIRST
family = LAGRANGE
block = '0 1'
[../]
[./u_z]
order = FIRST
family = LAGRANGE
block = '0 1'
[../]
[]
[AuxVariables]
######################################
##
## Auxiarilly variable definitions
## (can be intermediate variables
## or for postprocessed quantities)
##
######################################
######################################
##
## Global displacements
##
######################################
[./disp_x]
block = '0 1'
[../]
[./disp_y]
block = '0 1'
[../]
[./disp_z]
block = '0 1'
[../]
######################################
##
## Stress/strain tensor components
##
######################################
[./e00]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e01]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e10]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e11]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e12]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e22]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s00]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s01]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s10]
order = CONSTANT
family = MONOMIAL
[../]
[./s11]
order = CONSTANT
family = MONOMIAL
[../]
[./s12]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s22]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[]
[AuxKernels]
######################################
##
## Auxiarilly Kernel definitions
## (can be intermediate "operations"
## or for postprocessed quantities)
##
######################################
[./disp_x]
type = GlobalDisplacementAux
variable = disp_x
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 0
use_displaced_mesh = false
[../]
[./disp_y]
type = GlobalDisplacementAux
variable = disp_y
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 1
use_displaced_mesh = false
[../]
[./disp_z]
type = GlobalDisplacementAux
variable = disp_z
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 2
use_displaced_mesh = false
[../]
[./e00]
type = RankTwoAux
variable = e00
rank_two_tensor = total_strain
index_i = 0
index_j = 0
[../]
[./e01]
type = RankTwoAux
variable = e01
rank_two_tensor = total_strain
index_i = 0
index_j = 1
[../]
[./e10]
type = RankTwoAux
variable = e10
rank_two_tensor = total_strain
index_i = 1
index_j = 0
[../]
[./e12]
type = RankTwoAux
variable = e12
rank_two_tensor = total_strain
index_i = 1
index_j = 2
[../]
[./e11]
type = RankTwoAux
variable = e11
rank_two_tensor = total_strain
index_i = 1
index_j = 1
[../]
[./e22]
type = RankTwoAux
variable = e22
rank_two_tensor = total_strain
index_i = 2
index_j = 2
[../]
[./s00]
type = RankTwoAux
variable = s00
rank_two_tensor = stress
index_i = 0
index_j = 0
[../]
[./s01]
type = RankTwoAux
variable = s01
rank_two_tensor = stress
index_i = 0
index_j = 1
[../]
[./s10]
type = RankTwoAux
variable = s10
rank_two_tensor = stress
index_i = 1
index_j = 0
[../]
[./s12]
type = RankTwoAux
variable = s12
rank_two_tensor = stress
index_i = 1
index_j = 2
[../]
[./s11]
type = RankTwoAux
variable = s11
rank_two_tensor = stress
index_i = 1
index_j = 1
[../]
[./s22]
type = RankTwoAux
variable = s22
rank_two_tensor = stress
index_i = 2
index_j = 2
[../]
[]
[ScalarKernels]
######################################
##
## Necessary for PBC system
##
######################################
[./global_strain]
type = GlobalStrain
variable = global_strain
global_strain_uo = global_strain_uo
use_displaced_mesh = false
[../]
[]
[Materials]
#################################################
##
## add comments
##
## NOTE: there might be some Legendre transforms
## depending on what approach you use
## -i.e. inhomogeneous strain vs
## homogeneous strain [renormalized]
##=
##################################################
[./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_P_substr]
type = GenericConstantMaterial
prop_names = 'alpha1 alpha11 alpha12 alpha111 alpha112 alpha123 alpha1111 alpha1112 alpha1122 alpha1123'
prop_values = '10.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0'
block = '1'
[../]
############################################
##
## add comments
##
############################################
[./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'
[../]
[./mat_C_FE]
type = GenericConstantMaterial
prop_names = 'C11 C12 C44'
prop_values = '175.0 79.4 111.1'
block = '0'
[../]
[./mat_C_sub]
type = GenericConstantMaterial
prop_names = 'C11 C12 C44'
prop_values = '220.0 34.4 161.1'
block = '1'
[../]
##################################################
##=
## NOTE: Sign convention in Ferret for the
## electrostrictive coeff. is multiplied by
## an overall factor of (-1)
##
##################################################
[./mat_Q]
type = GenericConstantMaterial
prop_names = 'Q11 Q12 Q44'
prop_values = '-0.089 0.026 -0.03375'
block = '0 1'
[../]
[./mat_q]
type = GenericConstantMaterial
prop_names = 'q11 q12 q44'
prop_values = '-11.4 -0.01438 -7.5'
[../]
[./eigen_strain]
type = ComputeEigenstrain
# eigen_base = 'exx exy exz eyx eyy eyz ezx ezy ezz'
eigen_base = '1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0'
eigenstrain_name = eigenstrain
prefactor = 0.0
[../]
[./elasticity_tensor_1]
type = ComputeElasticityTensor
fill_method = symmetric9
###############################################
##
## symmetric9 fill_method is (default)
## C11 C12 C13 C22 C23 C33 C44 C55 C66
##
###############################################
C_ijkl = '175.0 79.4 79.4 175.0 79.4 175.0 111.1 111.1 111.1'
[../]
[./strain_1]
type = ComputeSmallStrain
global_strain = global_strain
eigenstrain_names = eigenstrain
[../]
[./stress_1]
type = ComputeLinearElasticStress
[../]
[./global_strain]
type = ComputeGlobalStrain
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
[../]
[./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)
##
###############################################
type = GenericConstantMaterial
prop_names = 'permittivity'
prop_values = '0.08854187'
[../]
[]
[Kernels]
###############################################
##
## Physical Kernel operators
## to enforce TDLGD evolution
##
###############################################
#Elastic problem
[./TensorMechanics]
use_displaced_mesh = false
eigenstrain_names = eigenstrain
[../]
[./bed_x]
type = BulkEnergyDerivativeEighth
variable = polar_x
component = 0
block = '0 1'
[../]
[./bed_y]
type = BulkEnergyDerivativeEighth
variable = polar_y
component = 1
block = '0 1'
[../]
[./bed_z]
type = BulkEnergyDerivativeEighth
variable = polar_z
component = 2
block = '0 1'
[../]
[./walled_x]
type = WallEnergyDerivative
variable = polar_x
component = 0
block = '0'
[../]
[./walled_y]
type = WallEnergyDerivative
variable = polar_y
component = 1
block = '0'
[../]
[./walled_z]
type = WallEnergyDerivative
variable = polar_z
component = 2
block = '0'
[../]
[./electrostr_ux]
type = ElectrostrictiveCouplingDispDerivative
variable = u_x
component = 0
block = '0 1'
[../]
[./electrostr_uy]
type = ElectrostrictiveCouplingDispDerivative
variable = u_y
component = 1
block = '0 1'
[../]
[./electrostr_uz]
type = ElectrostrictiveCouplingDispDerivative
variable = u_z
component = 2
block = '0 1'
[../]
[./electrostr_polar_coupled_x]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_x
component = 0
block = '0 1'
[../]
[./electrostr_polar_coupled_y]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_y
component = 1
block = '0 1'
[../]
[./electrostr_polar_coupled_z]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_z
component = 2
block = '0 1'
[../]
[./polar_x_electric_E]
type = PolarElectricEStrong
variable = potential_E_int
block = '0 1'
[../]
[./FE_E_int]
type = Electrostatics
variable = potential_E_int
block = '0 1'
[../]
[./polar_electric_px]
type = PolarElectricPStrong
variable = polar_x
component = 0
block = '0 1'
[../]
[./polar_electric_py]
type = PolarElectricPStrong
variable = polar_y
component = 1
block = '0 1'
[../]
[./polar_electric_pz]
type = PolarElectricPStrong
variable = polar_z
component = 2
block = '0 1'
[../]
[./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
[../]
[]
[BCs]
[./Periodic]
[./xy]
auto_direction = 'x y'
variable = 'u_x u_y u_z polar_x polar_y polar_z'
[../]
[../]
[./boundary_grounding]
type = DirichletBC
boundary = '52'
variable = potential_E_int
value = 0.0
[../]
# fix center point location
[./centerfix_x]
type = DirichletBC
boundary = 'back'
variable = u_x
value = 0
[../]
[./centerfix_y]
type = DirichletBC
boundary = 'back'
variable = u_y
value = 0
[../]
[./centerfix_z]
type = DirichletBC
boundary = 'back'
variable = u_z
value = 0
[../]
[]
[Postprocessors]
###############################################
##=
## Postprocessors (integrations over the
## computational domain) to calculate the total energy
## decomposed into linear combinations of the
## different physics.
##
###############################################
[./Fbulk]
type = BulkEnergyEighth
execute_on = 'timestep_end'
block = '0'
[../]
[./Fwall]
type = WallEnergy
execute_on = 'timestep_end'
block = '0'
[../]
[./Felastic]
type = ElasticEnergy
execute_on = 'timestep_end'
use_displaced_mesh = false
block = '0'
[../]
[./Fcoupled]
type = ElectrostrictiveCouplingEnergy
execute_on = 'timestep_end'
u_x = u_x
u_y = u_y
u_z = u_z
block = '0'
[../]
[./Felec]
type = ElectrostaticEnergy
execute_on = 'timestep_end'
block = '0'
[../]
[./Ftotal]
type = LinearCombinationPostprocessor
pp_names = 'Fbulk Fwall Fcoupled Felec'
pp_coefs = ' 1 1 1 1'
execute_on = 'timestep_end'
[../]
[./vol]
type = VolumePostprocessor
execute_on = 'timestep_end'
[../]
[./px]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 0
block = '0'
[../]
[./py]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 1
block = '0'
[../]
[./pz]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 2
block = '0'
[../]
[./perc_change]
type = PercentChangePostprocessor
postprocessor = Ftotal
execute_on = 'timestep_end'
[../]
[]
[UserObjects]
###############################################
##=
## GlobalStrain system to enforce periodicity
## in the anisotropic strain field
##
###############################################
[./global_strain_uo]
type = GlobalATiO3MaterialRVEUserObject
use_displaced_mesh = false
execute_on = 'Initial Linear Nonlinear'
applied_stress_tensor = '4.375 4.375 3.97 0.0 0.0 0.0'
block = '0 1'
[../]
###############################################
##
## terminator to end energy evolution when the energy difference
## between subsequent time steps is lower than 5e-6
##
## NOTE: can fail if the time step is small
##
###############################################
[./kill]
type = Terminator
expression = 'perc_change <= 1.0e-6'
[../]
[]
[Preconditioning]
###############################################
##
## Numerical preconditioning/solver options
##=
###############################################
[./smp]
type = SMP
full = true
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-ksp_gmres_restart -snes_atol -snes_rtol -ksp_rtol -pc_type -build_twosided'
petsc_options_value = ' 160 1e-8 1e-8 1e-5 bjacobi allreduce'
[../]
[]
[Executioner]
##########################################
##
## Time integration/solver options
##
##########################################
type = Transient
solve_type = 'PJFNK'
scheme = 'implicit-euler'
dtmin = 1e-13
dtmax = 0.6
l_max_its = 200
[./TimeStepper]
type = IterationAdaptiveDT
optimal_iterations = 6
growth_factor = 1.2
cutback_factor = 0.75
linear_iteration_ratio = 1000
dt = 0.6
[../]
verbose = true
nl_max_its = 20
num_steps = 5
[]
[Outputs]
###############################################
##===============
## Output options
##
###############################################
print_linear_residuals = false
perf_graph = false
[./out]
type = Exodus
file_base = out_pto_film_test
elemental_as_nodal = true
interval = 1
[../]
[]
(examples/domain_walls/BTO_wall_T298K.i)
[Mesh]
[gen]
############################################
##
## Type and dimension of the mesh
##
############################################
type = GeneratedMeshGenerator
dim = 3
nx = 60
ny = 2
nz = 2
#############################################
##
## Actual spatial coordinates of mesh.
## Jmax - Jmin = nJ/2 for J = x, y, z
## Units are in nanometers
##
#############################################
xmin = -30.0
xmax = 30.0
ymin = -1.0
ymax = 1.0
zmin = -1.0
zmax = 1.0
#############################################
##
## FE type/order (hexahedral, tetrahedral
##
#############################################
elem_type = HEX8
[]
[./cnode]
input = gen
############################################
##
## additional boundary sideset (one node)
## to zero one of the elastic displacement vectors
## vectors and eliminates rigid body translations
## from the degrees of freedom
##
## NOTE: This must conform with the about
## [Mesh] block settings
##
############################################
type = ExtraNodesetGenerator
coord = '0.0 0.0 0.0'
new_boundary = 100
[../]
[]
[GlobalParams]
len_scale = 1.0
polar_x = polar_x
polar_y = polar_y
polar_z = polar_z
potential_E_int = potential_E_int
vol = vol
displacements = 'u_x u_y u_z'
##############################################
##
## IMPORTANT(!): Units in Ferret are nm, kg,
## seconds, and attocoulombs
##
##############################################
u_x = u_x
u_y = u_y
u_z = u_z
[]
[Functions]
[./stripe1]
type = ParsedFunction
value = 0.1*cos(0.10471975512*(x+2)) #2pi/L = 0.10471975512
[../]
[]
[Variables]
#################################
##
## Variable definitions
## P, u, phi, e^global_ij
## and their initial conditions
##
#################################
[./global_strain]
order = SIXTH
family = SCALAR
[../]
[./polar_x]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.01e-4
max = 0.01e-4
[../]
[../]
[./polar_y]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.01e-4
max = 0.01e-4
[../]
[../]
[./polar_z]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = FunctionIC
function = 'stripe1'
[../]
[../]
[./potential_E_int]
order = FIRST
family = LAGRANGE
[../]
[./u_x]
order = FIRST
family = LAGRANGE
[../]
[./u_y]
order = FIRST
family = LAGRANGE
[../]
[./u_z]
order = FIRST
family = LAGRANGE
[../]
[]
[AuxVariables]
######################################
##
## Auxiarilly variable definitions
## (can be intermediate variables
## or for postprocessed quantities)
##
######################################
######################################
##
## Global displacements
##
######################################
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
######################################
##
## Stress/strain tensor components
##
######################################
[./e00]
order = CONSTANT
family = MONOMIAL
[../]
[./e01]
order = CONSTANT
family = MONOMIAL
[../]
[./e10]
order = CONSTANT
family = MONOMIAL
[../]
[./e11]
order = CONSTANT
family = MONOMIAL
[../]
[./e12]
order = CONSTANT
family = MONOMIAL
[../]
[./e22]
order = CONSTANT
family = MONOMIAL
[../]
[./s00]
order = CONSTANT
family = MONOMIAL
[../]
[./s01]
order = CONSTANT
family = MONOMIAL
[../]
[./s10]
order = CONSTANT
family = MONOMIAL
[../]
[./s11]
order = CONSTANT
family = MONOMIAL
[../]
[./s12]
order = CONSTANT
family = MONOMIAL
[../]
[./s22]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
######################################
##
## Auxiarilly Kernel definitions
## (can be intermediate "operations"
## or for postprocessed quantities)
##
######################################
[./disp_x]
type = GlobalDisplacementAux
variable = disp_x
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 0
use_displaced_mesh = false
[../]
[./disp_y]
type = GlobalDisplacementAux
variable = disp_y
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 1
use_displaced_mesh = false
[../]
[./disp_z]
type = GlobalDisplacementAux
variable = disp_z
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 2
use_displaced_mesh = false
[../]
[./e00]
type = RankTwoAux
variable = e00
rank_two_tensor = total_strain
index_i = 0
index_j = 0
[../]
[./e01]
type = RankTwoAux
variable = e01
rank_two_tensor = total_strain
index_i = 0
index_j = 1
[../]
[./e10]
type = RankTwoAux
variable = e10
rank_two_tensor = total_strain
index_i = 1
index_j = 0
[../]
[./e12]
type = RankTwoAux
variable = e12
rank_two_tensor = total_strain
index_i = 1
index_j = 2
[../]
[./e11]
type = RankTwoAux
variable = e11
rank_two_tensor = total_strain
index_i = 1
index_j = 1
[../]
[./e22]
type = RankTwoAux
variable = e22
rank_two_tensor = total_strain
index_i = 2
index_j = 2
[../]
[./s00]
type = RankTwoAux
variable = s00
rank_two_tensor = stress
index_i = 0
index_j = 0
[../]
[./s01]
type = RankTwoAux
variable = s01
rank_two_tensor = stress
index_i = 0
index_j = 1
[../]
[./s10]
type = RankTwoAux
variable = s10
rank_two_tensor = stress
index_i = 1
index_j = 0
[../]
[./s12]
type = RankTwoAux
variable = s12
rank_two_tensor = stress
index_i = 1
index_j = 2
[../]
[./s11]
type = RankTwoAux
variable = s11
rank_two_tensor = stress
index_i = 1
index_j = 1
[../]
[./s22]
type = RankTwoAux
variable = s22
rank_two_tensor = stress
index_i = 2
index_j = 2
[../]
[]
[ScalarKernels]
######################################
##
## Necessary for PBC system
##
######################################
[./global_strain]
type = GlobalStrain
variable = global_strain
global_strain_uo = global_strain_uo
use_displaced_mesh = false
[../]
[]
[Materials]
#################################################
##
## Bulk free energy and electrostrictive
## coefficients gleaned from
## Marton and Hlinka
## Phys. Rev. B. 74, 104014, (2006)
##
## NOTE: there might be some Legendre transforms
## depending on what approach you use
## -i.e. inhomogeneous strain vs
## homogeneous strain [renormalized]
##
##################################################
[./Landau_P]
type = GenericConstantMaterial
prop_names = 'alpha1 alpha11 alpha12 alpha111 alpha112 alpha123 alpha1111 alpha1112 alpha1122 alpha1123'
prop_values = '-0.027721 -0.64755 0.323 8.004 4.47 4.91 0.0 0.0 0.0 0.0'
[../]
############################################
##
## Gradient coefficients from
## Marton and Hlinka
## Phys. Rev. B. 74, 104014, (2006)
##
############################################
[./Landau_G]
type = GenericConstantMaterial
prop_names = 'G110 G11_G110 G12_G110 G44_G110 G44P_G110'
prop_values = '0.5 0.51 -0.02 0.02 0.0'
[../]
[./mat_C]
type = GenericConstantMaterial
prop_names = 'C11 C12 C44'
prop_values = '275.0 179.0 54.3'
[../]
##################################################
##
## NOTE: Sign convention in Ferret for the
## electrostrictive coeff. is multiplied by
## an overall factor of (-1)
##
##################################################
[./mat_Q]
type = GenericConstantMaterial
prop_names = 'Q11 Q12 Q44'
prop_values = '-0.11 0.045 -0.029'
[../]
[./mat_q]
type = GenericConstantMaterial
prop_names = 'q11 q12 q44'
prop_values = '-14.2 0.74 -1.57'
[../]
[./eigen_strain]
type = ComputeEigenstrain
eigen_base = '0. 0 0 0 0 0 0 0 0'
eigenstrain_name = eigenstrain
prefactor = 0.0
[../]
[./elasticity_tensor_1]
type = ComputeElasticityTensor
fill_method = symmetric9
###############################################
##
## symmetric9 fill_method is (default)
## C11 C12 C13 C22 C23 C33 C44 C55 C66
##
###############################################
C_ijkl = '275.0 179.0 179.0 275.0 179.0 275.0 54.3 54.3 54.3'
[../]
[./strain_1]
type = ComputeSmallStrain
global_strain = global_strain
eigenstrain_names = eigenstrain
[../]
[./stress_1]
type = ComputeLinearElasticStress
[../]
[./global_strain]
type = ComputeGlobalStrain
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
[../]
[./slab_ferroelectric]
type = ComputeElectrostrictiveTensor
Q_mnkl = '-0.11 0.045 0.045 -0.11 0.045 -0.11 -0.029 -0.029 -0.029'
C_ijkl = '275.0 179.0 179.0 275.0 179.0 275.0 54.3 54.3 54.3'
[../]
[./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)
##
###############################################
type = GenericConstantMaterial
prop_names = 'permittivity'
prop_values = '0.08854187'
[../]
[]
[Kernels]
###############################################
##
## Physical Kernel operators
## to enforce TDLGD evolution
##
###############################################
#Elastic problem
[./TensorMechanics]
use_displaced_mesh = false
[../]
[./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
[../]
[./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
[../]
[./electrostr_ux]
type = ElectrostrictiveCouplingDispDerivative
variable = u_x
component = 0
[../]
[./electrostr_uy]
type = ElectrostrictiveCouplingDispDerivative
variable = u_y
component = 1
[../]
[./electrostr_uz]
type = ElectrostrictiveCouplingDispDerivative
variable = u_z
component = 2
[../]
[./electrostr_polar_coupled_x]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_x
component = 0
[../]
[./electrostr_polar_coupled_y]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_y
component = 1
[../]
[./electrostr_polar_coupled_z]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_z
component = 2
[../]
[./polar_x_electric_E]
type = PolarElectricEStrong
variable = potential_E_int
[../]
[./FE_E_int]
type = Electrostatics
variable = potential_E_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
[../]
[./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
[../]
[]
[BCs]
[./Periodic]
[./xyz]
auto_direction = 'x y z'
variable = 'u_x u_y u_z polar_x polar_y polar_z potential_E_int'
[../]
[../]
# fix center point location
[./centerfix_x]
type = DirichletBC
boundary = 100
variable = u_x
value = 0
[../]
[./centerfix_y]
type = DirichletBC
boundary = 100
variable = u_y
value = 0
[../]
[./centerfix_z]
type = DirichletBC
boundary = 100
variable = u_z
value = 0
[../]
[]
[Postprocessors]
###############################################
##
## Postprocessors (integrations over the
## computational domain) to calculate the total energy
## decomposed into linear combinations of the
## different physics.
##
###############################################
[./Fb]
type = BulkEnergyEighth
execute_on = 'initial timestep_end'
[../]
[./Fw]
type = WallEnergy
execute_on = 'initial timestep_end'
[../]
[./Fela]
type = ElasticEnergy
execute_on = 'initial timestep_end'
use_displaced_mesh = false
[../]
[./Fc]
type = ElectrostrictiveCouplingEnergy
execute_on = 'initial timestep_end'
[../]
[./Fele]
type = ElectrostaticEnergy
execute_on = 'initial timestep_end'
[../]
[./Ftot]
type = LinearCombinationPostprocessor
pp_names = 'Fb Fw Fc Fele'
pp_coefs = ' 1 1 1 1'
execute_on = 'initial timestep_end'
[../]
[./vol]
type = VolumePostprocessor
execute_on = 'initial timestep_end'
[../]
[./px]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 0
[../]
[./py]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 1
[../]
[./pz]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 2
[../]
[./perc_change]
type = PercentChangePostprocessor
postprocessor = Ftot
execute_on = 'initial timestep_end'
[../]
[]
[UserObjects]
###############################################
##
## GlobalStrain system to enforce periodicity
## in the anisotropic strain field
##
###############################################
[./global_strain_uo]
type = GlobalATiO3MaterialRVEUserObject
use_displaced_mesh = false
execute_on = 'Initial Linear Nonlinear'
[../]
###############################################
##
## terminator to end energy evolution when the energy difference
## between subsequent time steps is lower than 5e-6
##
## NOTE: can fail if the time step is small
##
###############################################
[./kill]
type = Terminator
expression = 'perc_change <= 1.0e-6'
[../]
[]
[Preconditioning]
###############################################
##
## Numerical preconditioning/solver options
##
###############################################
[./smp]
type = SMP
full = true
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-ksp_gmres_restart -snes_atol -snes_rtol -ksp_rtol -pc_type -build_twosided'
petsc_options_value = ' 160 1e-10 1e-8 1e-6 bjacobi allreduce'
[../]
[]
[Executioner]
##########################################
##
## Time integration/solver options
##
##########################################
type = Transient
solve_type = 'PJFNK'
scheme = 'implicit-euler'
dtmin = 1e-13
###########################################
##
## dtmax is material dependent!
## for PTO is about 0.8 but BTO more like 3-10
##
###########################################
dtmax = 3.0
l_max_its = 200
[./TimeStepper]
type = IterationAdaptiveDT
optimal_iterations = 8
cutback_factor = 0.75
linear_iteration_ratio = 1000
dt = 0.3
[../]
verbose = true
[]
[Outputs]
###############################################
##
## Output options
##
###############################################
print_linear_residuals = false
perf_graph = false
[./out]
type = Exodus
file_base = out_bto_wall_T298K_2
elemental_as_nodal = true
[../]
[./console]
type = Console
time_precision = 3
[../]
[]
(examples/monodomain/BTO_monodomain_Tdef.i)
a1temp = a1def
a11temp = a11def
a111temp = a111def
[Mesh]
[gen]
############################################
##
## Type and dimension of the mesh
##
############################################
type = GeneratedMeshGenerator
dim = 3
#############################################
##
## Grid definition. Note that it should be
## nJ = 2*(Jmax-Jmin) for J = x, y, z
##
#############################################
nx = 8
ny = 8
nz = 8
#############################################
##
## Actual spatial coordinates of mesh.
## Jmax - Jmin = nJ/2 for J = x, y, z
## Units are in nanometers
##
#############################################
xmin = -2.0
xmax = 2.0
ymin = -2.0
ymax = 2.0
zmin = -2.0
zmax = 2.0
#############################################
##
## FE type/order (hexahedral, tetrahedral
##
#############################################
elem_type = HEX8
[]
[./cnode]
input = gen
############################################
##
## additional boundary sideset (one node)
## to zero one of the elastic displacement vectors
## vectors and eliminates rigid body translations
## from the degrees of freedom
##
## NOTE: This must conform with the about
## [Mesh] block settings
##
############################################
type = ExtraNodesetGenerator
coord = '0.0 0.0 0.0'
new_boundary = 100
[../]
[]
[GlobalParams]
len_scale = 1.0
polar_x = polar_x
polar_y = polar_y
polar_z = polar_z
potential_E_int = potential_E_int
vol = vol
displacements = 'u_x u_y u_z'
##############################################
##
## IMPORTANT(!): Units in Ferret are nm, kg,
## seconds, and attocoulombs
##
##############################################
u_x = u_x
u_y = u_y
u_z = u_z
[]
[Variables]
#################################
##
## Variable definitions
## P, u, phi, e^global_ij
## and their initial conditions
##
#################################
[./global_strain]
order = SIXTH
family = SCALAR
[../]
[./polar_x]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = 0.075
max = 0.1
[../]
[../]
[./polar_y]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = 0.075
max = 0.1
[../]
[../]
[./polar_z]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = 0.075
max = 0.1
[../]
[../]
[./potential_E_int]
order = FIRST
family = LAGRANGE
[../]
[./u_x]
order = FIRST
family = LAGRANGE
[../]
[./u_y]
order = FIRST
family = LAGRANGE
[../]
[./u_z]
order = FIRST
family = LAGRANGE
[../]
[]
[AuxVariables]
######################################
##
## Auxiarilly variable definitions
## (can be intermediate variables
## or for postprocessed quantities)
##
######################################
######################################
##
## Global displacements
##
######################################
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
######################################
##
## Stress/strain tensor components
##
######################################
[./e00]
order = CONSTANT
family = MONOMIAL
[../]
[./e01]
order = CONSTANT
family = MONOMIAL
[../]
[./e10]
order = CONSTANT
family = MONOMIAL
[../]
[./e11]
order = CONSTANT
family = MONOMIAL
[../]
[./e12]
order = CONSTANT
family = MONOMIAL
[../]
[./e22]
order = CONSTANT
family = MONOMIAL
[../]
[./s00]
order = CONSTANT
family = MONOMIAL
[../]
[./s01]
order = CONSTANT
family = MONOMIAL
[../]
[./s10]
order = CONSTANT
family = MONOMIAL
[../]
[./s11]
order = CONSTANT
family = MONOMIAL
[../]
[./s12]
order = CONSTANT
family = MONOMIAL
[../]
[./s22]
order = CONSTANT
family = MONOMIAL
[../]
[./eigs00]
order = CONSTANT
family = MONOMIAL
[../]
[./eigs11]
order = CONSTANT
family = MONOMIAL
[../]
[./eigs22]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
######################################
##
## Auxiarilly Kernel definitions
## (can be intermediate "operations"
## or for postprocessed quantities)
##
######################################
[./disp_x]
type = GlobalDisplacementAux
variable = disp_x
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 0
use_displaced_mesh = false
[../]
[./disp_y]
type = GlobalDisplacementAux
variable = disp_y
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 1
use_displaced_mesh = false
[../]
[./disp_z]
type = GlobalDisplacementAux
variable = disp_z
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 2
use_displaced_mesh = false
[../]
[./e00]
type = RankTwoAux
variable = e00
rank_two_tensor = total_strain
index_i = 0
index_j = 0
[../]
[./e01]
type = RankTwoAux
variable = e01
rank_two_tensor = total_strain
index_i = 0
index_j = 1
[../]
[./e10]
type = RankTwoAux
variable = e10
rank_two_tensor = total_strain
index_i = 1
index_j = 0
[../]
[./e12]
type = RankTwoAux
variable = e12
rank_two_tensor = total_strain
index_i = 1
index_j = 2
[../]
[./e11]
type = RankTwoAux
variable = e11
rank_two_tensor = total_strain
index_i = 1
index_j = 1
[../]
[./e22]
type = RankTwoAux
variable = e22
rank_two_tensor = total_strain
index_i = 2
index_j = 2
[../]
[./s00]
type = RankTwoAux
variable = s00
rank_two_tensor = stress
index_i = 0
index_j = 0
[../]
[./s01]
type = RankTwoAux
variable = s01
rank_two_tensor = stress
index_i = 0
index_j = 1
[../]
[./s10]
type = RankTwoAux
variable = s10
rank_two_tensor = stress
index_i = 1
index_j = 0
[../]
[./s12]
type = RankTwoAux
variable = s12
rank_two_tensor = stress
index_i = 1
index_j = 2
[../]
[./s11]
type = RankTwoAux
variable = s11
rank_two_tensor = stress
index_i = 1
index_j = 1
[../]
[./s22]
type = RankTwoAux
variable = s22
rank_two_tensor = stress
index_i = 2
index_j = 2
[../]
[./eigs00]
type = LocalABO3EigenstressAux
variable = eigs00
index_i = 0
index_j = 0
[../]
[./eigs11]
type = LocalABO3EigenstressAux
variable = eigs11
index_i = 1
index_j = 1
[../]
[./eigs22]
type = LocalABO3EigenstressAux
variable = eigs22
index_i = 2
index_j = 2
[../]
[]
[ScalarKernels]
######################################
##
## Necessary for PBC system
##
######################################
[./global_strain]
type = GlobalStrain
variable = global_strain
global_strain_uo = global_strain_uo
use_displaced_mesh = false
[../]
[]
[Materials]
#################################################
##
## Bulk free energy and electrostrictive
## coefficients gleaned from
## Marton and Hlinka
## Phys. Rev. B. 74, 104014, (2006)
##
## NOTE: there might be some Legendre transforms
## depending on what approach you use
## -i.e. inhomogeneous strain vs
## homogeneous strain [renormalized]
##
##################################################
[./Landau_P]
type = GenericConstantMaterial
prop_names = 'alpha1 alpha11 alpha12 alpha111 alpha112 alpha123 alpha1111 alpha1112 alpha1122 alpha1123'
prop_values = '${a1temp} ${a11temp} 0.323 ${a111temp} 4.47 4.91 0.0 0.0 0.0 0.0'
[../]
############################################
##
## Gradient coefficients from
## Marton and Hlinka
## Phys. Rev. B. 74, 104014, (2006)
##
############################################
[./Landau_G]
type = GenericConstantMaterial
prop_names = 'G110 G11_G110 G12_G110 G44_G110 G44P_G110'
prop_values = '0.5 0.51 -0.02 0.02 0.0'
[../]
[./mat_C]
type = GenericConstantMaterial
prop_names = 'C11 C12 C44'
prop_values = '275.0 179.0 54.3'
[../]
##################################################
##
## NOTE: Sign convention in Ferret for the
## electrostrictive coeff. is multiplied by
## an overall factor of (-1)
##
##################################################
[./mat_Q]
type = GenericConstantMaterial
prop_names = 'Q11 Q12 Q44'
prop_values = '-0.11 0.045 -0.029'
[../]
[./mat_q]
type = GenericConstantMaterial
prop_names = 'q11 q12 q44'
prop_values = '-14.2 0.74 -1.57'
[../]
[./eigen_strain]
type = ComputeEigenstrain
eigen_base = '1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0'
eigenstrain_name = eigenstrain
prefactor = 0.0
[../]
[./elasticity_tensor_1]
type = ComputeElasticityTensor
fill_method = symmetric9
###############################################
##
## symmetric9 fill_method is (default)
## C11 C12 C13 C22 C23 C33 C44 C55 C66
##
###############################################
C_ijkl = '275.0 179.0 179.0 275.0 179.0 275.0 54.3 54.3 54.3'
[../]
[./strain_1]
type = ComputeSmallStrain
global_strain = global_strain
eigenstrain_names = eigenstrain
[../]
[./stress_1]
type = ComputeLinearElasticStress
[../]
[./global_strain]
type = ComputeGlobalStrain
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
[../]
[./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)
##
###############################################
type = GenericConstantMaterial
prop_names = 'permittivity'
prop_values = '0.08854187'
[../]
[]
[Kernels]
###############################################
##
## Physical Kernel operators
## to enforce TDLGD evolution
##
###############################################
#Elastic problem
[./TensorMechanics]
use_displaced_mesh = false
eigenstrains_name = eigenstrain
[../]
[./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
[../]
[./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
[../]
[./electrostr_ux]
type = ElectrostrictiveCouplingDispDerivative
variable = u_x
component = 0
[../]
[./electrostr_uy]
type = ElectrostrictiveCouplingDispDerivative
variable = u_y
component = 1
[../]
[./electrostr_uz]
type = ElectrostrictiveCouplingDispDerivative
variable = u_z
component = 2
[../]
[./electrostr_polar_coupled_x]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_x
component = 0
[../]
[./electrostr_polar_coupled_y]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_y
component = 1
[../]
[./electrostr_polar_coupled_z]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_z
component = 2
[../]
[./polar_x_electric_E]
type = PolarElectricEStrong
variable = potential_E_int
[../]
[./FE_E_int]
type = Electrostatics
variable = potential_E_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
[../]
[./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
[../]
[]
[BCs]
[./Periodic]
[./xyz]
auto_direction = 'x y z'
variable = 'u_x u_y u_z polar_x polar_y polar_z potential_E_int'
[../]
[../]
# fix center point location
[./centerfix_x]
type = DirichletBC
boundary = 100
variable = u_x
value = 0
[../]
[./centerfix_y]
type = DirichletBC
boundary = 100
variable = u_y
value = 0
[../]
[./centerfix_z]
type = DirichletBC
boundary = 100
variable = u_z
value = 0
[../]
[]
[Postprocessors]
###############################################
##
## Postprocessors (integrations over the
## computational domain) to calculate the total energy
## decomposed into linear combinations of the
## different physics.
##
###############################################
[./avePx]
type = ElementAverageValue
variable = polar_x
execute_on = 'timestep_end'
[../]
[./avePy]
type = ElementAverageValue
variable = polar_y
execute_on = 'timestep_end'
[../]
[./avePz]
type = ElementAverageValue
variable = polar_z
execute_on = 'timestep_end'
[../]
[./ave_e00]
type = ElementAverageValue
variable = e00
execute_on = 'timestep_end'
[../]
[./ave_e11]
type = ElementAverageValue
variable = e11
execute_on = 'timestep_end'
[../]
[./ave_e22]
type = ElementAverageValue
variable = e22
execute_on = 'timestep_end'
[../]
[./Fb]
type = BulkEnergyEighth
execute_on = 'timestep_end'
[../]
[./Fw]
type = WallEnergy
execute_on = 'timestep_end'
[../]
[./Fela]
type = ElasticEnergy
execute_on = 'timestep_end'
use_displaced_mesh = false
[../]
[./Fc]
type = ElectrostrictiveCouplingEnergy
execute_on = 'timestep_end'
[../]
[./Fele]
type = ElectrostaticEnergy
execute_on = 'timestep_end'
[../]
[./Ftot]
type = LinearCombinationPostprocessor
pp_names = 'Fb Fw Fc Fele'
pp_coefs = ' 1 1 1 1'
execute_on = 'timestep_end'
[../]
[./vol]
type = VolumePostprocessor
execute_on = 'timestep_end'
[../]
[./px]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 0
[../]
[./py]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 1
[../]
[./pz]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 2
[../]
[./perc_change]
type = PercentChangePostprocessor
postprocessor = Ftot
execute_on = 'timestep_end'
[../]
[]
[UserObjects]
###############################################
##
## GlobalStrain system to enforce periodicity
## in the anisotropic strain field
##
###############################################
[./global_strain_uo]
type = GlobalATiO3MaterialRVEUserObject
use_displaced_mesh = false
execute_on = 'Initial Linear Nonlinear'
[../]
###############################################
##
## terminator to end energy evolution when the energy difference
## between subsequent time steps is lower than 5e-6
##
## NOTE: can fail if the time step is smallhotkey for tilde
##
###############################################
[./kill]
type = Terminator
expression = 'perc_change <= 1.0e-6'
[../]
[]
[Preconditioning]
###############################################
##
## Numerical preconditioning/solver options
##
###############################################
[./smp]
type = SMP
full = true
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-ksp_gmres_restart -snes_atol -snes_rtol -ksp_rtol -pc_type -build_twosided'
petsc_options_value = ' 160 1e-8 1e-6 1e-6 bjacobi allreduce'
[../]
[]
[Executioner]
##########################################
##
## Time integration/solver options
##
##########################################
type = Transient
solve_type = 'PJFNK'
scheme = 'implicit-euler'
dtmin = 1e-13
###########################################
##
## dtmax is material dependent!
##
###########################################
dtmax = 3.0
l_max_its = 200
[./TimeStepper]
type = IterationAdaptiveDT
optimal_iterations = 8
cutback_factor = 0.75
linear_iteration_ratio = 1000
dt = 0.3
[../]
verbose = true
[]
[Outputs]
###############################################
##
## Output options
##
###############################################
print_linear_residuals = false
perf_graph = false
[./outCSV]
type = CSV
file_base = out_bto_monodomain_Tdef
[../]
[]
(tutorial/film.i)
[Mesh]
[gen]
############################################
##
## Type and dimension of the mesh
##
############################################
type = GeneratedMeshGenerator
dim = 3
nx = 32
ny = 32
nz = 30
xmin = -16.0
xmax = 16.0
ymin = -16.0
ymax = 16.0
zmin = -10.0
zmax = 20.0
#############################################
##
## FE type/order (hexahedral, tetrahedral
##
#############################################
elem_type = HEX8
[]
[./cnode]
input = gen
############################################
##
## additional boundary sideset (one node)
## to zero one of the elastic displacement vectors
## vectors and eliminates rigid body translations
## from the degrees of freedom
##
## NOTE: This must conform with the about
## [Mesh] block settings
##
############################################
type = ExtraNodesetGenerator
coord = '-16.0 -16.0 -10.0'
new_boundary = 100
[../]
[subdomains]
type = SubdomainBoundingBoxGenerator
input = cnode
bottom_left = '-16.0 -16.0 -10.0'
block_id = 1
top_right = '16.0 16.0 0'
location = INSIDE
[]
[film_interface]
type = SideSetsBetweenSubdomainsGenerator
input = subdomains
primary_block = 0
paired_block = 1
new_boundary = 52
[]
[film_surface]
type = SideSetsFromNormalsGenerator
input = film_interface
normals = '0 0 1'
fixed_normal = true
new_boundary = '107'
[]
[substrate_bottom]
type = SideSetsFromNormalsGenerator
input = film_surface
normals = '0 0 -1'
fixed_normal = true
new_boundary = '108'
[]
[]
[GlobalParams]
len_scale = 1.0
polar_x = polar_x
polar_y = polar_y
polar_z = polar_z
potential_E_int = potential_E_int
displacements = 'u_x u_y u_z'
##############################################
##=
## IMPORTANT(!): Units in Ferret are nm, kg,
## seconds, and attocoulombs
##
##############################################
vol = vol
u_x = u_x
u_y = u_y
u_z = u_z
[]
[Variables]
#################################
##
## Variable definitions
## P, u, phi, e^global_ij
## and their initial conditions
##
#################################
[./global_strain]
order = SIXTH
family = SCALAR
[../]
[./polar_x]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -1e-2
max = 1e-2
[../]
block = '0'
[../]
[./polar_y]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -1e-2
max = 1e-2
[../]
block = '0'
[../]
[./polar_z]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -1e-2
max = 1e-2
[../]
block = '0'
[../]
[./potential_E_int]
order = FIRST
family = LAGRANGE
block = '0 1'
[../]
[./u_x]
order = FIRST
family = LAGRANGE
block = '0 1'
[../]
[./u_y]
order = FIRST
family = LAGRANGE
block = '0 1'
[../]
[./u_z]
order = FIRST
family = LAGRANGE
block = '0 1'
[../]
[]
[Kernels]
###############################################
##
## Physical Kernel operators
## to enforce TDLGD evolution
##
###############################################
#Elastic problem
[./TensorMechanics]
use_displaced_mesh = false
eigenstrain_names = eigenstrain
[../]
[./bed_x]
type = BulkEnergyDerivativeEighth
variable = polar_x
component = 0
block = '0'
[../]
[./bed_y]
type = BulkEnergyDerivativeEighth
variable = polar_y
component = 1
block = '0'
[../]
[./bed_z]
type = BulkEnergyDerivativeEighth
variable = polar_z
component = 2
block = '0'
[../]
[./walled_x]
type = WallEnergyDerivative
variable = polar_x
component = 0
block = '0'
[../]
[./walled_y]
type = WallEnergyDerivative
variable = polar_y
component = 1
block = '0'
[../]
[./walled_z]
type = WallEnergyDerivative
variable = polar_z
component = 2
block = '0'
[../]
[./walled2_x]
type = Wall2EnergyDerivative
variable = polar_x
component = 0
block = '0'
[../]
[./walled2_y]
type = Wall2EnergyDerivative
variable = polar_y
component = 1
block = '0'
[../]
[./walled2_z]
type = Wall2EnergyDerivative
variable = polar_z
component = 2
block = '0'
[../]
[./electrostr_ux]
type = ElectrostrictiveCouplingDispDerivative
variable = u_x
component = 0
block = '0'
[../]
[./electrostr_uy]
type = ElectrostrictiveCouplingDispDerivative
variable = u_y
component = 1
block = '0'
[../]
[./electrostr_uz]
type = ElectrostrictiveCouplingDispDerivative
variable = u_z
component = 2
block = '0'
[../]
[./electrostr_polar_coupled_x]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_x
component = 0
block = '0'
[../]
[./electrostr_polar_coupled_y]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_y
component = 1
block = '0'
[../]
[./electrostr_polar_coupled_z]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_z
component = 2
block = '0'
[../]
[./polar_x_electric_E]
type = PolarElectricEStrong
variable = potential_E_int
block = '0'
[../]
[./FE_E_int]
type = Electrostatics
variable = potential_E_int
block = '0 1'
[../]
[./polar_electric_px]
type = PolarElectricPStrong
variable = polar_x
component = 0
block = '0'
[../]
[./polar_electric_py]
type = PolarElectricPStrong
variable = polar_y
component = 1
block = '0'
[../]
[./polar_electric_pz]
type = PolarElectricPStrong
variable = polar_z
component = 2
block = '0'
[../]
[./polar_x_time]
type = TimeDerivativeScaled
variable=polar_x
time_scale = 1.0
block = '0'
[../]
[./polar_y_time]
type = TimeDerivativeScaled
variable = polar_y
time_scale = 1.0
block = '0'
[../]
[./polar_z_time]
type = TimeDerivativeScaled
variable = polar_z
time_scale = 1.0
block = '0'
[../]
[./u_x_time]
type = TimeDerivativeScaled
variable = u_x
time_scale = 1.0
[../]
[./u_y_time]
type = TimeDerivativeScaled
variable = u_y
time_scale = 1.0
[../]
[./u_z_time]
type = TimeDerivativeScaled
variable = u_z
time_scale = 1.0
[../]
[]
[AuxVariables]
######################################
##
## Auxiarilly variable definitions
## (can be intermediate variables
## or for postprocessed quantities)
##
######################################
######################################
##
## Global displacements
##
######################################
[./disp_x]
block = '0 1'
[../]
[./disp_y]
block = '0 1'
[../]
[./disp_z]
block = '0 1'
[../]
######################################
##
## Stress/strain tensor components
##
######################################
[./e00]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e01]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e10]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e11]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e12]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e22]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s00]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s01]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s10]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s11]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s12]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s22]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./E_x]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./E_y]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./E_z]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[]
[AuxKernels]
######################################
##
## Auxiarilly Kernel definitions
## (can be intermediate "operations"
## or for postprocessed quantities)
##
######################################
[./disp_x]
type = GlobalDisplacementAux
variable = disp_x
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 0
use_displaced_mesh = false
[../]
[./disp_y]
type = GlobalDisplacementAux
variable = disp_y
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 1
use_displaced_mesh = false
[../]
[./disp_z]
type = GlobalDisplacementAux
variable = disp_z
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 2
use_displaced_mesh = false
[../]
[./e00]
type = RankTwoAux
variable = e00
rank_two_tensor = total_strain
index_i = 0
index_j = 0
[../]
[./e01]
type = RankTwoAux
variable = e01
rank_two_tensor = total_strain
index_i = 0
index_j = 1
[../]
[./e10]
type = RankTwoAux
variable = e10
rank_two_tensor = total_strain
index_i = 1
index_j = 0
[../]
[./e12]
type = RankTwoAux
variable = e12
rank_two_tensor = total_strain
index_i = 1
index_j = 2
[../]
[./e11]
type = RankTwoAux
variable = e11
rank_two_tensor = total_strain
index_i = 1
index_j = 1
[../]
[./e22]
type = RankTwoAux
variable = e22
rank_two_tensor = total_strain
index_i = 2
index_j = 2
[../]
[./s00]
type = RankTwoAux
variable = s00
rank_two_tensor = stress
index_i = 0
index_j = 0
[../]
[./s01]
type = RankTwoAux
variable = s01
rank_two_tensor = stress
index_i = 0
index_j = 1
[../]
[./s10]
type = RankTwoAux
variable = s10
rank_two_tensor = stress
index_i = 1
index_j = 0
[../]
[./s12]
type = RankTwoAux
variable = s12
rank_two_tensor = stress
index_i = 1
index_j = 2
[../]
[./s11]
type = RankTwoAux
variable = s11
rank_two_tensor = stress
index_i = 1
index_j = 1
[../]
[./s22]
type = RankTwoAux
variable = s22
rank_two_tensor = stress
index_i = 2
index_j = 2
[../]
[]
[ScalarKernels]
######################################
##
## Necessary for PBC system
##
######################################
[./global_strain]
type = GlobalStrain
variable = global_strain
global_strain_uo = global_strain_uo
use_displaced_mesh = false
[../]
[]
[Materials]
#################################################
##
## Landau coefficients from Li et al (2001)
##
##################################################
[./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_P_substr]
type = GenericConstantMaterial
prop_names = 'alpha1 alpha11 alpha12 alpha111 alpha112 alpha123 alpha1111 alpha1112 alpha1122 alpha1123'
prop_values = '10.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0'
block = '1'
[../]
[./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'
[../]
[./mat_C_FE]
type = GenericConstantMaterial
prop_names = 'C11 C12 C44'
prop_values = '175.0 79.4 111.1'
block = '0'
[../]
[./mat_C_sub]
type = GenericConstantMaterial
prop_names = 'C11 C12 C44'
prop_values = '220.0 34.4 161.1'
block = '1'
[../]
##################################################
##=
## NOTE: Sign convention in Ferret for the
## electrostrictive coeff. is multiplied by
## an overall factor of (-1)
##
##################################################
[./mat_Q]
type = GenericConstantMaterial
prop_names = 'Q11 Q12 Q44'
prop_values = '-0.089 0.026 -0.03375'
block = '0 1'
[../]
[./mat_q]
type = GenericConstantMaterial
prop_names = 'q11 q12 q44'
prop_values = '-11.4 0.01438 -7.5'
[../]
[./eigen_strain]
type = ComputeEigenstrain
# eigen_base = 'exx exy exz eyx eyy eyz ezx ezy ezz'
eigen_base = '1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0'
eigenstrain_name = eigenstrain
prefactor = 0.0
[../]
[./elasticity_tensor_1]
type = ComputeElasticityTensor
fill_method = symmetric9
###############################################
##
## symmetric9 fill_method is (default)
## C11 C12 C13 C22 C23 C33 C44 C55 C66
##
###############################################
C_ijkl = '175.0 79.4 79.4 175.0 79.4 175.0 111.1 111.1 111.1'
[../]
[./strain_1]
type = ComputeSmallStrain
global_strain = global_strain
eigenstrain_names = eigenstrain
[../]
[./stress_1]
type = ComputeLinearElasticStress
[../]
[./global_strain]
type = ComputeGlobalStrain
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
[../]
[./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)
##
###############################################
type = GenericConstantMaterial
prop_names = 'permittivity'
prop_values = '0.08854187'
[../]
[]
[BCs]
[./Periodic]
[./xy]
auto_direction = 'x y'
variable = 'u_x u_y u_z polar_x polar_y polar_z potential_E_int'
[../]
[../]
[./boundary_interface_grounding]
type = DirichletBC
boundary = '52'
variable = potential_E_int
value = 0.0
[../]
# fix center point location
[./centerfix_x]
type = DirichletBC
boundary = '108'
variable = u_x
value = 0
[../]
[./centerfix_y]
type = DirichletBC
boundary = '108'
variable = u_y
value = 0
[../]
[./centerfix_z]
type = DirichletBC
boundary = '108'
variable = u_z
value = 0
[../]
[]
[Postprocessors]
###############################################
##=
## Postprocessors (integrations over the
## computational domain) to calculate the total
## energy decomposed into linear combinations of
## the different physics.
##
###############################################
[./Fbulk]
type = BulkEnergyEighth
execute_on = 'timestep_end'
block = '0'
[../]
[./Fwall]
type = WallEnergy
execute_on = 'timestep_end'
block = '0'
[../]
[./Felastic]
type = ElasticEnergy
execute_on = 'timestep_end'
use_displaced_mesh = false
block = '0'
[../]
[./Fcoupled]
type = ElectrostrictiveCouplingEnergy
execute_on = 'timestep_end'
block = '0'
[../]
[./Felec]
type = ElectrostaticEnergy
execute_on = 'timestep_end'
block = '0'
[../]
[./Ftotal]
type = LinearCombinationPostprocessor
pp_names = 'Fbulk Fwall Fcoupled Felec'
pp_coefs = '0.160218 0.160218 0.160218 0.160218' #converted to eV
execute_on = 'timestep_end'
[../]
[./vol]
type = VolumePostprocessor
execute_on = 'timestep_end'
[../]
[./px]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 0
block = '0'
[../]
[./py]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 1
block = '0'
[../]
[./pz]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 2
block = '0'
[../]
[./perc_change]
type = PercentChangePostprocessor
postprocessor = Ftotal
execute_on = 'timestep_end'
[../]
[./elapsed]
type = PerfGraphData
section_name = "Root" # for profiling the problem [on]
data_type = total
[../]
[]
[UserObjects]
###############################################
##
## GlobalStrain system to enforce periodicity
## in the anisotropic strain field
##
###############################################
[./global_strain_uo]
type = GlobalATiO3MaterialRVEUserObject
use_displaced_mesh = false
execute_on = 'Initial Linear Nonlinear'
applied_stress_tensor = '2.1 2.1 1.9056 0.0 0.0 0.0'
block = '0'
[../]
###############################################
##
## terminator to end energy evolution when the energy difference
## between subsequent time steps is lower than 5e-6
##
## NOTE: can fail if the time step is small
##
###############################################
[./kill]
type = Terminator
expression = 'perc_change <= 5.0e-4'
[../]
[]
[Preconditioning]
###############################################
##
## Numerical preconditioning/solver options
##
###############################################
[./smp]
type = SMP
full = true
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-ksp_gmres_restart -snes_atol -snes_rtol -ksp_rtol -pc_type -build_twosided'
petsc_options_value = ' 80 1e-8 1e-6 1e-5 bjacobi allreduce'
[../]
[]
[Executioner]
##########################################
##
## Time integ=ration/solver options
##
##########################################
type = Transient
solve_type = 'PJFNK'
scheme = 'bdf2'
dtmin = 1e-13
dtmax = 0.6
l_max_its = 200
[./TimeStepper]
type = IterationAdaptiveDT
optimal_iterations = 6
growth_factor = 1.2
cutback_factor = 0.75
linear_iteration_ratio = 1000
dt = 0.6
[../]
verbose = true
nl_max_its = 20
[]
[Outputs]
###############################################
##==
## Output options
##
###############################################
print_linear_residuals = false
perf_graph = false
[./out]
type = Exodus
file_base = out_PTOfilm_e12_T298K_E0_E0
elemental_as_nodal = true
interval = 1
[../]
[]
(examples/films/PZTfilm_e12_T298K_E0.i)
[Mesh]
[gen]
############################################
##
## Type and dimension of the mesh
##
############################################
type = GeneratedMeshGenerator
dim = 3
#############################################
##
## Grid definition. Note that it should be
## nJ = 2*(Jmax-Jmin) for J = x, y, z
##
#############################################
nx = 32
ny = 32
nz = 30
#############################################
##
## Actual spatial coordinates of mesh.
## Jmax - Jmin = nJ/2 for J = x, y, z
## Units are in nanometers
##
#############################################
xmin = -16.0
xmax = 16.0
ymin = -16.0
ymax = 16.0
zmin = -10.0
zmax = 20.0
#############################################
##
## FE type/order (hexahedral, tetrahedral
##
#############################################
elem_type = HEX8
[]
[./cnode]
input = gen
############################################
##
## additional boundary sideset (one node)
## to zero one of the elastic displacement vectors
## vectors and eliminates rigid body translations
## from the degrees of freedom
##
## NOTE: This must conform with the about
## [Mesh] block settings
##
############################################
type = ExtraNodesetGenerator
coord = '-16.0 -16.0 -10.0'
new_boundary = 100
[../]
[subdomains]
type = SubdomainBoundingBoxGenerator
input = cnode
bottom_left = '-16.0 -16.0 -10.0'
block_id = 1
top_right = '16.0 16.0 0'
location = INSIDE
[]
[film_interface]
type = SideSetsBetweenSubdomainsGenerator
input = subdomains
primary_block = 0
paired_block = 1
new_boundary = 52
[]
[film_surface]
type = SideSetsFromNormalsGenerator
input = film_interface
normals = '0 0 1'
fixed_normal = true
new_boundary = '107'
[]
[substrate_bottom]
type = SideSetsFromNormalsGenerator
input = film_surface
normals = '0 0 -1'
fixed_normal = true
new_boundary = '108'
[]
[]
[GlobalParams]
len_scale = 1.0
polar_x = polar_x
polar_y = polar_y
polar_z = polar_z
potential_E_int = potential_E_int
displacements = 'u_x u_y u_z'
##############################################
##=
## IMPORTANT(!): Units in Ferret are nm, kg,
## seconds, and attocoulombs
##
##############################################
vol = vol
u_x = u_x
u_y = u_y
u_z = u_z
[]
[Variables]
#################################
##
## Variable definitions
## P, u, phi, e^global_ij
## and their initial conditions
##
#################################
[./global_strain]
order = SIXTH
family = SCALAR
[../]
[./polar_x]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -1e-2
max = 1e-2
[../]
block = '0'
[../]
[./polar_y]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -1e-2
max = 1e-2
[../]
block = '0'
[../]
[./polar_z]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -1e-2
max = 1e-2
[../]
block = '0'
[../]
[./potential_E_int]
order = FIRST
family = LAGRANGE
block = '0 1'
[../]
[./u_x]
order = FIRST
family = LAGRANGE
block = '0 1'
[../]
[./u_y]
order = FIRST
family = LAGRANGE
block = '0 1'
[../]
[./u_z]
order = FIRST
family = LAGRANGE
block = '0 1'
[../]
[]
[AuxVariables]
######################################
##
## Auxiarilly variable definitions
## (can be intermediate variables
## or for postprocessed quantities)
##
######################################
######################################
##
## Global displacements
##
######################################
[./disp_x]
block = '0 1'
[../]
[./disp_y]
block = '0 1'
[../]
[./disp_z]
block = '0 1'
[../]
######################################
##
## Stress/strain tensor components
##
######################################
[./e00]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e01]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e10]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e11]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e12]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e22]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s00]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s01]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s10]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s11]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s12]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s22]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./E_x]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./E_y]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./E_z]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[]
[AuxKernels]
######################################
##
## Auxiarilly Kernel definitions
## (can be intermediate "operations"
## or for postprocessed quantities)
##
######################################
[./disp_x]
type = GlobalDisplacementAux
variable = disp_x
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 0
use_displaced_mesh = false
[../]
[./disp_y]
type = GlobalDisplacementAux
variable = disp_y
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 1
use_displaced_mesh = false
[../]
[./disp_z]
type = GlobalDisplacementAux
variable = disp_z
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 2
use_displaced_mesh = false
[../]
[./e00]
type = RankTwoAux
variable = e00
rank_two_tensor = total_strain
index_i = 0
index_j = 0
[../]
[./e01]
type = RankTwoAux
variable = e01
rank_two_tensor = total_strain
index_i = 0
index_j = 1
[../]
[./e10]
type = RankTwoAux
variable = e10
rank_two_tensor = total_strain
index_i = 1
index_j = 0
[../]
[./e12]
type = RankTwoAux
variable = e12
rank_two_tensor = total_strain
index_i = 1
index_j = 2
[../]
[./e11]
type = RankTwoAux
variable = e11
rank_two_tensor = total_strain
index_i = 1
index_j = 1
[../]
[./e22]
type = RankTwoAux
variable = e22
rank_two_tensor = total_strain
index_i = 2
index_j = 2
[../]
[./s00]
type = RankTwoAux
variable = s00
rank_two_tensor = stress
index_i = 0
index_j = 0
[../]
[./s01]
type = RankTwoAux
variable = s01
rank_two_tensor = stress
index_i = 0
index_j = 1
[../]
[./s10]
type = RankTwoAux
variable = s10
rank_two_tensor = stress
index_i = 1
index_j = 0
[../]
[./s12]
type = RankTwoAux
variable = s12
rank_two_tensor = stress
index_i = 1
index_j = 2
[../]
[./s11]
type = RankTwoAux
variable = s11
rank_two_tensor = stress
index_i = 1
index_j = 1
[../]
[./s22]
type = RankTwoAux
variable = s22
rank_two_tensor = stress
index_i = 2
index_j = 2
[../]
[./ex]
type = ElecFieldAux
variable = E_x
component = 0
[../]
[./ey]
type = ElecFieldAux
variable = E_y
component = 1
[../]
[./ez]
type = ElecFieldAux
variable = E_z
component = 2
[../]
[]
[ScalarKernels]
######################################
##
## Necessary for PBC system
##
######################################
[./global_strain]
type = GlobalStrain
variable = global_strain
global_strain_uo = global_strain_uo
use_displaced_mesh = false
[../]
[]
[Materials]
#################################################
##
## Landau, electrostrictive, elastic coefficients
## M. Mtebwa, A. K. Tagantsev, and N. Setter
## AIP Adv. 4, 127150 (2014).
##
##################################################
[./Landau_P_FE]
type = GenericConstantMaterial
prop_names = 'alpha1 alpha11 alpha12 alpha111 alpha112 alpha123 alpha1111 alpha1112 alpha1122 alpha1123'
prop_values = '-0.048889 0.04764 0.1336 0.1735 0.6128 -2.894 0.0 0.0 0.0 0.0'
block = '0'
[../]
[./Landau_P_substr]
type = GenericConstantMaterial
prop_names = 'alpha1 alpha11 alpha12 alpha111 alpha112 alpha123 alpha1111 alpha1112 alpha1122 alpha1123'
prop_values = '10.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0'
block = '1'
[../]
[./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'
[../]
[./mat_C_FE]
type = GenericConstantMaterial
prop_names = 'C11 C12 C44'
prop_values = '179.073 66.71 82.6446'
block = '0'
[../]
[./mat_C_sub]
type = GenericConstantMaterial
prop_names = 'C11 C12 C44'
prop_values = '220.0 34.4 161.1'
block = '1'
[../]
##################################################
##=
## NOTE: Sign convention in Ferret for the
## electrostrictive coeff. is multiplied by
## an overall factor of (-1)
##
##################################################
[./mat_Q]
type = GenericConstantMaterial
prop_names = 'Q11 Q12 Q44'
prop_values = '-0.0966 0.046 -0.0819'
block = '0 1'
[../]
[./mat_q]
type = GenericConstantMaterial
prop_names = 'q11 q12 q44'
prop_values = '-11.1608 4.86165 -6.76859'
[../]
[./eigen_strain]
type = ComputeEigenstrain
# eigen_base = 'exx exy exz eyx eyy eyz ezx ezy ezz'
eigen_base = '1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0'
eigenstrain_name = eigenstrain
prefactor = 0.0
[../]
[./elasticity_tensor_1]
type = ComputeElasticityTensor
fill_method = symmetric9
###############################################
##
## symmetric9 fill_method is (default)
## C11 C12 C13 C22 C23 C33 C44 C55 C66
##
###############################################
C_ijkl = '179.073 66.71 66.71 179.073 66.71 179.073 82.6446 82.6446 82.6446'
[../]
[./strain_1]
type = ComputeSmallStrain
global_strain = global_strain
eigenstrain_names = eigenstrain
[../]
[./stress_1]
type = ComputeLinearElasticStress
[../]
[./global_strain]
type = ComputeGlobalStrain
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
[../]
[./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)
##
###############################################
type = GenericConstantMaterial
prop_names = 'permittivity'
prop_values = '0.08854187'
[../]
[]
[Kernels]
###############################################
##
## Physical Kernel operators
## to enforce TDLGD evolution
##
###############################################
#Elastic problem
[./TensorMechanics]
use_displaced_mesh = false
eigenstrain_names = eigenstrain
[../]
[./bed_x]
type = BulkEnergyDerivativeEighth
variable = polar_x
component = 0
block = '0'
[../]
[./bed_y]
type = BulkEnergyDerivativeEighth
variable = polar_y
component = 1
block = '0'
[../]
[./bed_z]
type = BulkEnergyDerivativeEighth
variable = polar_z
component = 2
block = '0'
[../]
[./walled_x]
type = WallEnergyDerivative
variable = polar_x
component = 0
block = '0'
[../]
[./walled_y]
type = WallEnergyDerivative
variable = polar_y
component = 1
block = '0'
[../]
[./walled_z]
type = WallEnergyDerivative
variable = polar_z
component = 2
block = '0'
[../]
[./electrostr_ux]
type = ElectrostrictiveCouplingDispDerivative
variable = u_x
component = 0
block = '0'
[../]
[./electrostr_uy]
type = ElectrostrictiveCouplingDispDerivative
variable = u_y
component = 1
block = '0'
[../]
[./electrostr_uz]
type = ElectrostrictiveCouplingDispDerivative
variable = u_z
component = 2
block = '0'
[../]
[./electrostr_polar_coupled_x]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_x
component = 0
block = '0'
[../]
[./electrostr_polar_coupled_y]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_y
component = 1
block = '0'
[../]
[./electrostr_polar_coupled_z]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_z
component = 2
block = '0'
[../]
[./polar_x_electric_E]
type = PolarElectricEStrong
variable = potential_E_int
block = '0'
[../]
[./FE_E_int]
type = Electrostatics
variable = potential_E_int
block = '0 1'
[../]
[./polar_electric_px]
type = PolarElectricPStrong
variable = polar_x
component = 0
block = '0'
[../]
[./polar_electric_py]
type = PolarElectricPStrong
variable = polar_y
component = 1
block = '0'
[../]
[./polar_electric_pz]
type = PolarElectricPStrong
variable = polar_z
component = 2
block = '0'
[../]
[./polar_x_time]
type = TimeDerivativeScaled
variable=polar_x
time_scale = 1.0
block = '0'
[../]
[./polar_y_time]
type = TimeDerivativeScaled
variable = polar_y
time_scale = 1.0
block = '0'
[../]
[./polar_z_time]
type = TimeDerivativeScaled
variable = polar_z
time_scale = 1.0
block = '0'
[../]
[./u_x_time]
type = TimeDerivativeScaled
variable = u_x
time_scale = 1.0
[../]
[./u_y_time]
type = TimeDerivativeScaled
variable = u_y
time_scale = 1.0
[../]
[./u_z_time]
type = TimeDerivativeScaled
variable = u_z
time_scale = 1.0
[../]
[]
[BCs]
[./Periodic]
[./xy]
auto_direction = 'x y'
variable = 'u_x u_y u_z polar_x polar_y polar_z potential_E_int'
[../]
[../]
[./boundary_interface_grounding]
type = DirichletBC
boundary = '52'
variable = potential_E_int
value = 0.0
[../]
# fix center point location
[./centerfix_x]
type = DirichletBC
boundary = '108'
variable = u_x
value = 0
[../]
[./centerfix_y]
type = DirichletBC
boundary = '108'
variable = u_y
value = 0
[../]
[./centerfix_z]
type = DirichletBC
boundary = '108'
variable = u_z
value = 0
[../]
[]
[Postprocessors]
###############################################
##=
## Postprocessors (integrations over the
## computational domain) to calculate the total
## energy decomposed into linear combinations of
## the different physics.
##
###############################################
[./Fbulk]
type = BulkEnergyEighth
execute_on = 'timestep_end'
block = '0'
[../]
[./Fwall]
type = WallEnergy
execute_on = 'timestep_end'
block = '0'
[../]
[./Felastic]
type = ElasticEnergy
execute_on = 'timestep_end'
use_displaced_mesh = false
block = '0'
[../]
[./Fcoupled]
type = ElectrostrictiveCouplingEnergy
execute_on = 'timestep_end'
block = '0'
[../]
[./Felec]
type = ElectrostaticEnergy
execute_on = 'timestep_end'
block = '0'
[../]
[./Ftotal]
type = LinearCombinationPostprocessor
pp_names = 'Fbulk Fwall Fcoupled Felec'
pp_coefs = '0.160218 0.160218 0.160218 0.160218' #converted to eV
execute_on = 'timestep_end'
[../]
[./vol]
type = VolumePostprocessor
execute_on = 'timestep_end'
[../]
[./px]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 0
block = '0'
[../]
[./py]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 1
block = '0'
[../]
[./pz]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 2
block = '0'
[../]
[./perc_change]
type = PercentChangePostprocessor
postprocessor = Ftotal
execute_on = 'timestep_end'
[../]
[./elapsed]
type = PerfGraphData
section_name = "Root" # for profiling the problem [on]
data_type = total
[../]
[]
[UserObjects]
###############################################
##
## GlobalStrain system to enforce periodicity
## in the anisotropic strain field
##
###############################################
[./global_strain_uo]
type = GlobalATiO3MaterialRVEUserObject
use_displaced_mesh = false
execute_on = 'Initial Linear Nonlinear'
applied_stress_tensor = '2.1 2.1 1.9056 0.0 0.0 0.0'
block = '0'
[../]
###############################################
##
## terminator to end energy evolution when the energy difference
## between subsequent time steps is lower than 5e-6
##
## NOTE: can fail if the time step is small
##
###############################################
[./kill]
type = Terminator
expression = 'perc_change <= 5.0e-4'
[../]
[]
[Preconditioning]
###############################################
##
## Numerical preconditioning/solver options
##
###############################################
[./smp]
type = SMP
full = true
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-ksp_gmres_restart -snes_atol -snes_rtol -ksp_rtol -pc_type -build_twosided'
petsc_options_value = ' 80 1e-8 1e-6 1e-5 bjacobi allreduce'
[../]
[]
[Executioner]
##########################################
##
## Time integ=ration/solver options
##
##########################################
type = Transient
solve_type = 'PJFNK'
scheme = 'bdf2'
dtmin = 1e-13
dtmax = 0.6
l_max_its = 200
[./TimeStepper]
type = IterationAdaptiveDT
optimal_iterations = 6
growth_factor = 1.2
cutback_factor = 0.75
linear_iteration_ratio = 1000
dt = 0.6
[../]
verbose = true
nl_max_its = 20
[]
[Outputs]
###############################################
##==
## Output options
##
###############################################
print_linear_residuals = false
perf_graph = false
[./out]
type = Exodus
file_base = out_PTOfilm_e12_T298K_E0_E0
elemental_as_nodal = true
interval = 1
[../]
[]
(examples/monodomain/PZT_monodomain_Tdef.i)
a1temp = a1def
[Mesh]
[gen]
############################################
##
## Type and dimension of the mesh
##
############################################
type = GeneratedMeshGenerator
dim = 3
#############################################
##
## Grid definition. Note that it should be
## nJ = 2*(Jmax-Jmin) for J = x, y, z
##
#############################################
nx = 8
ny = 8
nz = 8
#############################################
##
## Actual spatial coordinates of mesh.
## Jmax - Jmin = nJ/2 for J = x, y, z
## Units are in nanometers
##
#############################################
xmin = -2.0
xmax = 2.0
ymin = -2.0
ymax = 2.0
zmin = -2.0
zmax = 2.0
#############################################
##
## FE type/order (hexahedral, tetrahedral
##
#############################################
elem_type = HEX8
[]
[./cnode]
input = gen
############################################
##
## additional boundary sideset (one node)
## to zero one of the elastic displacement vectors
## vectors and eliminates rigid body translations
## from the degrees of freedom
##
## NOTE: This must conform with the about
## [Mesh] block settings
##
############################################
type = ExtraNodesetGenerator
coord = '0.0 0.0 0.0'
new_boundary = 100
[../]
[]
[GlobalParams]
len_scale = 1.0
polar_x = polar_x
polar_y = polar_y
polar_z = polar_z
potential_E_int = potential_E_int
vol = vol
displacements = 'u_x u_y u_z'
##############################################
##
## IMPORTANT(!): Units in Ferret are nm, kg,
## seconds, and attocoulombs
##
##############################################
u_x = u_x
u_y = u_y
u_z = u_z
[]
[Variables]
#################################
##
## Variable definitions
## P, u, phi, e^global_ij
## and their initial conditions
##
#################################
[./global_strain]
order = SIXTH
family = SCALAR
[../]
[./polar_x]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.1e-6
max = 0.1e-6
[../]
[../]
[./polar_y]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.1e-6
max = 0.1e-6
[../]
[../]
[./polar_z]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = 0.05
max = 0.1
[../]
[../]
[./potential_E_int]
order = FIRST
family = LAGRANGE
[../]
[./u_x]
order = FIRST
family = LAGRANGE
[../]
[./u_y]
order = FIRST
family = LAGRANGE
[../]
[./u_z]
order = FIRST
family = LAGRANGE
[../]
[]
[AuxVariables]
######################################
##
## Auxiarilly variable definitions
## (can be intermediate variables
## or for postprocessed quantities)
##
######################################
######################################
##
## Global displacements
##
######################################
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
######################################
##
## Stress/strain tensor components
##
######################################
[./e00]
order = CONSTANT
family = MONOMIAL
[../]
[./e01]
order = CONSTANT
family = MONOMIAL
[../]
[./e10]
order = CONSTANT
family = MONOMIAL
[../]
[./e11]
order = CONSTANT
family = MONOMIAL
[../]
[./e12]
order = CONSTANT
family = MONOMIAL
[../]
[./e22]
order = CONSTANT
family = MONOMIAL
[../]
[./s00]
order = CONSTANT
family = MONOMIAL
[../]
[./s01]
order = CONSTANT
family = MONOMIAL
[../]
[./s10]
order = CONSTANT
family = MONOMIAL
[../]
[./s11]
order = CONSTANT
family = MONOMIAL
[../]
[./s12]
order = CONSTANT
family = MONOMIAL
[../]
[./s22]
order = CONSTANT
family = MONOMIAL
[../]
[./eigs00]
order = CONSTANT
family = MONOMIAL
[../]
[./eigs11]
order = CONSTANT
family = MONOMIAL
[../]
[./eigs22]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
######################################
##
## Auxiarilly Kernel definitions
## (can be intermediate "operations"
## or for postprocessed quantities)
##
######################################
[./disp_x]
type = GlobalDisplacementAux
variable = disp_x
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 0
use_displaced_mesh = false
[../]
[./disp_y]
type = GlobalDisplacementAux
variable = disp_y
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 1
use_displaced_mesh = false
[../]
[./disp_z]
type = GlobalDisplacementAux
variable = disp_z
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 2
use_displaced_mesh = false
[../]
[./e00]
type = RankTwoAux
variable = e00
rank_two_tensor = total_strain
index_i = 0
index_j = 0
[../]
[./e01]
type = RankTwoAux
variable = e01
rank_two_tensor = total_strain
index_i = 0
index_j = 1
[../]
[./e10]
type = RankTwoAux
variable = e10
rank_two_tensor = total_strain
index_i = 1
index_j = 0
[../]
[./e12]
type = RankTwoAux
variable = e12
rank_two_tensor = total_strain
index_i = 1
index_j = 2
[../]
[./e11]
type = RankTwoAux
variable = e11
rank_two_tensor = total_strain
index_i = 1
index_j = 1
[../]
[./e22]
type = RankTwoAux
variable = e22
rank_two_tensor = total_strain
index_i = 2
index_j = 2
[../]
[./s00]
type = RankTwoAux
variable = s00
rank_two_tensor = stress
index_i = 0
index_j = 0
[../]
[./s01]
type = RankTwoAux
variable = s01
rank_two_tensor = stress
index_i = 0
index_j = 1
[../]
[./s10]
type = RankTwoAux
variable = s10
rank_two_tensor = stress
index_i = 1
index_j = 0
[../]
[./s12]
type = RankTwoAux
variable = s12
rank_two_tensor = stress
index_i = 1
index_j = 2
[../]
[./s11]
type = RankTwoAux
variable = s11
rank_two_tensor = stress
index_i = 1
index_j = 1
[../]
[./s22]
type = RankTwoAux
variable = s22
rank_two_tensor = stress
index_i = 2
index_j = 2
[../]
[./eigs00]
type = LocalABO3EigenstressAux
variable = eigs00
index_i = 0
index_j = 0
[../]
[./eigs11]
type = LocalABO3EigenstressAux
variable = eigs11
index_i = 1
index_j = 1
[../]
[./eigs22]
type = LocalABO3EigenstressAux
variable = eigs22
index_i = 2
index_j = 2
[../]
[]
[ScalarKernels]
######################################
##
## Necessary for PBC system
##
######################################
[./global_strain]
type = GlobalStrain
variable = global_strain
global_strain_uo = global_strain_uo
use_displaced_mesh = false
[../]
[]
[Materials]
#################################################
##
## Landau, electrostrictive, elastic coefficients
## M. Mtebwa, A. K. Tagantsev, and N. Setter
## AIP Adv. 4, 127150 (2014).
##
##################################################
[./Landau_P_FE]
type = GenericConstantMaterial
prop_names = 'alpha1 alpha11 alpha12 alpha111 alpha112 alpha123 alpha1111 alpha1112 alpha1122 alpha1123'
prop_values = '${a1temp} 0.04764 0.1336 0.1735 0.6128 -2.894 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'
[../]
[./mat_C_FE]
type = GenericConstantMaterial
prop_names = 'C11 C12 C44'
prop_values = '179.073 66.71 82.6446'
block = '0'
[../]
##################################################
##
## NOTE: Sign convention in Ferret for the
## electrostrictive coeff. is multiplied by
## an overall factor of (-1)
##
##################################################
[./mat_Q]
type = GenericConstantMaterial
prop_names = 'Q11 Q12 Q44'
prop_values = '-0.0966 0.046 -0.0819'
block = '0'
[../]
[./mat_q]
type = GenericConstantMaterial
prop_names = 'q11 q12 q44'
prop_values = '-11.1608 4.86165 -6.76859'
[../]
[./eigen_strain]
type = ComputeEigenstrain
# eigen_base = 'exx exy exz eyx eyy eyz ezx ezy ezz'
eigen_base = '1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0'
eigenstrain_name = eigenstrain
prefactor = 0.0
[../]
[./elasticity_tensor_1]
type = ComputeElasticityTensor
fill_method = symmetric9
###############################################
##
## symmetric9 fill_method is (default)
## C11 C12 C13 C22 C23 C33 C44 C55 C66
##
###############################################
C_ijkl = '179.073 66.71 66.71 179.073 66.71 179.073 82.6446 82.6446 82.6446'
[../]
[./strain_1]
type = ComputeSmallStrain
global_strain = global_strain
eigenstrain_names = eigenstrain
[../]
[./stress_1]
type = ComputeLinearElasticStress
[../]
[./global_strain]
type = ComputeGlobalStrain
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
[../]
[./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)
##
###############################################
type = GenericConstantMaterial
prop_names = 'permittivity'
prop_values = '0.08854187'
[../]
[]
[Kernels]
###############################################
##
## Physical Kernel operators
## to enforce TDLGD evolution
##
###############################################
#Elastic problem
[./TensorMechanics]
use_displaced_mesh = false
eigenstrains_name = eigenstrain
[../]
[./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
[../]
[./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
[../]
[./electrostr_ux]
type = ElectrostrictiveCouplingDispDerivative
variable = u_x
component = 0
[../]
[./electrostr_uy]
type = ElectrostrictiveCouplingDispDerivative
variable = u_y
component = 1
[../]
[./electrostr_uz]
type = ElectrostrictiveCouplingDispDerivative
variable = u_z
component = 2
[../]
[./electrostr_polar_coupled_x]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_x
component = 0
[../]
[./electrostr_polar_coupled_y]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_y
component = 1
[../]
[./electrostr_polar_coupled_z]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_z
component = 2
[../]
[./polar_x_electric_E]
type = PolarElectricEStrong
variable = potential_E_int
[../]
[./FE_E_int]
type = Electrostatics
variable = potential_E_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
[../]
[./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
[../]
[]
[BCs]
[./Periodic]
[./xyz]
auto_direction = 'x y z'
variable = 'u_x u_y u_z polar_x polar_y polar_z'
[../]
[../]
[./boundary_grounding]
type = DirichletBC
boundary = '0 1 2 3 4 5'
variable = potential_E_int
value = 0.0
[../]
# fix center point location
[./centerfix_x]
type = DirichletBC
boundary = 100
variable = u_x
value = 0
[../]
[./centerfix_y]
type = DirichletBC
boundary = 100
variable = u_y
value = 0
[../]
[./centerfix_z]
type = DirichletBC
boundary = 100
variable = u_z
value = 0
[../]
[]
[Postprocessors]
###############################################
##
## Postprocessors (integrations over the
## computational domain) to calculate the total energy
## decomposed into linear combinations of the
## different physics.
##
###############################################
[./avePx]
type = ElementAverageValue
variable = polar_x
execute_on = 'timestep_end'
[../]
[./avePy]
type = ElementAverageValue
variable = polar_y
execute_on = 'timestep_end'
[../]
[./avePz]
type = ElementAverageValue
variable = polar_z
execute_on = 'timestep_end'
[../]
[./ave_e00]
type = ElementAverageValue
variable = e00
execute_on = 'timestep_end'
[../]
[./ave_e11]
type = ElementAverageValue
variable = e11
execute_on = 'timestep_end'
[../]
[./ave_e22]
type = ElementAverageValue
variable = e22
execute_on = 'timestep_end'
[../]
[./Fb]
type = BulkEnergyEighth
execute_on = 'timestep_end'
[../]
[./Fw]
type = WallEnergy
execute_on = 'timestep_end'
[../]
[./Fela]
type = ElasticEnergy
execute_on = 'timestep_end'
use_displaced_mesh = false
[../]
[./Fc]
type = ElectrostrictiveCouplingEnergy
execute_on = 'timestep_end'
[../]
[./Fele]
type = ElectrostaticEnergy
execute_on = 'timestep_end'
[../]
[./Ftot]
type = LinearCombinationPostprocessor
pp_names = 'Fb Fw Fc Fele'
pp_coefs = ' 1 1 1 1'
execute_on = 'timestep_end'
[../]
[./vol]
type = VolumePostprocessor
execute_on = 'timestep_end'
[../]
[./px]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 0
[../]
[./py]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 1
[../]
[./pz]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 2
[../]
[./perc_change]
type = PercentChangePostprocessor
postprocessor = Ftot
execute_on = 'timestep_end'
[../]
[]
[UserObjects]
###############################################
##
## GlobalStrain system to enforce periodicity
## in the anisotropic strain field
##
###############################################
[./global_strain_uo]
type = GlobalATiO3MaterialRVEUserObject
use_displaced_mesh = false
execute_on = 'Initial Linear Nonlinear'
[../]
###############################################
##
## terminator to end energy evolution when the energy difference
## between subsequent time steps is lower than 5e-6
##
## NOTE: can fail if the time step is smallhotkey for tilde
##
###############################################
[./kill]
type = Terminator
expression = 'perc_change <= 1.0e-6'
[../]
[]
[Preconditioning]
###############################################
##
## Numerical preconditioning/solver options
##
###############################################
[./smp]
type = SMP
full = true
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-ksp_gmres_restart -snes_atol -snes_rtol -ksp_rtol -pc_type -build_twosided'
petsc_options_value = ' 160 1e-8 1e-6 1e-6 bjacobi allreduce'
[../]
[]
[Executioner]
##########################################
##
## Time integration/solver options
##
##########################################
type = Transient
solve_type = 'PJFNK'
scheme = 'implicit-euler'
dtmin = 1e-13
###########################################
##
## dtmax is material dependent!
##
###########################################
dtmax = 1.0
l_max_its = 200
[./TimeStepper]
type = IterationAdaptiveDT
optimal_iterations = 8
cutback_factor = 0.75
linear_iteration_ratio = 1000
dt = 0.3
[../]
verbose = true
[]
[Outputs]
###############################################
##
## Output options
##
###############################################
print_linear_residuals = false
perf_graph = false
[./outCSV]
type = CSV
file_base = out_pzt_monodomain_Tdef
[../]
[]
(examples/films/PTOfilm_e12_T298K_E0.i)
[Mesh]
[gen]
############################################
##
## Type and dimension of the mesh
##
############################################
type = GeneratedMeshGenerator
dim = 3
#############################################
##
## Grid definition. Note that it should be
## nJ = 2*(Jmax-Jmin) for J = x, y, z
##
#############################################
nx = 32
ny = 32
nz = 30
#############################################
##
## Actual spatial coordinates of mesh.
## Jmax - Jmin = nJ/2 for J = x, y, z
## Units are in nanometers
##
#############################################
xmin = -16.0
xmax = 16.0
ymin = -16.0
ymax = 16.0
zmin = -10.0
zmax = 20.0
#############################################
##
## FE type/order (hexahedral, tetrahedral
##
#############################################
elem_type = HEX8
[]
[./cnode]
input = gen
############################################
##
## additional boundary sideset (one node)
## to zero one of the elastic displacement vectors
## vectors and eliminates rigid body translations
## from the degrees of freedom
##
## NOTE: This must conform with the about
## [Mesh] block settings
##
############################################
type = ExtraNodesetGenerator
coord = '-16.0 -16.0 -10.0'
new_boundary = 100
[../]
[subdomains]
type = SubdomainBoundingBoxGenerator
input = cnode
bottom_left = '-16.0 -16.0 -10.0'
block_id = 1
top_right = '16.0 16.0 0'
location = INSIDE
[]
[film_interface]
type = SideSetsBetweenSubdomainsGenerator
input = subdomains
primary_block = 0
paired_block = 1
new_boundary = 52
[]
[film_surface]
type = SideSetsFromNormalsGenerator
input = film_interface
normals = '0 0 1'
fixed_normal = true
new_boundary = '107'
[]
[substrate_bottom]
type = SideSetsFromNormalsGenerator
input = film_surface
normals = '0 0 -1'
fixed_normal = true
new_boundary = '108'
[]
[]
[GlobalParams]
len_scale = 1.0
polar_x = polar_x
polar_y = polar_y
polar_z = polar_z
potential_E_int = potential_E_int
displacements = 'u_x u_y u_z'
##############################################
##=
## IMPORTANT(!): Units in Ferret are nm, kg,
## seconds, and attocoulombs
##
##############################################
vol = vol
u_x = u_x
u_y = u_y
u_z = u_z
[]
[Variables]
#################################
##
## Variable definitions
## P, u, phi, e^global_ij
## and their initial conditions
##
#################################
[./global_strain]
order = SIXTH
family = SCALAR
[../]
[./polar_x]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -1e-2
max = 1e-2
[../]
block = '0'
[../]
[./polar_y]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -1e-2
max = 1e-2
[../]
block = '0'
[../]
[./polar_z]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -1e-2
max = 1e-2
[../]
block = '0'
[../]
[./potential_E_int]
order = FIRST
family = LAGRANGE
block = '0 1'
[../]
[./u_x]
order = FIRST
family = LAGRANGE
block = '0 1'
[../]
[./u_y]
order = FIRST
family = LAGRANGE
block = '0 1'
[../]
[./u_z]
order = FIRST
family = LAGRANGE
block = '0 1'
[../]
[]
[AuxVariables]
######################################
##
## Auxiarilly variable definitions
## (can be intermediate variables
## or for postprocessed quantities)
##
######################################
######################################
##
## Global displacements
##
######################################
[./disp_x]
block = '0 1'
[../]
[./disp_y]
block = '0 1'
[../]
[./disp_z]
block = '0 1'
[../]
######################################
##
## Stress/strain tensor components
##
######################################
[./e00]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e01]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e10]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e11]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e12]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e22]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s00]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s01]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s10]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s11]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s12]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s22]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./E_x]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./E_y]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./E_z]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[]
[AuxKernels]
######################################
##
## Auxiarilly Kernel definitions
## (can be intermediate "operations"
## or for postprocessed quantities)
##
######################################
[./disp_x]
type = GlobalDisplacementAux
variable = disp_x
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 0
use_displaced_mesh = false
[../]
[./disp_y]
type = GlobalDisplacementAux
variable = disp_y
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 1
use_displaced_mesh = false
[../]
[./disp_z]
type = GlobalDisplacementAux
variable = disp_z
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 2
use_displaced_mesh = false
[../]
[./e00]
type = RankTwoAux
variable = e00
rank_two_tensor = total_strain
index_i = 0
index_j = 0
[../]
[./e01]
type = RankTwoAux
variable = e01
rank_two_tensor = total_strain
index_i = 0
index_j = 1
[../]
[./e10]
type = RankTwoAux
variable = e10
rank_two_tensor = total_strain
index_i = 1
index_j = 0
[../]
[./e12]
type = RankTwoAux
variable = e12
rank_two_tensor = total_strain
index_i = 1
index_j = 2
[../]
[./e11]
type = RankTwoAux
variable = e11
rank_two_tensor = total_strain
index_i = 1
index_j = 1
[../]
[./e22]
type = RankTwoAux
variable = e22
rank_two_tensor = total_strain
index_i = 2
index_j = 2
[../]
[./s00]
type = RankTwoAux
variable = s00
rank_two_tensor = stress
index_i = 0
index_j = 0
[../]
[./s01]
type = RankTwoAux
variable = s01
rank_two_tensor = stress
index_i = 0
index_j = 1
[../]
[./s10]
type = RankTwoAux
variable = s10
rank_two_tensor = stress
index_i = 1
index_j = 0
[../]
[./s12]
type = RankTwoAux
variable = s12
rank_two_tensor = stress
index_i = 1
index_j = 2
[../]
[./s11]
type = RankTwoAux
variable = s11
rank_two_tensor = stress
index_i = 1
index_j = 1
[../]
[./s22]
type = RankTwoAux
variable = s22
rank_two_tensor = stress
index_i = 2
index_j = 2
[../]
[./ex]
type = ElecFieldAux
variable = E_x
component = 0
[../]
[./ey]
type = ElecFieldAux
variable = E_y
component = 1
[../]
[./ez]
type = ElecFieldAux
variable = E_z
component = 2
[../]
[]
[ScalarKernels]
######################################
##
## Necessary for PBC system
##
######################################
[./global_strain]
type = GlobalStrain
variable = global_strain
global_strain_uo = global_strain_uo
use_displaced_mesh = false
[../]
[]
[Materials]
#################################################
##
## Landau coefficients from Li et al (2001)
##
##################################################
[./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_P_substr]
type = GenericConstantMaterial
prop_names = 'alpha1 alpha11 alpha12 alpha111 alpha112 alpha123 alpha1111 alpha1112 alpha1122 alpha1123'
prop_values = '10.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0'
block = '1'
[../]
[./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'
[../]
[./mat_C_FE]
type = GenericConstantMaterial
prop_names = 'C11 C12 C44'
prop_values = '175.0 79.4 111.1'
block = '0'
[../]
[./mat_C_sub]
type = GenericConstantMaterial
prop_names = 'C11 C12 C44'
prop_values = '220.0 34.4 161.1'
block = '1'
[../]
##################################################
##=
## NOTE: Sign convention in Ferret for the
## electrostrictive coeff. is multiplied by
## an overall factor of (-1)
##
##################################################
[./mat_Q]
type = GenericConstantMaterial
prop_names = 'Q11 Q12 Q44'
prop_values = '-0.089 0.026 -0.03375'
block = '0 1'
[../]
[./mat_q]
type = GenericConstantMaterial
prop_names = 'q11 q12 q44'
prop_values = '-11.4 -0.01438 -7.5'
[../]
[./eigen_strain]
type = ComputeEigenstrain
# eigen_base = 'exx exy exz eyx eyy eyz ezx ezy ezz'
eigen_base = '1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0'
eigenstrain_name = eigenstrain
prefactor = 0.0
[../]
[./elasticity_tensor_1]
type = ComputeElasticityTensor
fill_method = symmetric9
###############################################
##
## symmetric9 fill_method is (default)
## C11 C12 C13 C22 C23 C33 C44 C55 C66
##
###############################################
C_ijkl = '175.0 79.4 79.4 175.0 79.4 175.0 111.1 111.1 111.1'
[../]
[./strain_1]
type = ComputeSmallStrain
global_strain = global_strain
eigenstrain_names = eigenstrain
[../]
[./stress_1]
type = ComputeLinearElasticStress
[../]
[./global_strain]
type = ComputeGlobalStrain
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
[../]
[./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)
##
###############################################
type = GenericConstantMaterial
prop_names = 'permittivity'
prop_values = '0.08854187'
[../]
[]
[Kernels]
###############################################
##
## Physical Kernel operators
## to enforce TDLGD evolution
##
###############################################
#Elastic problem
[./TensorMechanics]
use_displaced_mesh = false
eigenstrain_names = eigenstrain
[../]
[./bed_x]
type = BulkEnergyDerivativeEighth
variable = polar_x
component = 0
block = '0'
[../]
[./bed_y]
type = BulkEnergyDerivativeEighth
variable = polar_y
component = 1
block = '0'
[../]
[./bed_z]
type = BulkEnergyDerivativeEighth
variable = polar_z
component = 2
block = '0'
[../]
[./walled_x]
type = WallEnergyDerivative
variable = polar_x
component = 0
block = '0'
[../]
[./walled_y]
type = WallEnergyDerivative
variable = polar_y
component = 1
block = '0'
[../]
[./walled_z]
type = WallEnergyDerivative
variable = polar_z
component = 2
block = '0'
[../]
[./electrostr_ux]
type = ElectrostrictiveCouplingDispDerivative
variable = u_x
component = 0
block = '0'
[../]
[./electrostr_uy]
type = ElectrostrictiveCouplingDispDerivative
variable = u_y
component = 1
block = '0'
[../]
[./electrostr_uz]
type = ElectrostrictiveCouplingDispDerivative
variable = u_z
component = 2
block = '0'
[../]
[./electrostr_polar_coupled_x]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_x
component = 0
block = '0'
[../]
[./electrostr_polar_coupled_y]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_y
component = 1
block = '0'
[../]
[./electrostr_polar_coupled_z]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_z
component = 2
block = '0'
[../]
[./polar_x_electric_E]
type = PolarElectricEStrong
variable = potential_E_int
block = '0'
[../]
[./FE_E_int]
type = Electrostatics
variable = potential_E_int
block = '0 1'
[../]
[./polar_electric_px]
type = PolarElectricPStrong
variable = polar_x
component = 0
block = '0'
[../]
[./polar_electric_py]
type = PolarElectricPStrong
variable = polar_y
component = 1
block = '0'
[../]
[./polar_electric_pz]
type = PolarElectricPStrong
variable = polar_z
component = 2
block = '0'
[../]
[./polar_x_time]
type = TimeDerivativeScaled
variable=polar_x
time_scale = 1.0
block = '0'
[../]
[./polar_y_time]
type = TimeDerivativeScaled
variable = polar_y
time_scale = 1.0
block = '0'
[../]
[./polar_z_time]
type = TimeDerivativeScaled
variable = polar_z
time_scale = 1.0
block = '0'
[../]
[./u_x_time]
type = TimeDerivativeScaled
variable = u_x
time_scale = 1.0
[../]
[./u_y_time]
type = TimeDerivativeScaled
variable = u_y
time_scale = 1.0
[../]
[./u_z_time]
type = TimeDerivativeScaled
variable = u_z
time_scale = 1.0
[../]
[]
[BCs]
[./Periodic]
[./xy]
auto_direction = 'x y'
variable = 'u_x u_y u_z polar_x polar_y polar_z potential_E_int'
[../]
[../]
[./boundary_interface_grounding]
type = DirichletBC
boundary = '52'
variable = potential_E_int
value = 0.0
[../]
# fix center point location
[./centerfix_x]
type = DirichletBC
boundary = '108'
variable = u_x
value = 0
[../]
[./centerfix_y]
type = DirichletBC
boundary = '108'
variable = u_y
value = 0
[../]
[./centerfix_z]
type = DirichletBC
boundary = '108'
variable = u_z
value = 0
[../]
[]
[Postprocessors]
###############################################
##=
## Postprocessors (integrations over the
## computational domain) to calculate the total
## energy decomposed into linear combinations of
## the different physics.
##
###############################################
[./Fbulk]
type = BulkEnergyEighth
execute_on = 'timestep_end'
block = '0'
[../]
[./Fwall]
type = WallEnergy
execute_on = 'timestep_end'
block = '0'
[../]
[./Felastic]
type = ElasticEnergy
execute_on = 'timestep_end'
use_displaced_mesh = false
block = '0'
[../]
[./Fcoupled]
type = ElectrostrictiveCouplingEnergy
execute_on = 'timestep_end'
block = '0'
[../]
[./Felec]
type = ElectrostaticEnergy
execute_on = 'timestep_end'
block = '0'
[../]
[./Ftotal]
type = LinearCombinationPostprocessor
pp_names = 'Fbulk Fwall Fcoupled Felec'
pp_coefs = '0.160218 0.160218 0.160218 0.160218' #converted to eV
execute_on = 'timestep_end'
[../]
[./vol]
type = VolumePostprocessor
execute_on = 'timestep_end'
[../]
[./px]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 0
block = '0'
[../]
[./py]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 1
block = '0'
[../]
[./pz]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 2
block = '0'
[../]
[./perc_change]
type = PercentChangePostprocessor
postprocessor = Ftotal
execute_on = 'timestep_end'
[../]
[./elapsed]
type = PerfGraphData
section_name = "Root" # for profiling the problem [on]
data_type = total
[../]
[]
[UserObjects]
###############################################
##
## GlobalStrain system to enforce periodicity
## in the anisotropic strain field
##
###############################################
[./global_strain_uo]
type = GlobalATiO3MaterialRVEUserObject
use_displaced_mesh = false
execute_on = 'Initial Linear Nonlinear'
applied_stress_tensor = '2.1 2.1 1.9056 0.0 0.0 0.0'
block = '0'
[../]
###############################################
##
## terminator to end energy evolution when the energy difference
## between subsequent time steps is lower than 5e-6
##
## NOTE: can fail if the time step is small
##
###############################################
[./kill]
type = Terminator
expression = 'perc_change <= 5.0e-4'
[../]
[]
[Preconditioning]
###############################################
##
## Numerical preconditioning/solver options
##
###############################################
[./smp]
type = SMP
full = true
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-ksp_gmres_restart -snes_atol -snes_rtol -ksp_rtol -pc_type -build_twosided'
petsc_options_value = ' 80 1e-8 1e-6 1e-5 bjacobi allreduce'
[../]
[]
[Executioner]
##########################################
##
## Time integ=ration/solver options
##
##########################################
type = Transient
solve_type = 'PJFNK'
scheme = 'bdf2'
dtmin = 1e-13
dtmax = 0.6
l_max_its = 200
[./TimeStepper]
type = IterationAdaptiveDT
optimal_iterations = 6
growth_factor = 1.2
cutback_factor = 0.75
linear_iteration_ratio = 1000
dt = 0.6
[../]
verbose = true
nl_max_its = 20
[]
[Outputs]
###############################################
##==
## Output options
##
###############################################
print_linear_residuals = false
perf_graph = false
[./out]
type = Exodus
file_base = out_PTOfilm_e12_T298K_E0_E0
elemental_as_nodal = true
interval = 1
[../]
[]
(examples/monodomain/PTO_monodomain_Tdef.i)
a1temp = a1def
[Mesh]
[gen]
############################################
##
## Type and dimension of the mesh
##
############################################
type = GeneratedMeshGenerator
dim = 3
#############################################
##
## Grid definition. Note that it should be
## nJ = 2*(Jmax-Jmin) for J = x, y, z
##
#############################################
nx = 8
ny = 8
nz = 8
#############################################
##
## Actual spatial coordinates of mesh.
## Jmax - Jmin = nJ/2 for J = x, y, z
## Units are in nanometers
##
#############################################
xmin = -2.0
xmax = 2.0
ymin = -2.0
ymax = 2.0
zmin = -2.0
zmax = 2.0
#############################################
##
## FE type/order (hexahedral, tetrahedral
##
#############################################
elem_type = HEX8
[]
[./cnode]
input = gen
############################################
##
## additional boundary sideset (one node)
## to zero one of the elastic displacement vectors
## vectors and eliminates rigid body translations
## from the degrees of freedom
##
## NOTE: This must conform with the about
## [Mesh] block settings
##
############################################
type = ExtraNodesetGenerator
coord = '0.0 0.0 0.0'
new_boundary = 100
[../]
[]
[GlobalParams]
len_scale = 1.0
polar_x = polar_x
polar_y = polar_y
polar_z = polar_z
potential_E_int = potential_E_int
vol = vol
displacements = 'u_x u_y u_z'
##############################################
##
## IMPORTANT(!): Units in Ferret are nm, kg,
## seconds, and attocoulombs
##
##############################################
u_x = u_x
u_y = u_y
u_z = u_z
[]
[Variables]
#################################
##
## Variable definitions
## P, u, phi, e^global_ij
## and their initial conditions
##
#################################
[./global_strain]
order = SIXTH
family = SCALAR
[../]
[./polar_x]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.1e-6
max = 0.1e-6
[../]
[../]
[./polar_y]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.1e-6
max = 0.1e-6
[../]
[../]
[./polar_z]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = 0.05
max = 0.1
[../]
[../]
[./potential_E_int]
order = FIRST
family = LAGRANGE
[../]
[./u_x]
order = FIRST
family = LAGRANGE
[../]
[./u_y]
order = FIRST
family = LAGRANGE
[../]
[./u_z]
order = FIRST
family = LAGRANGE
[../]
[]
[AuxVariables]
######################################
##
## Auxiarilly variable definitions
## (can be intermediate variables
## or for postprocessed quantities)
##
######################################
######################################
##
## Global displacements
##
######################################
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
######################################
##
## Stress/strain tensor components
##
######################################
[./e00]
order = CONSTANT
family = MONOMIAL
[../]
[./e01]
order = CONSTANT
family = MONOMIAL
[../]
[./e10]
order = CONSTANT
family = MONOMIAL
[../]
[./e11]
order = CONSTANT
family = MONOMIAL
[../]
[./e12]
order = CONSTANT
family = MONOMIAL
[../]
[./e22]
order = CONSTANT
family = MONOMIAL
[../]
[./s00]
order = CONSTANT
family = MONOMIAL
[../]
[./s01]
order = CONSTANT
family = MONOMIAL
[../]
[./s10]
order = CONSTANT
family = MONOMIAL
[../]
[./s11]
order = CONSTANT
family = MONOMIAL
[../]
[./s12]
order = CONSTANT
family = MONOMIAL
[../]
[./s22]
order = CONSTANT
family = MONOMIAL
[../]
[./eigs00]
order = CONSTANT
family = MONOMIAL
[../]
[./eigs11]
order = CONSTANT
family = MONOMIAL
[../]
[./eigs22]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
######################################
##
## Auxiarilly Kernel definitions
## (can be intermediate "operations"
## or for postprocessed quantities)
##
######################################
[./disp_x]
type = GlobalDisplacementAux
variable = disp_x
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 0
use_displaced_mesh = false
[../]
[./disp_y]
type = GlobalDisplacementAux
variable = disp_y
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 1
use_displaced_mesh = false
[../]
[./disp_z]
type = GlobalDisplacementAux
variable = disp_z
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 2
use_displaced_mesh = false
[../]
[./e00]
type = RankTwoAux
variable = e00
rank_two_tensor = total_strain
index_i = 0
index_j = 0
[../]
[./e01]
type = RankTwoAux
variable = e01
rank_two_tensor = total_strain
index_i = 0
index_j = 1
[../]
[./e10]
type = RankTwoAux
variable = e10
rank_two_tensor = total_strain
index_i = 1
index_j = 0
[../]
[./e12]
type = RankTwoAux
variable = e12
rank_two_tensor = total_strain
index_i = 1
index_j = 2
[../]
[./e11]
type = RankTwoAux
variable = e11
rank_two_tensor = total_strain
index_i = 1
index_j = 1
[../]
[./e22]
type = RankTwoAux
variable = e22
rank_two_tensor = total_strain
index_i = 2
index_j = 2
[../]
[./s00]
type = RankTwoAux
variable = s00
rank_two_tensor = stress
index_i = 0
index_j = 0
[../]
[./s01]
type = RankTwoAux
variable = s01
rank_two_tensor = stress
index_i = 0
index_j = 1
[../]
[./s10]
type = RankTwoAux
variable = s10
rank_two_tensor = stress
index_i = 1
index_j = 0
[../]
[./s12]
type = RankTwoAux
variable = s12
rank_two_tensor = stress
index_i = 1
index_j = 2
[../]
[./s11]
type = RankTwoAux
variable = s11
rank_two_tensor = stress
index_i = 1
index_j = 1
[../]
[./s22]
type = RankTwoAux
variable = s22
rank_two_tensor = stress
index_i = 2
index_j = 2
[../]
[./eigs00]
type = LocalABO3EigenstressAux
variable = eigs00
index_i = 0
index_j = 0
[../]
[./eigs11]
type = LocalABO3EigenstressAux
variable = eigs11
index_i = 1
index_j = 1
[../]
[./eigs22]
type = LocalABO3EigenstressAux
variable = eigs22
index_i = 2
index_j = 2
[../]
[]
[ScalarKernels]
######################################
##
## Necessary for PBC system
##
######################################
[./global_strain]
type = GlobalStrain
variable = global_strain
global_strain_uo = global_strain_uo
use_displaced_mesh = false
[../]
[]
[Materials]
#################################################
##
## Bulk free energy and electrostrictive
## coefficients gleaned from
## Marton and Hlinka
## Phys. Rev. B. 74, 104014, (2006)
##
## NOTE: there might be some Legendre transforms
## depending on what approach you use
## -i.e. inhomogeneous strain vs
## homogeneous strain [renormalized]
##
##################################################
[./Landau_P]
type = GenericConstantMaterial
prop_names = 'alpha1 alpha11 alpha12 alpha111 alpha112 alpha123 alpha1111 alpha1112 alpha1122 alpha1123'
prop_values = '${a1temp} -0.073 0.75 0.26 0.61 -3.67 0.0 0.0 0.0 0.0'
[../]
############################################
##
## Gradient coefficients from
## Marton and Hlinka
## Phys. Rev. B. 74, 104014, (2006)
##
############################################
[./Landau_G]
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'
[../]
[./mat_C]
type = GenericConstantMaterial
prop_names = 'C11 C12 C44'
prop_values = '175.0 79.4 111.1'
[../]
##################################################
##
## NOTE: Sign convention in Ferret for the
## electrostrictive coeff. is multiplied by
## an overall factor of (-1)
##
##################################################
[./mat_Q]
type = GenericConstantMaterial
prop_names = 'Q11 Q12 Q44'
prop_values = '-0.089 0.026 -0.03375'
[../]
[./mat_q]
type = GenericConstantMaterial
prop_names = 'q11 q12 q44'
prop_values = '-11.4 -0.01438 -7.5'
[../]
[./eigen_strain]
type = ComputeEigenstrain
eigen_base = '1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0'
eigenstrain_name = eigenstrain
prefactor = 0.0
[../]
[./elasticity_tensor_1]
type = ComputeElasticityTensor
fill_method = symmetric9
###############################################
##
## symmetric9 fill_method is (default)
## C11 C12 C13 C22 C23 C33 C44 C55 C66
##
###############################################
C_ijkl = '175.0 79.4 79.4 175.0 79.4 175.0 111.1 111.1 111.1'
[../]
[./strain_1]
type = ComputeSmallStrain
global_strain = global_strain
eigenstrain_names = eigenstrain
[../]
[./stress_1]
type = ComputeLinearElasticStress
[../]
[./global_strain]
type = ComputeGlobalStrain
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
[../]
[./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)
##
###############################################
type = GenericConstantMaterial
prop_names = 'permittivity'
prop_values = '0.08854187'
[../]
[]
[Kernels]
###############################################
##
## Physical Kernel operators
## to enforce TDLGD evolution
##
###############################################
#Elastic problem
[./TensorMechanics]
use_displaced_mesh = false
eigenstrains_name = eigenstrain
[../]
[./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
[../]
[./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
[../]
[./electrostr_ux]
type = ElectrostrictiveCouplingDispDerivative
variable = u_x
component = 0
[../]
[./electrostr_uy]
type = ElectrostrictiveCouplingDispDerivative
variable = u_y
component = 1
[../]
[./electrostr_uz]
type = ElectrostrictiveCouplingDispDerivative
variable = u_z
component = 2
[../]
[./electrostr_polar_coupled_x]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_x
component = 0
[../]
[./electrostr_polar_coupled_y]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_y
component = 1
[../]
[./electrostr_polar_coupled_z]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_z
component = 2
[../]
[./polar_x_electric_E]
type = PolarElectricEStrong
variable = potential_E_int
[../]
[./FE_E_int]
type = Electrostatics
variable = potential_E_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
[../]
[./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
[../]
[]
[BCs]
[./Periodic]
[./xyz]
auto_direction = 'x y z'
variable = 'u_x u_y u_z polar_x polar_y polar_z'
[../]
[../]
[./boundary_grounding]
type = DirichletBC
boundary = '0 1 2 3 4 5'
variable = potential_E_int
value = 0.0
[../]
# fix center point location
[./centerfix_x]
type = DirichletBC
boundary = 100
variable = u_x
value = 0
[../]
[./centerfix_y]
type = DirichletBC
boundary = 100
variable = u_y
value = 0
[../]
[./centerfix_z]
type = DirichletBC
boundary = 100
variable = u_z
value = 0
[../]
[]
[Postprocessors]
###############################################
##
## Postprocessors (integrations over the
## computational domain) to calculate the total energy
## decomposed into linear combinations of the
## different physics.
##
###############################################
[./avePx]
type = ElementAverageValue
variable = polar_x
execute_on = 'timestep_end'
[../]
[./avePy]
type = ElementAverageValue
variable = polar_y
execute_on = 'timestep_end'
[../]
[./avePz]
type = ElementAverageValue
variable = polar_z
execute_on = 'timestep_end'
[../]
[./ave_e00]
type = ElementAverageValue
variable = e00
execute_on = 'timestep_end'
[../]
[./ave_e11]
type = ElementAverageValue
variable = e11
execute_on = 'timestep_end'
[../]
[./ave_e22]
type = ElementAverageValue
variable = e22
execute_on = 'timestep_end'
[../]
[./Fb]
type = BulkEnergyEighth
execute_on = 'timestep_end'
[../]
[./Fw]
type = WallEnergy
execute_on = 'timestep_end'
[../]
[./Fela]
type = ElasticEnergy
execute_on = 'timestep_end'
use_displaced_mesh = false
[../]
[./Fc]
type = ElectrostrictiveCouplingEnergy
execute_on = 'timestep_end'
[../]
[./Fele]
type = ElectrostaticEnergy
execute_on = 'timestep_end'
[../]
[./Ftot]
type = LinearCombinationPostprocessor
pp_names = 'Fb Fw Fc Fele'
pp_coefs = ' 1 1 1 1'
execute_on = 'timestep_end'
[../]
[./vol]
type = VolumePostprocessor
execute_on = 'timestep_end'
[../]
[./px]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 0
[../]
[./py]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 1
[../]
[./pz]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 2
[../]
[./perc_change]
type = PercentChangePostprocessor
postprocessor = Ftot
execute_on = 'timestep_end'
[../]
[]
[UserObjects]
###############################################
##
## GlobalStrain system to enforce periodicity
## in the anisotropic strain field
##
###############################################
[./global_strain_uo]
type = GlobalATiO3MaterialRVEUserObject
use_displaced_mesh = false
execute_on = 'Initial Linear Nonlinear'
[../]
###############################################
##
## terminator to end energy evolution when the energy difference
## between subsequent time steps is lower than 5e-6
##
## NOTE: can fail if the time step is smallhotkey for tilde
##
###############################################
[./kill]
type = Terminator
expression = 'perc_change <= 1.0e-6'
[../]
[]
[Preconditioning]
###############################################
##
## Numerical preconditioning/solver options
##
###############################################
[./smp]
type = SMP
full = true
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-ksp_gmres_restart -snes_atol -snes_rtol -ksp_rtol -pc_type -build_twosided'
petsc_options_value = ' 160 1e-8 1e-6 1e-6 bjacobi allreduce'
[../]
[]
[Executioner]
##########################################
##
## Time integration/solver options
##
##########################################
type = Transient
solve_type = 'PJFNK'
scheme = 'implicit-euler'
dtmin = 1e-13
###########################################
##
## dtmax is material dependent!
##
###########################################
dtmax = 1.0
l_max_its = 200
[./TimeStepper]
type = IterationAdaptiveDT
optimal_iterations = 8
cutback_factor = 0.75
linear_iteration_ratio = 1000
dt = 0.3
[../]
verbose = true
[]
[Outputs]
###############################################
##
## Output options
##
###############################################
print_linear_residuals = false
perf_graph = false
[./outCSV]
type = CSV
file_base = out_pto_monodomain_Tdef
[../]
[]
(test/tests/film/PTO_film3_T298K_noP.i)
[Mesh]
[gen]
############################################
##=
## Type and dimension of the mesh
##=
############################################
type = GeneratedMeshGenerator
dim = 3
#############################################
##
## Grid definition. Note that it should be
## nJ = 2*(Jmax-Jmin) for J = x, y, z
##=
#############################################
nx = 4
ny = 4
nz = 18 #need approximately 90 processors for this since dof is 953133 and we want 10k dof/proc
#############################################
##
## Actual spatial coordinates of mesh.
## Jmax - Jmin = nJ/2 for J = x, y, z
## Units are in nanometers
##
#############################################
xmin = -2.0
xmax = 2.0
ymin = -2.0
ymax = 2.0
zmin = -10.0
zmax = 8.0
#############################################
##
## FE type/order (hexahedral, tetrahedral
##
#############################################
elem_type = HEX8
[]
[./cnode]
input = gen
############################################
##
## additional boundary sideset (one node)
## to zero one of the elastic displacement vectors
## vectors and eliminates rigid body translations
## from the degrees of freedom
##
## NOTE: This must conform with the about
## [Mesh] block settings
##=
############################################
type = ExtraNodesetGenerator
coord = '-2.0 -2.0 -10.0'
new_boundary = 100
[../]
[subdomains]
type = SubdomainBoundingBoxGenerator
input = cnode
bottom_left = '-2.0 -2.0 -10.0'
block_id = 1
top_right = '2.0 2.0 0'
location = INSIDE
[]
[film_interface]
type = SideSetsBetweenSubdomainsGenerator
input = subdomains
primary_block = 0
paired_block = 1
new_boundary = 52
[]
[]
[GlobalParams]
len_scale = 1.0
polar_x = polar_x
polar_y = polar_y
polar_z = polar_z
potential_E_int = potential_E_int
displacements = 'u_x u_y u_z'
vol = vol
##############################################
##=
## IMPORTANT(!): Units in Ferret are nm, kg,
## seconds, and attocoulombs
##
##############################################
u_x = u_x
u_y = u_y
u_z = u_z
[]
[Variables]
#################################
##
## Variable definitions
## P, u, phi, e^global_ij
## and their initial conditions
##
#################################
[./global_strain]
order = SIXTH
family = SCALAR
[../]
[./polar_x]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.1e-5
max = 0.1e-5
[../]
block = '0'
[../]
[./polar_y]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.1e-5
max = 0.1e-5
[../]
block = '0'
[../]
[./polar_z]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.1e-5
max = 0.1e-5
[../]
block = '0'
[../]
[./potential_E_int]
order = FIRST
family = LAGRANGE
block = '0 1'
[../]
[./u_x]
order = FIRST
family = LAGRANGE
block = '0 1'
[../]
[./u_y]
order = FIRST
family = LAGRANGE
block = '0 1'
[../]
[./u_z]
order = FIRST
family = LAGRANGE
block = '0 1'
[../]
[]
[AuxVariables]
######################################
##
## Auxiarilly variable definitions
## (can be intermediate variables
## or for postprocessed quantities)
##
######################################
######################################
##
## Global displacements
##
######################################
[./disp_x]
block = '0 1'
[../]
[./disp_y]
block = '0 1'
[../]
[./disp_z]
block = '0 1'
[../]
######################################
##
## Stress/strain tensor components
##
######################################
[./e00]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e01]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e10]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e11]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e12]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./e22]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s00]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s01]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s10]
order = CONSTANT
family = MONOMIAL
[../]
[./s11]
order = CONSTANT
family = MONOMIAL
[../]
[./s12]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[./s22]
order = CONSTANT
family = MONOMIAL
block = '0 1'
[../]
[]
[AuxKernels]
######################################
##
## Auxiarilly Kernel definitions
## (can be intermediate "operations"
## or for postprocessed quantities)
##
######################################
[./disp_x]
type = GlobalDisplacementAux
variable = disp_x
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 0
use_displaced_mesh = false
[../]
[./disp_y]
type = GlobalDisplacementAux
variable = disp_y
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 1
use_displaced_mesh = false
[../]
[./disp_z]
type = GlobalDisplacementAux
variable = disp_z
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
component = 2
use_displaced_mesh = false
[../]
[./e00]
type = RankTwoAux
variable = e00
rank_two_tensor = total_strain
index_i = 0
index_j = 0
[../]
[./e01]
type = RankTwoAux
variable = e01
rank_two_tensor = total_strain
index_i = 0
index_j = 1
[../]
[./e10]
type = RankTwoAux
variable = e10
rank_two_tensor = total_strain
index_i = 1
index_j = 0
[../]
[./e12]
type = RankTwoAux
variable = e12
rank_two_tensor = total_strain
index_i = 1
index_j = 2
[../]
[./e11]
type = RankTwoAux
variable = e11
rank_two_tensor = total_strain
index_i = 1
index_j = 1
[../]
[./e22]
type = RankTwoAux
variable = e22
rank_two_tensor = total_strain
index_i = 2
index_j = 2
[../]
[./s00]
type = RankTwoAux
variable = s00
rank_two_tensor = stress
index_i = 0
index_j = 0
[../]
[./s01]
type = RankTwoAux
variable = s01
rank_two_tensor = stress
index_i = 0
index_j = 1
[../]
[./s10]
type = RankTwoAux
variable = s10
rank_two_tensor = stress
index_i = 1
index_j = 0
[../]
[./s12]
type = RankTwoAux
variable = s12
rank_two_tensor = stress
index_i = 1
index_j = 2
[../]
[./s11]
type = RankTwoAux
variable = s11
rank_two_tensor = stress
index_i = 1
index_j = 1
[../]
[./s22]
type = RankTwoAux
variable = s22
rank_two_tensor = stress
index_i = 2
index_j = 2
[../]
[]
[ScalarKernels]
######################################
##
## Necessary for PBC system
##
######################################
[./global_strain]
type = GlobalStrain
variable = global_strain
global_strain_uo = global_strain_uo
use_displaced_mesh = false
[../]
[]
[Materials]
#################################################
##
## add comments
##
## NOTE: there might be some Legendre transforms
## depending on what approach you use
## -i.e. inhomogeneous strain vs
## homogeneous strain [renormalized]
##=
##################################################
[./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_P_substr]
type = GenericConstantMaterial
prop_names = 'alpha1 alpha11 alpha12 alpha111 alpha112 alpha123 alpha1111 alpha1112 alpha1122 alpha1123'
prop_values = '10.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0'
block = '1'
[../]
############################################
##
## add comments
##
############################################
[./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'
[../]
[./mat_C_FE]
type = GenericConstantMaterial
prop_names = 'C11 C12 C44'
prop_values = '175.0 79.4 111.1'
block = '0'
[../]
[./mat_C_sub]
type = GenericConstantMaterial
prop_names = 'C11 C12 C44'
prop_values = '220.0 34.4 161.1'
block = '1'
[../]
##################################################
##=
## NOTE: Sign convention in Ferret for the
## electrostrictive coeff. is multiplied by
## an overall factor of (-1)
##
##################################################
[./mat_Q]
type = GenericConstantMaterial
prop_names = 'Q11 Q12 Q44'
prop_values = '-0.089 0.026 -0.03375'
block = '0 1'
[../]
[./mat_q]
type = GenericConstantMaterial
prop_names = 'q11 q12 q44'
prop_values = '-11.4 -0.01438 -7.5'
[../]
[./eigen_strain]
type = ComputeEigenstrain
# eigen_base = 'exx exy exz eyx eyy eyz ezx ezy ezz'
eigen_base = '1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0'
eigenstrain_name = eigenstrain
prefactor = 0.0
[../]
[./elasticity_tensor_1]
type = ComputeElasticityTensor
fill_method = symmetric9
###############################################
##
## symmetric9 fill_method is (default)
## C11 C12 C13 C22 C23 C33 C44 C55 C66
##
###############################################
C_ijkl = '175.0 79.4 79.4 175.0 79.4 175.0 111.1 111.1 111.1'
[../]
[./strain_1]
type = ComputeSmallStrain
global_strain = global_strain
eigenstrain_names = eigenstrain
[../]
[./stress_1]
type = ComputeLinearElasticStress
[../]
[./global_strain]
type = ComputeGlobalStrain
scalar_global_strain = global_strain
global_strain_uo = global_strain_uo
[../]
[./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)
##
###############################################
type = GenericConstantMaterial
prop_names = 'permittivity'
prop_values = '0.08854187'
[../]
[]
[Kernels]
###############################################
##
## Physical Kernel operators
## to enforce TDLGD evolution
##
###############################################
#Elastic problem
[./TensorMechanics]
use_displaced_mesh = false
eigenstrain_names = eigenstrain
[../]
[./bed_x]
type = BulkEnergyDerivativeEighth
variable = polar_x
component = 0
block = '0'
[../]
[./bed_y]
type = BulkEnergyDerivativeEighth
variable = polar_y
component = 1
block = '0'
[../]
[./bed_z]
type = BulkEnergyDerivativeEighth
variable = polar_z
component = 2
block = '0'
[../]
[./walled_x]
type = WallEnergyDerivative
variable = polar_x
component = 0
block = '0'
[../]
[./walled_y]
type = WallEnergyDerivative
variable = polar_y
component = 1
block = '0'
[../]
[./walled_z]
type = WallEnergyDerivative
variable = polar_z
component = 2
block = '0'
[../]
[./electrostr_ux]
type = ElectrostrictiveCouplingDispDerivative
variable = u_x
component = 0
block = '0'
[../]
[./electrostr_uy]
type = ElectrostrictiveCouplingDispDerivative
variable = u_y
component = 1
block = '0'
[../]
[./electrostr_uz]
type = ElectrostrictiveCouplingDispDerivative
variable = u_z
component = 2
block = '0'
[../]
[./electrostr_polar_coupled_x]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_x
component = 0
block = '0'
[../]
[./electrostr_polar_coupled_y]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_y
component = 1
block = '0'
[../]
[./electrostr_polar_coupled_z]
type = ElectrostrictiveCouplingPolarDerivative
variable = polar_z
component = 2
block = '0'
[../]
[./polar_x_electric_E]
type = PolarElectricEStrong
variable = potential_E_int
block = '0'
[../]
[./FE_E_int]
type = Electrostatics
variable = potential_E_int
block = '0 1'
[../]
[./polar_electric_px]
type = PolarElectricPStrong
variable = polar_x
component = 0
block = '0'
[../]
[./polar_electric_py]
type = PolarElectricPStrong
variable = polar_y
component = 1
block = '0'
[../]
[./polar_electric_pz]
type = PolarElectricPStrong
variable = polar_z
component = 2
block = '0'
[../]
[./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
[../]
[]
[BCs]
[./Periodic]
[./xy]
auto_direction = 'x y'
variable = 'u_x u_y u_z polar_x polar_y polar_z'
[../]
[../]
[./boundary_grounding]
type = DirichletBC
boundary = '52'
variable = potential_E_int
value = 0.0
[../]
# fix center point location
[./centerfix_x]
type = DirichletBC
boundary = 'back'
variable = u_x
value = 0
[../]
[./centerfix_y]
type = DirichletBC
boundary = 'back'
variable = u_y
value = 0
[../]
[./centerfix_z]
type = DirichletBC
boundary = 'back'
variable = u_z
value = 0
[../]
[]
[Postprocessors]
###############################################
##=
## Postprocessors (integrations over the
## computational domain) to calculate the total energy
## decomposed into linear combinations of the
## different physics.
##
###############################################
[./Fbulk]
type = BulkEnergyEighth
execute_on = 'timestep_end'
block = '0'
[../]
[./Fwall]
type = WallEnergy
execute_on = 'timestep_end'
block = '0'
[../]
[./Felastic]
type = ElasticEnergy
execute_on = 'timestep_end'
use_displaced_mesh = false
block = '0'
[../]
[./Fcoupled]
type = ElectrostrictiveCouplingEnergy
execute_on = 'timestep_end'
u_x = u_x
u_y = u_y
u_z = u_z
block = '0'
[../]
[./Felec]
type = ElectrostaticEnergy
execute_on = 'timestep_end'
block = '0'
[../]
[./Ftotal]
type = LinearCombinationPostprocessor
pp_names = 'Fbulk Fwall Fcoupled Felec'
pp_coefs = ' 1 1 1 1'
execute_on = 'timestep_end'
[../]
[./vol]
type = VolumePostprocessor
execute_on = 'timestep_end'
[../]
[./px]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 0
block = '0'
[../]
[./py]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 1
block = '0'
[../]
[./pz]
type = DomainVariantPopulation
execute_on = 'timestep_end'
component = 2
block = '0'
[../]
[./perc_change]
type = PercentChangePostprocessor
postprocessor = Ftotal
execute_on = 'timestep_end'
[../]
[]
[UserObjects]
###############################################
##=
## GlobalStrain system to enforce periodicity
## in the anisotropic strain field
##
###############################################
[./global_strain_uo]
type = GlobalATiO3MaterialRVEUserObject
use_displaced_mesh = false
execute_on = 'Initial Linear Nonlinear'
applied_stress_tensor = '4.375 4.375 3.97 0.0 0.0 0.0'
block = '0'
[../]
###############################################
##
## terminator to end energy evolution when the energy difference
## between subsequent time steps is lower than 5e-6
##
## NOTE: can fail if the time step is small
##
###############################################
[./kill]
type = Terminator
expression = 'perc_change <= 1.0e-6'
[../]
[]
[Preconditioning]
###############################################
##
## Numerical preconditioning/solver options
##=
###############################################
[./smp]
type = SMP
full = true
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-ksp_gmres_restart -snes_atol -snes_rtol -ksp_rtol -pc_type -build_twosided'
petsc_options_value = ' 160 1e-8 1e-8 1e-5 bjacobi allreduce'
[../]
[]
[Executioner]
##########################################
##
## Time integration/solver options
##
##########################################
type = Transient
solve_type = 'PJFNK'
scheme = 'implicit-euler'
dtmin = 1e-13
dtmax = 0.6
l_max_its = 200
[./TimeStepper]
type = IterationAdaptiveDT
optimal_iterations = 6
growth_factor = 1.2
cutback_factor = 0.75
linear_iteration_ratio = 1000
dt = 0.6
[../]
verbose = true
nl_max_its = 20
num_steps = 5
[]
[Outputs]
###############################################
##===============
## Output options
##
###############################################
print_linear_residuals = false
perf_graph = false
[./out]
type = Exodus
file_base = out_pto_film_noP_test
elemental_as_nodal = true
interval = 1
[../]
[]