- C_ijklelastic stiffness tensor for material
C++ Type:std::vector<double>
Controllable:No
Description:elastic stiffness tensor for material
- Q_mnklelectrostrictive tensor for material
C++ Type:std::vector<double>
Controllable:No
Description:electrostrictive tensor for material
ComputeElectrostrictiveTensor
buildconstruction:Undocumented Class
The ComputeElectrostrictiveTensor has not been documented. The content listed below should be used as a starting point for documenting the class, which includes the typical automatic documentation associated with a MooseObject; however, what is contained is ultimately determined by what is necessary to make the documentation clear for users.
Compute an electrostrictive tensor.
Overview
Calculates the electrostrictive tensor components . Here,
with and are the elastic stiffness and electrostictive coefficients respectively. The tensor may be rotated via,
with the internal RotationTensor operation in MOOSE utils. The rotation operator accepts Euler angles in the standard Bunge sequence ().
Example Input File Syntax
Input Parameters
- base_nameOptional parameter that allows the user to define multiple mechanics material systems on the same block, i.e. for multiple phases
C++ Type:std::string
Controllable:No
Description:Optional parameter that allows the user to define multiple mechanics material systems on the same block, i.e. for multiple phases
- 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
- boundaryThe list of boundaries (ids or names) from the mesh where this object applies
C++ Type:std::vector<BoundaryName>
Controllable:No
Description:The list of boundaries (ids or names) from the mesh where this object applies
- computeTrueWhen false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
Default:True
C++ Type:bool
Controllable:No
Description:When false, MOOSE will not call compute methods on this material. The user must call computeProperties() after retrieving the MaterialBase via MaterialBasePropertyInterface::getMaterialBase(). Non-computed MaterialBases are not sorted for dependencies.
- compute_electrostrictive_coeffFalsecompute the electrostrictive coefficients Q_mnkl
Default:False
C++ Type:bool
Controllable:No
Description:compute the electrostrictive coefficients Q_mnkl
- constant_onNONEWhen ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
Default:NONE
C++ Type:MooseEnum
Options:NONE, ELEMENT, SUBDOMAIN
Controllable:No
Description:When ELEMENT, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps.When SUBDOMAIN, MOOSE will only call computeQpProperties() for the 0th quadrature point, and then copy that value to the other qps. Evaluations on element qps will be skipped
- declare_suffixAn optional suffix parameter that can be appended to any declared 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 declared properties. The suffix will be prepended with a '_' character.
- euler_angle_10Euler angle in direction 1
Default:0
C++ Type:double
Controllable:No
Description:Euler angle in direction 1
- euler_angle_20Euler angle in direction 2
Default:0
C++ Type:double
Controllable:No
Description:Euler angle in direction 2
- euler_angle_30Euler angle in direction 3
Default:0
C++ Type:double
Controllable:No
Description:Euler angle in direction 3
- fill_methodsymmetric9The fill method
Default:symmetric9
C++ Type:MooseEnum
Options:antisymmetric, symmetric9, symmetric21, general_isotropic, symmetric_isotropic, symmetric_isotropic_E_nu, antisymmetric_isotropic, axisymmetric_rz, general, principal, orthotropic
Controllable:No
Description:The fill method
- 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
- 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.
- 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
- 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
- output_propertiesList of material properties, from this material, to output (outputs must also be defined to an output type)
C++ Type:std::vector<std::string>
Controllable:No
Description:List of material properties, from this material, to output (outputs must also be defined to an output type)
- outputsnone Vector of output names where you would like to restrict the output of variables(s) associated with this object
Default:none
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
Outputs Parameters
Input Files
- (examples/domain_walls/BTO_90wall_T298K.i)
- (test/tests/electrooptics/BTO_monodomain_T298K_REF.i)
- (examples/domain_walls/BTO_wall_T298K.i)
- (test/tests/domain_wall/test_BTO_domain_wall.i)
- (examples/other/PZT_nanowire_test.i)
- (tutorial/ferroelectric_domain_wall.i)
- (test/tests/ics/PTOtest_fluc.i)
- (test/tests/coupled/PTOtest_3D.i)
- (test/tests/coupled/PTOtest_2D.i)
- (test/tests/electrooptics/BTO_monodomain_T298K_REFnoEO.i)
- (test/tests/pbc/pbc.i)
- (test/tests/auxkernels/microforce.i)
(examples/domain_walls/BTO_90wall_T298K.i)
[Mesh]
[gen]
############################################
##
## Type and dimension of the mesh
##
############################################
type = GeneratedMeshGenerator
dim = 3
nx = 60
ny = 60
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 = -30.0
ymax = 30.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.10471975511965977*((x+5)+(y+5)))+0.1 # 0.10471975511965977 = pi/(xmax-xmin)
[../]
[./stripe2]
type = ParsedFunction
value = -0.1*cos(0.10471975511965977*((x+5)+(y+5)))+0.1
[../]
[]
[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 = FunctionIC
function = 'stripe1'
[../]
[../]
[./polar_y]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = FunctionIC
function = 'stripe2'
[../]
[../]
[./polar_z]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.01e-4
max = 0.01e-4
[../]
[../]
[./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
[../]
[./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]
[./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'
[../]
[./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-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!
## 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_90wall_T298K_2
elemental_as_nodal = true
[../]
[./console]
type = Console
time_precision = 3
[../]
[]
(test/tests/electrooptics/BTO_monodomain_T298K_REF.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 = 2
ny = 2
nz = 2
#############################################
##
## Actual spatial coordinates of mesh.
## Jmax - Jmin = nJ/2 for J = x, y, z
## Units are in nanometers
##
#############################################
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
#############################################
##
## 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.5 -0.5 -0.5'
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
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.01*cos(0.08975979010256552*(x)) # 0.08975979010256552 = 2pi/(xmax-xmin)
[../]
[]
[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 = RandomIC
min = 0.01
max = 0.02
[../]
[../]
[./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
##
######################################
[./stress_xx_elastic]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yy_elastic]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xy_elastic]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xz_elastic]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_zz_elastic]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yz_elastic]
order = CONSTANT
family = MONOMIAL
[../]
[./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
[../]
######################################
##
## Principle change in refractive and
## refractive index components
##
######################################
[./dpn_1]
order = CONSTANT
family = MONOMIAL
[../]
[./dpn_2]
order = CONSTANT
family = MONOMIAL
[../]
[./dpn_3]
order = CONSTANT
family = MONOMIAL
[../]
[./den_1]
order = CONSTANT
family = MONOMIAL
[../]
[./den_2]
order = CONSTANT
family = MONOMIAL
[../]
[./den_3]
order = CONSTANT
family = MONOMIAL
[../]
[./n_1]
order = CONSTANT
family = MONOMIAL
[../]
[./n_2]
order = CONSTANT
family = MONOMIAL
[../]
[./n_3]
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
[../]
[./dn_p1]
type = PolarOpticChangeInRefractiveIndex
component = 0
variable = dpn_1
execute_on = 'timestep_end'
[../]
[./dn_p2]
type = PolarOpticChangeInRefractiveIndex
component = 1
variable = dpn_2
execute_on = 'timestep_end'
[../]
[./dn_p3]
type = PolarOpticChangeInRefractiveIndex
component = 2
variable = dpn_3
execute_on = 'timestep_end'
[../]
[./dn_e1]
type = ElastoChangeInRefractiveIndex
component = 0
variable = den_1
u_x = u_x
u_y = u_y
u_z = u_z
execute_on = 'timestep_end'
[../]
[./dn_e2]
type = ElastoChangeInRefractiveIndex
component = 1
variable = den_2
u_x = u_x
u_y = u_y
u_z = u_z
execute_on = 'timestep_end'
[../]
[./dn_e3]
type = ElastoChangeInRefractiveIndex
component = 2
variable = den_3
u_x = u_x
u_y = u_y
u_z = u_z
execute_on = 'timestep_end'
[../]
[./n_1_c]
type = ReworkedRefractiveIndex
variable = n_1
elasto = true
electro = false
polar = true
component = 0
var1 = dpn_1
var2 = den_1
execute_on = 'timestep_end'
[../]
[./n_2_c]
type = ReworkedRefractiveIndex
variable = n_2
elasto = true
electro = false
polar = true
component = 1
var1 = dpn_2
var2 = den_2
execute_on = 'timestep_end'
[../]
[./n_3_c]
type = ReworkedRefractiveIndex
variable = n_3
elasto = true
electro = false
polar = true
component = 2
var1 = dpn_3
var2 = den_3
execute_on = 'timestep_end'
[../]
[]
[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'
[../]
[./ref_index]
############################################
##
## Cubic refractive index BTO
## (need correct value)
##
############################################
type = GenericConstantMaterial
prop_names = 'n1 n2 n3 n4 n5 n6'
prop_values = '2.4272 2.4272 2.4272 0.0 0.0 0.0'
[../]
[./po_tensor]
############################################
##
## PO tensor coefficients from
## Bernasconi and Günter
## J. Appl. Phys., 78, 2651 (1995)
##
############################################
type = GenericConstantMaterial
prop_names = 'g1111 g1122 g1212'
prop_values = '0.15 0.038 0.07'
[../]
[./eo_tensor]
############################################
##
## Elastoptic tensor coefficients from
## Bernasconi and Günter
## J. Appl. Phys., 78, 2651 (1995)
##
############################################
type = GenericConstantMaterial
prop_names = 'p1111 p1122 p1212'
prop_values = '0.37 0.11 -0.30'
[../]
[]
[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
###############################################
##
## 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)
##
###############################################
permittivity = 0.08854187
[../]
[./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.
##
###############################################
[./Fbulk]
type = BulkEnergyEighth
execute_on = 'initial timestep_end'
[../]
[./Fwall]
type = WallEnergy
execute_on = 'initial timestep_end'
[../]
[./Felastic]
type = ElasticEnergy
execute_on = 'initial timestep_end'
use_displaced_mesh = false
[../]
[./Fcoupled]
type = ElectrostrictiveCouplingEnergy
execute_on = 'initial timestep_end'
[../]
[./Felec]
type = ElectrostaticEnergy
execute_on = 'initial timestep_end'
[../]
[./Ftotal]
type = LinearCombinationPostprocessor
pp_names = 'Fbulk Fwall Fcoupled Felec'
pp_coefs = ' 1 1 1 1'
execute_on = 'initial timestep_end'
[../]
[./perc_change]
type = PercentChangePostprocessor
postprocessor = Ftotal
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
num_steps = 5
[]
[Outputs]
###############################################
##
## Output options
##
###############################################
print_linear_residuals = false
perf_graph = false
[./out]
type = Exodus
file_base = test_bire
elemental_as_nodal = true
[../]
[]
(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
[../]
[]
(test/tests/domain_wall/test_BTO_domain_wall.i)
############################################
##
## to use this file, scale up size and add
## the functionIC to the polar field for the
## appropriate wall
## current Nov 2020 tests benchmark this
## within 10% of FFT calculations using the
## same functional
##
############################################
[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 = 2
ny = 2
nz = 2
#############################################
##
## Actual spatial coordinates of mesh.
## Jmax - Jmin = nJ/2 for J = x, y, z
## Units are in nanometers
##
#############################################
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
#############################################
##
## 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.5 -0.5 -0.5'
new_boundary = 100
[../]
[]
[GlobalParams]
len_scale = 1.0
############################################
##
## Gradient coefficients from
## Marton and Hlinka
## Phys. Rev. B. 74, 104014, (2006)
##
############################################
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
##
##############################################
u_x = u_x
u_y = u_y
u_z = u_z
[]
[Functions]
[./stripe1]
type = ParsedFunction
value = 0.01*cos(0.08975979010256552*(x))
[../]
[]
[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 = RandomIC
min = -0.01e-4
max = 0.01e-4
[../]
[../]
[./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)
##
######################################
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[./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
[../]
[]
[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
[../]
[]
[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)
##
## NOTE1: there might be some Legendre transforms
## depending on what approach you use
## -i.e. inhomogeneous strain vs
## homogeneous strain [renormalized]
##
## NOTE2: Sign convention in Ferret for the
## electrostrictive coeff. is multiplied by
## an overall factor of (-1)
##
##################################################
[./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'
[../]
[./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'
[../]
[./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
[../]
[./wall2ed_x]
type = Wall2EnergyDerivative
variable = polar_x
component = 0
[../]
[./wall2ed_y]
type = Wall2EnergyDerivative
variable = polar_y
component = 1
[../]
[./wall2ed_z]
type = Wall2EnergyDerivative
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.
##
###############################################
[./Fbulk]
type = BulkEnergyEighth
execute_on = 'initial timestep_end'
[../]
[./Fwall]
type = WallEnergy
execute_on = 'initial timestep_end'
[../]
[./Felastic]
type = ElasticEnergy
execute_on = 'initial timestep_end'
use_displaced_mesh = false
[../]
[./Fcoupled]
type = ElectrostrictiveCouplingEnergy
execute_on = 'initial timestep_end'
[../]
[./Felec]
type = ElectrostaticEnergy
execute_on = 'initial timestep_end'
[../]
[./Ftotal]
type = LinearCombinationPostprocessor
pp_names = 'Fbulk Fwall Fcoupled Felec'
pp_coefs = ' 1 1 1 1'
execute_on = 'initial timestep_end'
[../]
[./perc_change]
type = PercentChangePostprocessor
postprocessor = Ftotal
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-8 1e-8 1e-8 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
num_steps = 4
[]
[Outputs]
###############################################
##
## Output options
##
###############################################
print_linear_residuals = false
perf_graph = false
[./out]
type = Exodus
file_base = test
elemental_as_nodal = true
[../]
[]
(examples/other/PZT_nanowire_test.i)
[Mesh]
file = exodus_nanowire_rad20_h60.e
[]
[GlobalParams]
len_scale = 1.0
alpha1 = -0.14937
alpha11 = -0.0305
alpha111 = 0.2475
alpha12 = 0.632
alpha112 = 0.9684
alpha123 = -4.901
G110 = 0.173
G11_G110 = 2.0
G12_G110 = 0
G44_G110 = 1.0
G44P_G110 = 1.0
polar_x = polar_x
polar_y = polar_y
polar_z = polar_z
potential_int = potential_int
disp_x = disp_x
disp_y = disp_y
disp_z = disp_z
displacements = 'disp_x disp_y disp_z'
prefactor = 0.00 #negative = tension, positive = compression
[]
[Variables]
[./polar_x]
block = '1'
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.5e-5
max = 0.5e-5
seed = 1
[../]
[../]
[./polar_y]
block = '1'
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.5e-5
max = 0.5e-5
seed = 1
[../]
[../]
[./polar_z]
block = '1'
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.5e-5
max = 0.5e-5
seed = 1
[../]
[../]
[./potential_int]
order = FIRST
family = LAGRANGE
block = '1 2'
[./InitialCondition]
type = RandomIC
min = -0.5e-5
max = 0.5e-5
seed = 1
[../]
[../]
[./disp_x]
order = FIRST
family = LAGRANGE
block = '1 2'
[../]
[./disp_y]
order = FIRST
family = LAGRANGE
block = '1 2'
[../]
[./disp_z]
order = FIRST
family = LAGRANGE
block = '1 2'
[../]
[]
[AuxVariables]
[./stress_xx_elastic]
order = CONSTANT
family = MONOMIAL
#initial_from_file_var = stress_xx
[../]
[./stress_yy_elastic]
order = CONSTANT
family = MONOMIAL
#initial_from_file_var = stress_yy
[../]
[./stress_xy_elastic]
order = CONSTANT
family = MONOMIAL
#initial_from_file_var = stress_xy
[../]
[./stress_xz_elastic]
order = CONSTANT
family = MONOMIAL
#initial_from_file_var = stress_xz
[../]
[./stress_zz_elastic]
order = CONSTANT
family = MONOMIAL
#initial_from_file_var = stress_zz
[../]
[./stress_yz_elastic]
order = CONSTANT
family = MONOMIAL
#initial_from_file_var = stress_yz
[../]
[./strain_xx_elastic]
order = CONSTANT
family = MONOMIAL
#initial_from_file_var = strain_xx
[../]
[./strain_yy_elastic]
order = CONSTANT
family = MONOMIAL
#initial_from_file_var = strain_yy
[../]
[./strain_xy_elastic]
order = CONSTANT
family = MONOMIAL
#initial_from_file_var = strain_xy
[../]
[./strain_xz_elastic]
order = CONSTANT
family = MONOMIAL
#initial_from_file_var = strain_xz
[../]
[./strain_zz_elastic]
order = CONSTANT
family = MONOMIAL
#initial_from_file_var = strain_zz
[../]
[./strain_yz_elastic]
order = CONSTANT
family = MONOMIAL
#initial_from_file_var = strain_yz
[../]
[./chern]
order = CONSTANT
family = MONOMIAL
[../]
[./chernMag]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./cherndens]
type = ChernSimonsDensity
variable = chern
[../]
[./cherndensMag]
type = ChernSimonsDensityMag
block = '1'
variable = chernMag
[../]
[./matl_e11]
type = RankTwoAux
rank_two_tensor = elastic_strain
index_i = 0
index_j = 0
variable = strain_xx_elastic
execute_on = 'timestep_end'
[../]
[./matl_e12]
type = RankTwoAux
rank_two_tensor = elastic_strain
index_i = 0
index_j = 1
variable = strain_xy_elastic
execute_on = 'timestep_end'
[../]
[./matl_e13]
type = RankTwoAux
rank_two_tensor = elastic_strain
index_i = 0
index_j = 2
variable = strain_xz_elastic
execute_on = 'timestep_end'
[../]
[./matl_e22]
type = RankTwoAux
rank_two_tensor = elastic_strain
index_i = 1
index_j = 1
variable = strain_yy_elastic
execute_on = 'timestep_end'
[../]
[./matl_e23]
type = RankTwoAux
rank_two_tensor = elastic_strain
index_i = 1
index_j = 2
variable = strain_yz_elastic
execute_on = 'timestep_end'
[../]
[./matl_e33]
type = RankTwoAux
rank_two_tensor = elastic_strain
index_i = 2
index_j = 2
variable = strain_zz_elastic
execute_on = 'timestep_end'
[../]
[./matl_s11]
type = RankTwoAux
rank_two_tensor = stress
index_i = 0
index_j = 0
variable = stress_xx_elastic
execute_on = 'timestep_end'
[../]
[./matl_s12]
type = RankTwoAux
rank_two_tensor = stress
index_i = 0
index_j = 1
variable = stress_xy_elastic
execute_on = 'timestep_end'
[../]
[./matl_s13]
type = RankTwoAux
rank_two_tensor = stress
index_i = 0
index_j = 2
variable = stress_xz_elastic
execute_on = 'timestep_end'
[../]
[./matl_s22]
type = RankTwoAux
rank_two_tensor = stress
index_i = 1
index_j = 1
variable = stress_yy_elastic
execute_on = 'timestep_end'
[../]
[./matl_s23]
type = RankTwoAux
rank_two_tensor = stress
index_i = 1
index_j = 2
variable = stress_yz_elastic
execute_on = 'timestep_end'
[../]
[./matl_s33]
type = RankTwoAux
rank_two_tensor = stress
index_i = 2
index_j = 2
variable = stress_zz_elastic
execute_on = 'timestep_end'
[../]
[]
[Materials]
[./eigen_strain_zz] #Use for stress-free strain (ie epitaxial)
type = ComputeEigenstrain
block = '1 2'
# eigen_base = 'exx exy exz eyx eyy eyz ezx ezy ezz'
eigen_base = '1 0 0 0 1 0 0 0 0'
eigenstrain_name = eigenstrain
[../]
[./elasticity_tensor_1]
type = ComputeElasticityTensor
fill_method = symmetric9
#from MaterialsProject
C_ijkl = '281 115.74 115.74 281 115.74 281 97.18 97.18 97.18'
block = '1'
[../]
[./strain_1]
type = ComputeSmallStrain
block = '1'
eigenstrain_names = eigenstrain
[../]
[./stress_1]
type = ComputeLinearElasticStress
block = '1'
[../]
[./slab_ferroelectric]
block = '1'
type = ComputeElectrostrictiveTensor
Q_mnkl = '0.0842 -0.02446 -0.02446 0.0842 -0.02446 0.0842 0.06417 0.06417 0.06417'
C_ijkl = '281 115.74 115.74 281 115.74 281 97.18 97.18 97.18'
[../]
[./elasticity_tensor_2]
type = ComputeElasticityTensor
#averaged from BulkMod
C_ijkl = '319 99.6 99.6 319 99.6 319 109.53 109.53 109.53'
fill_method = symmetric9
block = '2'
[../]
[./strain_2]
type = ComputeSmallStrain
block = '2'
eigenstrain_names = eigenstrain
[../]
[./stress_2]
type = ComputeLinearElasticStress
block = '2'
[../]
[]
[Kernels]
#Elastic problem
[./TensorMechanics]
#This is an action block
[../]
#Bulk energy density
[./bed_x]
type = BulkEnergyDerivativeSixth
variable = polar_x
component = 0
[../]
[./bed_y]
type = BulkEnergyDerivativeSixth
variable = polar_y
component = 1
[../]
[./bed_z]
type = BulkEnergyDerivativeSixth
variable = polar_z
component = 2
[../]
##Wall energy penalty
[./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
[../]
##Polarization-strain coupling
[./ferroelectriccouplingp_xx]
type = FerroelectricCouplingP
variable = polar_x
component = 0
[../]
[./ferroelectriccouplingp_yy]
type = FerroelectricCouplingP
variable = polar_y
component = 1
[../]
[./ferroelectriccouplingp_zz]
type = FerroelectricCouplingP
variable = polar_z
component = 2
[../]
[./ferroelectriccouplingX_xx]
type = FerroelectricCouplingX
block = '1'
variable = disp_x
component = 0
[../]
[./ferroelectriccouplingX_yy]
type = FerroelectricCouplingX
block = '1'
variable = disp_y
component = 1
[../]
[./ferroelectriccouplingX_zz]
type = FerroelectricCouplingX
block = '1'
variable = disp_z
component = 2
[../]
##Electrostatics
[./polar_x_electric_E]
type = PolarElectricEStrong
variable = potential_int
permittivity = 0.08854187
[../]
[./FE_E_int]
type = Electrostatics
variable = potential_int
block = '1'
permittivity = 0.08854187
[../]
[./DIE_E_int]
type = Electrostatics
variable = potential_int
block = '2'
permittivity = 2.6562561
[../]
[./polar_electric_px]
type = PolarElectricPStrong
variable = polar_x
component = 0
[../]
[./polar_electric_py]
type = PolarElectricPStrong
variable = polar_y
component = 1
[../]
[./polar_electric_pz]
type = PolarElectricPStrong
variable = polar_z
component = 2
[../]
##Time dependence
[./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]
[./potential_int_1]
type = DirichletBC
variable = potential_int
boundary = '5'
value = -0.00001
[../]
[./potential_int_2]
type = DirichletBC
variable = potential_int
boundary = '6'
value = -0.00001
[../]
[./disp_x]
type = DirichletBC
variable = disp_x
boundary = '1 2 5 6'
value = 0
[../]
[./disp_y]
type = DirichletBC
variable = disp_y
boundary = '1 2 5 6'
value = 0
[../]
[./disp_z]
type = DirichletBC
variable = disp_z
boundary = '1 2 5 6'
value = 0
[../]
[]
[Postprocessors]
[./avgChern]
block = '1'
type = ElementAverageValue
variable = chern
[../]
[./avgchernMag]
block = '1'
type = ElementAverageValue
variable = chernMag
[../]
[./Fbulk]
type = BulkEnergy
block = '1'
execute_on = 'timestep_end'
[../]
[./Fwall]
type = WallEnergy
block = '1'
execute_on = 'timestep_end'
[../]
[./Felastic]
type = ElasticEnergy
block = '1 2'
execute_on = 'timestep_end'
[../]
[./Fcoupled]
block = '1'
type = CoupledEnergy
execute_on = 'timestep_end'
[../]
[./Felec]
block = '1'
type = ElectrostaticEnergy
permittivity = 0.08854187
execute_on = 'timestep_end'
[../]
[./Ftotal]
type = TotalEnergyFlow
Fbulk = Fbulk
Fwall = Fwall
Fcoupled = Fcoupled
Felec = Felec
execute_on = 'timestep_end'
[../]
[./perc_change]
type = PercentChangePostprocessor
postprocessor = Ftotal
[../]
[]
[UserObjects]
[./kill]
type = Terminator
expression = 'perc_change <= 2.0e-3'
[../]
[]
[Preconditioning]
[./smp]
type = SMP
full = true
petsc_options = '-snes_view -snes_linesearch_monitor -snes_converged_reason -ksp_converged_reason'
petsc_options_iname = '-ksp_gmres_restart -snes_atol -snes_rtol -ksp_rtol -pc_type'
petsc_options_value = ' 121 1e-10 1e-6 1e-6 bjacobi'
[../]
[]
[Executioner]
type = Transient
[./TimeStepper]
type = IterationAdaptiveDT
dt = 0.7
#iteration_window = 3
optimal_iterations = 6 #should be 5 probably
growth_factor = 1.4
linear_iteration_ratio = 1000
cutback_factor = 0.8
[../]
solve_type = 'NEWTON' #"PJFNK, JFNK, NEWTON"
scheme = 'implicit-euler' #"implicit-euler, explicit-euler, crank-nicolson, bdf2, rk-2"
dtmin = 1e-13
dtmax = 0.7
[]
[Outputs]
print_linear_residuals = true
print_perf_log = true
[./out]
type = Exodus
file_base = out_PZT_nanowire_inSTO_test
elemental_as_nodal = true
interval = 8
[../]
[./outcsv]
type = CSV
file_base = out_PZT_nanowire_inSTO_test
[../]
[]
(tutorial/ferroelectric_domain_wall.i)
[Mesh]
[gen]
############################################
##
## Type and dimension of the mesh
##
############################################
type = GeneratedMeshGenerator
dim = 3
nx = 480
ny = 1
nz = 1
xmin = -60.0
xmax = 60.0
ymin = -0.25
ymax = 0.25
zmin = -0.25
zmax = 0.25
#############################################
##
## 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 = '-60.0 -0.25 -0.25'
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
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.05235987755*(x)) #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
[../]
[]
[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
[../]
[./walled2_x]
type = Wall2EnergyDerivative
variable = polar_x
component = 0
[../]
[./walled2_y]
type = Wall2EnergyDerivative
variable = polar_y
component = 1
[../]
[./walled2_z]
type = Wall2EnergyDerivative
variable = polar_z
component = 2
[../]
[./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
[../]
[]
[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'
[../]
[]
[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'
[../]
[./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 1e-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
growth_factor = 1.2
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_ferroelectric_domain_wall
elemental_as_nodal = true
execute_on = 'timestep_end'
interval = 2
[../]
[]
(test/tests/ics/PTOtest_fluc.i)
[Mesh]
type = GeneratedMesh
dim = 3
nx = 10
ny = 10
nz = 7
xmin = -4
xmax = 4
ymin = -4
ymax = 4
zmin = -3
zmax = 3
elem_type = HEX8
[]
[GlobalParams]
len_scale = 1.0
alpha1 = -0.1722883
alpha11 = -0.07253
alpha111 = 0.26
alpha12 = 0.75
alpha112 = 0.61
alpha123 = -3.67
polar_x = polar_x
polar_y = polar_y
polar_z = polar_z
potential_E_int = potential_E_int
disp_x = disp_x
disp_y = disp_y
disp_z = disp_z
displacements = 'disp_x disp_y disp_z'
prefactor = 0.01 #negative = tension, positive = compression
[]
[Variables]
[./polar_x]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = FluctuationsIC
epsilon = 1.0e-5
q1 = '354 1005 645'
q2 = '715 1065 1132'
q3 = '391 305 1106'
q4 = '1053 1116 627'
h = 0.22
[../]
[../]
[./polar_y]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = FluctuationsIC
epsilon = 1.0e-5
q1 = '354 850 10'
q2 = '715 28 5'
q3 = '391 305 1106'
q4 = '653 1116 627'
h = 0.15
[../]
[../]
[./polar_z]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = FluctuationsIC
epsilon = 1.0e-5
q1 = '860 165 645'
q2 = '715 665 1332'
q3 = '361 15 706'
q4 = '253 1116 627'
h = 0.05
[../]
[../]
[./potential_E_int]
order = FIRST
family = LAGRANGE
[../]
[./disp_x]
order = FIRST
family = LAGRANGE
[../]
[./disp_y]
order = FIRST
family = LAGRANGE
[../]
[./disp_z]
order = FIRST
family = LAGRANGE
[../]
[]
[Materials]
[./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'
[../]
[./eigen_strain_zz] #Use for stress-free strain (ie epitaxial)
type = ComputeEigenstrain
block = '0'
# eigen_base = 'exx exy exz eyx eyy eyz ezx ezy ezz'
eigen_base = '1 0 0 0 1 0 0 0 0'
eigenstrain_name = eigenstrain
[../]
[./elasticity_tensor_1]
type = ComputeElasticityTensor
fill_method = symmetric9
C_ijkl = '380. 150. 150. 380. 150. 380. 110. 110. 110.'
[../]
[./strain_1]
type = ComputeSmallStrain
eigenstrain_names = eigenstrain
[../]
[./stress_1]
type = ComputeLinearElasticStress
[../]
[./slab_ferroelectric]
type = ComputeElectrostrictiveTensor
Q_mnkl = '-0.089 0.026 0.026 -0.089 0.026 -0.089 -0.03375 -0.03375 -0.03375'
C_ijkl = '380. 150. 150. 380. 150. 380. 110. 110. 110.'
[../]
[./permitivitty_1]
type = GenericConstantMaterial
prop_names = 'permittivity'
prop_values = '0.08854187'
[../]
[]
[Kernels]
[./TensorMechanics]
[../]
[./bed_x]
type = BulkEnergyDerivativeSixth
variable = polar_x
component = 0
[../]
[./bed_y]
type = BulkEnergyDerivativeSixth
variable = polar_y
component = 1
[../]
[./bed_z]
type = BulkEnergyDerivativeSixth
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
[../]
[./ferroelectriccouplingp_xx]
type = FerroelectricCouplingP
variable = polar_x
component = 0
[../]
[./ferroelectriccouplingp_yy]
type = FerroelectricCouplingP
variable = polar_y
component = 1
[../]
[./ferroelectriccouplingp_zz]
type = FerroelectricCouplingP
variable = polar_z
component = 2
[../]
[./ferroelectriccouplingX_xx]
type = FerroelectricCouplingX
variable = disp_x
component = 0
[../]
[./ferroelectriccouplingX_yy]
type = FerroelectricCouplingX
variable = disp_y
component = 1
[../]
[./ferroelectriccouplingX_zz]
type = FerroelectricCouplingX
variable = disp_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]
[./disp_x_front]
type = DirichletBC
boundary = 'front'
value = 0.0
variable = disp_x
[../]
[./disp_y_front]
type = DirichletBC
boundary = 'front'
value = 0.0
variable = disp_y
[../]
[./disp_z_front]
type = DirichletBC
boundary = 'front'
value = 0.0
variable = disp_z
[../]
[./potential_E_int_front]
type = DirichletBC
boundary = 'front'
value = 0.0001
variable = potential_E_int
[../]
[./potential_E_int_back]
type = DirichletBC
boundary = 'back'
value = 0.0001
variable = potential_E_int
[../]
[./disp_x_back]
type = DirichletBC
boundary = 'back'
value = 0.0
variable = disp_x
[../]
[./disp_y_back]
type = DirichletBC
boundary = 'back'
value = 0.0
variable = disp_y
[../]
[./disp_z_back]
type = DirichletBC
boundary = 'back'
value = 0.0
variable = disp_z
[../]
[]
[Postprocessors]
[./Fbulk]
type = BulkEnergy
execute_on = 'initial timestep_end'
[../]
[./Fwall]
type = WallEnergy
execute_on = 'initial timestep_end'
[../]
[./Felastic]
type = ElasticEnergy
execute_on = 'initial timestep_end'
[../]
[./Fcoupled]
type = ElectrostrictiveEnergy
execute_on = 'initial timestep_end'
[../]
[./Felec]
type = ElectrostaticEnergy
execute_on = 'initial timestep_end'
[../]
[./Ftotal]
type = LinearCombinationPostprocessor
pp_names = 'Fbulk Fwall Fcoupled Felec'
pp_coefs = ' 1 1 1 1'
execute_on = 'initial timestep_end'
[../]
[./perc_change]
type = PercentChangePostprocessor
postprocessor = Ftotal
execute_on = 'initial timestep_end'
[../]
[]
[Preconditioning]
[./smp]
type = SMP
full = true
petsc_options_iname = '-ksp_gmres_restart -snes_rtol -ksp_rtol -pc_type'
petsc_options_value = ' 121 1e-6 1e-8 bjacobi'
[../]
[]
[Executioner]
type = Transient
solve_type = 'NEWTON'
scheme = 'implicit-euler'
dtmin = 1e-13
dtmax = 0.1
num_steps = 5
[]
[Outputs]
print_linear_residuals = false
[./out]
type = Exodus
file_base = outPTO_test_fluct
elemental_as_nodal = true
interval = 1
[../]
[]
(test/tests/coupled/PTOtest_3D.i)
[Mesh]
type = GeneratedMesh
dim = 3
nx = 6
ny = 6
nz = 4
xmin = -1.5
xmax = 1.5
ymin = -1.5
ymax = 1.5
zmin = -1
zmax = 1
elem_type = HEX8
[]
[GlobalParams]
len_scale = 1.0
alpha1 = -0.1722883 # (3.766(T-765.1)*10^5) C^{-2} nm^2 (T = 293 K)
alpha11 = -0.07253
alpha111 = 0.26
alpha12 = 0.75
alpha112 = 0.61
alpha123 = -3.67
polar_x = polar_x
polar_y = polar_y
polar_z = polar_z
potential_E_int = potential_E_int
disp_x = disp_x
disp_y = disp_y
disp_z = disp_z
displacements = 'disp_x disp_y disp_z'
prefactor = 0.01 #negative = tension,positive = compression
[]
[Variables]
[./polar_x]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.01e-5
max = 0.01e-5
seed = 5
legacy_generator = true
[../]
[../]
[./polar_y]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.01e-5
max = 0.01e-5
seed = 5
legacy_generator = true
[../]
[../]
[./polar_z]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.01e-5
max = 0.01e-5
seed = 5
legacy_generator = true
[../]
[../]
[./potential_E_int]
order = FIRST
family = LAGRANGE
[../]
[./disp_x]
order = FIRST
family = LAGRANGE
[../]
[./disp_y]
order = FIRST
family = LAGRANGE
[../]
[./disp_z]
order = FIRST
family = LAGRANGE
[../]
[]
[Materials]
[./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'
[../]
[./eigen_strain_zz] #Use for stress-free strain (ie epitaxial)
type = ComputeEigenstrain
block = '0'
# eigen_base = 'exx exy exz eyx eyy eyz ezx ezy ezz'
eigen_base = '1 0 0 0 1 0 0 0 0'
eigenstrain_name = eigenstrain
[../]
[./elasticity_tensor_1]
type = ComputeElasticityTensor
fill_method = symmetric9
C_ijkl = '380. 150. 150. 380. 150. 380. 110. 110. 110.'
[../]
[./strain_1]
type = ComputeSmallStrain
eigenstrain_names = eigenstrain
[../]
[./stress_1]
type = ComputeLinearElasticStress
[../]
[./slab_ferroelectric]
type = ComputeElectrostrictiveTensor
Q_mnkl = '-0.089 0.026 0.026 -0.089 0.026 -0.089 -0.03375 -0.03375 -0.03375'
C_ijkl = '380. 150. 150. 380. 150. 380. 110. 110. 110.'
[../]
[./permitivitty_1]
type = GenericConstantMaterial
prop_names = 'permittivity'
prop_values = '0.08854187'
[../]
[]
[Kernels]
#Elastic problem
[./TensorMechanics]
[../]
[./bed_x]
type = BulkEnergyDerivativeSixth
variable = polar_x
component = 0
[../]
[./bed_y]
type = BulkEnergyDerivativeSixth
variable = polar_y
component = 1
[../]
[./bed_z]
type = BulkEnergyDerivativeSixth
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
[../]
[./ferroelectriccouplingp_xx]
type = FerroelectricCouplingP
variable = polar_x
component = 0
[../]
[./ferroelectriccouplingp_yy]
type = FerroelectricCouplingP
variable = polar_y
component = 1
[../]
[./ferroelectriccouplingp_zz]
type = FerroelectricCouplingP
variable = polar_z
component = 2
[../]
[./ferroelectriccouplingX_xx]
type = FerroelectricCouplingX
variable = disp_x
component = 0
[../]
[./ferroelectriccouplingX_yy]
type = FerroelectricCouplingX
variable = disp_y
component = 1
[../]
[./ferroelectriccouplingX_zz]
type = FerroelectricCouplingX
variable = disp_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]
[./potential_E_int_front]
type = DirichletBC
boundary = 'front'
value = 0.0001
variable = potential_E_int
[../]
[./potential_E_int_back]
type = DirichletBC
boundary = 'back'
value = 0.0001
variable = potential_E_int
[../]
[]
[Postprocessors]
[./Fbulk]
type = BulkEnergy
execute_on = 'initial timestep_end'
[../]
[./Fwall]
type = WallEnergy
execute_on = 'initial timestep_end'
[../]
[./Felastic]
type = ElasticEnergy
execute_on = 'initial timestep_end'
[../]
[./Fcoupled]
type = ElectrostrictiveEnergy
execute_on = 'initial timestep_end'
[../]
[./Felec]
type = ElectrostaticEnergy
execute_on = 'initial timestep_end'
[../]
[./Ftotal]
type = LinearCombinationPostprocessor
pp_names = 'Fbulk Fwall Fcoupled Felec'
pp_coefs = ' 1 1 1 1'
execute_on = 'initial timestep_end'
[../]
[./perc_change]
type = PercentChangePostprocessor
postprocessor = Ftotal
execute_on = 'initial timestep_end'
[../]
[]
[Preconditioning]
[./smp]
type = SMP
full = true
petsc_options_iname = '-ksp_gmres_restart -snes_atol -snes_rtol -ksp_rtol -pc_type'
petsc_options_value = ' 121 1e-10 1e-8 1e-6 bjacobi'
[../]
[]
[Executioner]
type = Transient
solve_type = 'NEWTON'
scheme = 'implicit-euler'
dtmin = 1e-13
dtmax = 0.35
num_steps = 15
[]
[Outputs]
print_linear_residuals = false
[./out]
type = Exodus
file_base = outPTOchunk_test
elemental_as_nodal = true
interval = 1
[../]
[]
(test/tests/coupled/PTOtest_2D.i)
[Mesh]
type = GeneratedMesh
dim = 2
nx = 32
ny = 32
xmin = -8
xmax = 8
ymin = -8
ymax = 8
[]
[GlobalParams]
len_scale = 1.0
alpha1 = -0.1722883 # (3.766(T-765.1)*10^5) C^{-2} nm^2 (T = 293 K)
alpha11 = -0.07253
alpha111 = 0.26
alpha12 = 0.75
alpha112 = 0.61
alpha123 = -3.67
polar_x = polar_x
polar_y = polar_y
potential_E_int = potential_E_int
disp_x = disp_x
disp_y = disp_y
displacements = 'disp_x disp_y'
prefactor = 0.0
[]
[Variables]
[./polar_x]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.01e-4
max = 0.01e-4
seed = 5
legacy_generator = true
[../]
[../]
[./polar_y]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.01e-4
max = 0.01e-4
seed = 5
legacy_generator = true
[../]
[../]
[./potential_E_int]
order = FIRST
family = LAGRANGE
[../]
[./disp_x]
order = FIRST
family = LAGRANGE
[../]
[./disp_y]
order = FIRST
family = LAGRANGE
[../]
[]
[Materials]
[./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'
[../]
[./eigen_strain_zz] #Use for stress-free strain (ie epitaxial)
type = ComputeEigenstrain
block = '0'
# eigen_base = 'exx exy exz eyx eyy eyz ezx ezy ezz'
eigen_base = '1 0 0 0 1 0 0 0 0'
eigenstrain_name = eigenstrain
[../]
[./elasticity_tensor_1]
type = ComputeElasticityTensor
fill_method = symmetric9
C_ijkl = '380. 150. 150. 380. 150. 380. 110. 110. 110.'
[../]
[./strain_1]
type = ComputeSmallStrain
eigenstrain_names = eigenstrain
[../]
[./stress_1]
type = ComputeLinearElasticStress
[../]
[./slab_ferroelectric]
type = ComputeElectrostrictiveTensor
Q_mnkl = '-0.089 0.026 0.026 -0.089 0.026 -0.089 -0.03375 -0.03375 -0.03375'
C_ijkl = '380. 150. 150. 380. 150. 380. 110. 110. 110.'
[../]
[./permitivitty_1]
type = GenericConstantMaterial
prop_names = 'permittivity'
prop_values = '0.08854187'
[../]
[]
[Kernels]
#Elastic problem
[./TensorMechanics]
[../]
[./bed_x]
type = BulkEnergyDerivativeSixth
variable = polar_x
component = 0
[../]
[./bed_y]
type = BulkEnergyDerivativeSixth
variable = polar_y
component = 1
[../]
[./walled_x]
type = WallEnergyDerivative
variable = polar_x
component = 0
[../]
[./walled_y]
type = WallEnergyDerivative
variable = polar_y
component = 1
[../]
[./ferroelectriccouplingp_xx]
type = FerroelectricCouplingP
variable = polar_x
component = 0
[../]
[./ferroelectriccouplingp_yy]
type = FerroelectricCouplingP
variable = polar_y
component = 1
[../]
[./ferroelectriccouplingX_xx]
type = FerroelectricCouplingX
variable = disp_x
component = 0
[../]
[./ferroelectriccouplingX_yy]
type = FerroelectricCouplingX
variable = disp_y
component = 1
[../]
[./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_x_time]
type = TimeDerivativeScaled
variable=polar_x
time_scale = 1.0
[../]
[./polar_y_time]
type = TimeDerivativeScaled
variable=polar_y
time_scale = 1.0
[../]
[]
[BCs]
[./potential_E_int_front]
type = DirichletBC
boundary = 'top'
value = 0.0001
variable = potential_E_int
[../]
[./potential_E_int_back]
type = DirichletBC
boundary = 'bottom'
value = 0.0001
variable = potential_E_int
[../]
[./disp_x_back]
type = DirichletBC
boundary = 'top'
value = 0.0
variable = disp_x
[../]
[./disp_y_back]
type = DirichletBC
boundary = 'top'
value = 0.0
variable = disp_y
[../]
[]
[Postprocessors]
[./Fbulk]
type = BulkEnergy
execute_on = 'initial timestep_end'
[../]
[./Fwall]
type = WallEnergy
execute_on = 'initial timestep_end'
[../]
[./Felastic]
type = ElasticEnergy
execute_on = 'initial timestep_end'
[../]
[./Fcoupled]
type = ElectrostrictiveEnergy
execute_on = 'initial timestep_end'
[../]
[./Felec]
type = ElectrostaticEnergy
execute_on = 'initial timestep_end'
[../]
[./Ftotal]
type = LinearCombinationPostprocessor
pp_names = 'Fbulk Fwall Fcoupled Felec'
pp_coefs = ' 1 1 1 1'
execute_on = 'initial timestep_end'
[../]
[./perc_change]
type = PercentChangePostprocessor
postprocessor = Ftotal
execute_on = 'initial timestep_end'
[../]
[]
[Preconditioning]
[./smp]
type = SMP
full = true
petsc_options_iname = '-ksp_gmres_restart -snes_atol -snes_rtol -ksp_rtol -pc_type'
petsc_options_value = ' 121 1e-10 1e-8 1e-6 bjacobi'
[../]
[]
[Executioner]
type = Transient
solve_type = 'NEWTON'
scheme = 'implicit-euler'
#dt = 0.5
dtmin = 1e-13
dtmax = 0.1
num_steps = 10
[]
[Outputs]
print_linear_residuals = false
[./out]
type = Exodus
file_base = outPTOchunk_test_2D
elemental_as_nodal = true
interval = 1
[../]
[]
(test/tests/electrooptics/BTO_monodomain_T298K_REFnoEO.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 = 2
ny = 2
nz = 2
#############################################
##
## Actual spatial coordinates of mesh.
## Jmax - Jmin = nJ/2 for J = x, y, z
## Units are in nanometers
##
#############################################
xmin = -0.5
xmax = 0.5
ymin = -0.5
ymax = 0.5
zmin = -0.5
zmax = 0.5
#############################################
##
## 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.5 -0.5 -0.5'
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
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.01*cos(0.08975979010256552*(x)) # 0.08975979010256552 = 2pi/(xmax-xmin)
[../]
[]
[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.01
max = 0.02
[../]
[../]
[./polar_z]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.01e-4
max = 0.01e-4
[../]
[../]
[./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
##
######################################
[./stress_xx_elastic]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yy_elastic]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xy_elastic]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_xz_elastic]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_zz_elastic]
order = CONSTANT
family = MONOMIAL
[../]
[./stress_yz_elastic]
order = CONSTANT
family = MONOMIAL
[../]
[./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
[../]
######################################
##
## Principle change in refractive and
## refractive index components
##
######################################
[./dpn_1]
order = CONSTANT
family = MONOMIAL
[../]
[./dpn_2]
order = CONSTANT
family = MONOMIAL
[../]
[./dpn_3]
order = CONSTANT
family = MONOMIAL
[../]
[./den_1]
order = CONSTANT
family = MONOMIAL
[../]
[./den_2]
order = CONSTANT
family = MONOMIAL
[../]
[./den_3]
order = CONSTANT
family = MONOMIAL
[../]
[./n_1]
order = CONSTANT
family = MONOMIAL
[../]
[./n_2]
order = CONSTANT
family = MONOMIAL
[../]
[./n_3]
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
[../]
[./dn_p1]
type = PolarOpticChangeInRefractiveIndex
component = 0
variable = dpn_1
execute_on = 'timestep_end'
[../]
[./dn_p2]
type = PolarOpticChangeInRefractiveIndex
component = 1
variable = dpn_2
execute_on = 'timestep_end'
[../]
[./dn_p3]
type = PolarOpticChangeInRefractiveIndex
component = 2
variable = dpn_3
execute_on = 'timestep_end'
[../]
[./dn_e1]
type = ElastoChangeInRefractiveIndex
component = 0
variable = den_1
u_x = u_x
u_y = u_y
u_z = u_z
execute_on = 'timestep_end'
[../]
[./dn_e2]
type = ElastoChangeInRefractiveIndex
component = 1
variable = den_2
u_x = u_x
u_y = u_y
u_z = u_z
execute_on = 'timestep_end'
[../]
[./dn_e3]
type = ElastoChangeInRefractiveIndex
component = 2
variable = den_3
u_x = u_x
u_y = u_y
u_z = u_z
execute_on = 'timestep_end'
[../]
[./n_1_c]
type = ReworkedRefractiveIndex
variable = n_1
elasto = false
electro = false
polar = true
component = 0
var1 = dpn_1
execute_on = 'timestep_end'
[../]
[./n_2_c]
type = ReworkedRefractiveIndex
variable = n_2
elasto = false
electro = false
polar = true
component = 1
var1 = dpn_2
execute_on = 'timestep_end'
[../]
[./n_3_c]
type = ReworkedRefractiveIndex
variable = n_3
elasto = false
electro = false
polar = true
component = 2
var1 = dpn_3
execute_on = 'timestep_end'
[../]
[]
[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'
[../]
[./ref_index]
############################################
##
## Cubic refractive index BTO
## (need correct value)
##
############################################
type = GenericConstantMaterial
prop_names = 'n1 n2 n3 n4 n5 n6'
prop_values = '2.4272 2.4272 2.4272 0.0 0.0 0.0'
[../]
[./po_tensor]
############################################
##
## PO tensor coefficients from
## Bernasconi and Günter
## J. Appl. Phys., 78, 2651 (1995)
##
############################################
type = GenericConstantMaterial
prop_names = 'g1111 g1122 g1212'
prop_values = '0.15 0.038 0.07'
[../]
[./eo_tensor]
############################################
##
## Elastoptic tensor coefficients from
## Bernasconi and Günter
## J. Appl. Phys., 78, 2651 (1995)
##
############################################
type = GenericConstantMaterial
prop_names = 'p1111 p1122 p1212'
prop_values = '0.37 0.11 -0.30'
[../]
[./eps]
###############################################
##
## 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'
[../]
[../]
[./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.
##
###############################################
[./Fbulk]
type = BulkEnergyEighth
execute_on = 'initial timestep_end'
[../]
[./Fwall]
type = WallEnergy
execute_on = 'initial timestep_end'
[../]
[./Felastic]
type = ElasticEnergy
execute_on = 'initial timestep_end'
use_displaced_mesh = false
[../]
[./Fcoupled]
type = ElectrostrictiveCouplingEnergy
execute_on = 'initial timestep_end'
[../]
[./Felec]
type = ElectrostaticEnergy
execute_on = 'initial timestep_end'
[../]
[./Ftotal]
type = LinearCombinationPostprocessor
pp_names = 'Fbulk Fwall Fcoupled Felec'
pp_coefs = ' 1 1 1 1'
execute_on = 'initial timestep_end'
[../]
[./perc_change]
type = PercentChangePostprocessor
postprocessor = Ftotal
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
num_steps = 5
[]
[Outputs]
###############################################
##
## Output options
##
###############################################
print_linear_residuals = false
perf_graph = false
[./out]
type = Exodus
file_base = test_bire_no_EO
elemental_as_nodal = true
[../]
[]
(test/tests/pbc/pbc.i)
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 3
nx = 6
ny = 6
nz = 2
xmin = -6.0
xmax = 6.0
ymin = -6.0
ymax = 6.0
zmin = -2.0
zmax = 2.0
elem_type = HEX8
[]
[./cnode]
input = gen
type = ExtraNodesetGenerator
coord = '-6.0 -6.0 -2.0'
new_boundary = 100
[../]
# additional boundary sideset (one node) to zero one of the elastic displacement vectors - eliminates rigid body translations from the degrees of freedom
[]
[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'
u_x = u_x
u_y = u_y
u_z = u_z
[]
[Variables]
[./global_strain]
order = SIXTH
family = SCALAR
[../]
[./polar_x]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.01e-4
max = 0.01e-4
seed = 6
[../]
[../]
[./polar_y]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.01e-4
max = 0.01e-4
seed = 6
[../]
[../]
[./polar_z]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.01e-4
max = 0.01e-4
seed = 6
[../]
[../]
[./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]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[./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
[../]
[]
[AuxKernels]
[./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
[../]
[]
[ScalarKernels]
[./global_strain]
type = GlobalStrain
variable = global_strain
global_strain_uo = global_strain_uo
use_displaced_mesh = false
[../]
[]
[Materials]
[./Landau_P]
type = GenericConstantMaterial
prop_names = 'alpha1 alpha11 alpha12 alpha111 alpha112 alpha123 alpha1111 alpha1112 alpha1122 alpha1123'
prop_values = '-0.1722883 0.42 0.735 0.26 0.61 -3.67 0 0 0 0'
[../]
[./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 79.4 111.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 = '0. 0 0 0 0 0 0 0 0'
eigenstrain_name = eigenstrain
prefactor = 0.0
[../]
[./elasticity_tensor_1]
type = ComputeElasticityTensor
fill_method = symmetric9
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
[../]
[./slab_ferroelectric]
type = ComputeElectrostrictiveTensor
Q_mnkl = '-0.089 0.026 0.026 -0.089 0.026 -0.089 -0.03375 -0.03375 -0.03375'
C_ijkl = '175.0 79.4 79.4 175.0 79.4 175.0 111.1 111.1 111.1'
[../]
[./permitivitty_1]
type = GenericConstantMaterial
prop_names = 'permittivity'
prop_values = '0.08854187'
[../]
[]
[Kernels]
#note below we use a strain-renormalized functional for lead titanate (this is different than the stress-free functionals typically used)
# they are related by a Legendre transformation
#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]
[./Fbulk]
type = BulkEnergyEighth
execute_on = 'initial timestep_end'
[../]
[./Fwall]
type = WallEnergy
execute_on = 'initial timestep_end'
[../]
[./Felastic]
type = ElasticEnergy
execute_on = 'initial timestep_end'
use_displaced_mesh = false
[../]
[./Fcoupled]
type = ElectrostrictiveCouplingEnergy
execute_on = 'initial timestep_end'
[../]
[./Felec]
type = ElectrostaticEnergy
execute_on = 'initial timestep_end'
[../]
[./Ftotal]
type = LinearCombinationPostprocessor
pp_names = 'Fbulk Fwall Fcoupled Felec'
pp_coefs = ' 1 1 1 1'
execute_on = 'initial timestep_end'
[../]
[./perc_change]
type = PercentChangePostprocessor
postprocessor = Ftotal
execute_on = 'initial timestep_end'
[../]
#[./elapsed]
# type = PerfGraphData
# section_name = "Root"
# data_type = total
#[../]
#[./nodes]
# type = NumNodes
#[../]
[]
[UserObjects]
[./global_strain_uo]
type = GlobalATiO3MaterialRVEUserObject
use_displaced_mesh = false
execute_on = 'Initial Linear Nonlinear'
[../]
[./kill]
type = Terminator
expression = 'perc_change <= 1.0e-6'
[../]
[]
[Preconditioning]
[./smp]
type = SMP
full = true
petsc_options = '-snes_ksp_ew'
petsc_options_iname = '-ksp_gmres_restart -snes_atol -snes_rtol -ksp_rtol -pc_type -build_twosided'
petsc_options_value = ' 160 1e-8 1e-8 1e-8 bjacobi allreduce'
[../]
[]
#[Debug]
# show_var_residual_norms = true
#[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
scheme = 'implicit-euler'
dtmin = 1e-13
dtmax = 0.8
l_max_its = 200
[./TimeStepper]
type = IterationAdaptiveDT
optimal_iterations = 8
cutback_factor = 0.75
linear_iteration_ratio = 1000
dt = 0.3
[../]
num_steps = 10
#[./Adaptivity]
# coarsen_fraction = 0.1
# refine_fraction = 0.7
# max_h_level = 1
#[../]
[]
[Outputs]
print_linear_residuals = false
perf_graph = true
[./out]
type = Exodus
file_base = out_pbc
elemental_as_nodal = true
interval = 5
[../]
[./outCSV]
type = CSV
file_base = out_pbc
[../]
[]
(test/tests/auxkernels/microforce.i)
[Mesh]
[gen]
type = GeneratedMeshGenerator
dim = 3
nx = 4
ny = 4
nz = 4
xmin = -6.0
xmax = 6.0
ymin = -6.0
ymax = 6.0
zmin = -2.0
zmax = 2.0
elem_type = HEX8
[]
[./cnode]
input = gen
type = ExtraNodesetGenerator
coord = '-6.0 -6.0 -2.0'
new_boundary = 100
[../]
# additional boundary sideset (one node) to zero one of the elastic displacement vectors - eliminates rigid body translations from the degrees of freedom
[]
[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'
u_x = u_x
u_y = u_y
u_z = u_z
[]
[Variables]
[./global_strain]
order = SIXTH
family = SCALAR
[../]
[./polar_x]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.01e-4
max = 0.01e-4
seed = 6
[../]
[../]
[./polar_y]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.01e-4
max = 0.01e-4
seed = 6
[../]
[../]
[./polar_z]
order = FIRST
family = LAGRANGE
[./InitialCondition]
type = RandomIC
min = -0.01e-4
max = 0.01e-4
seed = 6
[../]
[../]
[./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]
[./disp_x]
[../]
[./disp_y]
[../]
[./disp_z]
[../]
[./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
[../]
[./Fb_x]
order = CONSTANT
family = MONOMIAL
[../]
[./Fb_y]
order = CONSTANT
family = MONOMIAL
[../]
[./Fb_z]
order = CONSTANT
family = MONOMIAL
[../]
[./Felstr_x]
order = CONSTANT
family = MONOMIAL
[../]
[./Felstr_y]
order = CONSTANT
family = MONOMIAL
[../]
[./Felstr_z]
order = CONSTANT
family = MONOMIAL
[../]
[./Felec_x]
order = CONSTANT
family = MONOMIAL
[../]
[./Felec_y]
order = CONSTANT
family = MONOMIAL
[../]
[./Felec_z]
order = CONSTANT
family = MONOMIAL
[../]
[./Fw_x]
order = CONSTANT
family = MONOMIAL
[../]
[./Fw_y]
order = CONSTANT
family = MONOMIAL
[../]
[./Fw_z]
order = CONSTANT
family = MONOMIAL
[../]
[]
[AuxKernels]
[./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
[../]
[./fbx]
type = MicroforceBulkEnergy
variable = Fb_x
component = 0
[../]
[./fby]
type = MicroforceBulkEnergy
variable = Fb_y
component = 1
[../]
[./fbz]
type = MicroforceBulkEnergy
variable = Fb_z
component = 2
[../]
[./felstrx]
type = MicroforceElectrostrictiveCouplingEnergy
variable = Felstr_x
component = 0
[../]
[./felstry]
type = MicroforceElectrostrictiveCouplingEnergy
variable = Felstr_y
component = 1
[../]
[./felstrz]
type = MicroforceElectrostrictiveCouplingEnergy
variable = Felstr_z
component = 2
[../]
[./felecx]
type = MicroforceElectrostaticEnergy
variable = Felec_x
component = 0
[../]
[./felecy]
type = MicroforceElectrostaticEnergy
variable = Felec_y
component = 1
[../]
[./felecz]
type = MicroforceElectrostaticEnergy
variable = Felec_z
component = 2
[../]
[./fwx]
type = MicroforceWallEnergy
variable = Fw_x
component = 0
[../]
[./fwy]
type = MicroforceWallEnergy
variable = Fw_y
component = 1
[../]
[./fwz]
type = MicroforceWallEnergy
variable = Fw_z
component = 2
[../]
[]
[ScalarKernels]
[./global_strain]
type = GlobalStrain
variable = global_strain
global_strain_uo = global_strain_uo
use_displaced_mesh = false
[../]
[]
[Materials]
[./Landau_P]
type = GenericConstantMaterial
prop_names = 'alpha1 alpha11 alpha12 alpha111 alpha112 alpha123 alpha1111 alpha1112 alpha1122 alpha1123'
prop_values = '-0.1722883 0.42 0.735 0.26 0.61 -3.67 0 0 0 0'
[../]
[./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 79.4 111.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 = '0. 0 0 0 0 0 0 0 0'
eigenstrain_name = eigenstrain
prefactor = 0.0
[../]
[./elasticity_tensor_1]
type = ComputeElasticityTensor
fill_method = symmetric9
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
[../]
[./slab_ferroelectric]
type = ComputeElectrostrictiveTensor
Q_mnkl = '-0.089 0.026 0.026 -0.089 0.026 -0.089 -0.03375 -0.03375 -0.03375'
C_ijkl = '175.0 79.4 79.4 175.0 79.4 175.0 111.1 111.1 111.1'
[../]
[./permitivitty_1]
type = GenericConstantMaterial
prop_names = 'permittivity'
prop_values = '0.08854187'
[../]
[]
[Kernels]
#note below we use a strain-renormalized functional for lead titanate (this is different than the stress-free functionals typically used)
# they are related by a Legendre transformation
#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]
[./Fbulk]
type = BulkEnergyEighth
execute_on = 'initial timestep_end'
[../]
[./Fwall]
type = WallEnergy
execute_on = 'initial timestep_end'
[../]
[./Felastic]
type = ElasticEnergy
execute_on = 'initial timestep_end'
use_displaced_mesh = false
[../]
[./Fcoupled]
type = ElectrostrictiveCouplingEnergy
execute_on = 'initial timestep_end'
[../]
[./Felec]
type = ElectrostaticEnergy
execute_on = 'initial timestep_end'
[../]
[./Ftotal]
type = LinearCombinationPostprocessor
pp_names = 'Fbulk Fwall Fcoupled Felec'
pp_coefs = ' 1 1 1 1'
execute_on = 'initial timestep_end'
[../]
[./perc_change]
type = PercentChangePostprocessor
postprocessor = Ftotal
execute_on = 'initial timestep_end'
[../]
#[./elapsed]
# type = PerfGraphData
# section_name = "Root"
# data_type = total
#[../]
#[./nodes]
# type = NumNodes
#[../]
[]
[UserObjects]
[./global_strain_uo]
type = GlobalATiO3MaterialRVEUserObject
use_displaced_mesh = false
execute_on = 'Initial Linear Nonlinear'
[../]
[./kill]
type = Terminator
expression = 'perc_change <= 1.0e-6'
[../]
[]
[Preconditioning]
[./smp]
type = SMP
full = true
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-8 bjacobi allreduce'
[../]
[]
#[Debug]
# show_var_residual_norms = true
#[]
[Executioner]
type = Transient
solve_type = 'PJFNK'
scheme = 'implicit-euler'
dtmin = 1e-13
dtmax = 0.8
l_max_its = 200
[./TimeStepper]
type = IterationAdaptiveDT
optimal_iterations = 8
cutback_factor = 0.75
linear_iteration_ratio = 1000
dt = 0.3
[../]
num_steps = 10
[]
[Outputs]
print_linear_residuals = false
perf_graph = true
[./out]
type = Exodus
file_base = out_microforce_test
elemental_as_nodal = true
interval = 1
[../]
[./outCSV]
type = CSV
file_base = out_microforce_test
[../]
[]