Updated docs (especially the dendrite app) (#91)

* updates to the dendritic solidification app, cleaning up the derivation. still needs some more explanation

* finished overhaul of the dendriticSolidification documentation, incl notational changes to the equations.h file

* update the user guide for v2.0.1, fix duplicate documentation file for precipitateEvolution

* minor tweaks to the dendritic solidification documentation

applications/dendriticSolidification/ICs_and_BCs.h

@@ -24,8 +24,8 @@ double InitialCondition<dim>::value (const Point<dim> &p, const unsigned int com
 
 	  // Initial condition for the concentration field
 	  if (index == 0){
-		  double deltaV = userInputs.get_model_constant_double("deltaV");
-		  scalar_IC = -deltaV;
+		  double delta = userInputs.get_model_constant_double("delta");
+		  scalar_IC = -delta;
 	  }
 	  // Initial condition for the order parameter field
 	  else if (index == 1) {

applications/dendriticSolidification/customPDE.h

@@ -40,9 +40,9 @@ class customPDE: public MatrixFreePDE<dim,degree>
 	// Model constants specific to this subclass
 	// ================================================================
 
-	double DV = userInputs.get_model_constant_double("DV");
+	double D = userInputs.get_model_constant_double("D");
 	double W0 = userInputs.get_model_constant_double("W0");
-	double deltaV = userInputs.get_model_constant_double("deltaV");
+	double delta = userInputs.get_model_constant_double("delta");
 	double epsilonM = userInputs.get_model_constant_double("epsilonM");
 	double theta0 = userInputs.get_model_constant_double("theta0");
 	double mult = userInputs.get_model_constant_double("mult");

applications/dendriticSolidification/equations.h

@@ -5,7 +5,7 @@
 // =================================================================================
 void variableAttributeLoader::loadVariableAttributes(){
 	// Variable 0
-	set_variable_name				(0,"T");
+	set_variable_name				(0,"u");
 	set_variable_type				(0,SCALAR);
 	set_variable_equation_type		(0,PARABOLIC);
 
@@ -17,7 +17,7 @@ void variableAttributeLoader::loadVariableAttributes(){
 	set_need_gradient_residual_term	(0,true);
 
     // Variable 1
-	set_variable_name				(1,"n");
+	set_variable_name				(1,"phi");
 	set_variable_type				(1,SCALAR);
 	set_variable_equation_type		(1,PARABOLIC);
 
@@ -50,9 +50,6 @@ void variableAttributeLoader::loadVariableAttributes(){
 // equations (or parts of residual equations) can be written below in "residualRHS".
 
 
-// Free energy expression (or rather, its derivative)
-
-
 
 // =================================================================================
 // residualRHS
@@ -68,54 +65,54 @@ template <int dim, int degree>
 void customPDE<dim,degree>::residualRHS(variableContainer<dim,degree,dealii::VectorizedArray<double> > & variable_list,
 				 dealii::Point<dim, dealii::VectorizedArray<double> > q_point_loc) const {
 
-// The temperature and its derivatives 
-scalarvalueType T = variable_list.get_scalar_value(0);
-scalargradType Tx = variable_list.get_scalar_gradient(0);
+// The temperature and its derivatives
+scalarvalueType u = variable_list.get_scalar_value(0);
+scalargradType ux = variable_list.get_scalar_gradient(0);
 
-// The order parameter and its derivatives 
-scalarvalueType n = variable_list.get_scalar_value(1);
-scalargradType nx = variable_list.get_scalar_gradient(1);
+// The order parameter and its derivatives
+scalarvalueType phi = variable_list.get_scalar_value(1);
+scalargradType phix = variable_list.get_scalar_gradient(1);
 
-// The order parameter chemical potential and its derivatives 
+// The order parameter chemical potential and its derivatives
 scalarvalueType mu = variable_list.get_scalar_value(2);
 
-double lambdaV = (DV/0.6267/W0/W0);
+// The coupling constant, determined from solvability theory
+double lambda = (D/0.6267/W0/W0);
 
-// Derivative of the free energy density with respect to n
-scalarvalueType fnV = (-(n-constV(lambdaV)*T*(constV(1.0)-n*n))*(constV(1.0)-n*n));
+// Derivative of the free energy density with respect to phi
+scalarvalueType f_phi = -(phi-constV(lambda)*u*(constV(1.0)-phi*phi))*(constV(1.0)-phi*phi);
 
+// The azimuthal angle
 scalarvalueType theta;
-
-for (unsigned i=0; i< n.n_array_elements;i++){
-	theta[i] = std::atan2(nx[1][i],nx[0][i]);
+for (unsigned i=0; i< phi.n_array_elements;i++){
+	theta[i] = std::atan2(phix[1][i],phix[0][i]);
 }
 
 // Anisotropic gradient energy coefficient, its derivative and square
-scalarvalueType WV = (constV(W0)*(constV(1.0)+constV(epsilonM)*std::cos(constV(mult)*(theta-constV(theta0)))));
-scalarvalueType WTV = (constV(W0)*(-constV(epsilonM)*constV(mult)*std::sin(constV(mult)*(theta-constV(theta0)))));
-scalarvalueType tauV = (WV*WV);
-
-
+scalarvalueType W = constV(W0)*(constV(1.0)+constV(epsilonM)*std::cos(constV(mult)*(theta-constV(theta0))));
+scalarvalueType W_theta = constV(-W0)*(constV(epsilonM)*constV(mult)*std::sin(constV(mult)*(theta-constV(theta0))));
+scalarvalueType tau = W/constV(W0);
 
+// The anisotropy term that enters in to the residual equation for mu
 scalargradType aniso;
-aniso[0] = -WV*WV*nx[0]+WV*WTV*nx[1];
-aniso[1] = -WV*WV*nx[1]-WV*WTV*nx[0];
+aniso[0] = W*W*phix[0]-W*W_theta*phix[1];
+aniso[1] = W*W*phix[1]+W*W_theta*phix[0];
 
 // Define required residuals
-scalarvalueType rTV = (T-constV(0.5)*mu*constV(userInputs.dtValue)/tauV);
-scalargradType rTxV = (constV(-DV*userInputs.dtValue)*Tx);
-scalarvalueType rnV = (n-constV(userInputs.dtValue)*mu/tauV);
-scalarvalueType rmuV = (fnV);
+scalarvalueType ruV = (u+constV(0.5)*mu*constV(userInputs.dtValue)/tau);
+scalargradType ruxV = (constV(-D*userInputs.dtValue)*ux);
+scalarvalueType rphiV = (phi+constV(userInputs.dtValue)*mu/tau);
+scalarvalueType rmuV = (-f_phi);
 scalargradType rmuxV = (-aniso);
 
-// Residuals for the equation to evolve the concentration 
-variable_list.set_scalar_value_residual_term(0,rTV);
-variable_list.set_scalar_gradient_residual_term(0,rTxV);
+// Residuals for the equation to evolve the concentration
+variable_list.set_scalar_value_residual_term(0,ruV);
+variable_list.set_scalar_gradient_residual_term(0,ruxV);
 
-// Residuals for the equation to evolve the order parameter 
-variable_list.set_scalar_value_residual_term(1,rnV);
+// Residuals for the equation to evolve the order parameter
+variable_list.set_scalar_value_residual_term(1,rphiV);
 
-// Residuals for the equation to evolve the order parameter chemical potential 
+// Residuals for the equation to evolve the order parameter chemical potential
 variable_list.set_scalar_value_residual_term(2,rmuV);
 variable_list.set_scalar_gradient_residual_term(2,rmuxV);
 

applications/dendriticSolidification/parameters.in

@@ -1,4 +1,4 @@
-# Parameter list for the Allen-Cahn example application
+# Parameter list for the dendriticSolidification example application
 # Refer to the PRISMS-PF manual for use of these parameters in the source code.
 
 # =================================================================================
@@ -44,7 +44,7 @@ set Max refinement level = 6
 set Min refinement level = 0
 
 # Set the fields used to determine the refinement using their index.
-set Refinement criteria fields = n
+set Refinement criteria fields = phi
 
 # Set the maximum and minimum value of the fields where the mesh should be refined
 set Refinement window max = 0.9999
@@ -83,8 +83,8 @@ set Number of time steps = 25000
 # 1.5 on the top and bottom for variable 'n' in 2D
 # set Boundary condition for variable n = NATURAL, NATURAL, DIRICHLET: 1.5, DIRICHLET: 1.5
 
-set Boundary condition for variable T = DIRICHLET: -0.75 # equal to -deltaV
-set Boundary condition for variable n = NATURAL
+set Boundary condition for variable u = DIRICHLET: -0.75 # equal to -deltaV
+set Boundary condition for variable phi = NATURAL
 set Boundary condition for variable mu = NATURAL
 
 # =================================================================================
@@ -97,13 +97,13 @@ set Boundary condition for variable mu = NATURAL
 # where [symmetry] is ISOTROPIC, TRANSVERSE, ORTHOTROPIC, or ANISOTROPIC.
 
 # Heat diffusion constant
-set Model constant DV = 1.0, DOUBLE
+set Model constant D = 1.0, DOUBLE
 
 # Isotropic gradient energy coefficient
 set Model constant W0 = 1.0, DOUBLE
 
 # Undercooling
-set Model constant deltaV = 0.75, DOUBLE
+set Model constant delta = 0.75, DOUBLE
 
 # Anisotropy paramter
 set Model constant epsilonM = 0.05, DOUBLE
@@ -126,12 +126,18 @@ set Output condition = EQUAL_SPACING
 set Number of outputs = 10
 
 # The number of time steps between updates being printed to the screen
-set Skip print steps = 1000
+set Skip print steps = 1 #1000
 
 # =================================================================================
 # Set the checkpoint/restart parameters
 # =================================================================================
-
+# Whether to start this simulation from the checkpoint of a previous simulation
 set Load from a checkpoint = false
+
+# Type of spacing between checkpoints ("EQUAL_SPACING", "LOG_SPACING", "N_PER_DECADE",
+# or "LIST")
 set Checkpoint condition = EQUAL_SPACING
-set Number of checkpoints = 10
+
+# Number of times the creates checkpoints (total number for "EQUAL_SPACING"
+# and "LOG_SPACING", number per decade for "N_PER_DECADE", ignored for "LIST")
+set Number of checkpoints = 2

applications/dendriticSolidification/tex_files/dendriticSolidification.tex

@@ -183,6 +183,8 @@
 This example application implements a simple model of dendritic solidification based on the CHiMaD Benchmark Problem 3, itself based on the model given in the following article: \\
 ``Multiscale Finite-Difference-Diffusion-Monte-Carlo Method for Simulating Dendritic Solidification'' by M. Plapp and A. Karma, \emph{Journal of Computational Physics}, 165, 592-619 (2000)
 
+This example application examines the non-isothermal solidification of a pure substance. The simulation starts with a circular solid seed in a uniformly undercooled liquid. As this seed grows, two variables are tracked, an order parameter, $\phi$, that denotes whether the material a liquid or solid and a nondimensional temperature,  $u$. The crystal structure of the solid is offset from the simulation frame for generality and to expose more readily any effects of the mesh on the dendrite shape.
+
 \section{Governing Equations}
 
 Consider a free energy density given by:
@@ -199,14 +201,14 @@ \section{Governing Equations}
 \end{equation}
 where $\lambda$ is a dimensionless coupling constant. The gradient energy coefficient, $W$, is given by 
 \begin{equation}
-W = W_0 [1+\epsilon_m \cos[m(\theta-\theta_0)]]
+W(\theta) = W_0 [1+\epsilon_m \cos[m(\theta-\theta_0)]]
 \end{equation}
-where, $W_0$, $\epsilon_m$, and $\theta_0$ are constants and $\theta$ is the in-plane azimuthal angle, where $\tan(\theta) = \frac{n_y}{n_x}$, where $n_y$ and $n_x$ are components of the normal vector $\hat{n} = \frac{\nabla \phi}{|\nabla \phi|}$.
+where, $W_0$, $\epsilon_m$, and $\theta_0$ are constants and $\theta$ is the in-plane azimuthal angle, where $\tan(\theta) = \frac{\partial \phi}{\partial y} / \frac{\partial \phi}{\partial x}$.
 
 The evolution equations are:
 \begin{gather}
-\frac{\partial u}{\partial t} = - \frac{\delta \Pi}{\delta \phi} = D \nabla^2 u + \frac{1}{2} \tau(\hat{n}) \frac{\partial \phi}{\partial t} \\
-\tau(\hat{n}) \frac{\partial \phi}{\partial t} = \frac{\partial f}{\partial \phi} + \frac{\partial}{\partial x} \left[ |\nabla \phi|^2 W(\hat{n}) \frac{\partial W(\hat{n})}{\partial \left( \frac{\partial \phi}{\partial x} \right)} \right] \hat{x} + \frac{\partial}{\partial y} \left[ |\nabla \phi|^2 W(\hat{n}) \frac{\partial W(\hat{n})}{\partial \left( \frac{\partial \phi}{\partial y} \right)} \right] \hat{y} 
+\frac{\partial u}{\partial t} = D \nabla^2 u + \frac{1}{2}  \frac{\partial \phi}{\partial t} \\
+\tau(\hat{n}) \frac{\partial \phi}{\partial t} = -\frac{\partial f}{\partial \phi} + \nabla \cdot \left[W^2(\theta) \nabla \phi \right]+  \frac{\partial}{\partial x} \left[ |\nabla \phi|^2 W(\theta) \frac{\partial W(\theta)}{\partial \left( \frac{\partial \phi}{\partial x} \right)} \right] + \frac{\partial}{\partial y} \left[ |\nabla \phi|^2 W(\theta) \frac{\partial W(\theta)}{\partial \left( \frac{\partial \phi}{\partial y} \right)} \right] 
 \end{gather}
 where
 \begin{gather}
@@ -216,35 +218,48 @@ \section{Governing Equations}
 
 The governing equations can be written more compactly using the variable $\mu$, the driving force for the phase transformation:
 \begin{gather}
-\frac{\partial u}{\partial t} = D \nabla^2 u - \frac{\mu}{\tau} \\
+\frac{\partial u}{\partial t} = D \nabla^2 u + \frac{\mu}{2 \tau} \\
 \tau(\hat{n}) \frac{\partial \phi}{\partial t} = \mu \\
-\mu = \left(\frac{\partial f}{\partial \phi}\right) - \nabla \cdot \left[\left(W^2 \frac{\partial \phi}{\partial x} - W \frac{\partial W}{\partial u} \frac{\partial \phi}{\partial y}\right)\hat{x}  + \left(W^2 \frac{\partial \phi}{\partial y} + W \frac{\partial W}{\partial u} \frac{\partial \phi}{\partial x}\right) \hat{y} \right]
+\mu = -\frac{\partial f}{\partial \phi} + \nabla \cdot \left[W^2(\theta) \nabla \phi \right]+  \frac{\partial}{\partial x} \left[ |\nabla \phi|^2 W(\theta) \frac{\partial W(\theta)}{\partial \left( \frac{\partial \phi}{\partial x} \right)} \right] + \frac{\partial}{\partial y} \left[ |\nabla \phi|^2 W\theta) \frac{\partial W(\theta)}{\partial \left( \frac{\partial \phi}{\partial y} \right)} \right] 
 \end{gather}
 
+The  $\frac{\partial W(\theta)}{\partial \left( \frac{\partial \phi}{\partial x} \right)}$ and $\frac{\partial W(\theta)}{\partial \left( \frac{\partial \phi}{\partial y} \right)}$ expressions can be evaluated using the chain rule, using $\theta$ as an intermediary (i.e. $\frac{\partial W(\theta)}{\partial \left( \frac{\partial \phi}{\partial x} \right)}=\frac{\partial W(\theta)}{\partial \theta} \frac{\partial \theta}{\partial \left( \frac{\partial \phi}{\partial x} \right)}$  and $\frac{\partial W(\theta)}{\partial \left( \frac{\partial \phi}{\partial y} \right)}=\frac{\partial W(\theta)}{\partial \theta} \frac{\partial \theta}{\partial \left( \frac{\partial \phi}{\partial y} \right)}$). Also, the last two terms can be expressed using a divergence operator, allowing them to be grouped with the second term, which will simplify matters later. Carrying out these transformations yields:
+\begin{multline}
+\mu = \left[ \phi - \lambda u \left(1 - \phi^2 \right) \right] \left(1-\phi^2\right) + \nabla \cdot \bigg[\left(W^2 \frac{\partial \phi}{\partial x} + W_0 \epsilon_m m W(\theta) \sin \left[ m \left(\theta - \theta_0 \right) \right] \frac{\partial \phi}{\partial y}\right)\hat{x}  \\
++ \left(W^2 \frac{\partial \phi}{\partial y} -W_0 \epsilon_m m W(\theta) \sin \left[ m \left(\theta - \theta_0 \right) \right] \frac{\partial \phi}{\partial x}\right) \hat{y} \bigg]
+\end{multline}
+
 \section{Model Constants}
 $W_0$: Controls the interfacial thickness, default value of 1.0. \\
 $\tau_0$: Controls the phase transformation kinetics, default value of 1.0. \\
 $\epsilon_m$: T the strength of the anisotropy, default value of 0.05. \\
 $D$: The thermal diffusion constant, default value of 1.0. \\
-$\Delta = \frac{T_m-T_0}{L/c_p}$: The level of undercooling, default value of 0.75. \\
-$\theta_0$ = The rotation angle of the anisotropy, default value of 0.125 ($\sim$7.16$^\circ$).
+$\Delta: \frac{T_m-T_0}{L/c_p}$: The level of undercooling, default value of 0.75. \\
+$\theta_0$: The rotation angle of the anisotropy with respect to the simulation frame, default value of 0.125 ($\sim$7.2$^\circ$).
 
 \section{Time Discretization}
 Considering forward Euler explicit time stepping, we have the time discretized kinetics equation:
 \begin{gather}
-u^{n+1} = \left(u^{n} - \frac{\mu^n \Delta t}{\tau}\right) + D \Delta t \nabla^2 u^n \\
-\phi^{n+1} = \left(\phi^n - \frac{\Delta t \mu^n}{\tau}\right) \\
- \mu^{n+1} = \left(\frac{\partial f^n}{\partial \phi}\right) - \nabla \cdot \left[\left(W^2 \frac{\partial \phi^n}{\partial x} - W \frac{\partial W}{\partial u^n} \frac{\partial \phi^n}{\partial y}\right)\hat{x}  + \left(W^2 \frac{\partial \phi^n}{\partial y} + W \frac{\partial W}{\partial u} \frac{\partial \phi^n}{\partial x}\right) \hat{y} \right]
+u^{n+1} = u^{n} + \Delta t \left( D  \nabla^2 u^n + \frac{\mu^n}{2 \tau} \right) \\
+\phi^{n+1} = \phi^n + \frac{\Delta t \mu^n}{\tau}
 \end{gather}
+\begin{multline}
+\mu^{n+1} =  \left[ \phi^n - \lambda u \left(1 - (\phi^n)^2 \right) \right] \left(1-(\phi^n)^2\right) + \nabla \cdot \bigg[\left(W^2 \frac{\partial \phi^n}{\partial x} + W_0 \epsilon_m m W(\theta^n) \sin \left[ m \left(\theta^n - \theta_0 \right) \right] \frac{\partial \phi^n}{\partial y}\right)\hat{x}  \\
++ \left(W^2 \frac{\partial \phi^n}{\partial y} -W_0 \epsilon_m m W(\theta^n) \sin \left[ m \left(\theta^n - \theta_0 \right) \right] \frac{\partial \phi^n}{\partial x}\right) \hat{y} \bigg]
+\end{multline}
 
 \section{Weak Formulation}
 
 \begin{gather}
-\int_{\Omega}   w  u^{n+1}  ~dV = \int_{\Omega}   w \underbrace{\left(u^{n} - \frac{\mu^n \Delta t}{\tau}\right)}_{r_u} + \nabla w \underbrace{(-D \Delta t \nabla u^n)}_{r_{ux}} ~dV \\
-\int_{\Omega}   w  \phi^{n+1}  ~dV = \int_{\Omega}   w \underbrace{\left(\phi^n - \frac{\Delta t \mu^n}{\tau}\right)}_{r_{\phi}} ~dV \\
-\int_{\Omega}   w  \mu^{n+1}  ~dV = \int_{\Omega}   w \underbrace{\left(\frac{\partial f^n}{\partial \phi}\right)}_{r_{\mu}} + \nabla w \cdot \underbrace{\left[\left(W^2 \frac{\partial \phi^n}{\partial x} - W \frac{\partial W}{\partial u^n} \frac{\partial \phi^n}{\partial y}\right)\hat{x}  + \left(W^2 \frac{\partial \phi^n}{\partial y} + W \frac{\partial W}{\partial u} \frac{\partial \phi^n}{\partial x}\right) \hat{y} \right]}_{r_{\phi x}}  ~dV 
-\end{gather}
-
+\int_{\Omega}   w  u^{n+1}  ~dV = \int_{\Omega}   w \underbrace{\left(u^{n} + \frac{\mu^n \Delta t}{2 \tau}\right)}_{r_u} + \nabla w \cdot \underbrace{(-D \Delta t \nabla u^n)}_{r_{ux}} ~dV \\
+\int_{\Omega}   w  \phi^{n+1}  ~dV = \int_{\Omega}   w \underbrace{\left(\phi^n + \frac{\Delta t \mu^n}{\tau}\right)}_{r_{\phi}} ~dV 
+\end{gather} \small
+\begin{multline}
+\int_{\Omega}   w  \mu^{n+1}  ~dV = \int_{\Omega}   w \underbrace{\left[ \phi^n - \lambda u \left(1 - (\phi^n)^2 \right) \right] \left(1-(\phi^n)^2\right)}_{r_{\mu}} \\
++ \nabla w \cdot \underbrace{\bigg[-\left(W^2 \frac{\partial \phi^n}{\partial x} + W_0 \epsilon_m m W(\theta^n) \sin \left[ m \left(\theta^n - \theta_0 \right) \right] \frac{\partial \phi^n}{\partial y}\right)\hat{x}
+- \left(W^2 \frac{\partial \phi^n}{\partial y} -W_0 \epsilon_m m W(\theta^n) \sin \left[ m \left(\theta^n - \theta_0 \right) \right] \frac{\partial \phi^n}{\partial x}\right) \hat{y} \bigg]}_{r_{\phi x}}  ~dV 
+\end{multline}
+\normalsize
 
 
 

applications/precipitateEvolution/parameters.in

@@ -95,10 +95,16 @@ set Skip print steps = 1000
 # =================================================================================
 # Set the checkpoint/restart parameters
 # =================================================================================
-
+# Whether to start this simulation from the checkpoint of a previous simulation
 set Load from a checkpoint = false
+
+# Type of spacing between checkpoints ("EQUAL_SPACING", "LOG_SPACING", "N_PER_DECADE",
+# or "LIST")
 set Checkpoint condition = EQUAL_SPACING
-set Number of checkpoints = 10
+
+# Number of times the creates checkpoints (total number for "EQUAL_SPACING"
+# and "LOG_SPACING", number per decade for "N_PER_DECADE", ignored for "LIST")
+set Number of checkpoints = 2
 
 # =================================================================================
 # Set the boundary conditions