Adding tests for tau

python/KaempferPlapp2009.py

@@ -42,7 +42,9 @@
 a = 3.19e-10
 k = 1.38e-23
 a_1 = (5./8.)*sqrt(2.)
-W = 1e-6
+W = 1.e-6
+alpha = 10.e-2
+m = 2.99e-26
 
 # Eq. (2)
 P_vs = exp(K_fit[0]*T**(-1) + K_fit[1]*T**0 + K_fit[2]*T**1 + K_fit[3]*T**2 + K_fit[4]*T**3 + K_fit[5]*log(T))
@@ -65,6 +67,18 @@
 lbda = a_1 * W / d_0_prime
 lbda_fixed = a_1 * W / d_0_prime_fixed
 
+# beta0, beta0' Eq. (26)
+beta_0 = 1/alpha * sqrt(2*pi*m /(k * T))
+beta_0_prime = (rho_vs/rho_i) * beta_0
+
+# Make fixed version for testing user prescribed d0
+beta_0_fixed = 5.5e5
+beta_0_prime_fixed = (rho_vs/rho_i) * beta_0_fixed
+
+# Tau Eq. (38)
+tau = beta_0_prime * W * lbda / a_1
+tau_fixed = beta_0_prime_fixed * W * lbda / a_1
+
 # Concentration equilibrium
 u_eq = (rho_vs - rho_vs.subs(T, 263.15)) / rho_i
 
@@ -99,3 +113,21 @@
 print "lambda_fixed(263.15) = ", lbda_fixed.evalf(subs={T: 263.15})
 print "lambda_fixed(268.15) = ", lbda_fixed.evalf(subs={T: 268.15})
 print ""
+print "beta0(263.15) = ", beta_0.evalf(subs={T: 263.15})
+print "beta0(268.15) = ", beta_0.evalf(subs={T: 268.15})
+print ""
+print "beta0'(263.15) = ", beta_0_prime.evalf(subs={T: 263.15})
+print "beta0'(268.15) = ", beta_0_prime.evalf(subs={T: 268.15})
+print ""
+print "tau(263.15) = ", tau.evalf(subs={T: 263.15})
+print "tau(268.15) = ", tau.evalf(subs={T: 268.15})
+print ""
+print "beta0_fixed(263.15) = ", beta_0_fixed
+print "beta0_fixed(268.15) = ", beta_0_fixed
+print ""
+print "beta0_fixed'(263.15) = ", beta_0_prime_fixed.evalf(subs={T: 263.15})
+print "beta0_fixed'(268.15) = ", beta_0_prime_fixed.evalf(subs={T: 268.15})
+print ""
+print "tau_fixed(263.15) = ", tau_fixed.evalf(subs={T: 263.15})
+print "tau_fixed(268.15) = ", tau_fixed.evalf(subs={T: 268.15})
+print ""

src/materials/PikaMaterial.C

@@ -94,6 +94,10 @@ PikaMaterial::computeQpProperties()
   const Real _beta_0_prime = (rho_vs / _density_ice) * interface_kinetic_coefficient;
   _tau[_qp] = _beta_0_prime * _interface_thickness * _lambda[_qp] / _a_1;
 
+  std::cout << "W = " << _interface_thickness << std::endl;
+  std::cout << "lambda = " << _lambda[_qp] << std::endl;
+  std::cout << "a_1 = " << _a_1 << std::endl;
+
   // u_eq; Eq. (33)
   _equilibrium_chemical_potential[_qp] = (rho_vs- _rho_vs_T_0) / _density_ice;
 

tests/materials/tau/gold/tau_data.csv

@@ -0,0 +1,2 @@
+time,T,beta_0,beta_0_prime,tau
+1,263.15,0.0719257332409219,1.67821307631168e-7,7.38190091686490e-5

tests/materials/tau/gold/tau_fixed_data.csv

@@ -0,0 +1,2 @@
+time,T,d_0,d_0_prime,lambda
+1,263.15,1.3e-09,3.03323567366001e-15,291399538.835259

tests/materials/tau/tau.i

@@ -0,0 +1,71 @@
+[Mesh]
+  type = GeneratedMesh
+  dim = 2
+[]
+
+[Variables]
+  [./T]
+    initial_condition = 263.15
+  [../]
+[]
+
+[PikaMaterials]
+  conductivity_air = 0.1
+  phase = -1
+  interface_thickness = 1e-6
+  temperature = T
+  output_properties = 'relaxation_time interface_kinetic_coefficient interface_kinetic_coefficient_prime'
+  outputs = all
+  debug = true
+[]
+
+[Postprocessors]
+  [./T]
+    type = PointValue
+    variable = T
+    point = '0.5 0.5 0'
+  [../]
+  [./beta_0]
+    type = PointValue
+    variable = interface_kinetic_coefficient
+    point = '0.5 0.5 0'
+  [../]
+  [./beta_0_prime]
+    type = PointValue
+    variable = interface_kinetic_coefficient_prime
+    point = '0.5 0.5 0'
+  [../]
+  [./tau]
+    type = PointValue
+    variable = relaxation_time
+    point = '0.5 0.5 0'
+  [../]
+[]
+
+[Problem]
+  type = FEProblem
+  solve = false
+  kernel_coverage_check = false
+[]
+
+[Executioner]
+  # Preconditioned JFNK (default)
+  type = Transient
+  num_steps = 1
+  solve_type = PJFNK
+  petsc_options_iname = '-pc_type -pc_hypre_type'
+  petsc_options_value = 'hypre boomeramg'
+[]
+
+[Outputs]
+  exodus = true
+  [./console]
+    type = Console
+    perf_log = true
+    nonlinear_residuals = true
+    linear_residuals = true
+  [../]
+  [./data]
+    type = CSV
+  [../]
+[]

tests/materials/tau/tau_fixed.i

@@ -0,0 +1,73 @@
+[Mesh]
+  type = GeneratedMesh
+  dim = 2
+[]
+
+[Variables]
+  [./T]
+    initial_condition = 263.15
+  [../]
+[]
+
+[PikaMaterials]
+  conductivity_air = 0.1
+  interface_thickness = 1e-6
+  capillary_length = 1.3e-9
+  phase = -1
+  temperature = T
+  output_properties = 'capillary_length capillary_length_prime phase_field_coupling_constant'
+  outputs = all
+  debug = true
+[]
+
+[Postprocessors]
+  [./T]
+    type = PointValue
+    variable = T
+    point = '0.5 0.5 0'
+
+  [../]
+  [./lambda]
+    type = PointValue
+    variable = phase_field_coupling_constant
+    point = '0.5 0.5 0'
+  [../]
+  [./d_0]
+    type = PointValue
+    variable = capillary_length
+    point = '0.5 0.5 0'
+  [../]
+  [./d_0_prime]
+    type = PointValue
+    variable = capillary_length_prime
+    point = '0.5 0.5 0'
+  [../]
+[]
+
+[Problem]
+  type = FEProblem
+  solve = false
+  kernel_coverage_check = false
+[]
+
+[Executioner]
+  # Preconditioned JFNK (default)
+  type = Transient
+  num_steps = 1
+  solve_type = PJFNK
+  petsc_options_iname = '-pc_type -pc_hypre_type'
+  petsc_options_value = 'hypre boomeramg'
+[]
+
+[Outputs]
+  exodus = true
+  [./console]
+    type = Console
+    perf_log = true
+    nonlinear_residuals = true
+    linear_residuals = true
+  [../]
+  [./data]
+    type = CSV
+  [../]
+[]

tests/materials/tau/tests

@@ -0,0 +1,14 @@
+[Tests]
+  [./tai]
+    # Gold for this test comes python/KaempferPlapp2009.py
+    type = 'CSVDiff'
+    input = 'tau.i'
+    csvdiff = 'tau_data.csv'
+  [../]
+  [./tau_beta0_fixed]
+    # Gold for this test comes python/KaempferPlapp2009.py
+    type = 'CSVDiff'
+    input = 'tau_fixed.i'
+    csvdiff = 'tau_fixed_data.csv'
+  [../]
+[]