From 2b55138e9fa959fc97ab8c43495885c5e7d3cb79 Mon Sep 17 00:00:00 2001
From: Jean-Luc Fattebert <fattebertj@ornl.gov>
Date: Tue, 12 Mar 2019 15:22:45 -0400
Subject: [PATCH] Merged some changes from master branch

---
 source/BeckermannCompositionRHSStrategy.h |  2 +-
 source/QuatModel.C                        |  3 +-
 source/QuatModelParameters.C              |  3 ++
 source/QuatModelParameters.h              |  9 +++-
 source/fortran/functions.m4               | 55 ++++++++++++++++++++---
 5 files changed, 63 insertions(+), 9 deletions(-)

diff --git a/source/BeckermannCompositionRHSStrategy.h b/source/BeckermannCompositionRHSStrategy.h
index 439e045f..8648962c 100644
--- a/source/BeckermannCompositionRHSStrategy.h
+++ b/source/BeckermannCompositionRHSStrategy.h
@@ -34,7 +34,7 @@
 // POSSIBILITY OF SUCH DAMAGE.
 // 
 #ifndef included_BeckermannCompositionRHSStrategy
-#define included_BeckermannCompositionRHSStrategy 
+#define included_BeckermannCompositionRHSStrategy
 
 #include "CompositionRHSStrategy.h"
 
diff --git a/source/QuatModel.C b/source/QuatModel.C
index 0521b43e..b3cc1a13 100644
--- a/source/QuatModel.C
+++ b/source/QuatModel.C
@@ -2772,7 +2772,8 @@ void QuatModel::preRunDiagnostics( void )
          if( found_ceq && !d_is_from_restart)
          {
             setRefPhaseConcentrationsToEquilibrium(ceq);
-            setPhaseConcentrationsToEquilibrium(ceq);
+            if( d_model_parameters.initPhaseConcAtEq() )
+               setPhaseConcentrationsToEquilibrium(ceq);
          }
          
       } // with_concentration
diff --git a/source/QuatModelParameters.C b/source/QuatModelParameters.C
index 0a77c21b..8eab2fe9 100644
--- a/source/QuatModelParameters.C
+++ b/source/QuatModelParameters.C
@@ -358,6 +358,9 @@ void QuatModelParameters::readConcDB(boost::shared_ptr<tbox::Database> conc_db)
       d_initc_in_phase.resize(nterms);
       conc_db->getDoubleArray("initc_in_phase",&d_initc_in_phase[0],nterms);
    }
+
+   d_init_phase_conc_eq =
+      conc_db->getBoolWithDefault("init_phase_conc_eq", true);
 }
 
 //=======================================================================
diff --git a/source/QuatModelParameters.h b/source/QuatModelParameters.h
index eb580f8b..4e11c9fd 100644
--- a/source/QuatModelParameters.h
+++ b/source/QuatModelParameters.h
@@ -356,7 +356,12 @@ class QuatModelParameters
    {
       return d_epsilon_phase/sqrt(32.*d_phase_well_scale);
    }
- 
+
+   bool initPhaseConcAtEq()const
+   {
+      return d_init_phase_conc_eq;
+   }
+
 private:
    void readNumberSpecies(boost::shared_ptr<tbox::Database> conc_db);
 
@@ -509,6 +514,8 @@ class QuatModelParameters
    bool d_visit_energy_output;
    bool d_visit_grain_output;
 
+   bool d_init_phase_conc_eq;
+
    void readMolarVolumes(boost::shared_ptr<tbox::Database> db);
 };
 
diff --git a/source/fortran/functions.m4 b/source/fortran/functions.m4
index 38d00ce6..45b4df4e 100644
--- a/source/fortran/functions.m4
+++ b/source/fortran/functions.m4
@@ -91,6 +91,19 @@ c        Valid values for type are "quadratic" and "pbg"
          phit = min( 1.d0, phit )
          interp_func = phit
 
+      else if ( type(1:1) .eq. 'a' ) then
+
+         if( phi < 0.08333333333333333d0 )then
+            phit = max( 0.d0, phi )
+            interp_func = 54.d0*phit*phit*phit
+         else if( phi>0.9166666666666666d0 )then
+            phit = min( 1.d0, phi )
+            phit = 1.d0-phit
+            interp_func = 1.d0-54.d0*phit*phit*phit
+         else
+            interp_func = 1.125*phi-0.0625
+         endif
+
       else if ( type(1:1) .eq. '3'  ) then
 
          phit = max( 0.d0, phi )
@@ -173,6 +186,19 @@ c-----------------------------------------------------------------------
              deriv_interp_func = 0.d0
          endif
 
+      else if ( type(1:1) .eq. 'a' ) then
+
+         if( phi < 0.08333333333333333d0 )then
+            phit = max( 0.d0, phi )
+            deriv_interp_func = 162.d0*phit*phit
+         else if( phi>0.9166666666666666d0 )then
+            phit = min( 1.d0, phi )
+            phit = 1.d0-phit
+            deriv_interp_func = 162.d0*phit*phit
+         else
+            deriv_interp_func = 1.125
+         endif
+
       else if ( type(1:1) .eq. '3' ) then
 
          phit = max( 0.d0, phi )
@@ -236,6 +262,19 @@ c=======================================================================
 
          second_deriv_interp_func = 0.d0
 
+      else if ( type(1:1) .eq. 'a' ) then
+
+         if( phi < 0.08333333333333333d0 )then
+            phit = max( 0.d0, phi )
+            second_deriv_interp_func = 324.d0*phit
+         else if( phi>0.9166666666666666d0 )then
+            phit = min( 1.d0, phi )
+            phit = 1.d0-phit
+            second_deriv_interp_func = -324.d0*phit
+         else
+            second_deriv_interp_func = 0.d0
+         endif
+
       else if ( type(1:1) .eq. 'c' ) then
 
          second_deriv_interp_func = 0.d0
@@ -538,13 +577,13 @@ c-----------------------------------------------------------------------
       double precision phi, phit
       character*(*) type1, type2
 
-      if ( type1(1:1) .eq. 'p' )then
+      if ( type1(1:1) .eq. type2(1:1) )then
 
-         if ( type2(1:1) .eq. 'p' )then 
+         interp_ratio = 1.d0
 
-            interp_ratio = 1.d0
+      else if ( type1(1:1) .eq. 'p' )then
 
-         else if ( type2(1:1) .eq. 'l' )then
+         if ( type2(1:1) .eq. 'l' )then
 
             phit = max( 0.d0, min( 1.d0, phi ) )
 
@@ -577,7 +616,11 @@ c-----------------------------------------------------------------------
       double precision phi, phit
       character*(*) type1, type2
 
-      if ( type1(1:1) .eq. 'p' )then
+      if ( type1(1:1) .eq. type2(1:1) )then
+
+         compl_interp_ratio = 1.d0
+
+      else if ( type1(1:1) .eq. 'p' )then
 
          if ( type2(1:1) .eq. 'p' )then
 
@@ -592,7 +635,7 @@ c-----------------------------------------------------------------------
 
          else
 
-            print *, "Error, interp_ratio: unknown/incompatible types"
+            print *, "Error, compl_interp_ratio: unknown/incomp. types"
             stop
 
          endif