Merge pull request #4 from LLNL/KKSdilute

KKS dilute

examples/DiluteAlCu/2d.input

@@ -0,0 +1,110 @@
+// Units are:
+//   length: um = 10^-6 m
+//   energy: pJ = 10^-12 J
+//
+model_type = "KWCcomplex"
+
+end_time = 1.e-2  // required
+max_timesteps = 5000  // optional, default is INT_MAX (alias: max_cycles)
+
+Verbosity {  // optional, block
+   level = 1  // optional, default=1
+}
+
+Visit {
+   interval = 1.e-4
+   interval_type = "time"
+}
+
+TimerManager {
+    timer_list = "solv::*::*","xfer::*::*","AMPE::*::*"
+    print_threshold = 0.0
+}
+
+ModelParameters {  // required block
+
+   epsilon_anisotropy = 0.03 //delta
+
+   H_parameter = 0. // we need quaternions to define anisotropy
+
+   epsilon_phi = 0.2   // [(pJ/um)^1/2]
+
+   phi_mobility_type = "kim"
+   //phi_mobility = 1.e3
+
+   orient_mobility = 1.
+   epsilon_orient = 1.e3 // to avoid negative diffusion for q
+
+   // free energy parameters:
+   // f(phi) = scale_energy_well*g(phi)
+   // where g is a double well potential
+   phi_well_scale = 0.3 // [pJ/um^3]
+   phi_well_func_type = "double"
+
+   Temperature{
+      temperature = 915.
+   }
+
+   ConcentrationModel {
+      model = "dilute"
+      //antitrapping = TRUE
+      rhs_form = "ebs"
+      molar_volume = 1.095e-5
+
+      liquidus_slope = -640 // K
+      meltingT = 933. // K
+      keq      = 0.14
+
+      diffusion_type = "temperature_dependent"
+      D_solid  =  3.e-1 //um**2/s
+      D_liquid  = 3.e3  //um**2/s
+      Q0_solid  = 0.    // [J/mol]
+      Q0_liquid = 0.    // [J/mol]
+   }
+
+   phi_interp_func_type ="pbg"
+   conc_interp_func_type = "l"
+   diffusion_interp_func_type = "l"
+
+   BoundaryConditions {
+      Phase{
+         boundary_0 = "slope", "0"
+         boundary_1 = "slope", "0"
+         boundary_2 = "slope", "0"
+         boundary_3 = "slope", "0"
+      }
+      Conc{
+         boundary_0 = "slope", "0"
+         boundary_1 = "slope", "0"
+         boundary_2 = "slope", "0"
+         boundary_3 = "slope", "0"
+      }
+      Quat{
+         boundary_0 = "slope", "0"
+         boundary_1 = "slope", "0"
+         boundary_2 = "slope", "0"
+         boundary_3 = "slope", "0"
+      }
+   }
+}
+
+InitialConditions {
+   filename = "512x512.nc"  // required
+   init_q   = 1., 0.
+}
+
+ScalarDiagnostics {
+   interval = 1.e-4
+   interval_type = "time"
+}   
+
+Integrator {
+   atol = 1.e-5
+}
+
+Geometry{
+   periodic_dimension = 0, 0
+   coarsest_level_resolution = 512, 512  // required
+   x_lo         = 0., 0. // lower end of computational domain.
+   x_up         = 32.0, 32.0 // upper end of computational domain.
+}

examples/DiluteAlCu/README

@@ -0,0 +1,5 @@
+python ~/GIT/AMPE/utils/make_nuclei.py -x 512 -y 512 -z 1 -r 15 \
+   --ngrains 1 --concentration-in 0.003 --concentration-out 0.02 \
+   512x512.nc 
+
+

source/CALPHADConcSolverBinary.C

@@ -35,6 +35,7 @@
 // 
 #include "CALPHADConcSolverBinary.h"
 #include "CALPHADFunctions.h"
+#include "xlogx.h"
 
 #include <iostream>
 #include <cmath>

source/CALPHADFunctions.C

@@ -34,63 +34,13 @@
 // POSSIBILITY OF SUCH DAMAGE.
 // 
 #include "CALPHADFunctions.h"
+#include "xlogx.h"
 
 #include <cmath>
 #include <cassert>
 #include <stddef.h>
 using namespace std;
 
-static const double s_smallx = 1.0e-8;
-static const double s_inv_smallx = 1./s_smallx;
-static const double s_log_smallx = log( s_smallx );
-static const double s_smallx_log_smallx = s_smallx * s_log_smallx;
-
-//-----------------------------------------------------------------------
-// C2 extension of x(log(x) function for x<=smallx
-//
-double xlogx( const double x )
-{
-   double r;
-
-   if ( x > s_smallx ) {
-      r = x * log( x );
-   }
-   else {
-      r = s_smallx_log_smallx +
-         ( x - s_smallx ) * s_log_smallx +
-         0.5 * ( x * x * s_inv_smallx - s_smallx );
-   }
-
-   return r;
-}
-
-double xlogx_deriv( const double x )
-{
-   double r;
-
-   if ( x > s_smallx ) {
-      r = log( x ) + 1.0;
-   }
-   else {
-      r = s_log_smallx + x * s_inv_smallx;
-   }
-
-   return r;
-}
-
-double xlogx_deriv2( const double x )
-{
-   double r;
-
-   if ( x > s_smallx ) {
-      r = 1./x;
-   }
-   else {
-      r = s_inv_smallx;
-   }
-
-   return r;
-}
 
 double CALPHADcomputeFMixBinary(
    const double l0,

source/CALPHADFunctions.h

@@ -38,10 +38,6 @@
 
 #include <vector>
 
-double xlogx( const double x );
-double xlogx_deriv( const double x );
-double xlogx_deriv2( const double x );
-
 double CALPHADcomputeFMixBinary(
    const double l0,
    const double l1,

source/CMakeLists.txt

@@ -19,6 +19,11 @@ MACRO(m4tof filename dir)
 ENDMACRO(m4tof)
 
 set(SOURCES ${CMAKE_SOURCE_DIR}/source/main.C
+            ${CMAKE_SOURCE_DIR}/source/KKSdiluteEquilibriumPhaseConcentrationsStrategy.C
+            ${CMAKE_SOURCE_DIR}/source/xlogx.C
+            ${CMAKE_SOURCE_DIR}/source/KKSdiluteBinaryConcentrationSolver.C
+            ${CMAKE_SOURCE_DIR}/source/KKSFreeEnergyFunctionDiluteBinary.C
+            ${CMAKE_SOURCE_DIR}/source/KKSdiluteBinary.C
             ${CMAKE_SOURCE_DIR}/source/CALPHADFunctions.C
             ${CMAKE_SOURCE_DIR}/source/NewtonSolver.C
             ${CMAKE_SOURCE_DIR}/source/DampedNewtonSolver.C
@@ -187,6 +192,7 @@ add_executable(plotCALPHADbinary plotCALPHADbinary.C
                                  CALPHADSpeciesPhaseGibbsEnergy.C
                                  CALPHADSpeciesPhaseGibbsEnergyExpansion.C
                                  CALPHADFunctions.C
+                                 xlogx.C
                                  CALPHADConcSolverBinary.C
                                  DampedNewtonSolver.C
                                  NewtonSolver.C
@@ -201,6 +207,7 @@ add_executable(plotEquilibriumCompositions
                   CALPHADSpeciesPhaseGibbsEnergy.C
                   CALPHADSpeciesPhaseGibbsEnergyExpansion.C
                   CALPHADFunctions.C
+                  xlogx.C
                   CALPHADConcSolverBinary.C
                   DampedNewtonSolver.C
                   NewtonSolver.C
@@ -217,6 +224,7 @@ add_executable(plotEquilibriumCompositionsTernary
                   CALPHADSpeciesPhaseGibbsEnergy.C
                   CALPHADSpeciesPhaseGibbsEnergyExpansion.C
                   CALPHADFunctions.C
+                  xlogx.C
                   CALPHADConcSolverTernary.C
                   DampedNewtonSolver.C
                   NewtonSolver.C
@@ -227,13 +235,15 @@ target_link_libraries(plotEquilibriumCompositionsTernary
                          ${PROJECT_LINK_LIBS_HDF5})
 
 add_executable(computeCALPHADbinaryEquilibrium
+               ${CMAKE_SOURCE_DIR}/source/xlogx.C
                computeCALPHADbinaryEquilibrium.C
                ${CMAKE_SOURCE_DIR}/source/CALPHADFreeEnergyFunctionsBinary.C
                ${CMAKE_SOURCE_DIR}/source/CALPHADEqConcSolverBinary.C
                ${CMAKE_SOURCE_DIR}/source/CALPHADConcSolverBinary.C
                ${CMAKE_SOURCE_DIR}/source/CALPHADSpeciesPhaseGibbsEnergy.C
                ${CMAKE_SOURCE_DIR}/source/CALPHADSpeciesPhaseGibbsEnergyExpansion.C
                ${CMAKE_SOURCE_DIR}/source/CALPHADFunctions.C
+               ${CMAKE_SOURCE_DIR}/source/xlogx.C
                ${CMAKE_SOURCE_DIR}/source/NewtonSolver.C
                ${CMAKE_SOURCE_DIR}/source/DampedNewtonSolver.C
                ${CMAKE_SOURCE_DIR}/source/math_utilities.C

source/DampedNewtonSolver.C

@@ -88,13 +88,15 @@ void DampedNewtonSolver::UpdateSolution(
          mwork[ii][jj] = fvec[ii];
       }
 
-      del[jj] = D_inv * Determinant( mwork );
+      const double Dmwork = Determinant( mwork );
+      //cout << "nn="<<nn<<", Dmwork="<<Dmwork <<endl;
+      del[jj] = D_inv * Dmwork;
 
       const double maxdel=0.25;
       if( fabs(del[jj])>maxdel)
          del[jj] = del[jj]>0 ? maxdel : -maxdel;
 
-      //cout << "del_c[" << jj << "] = " << del_c[jj] << endl;
+      //cout << "del[" << jj << "] = " << del[jj] << endl;
    }
 
    double w = d_alpha;

source/FieldsInitializer.C

@@ -138,6 +138,7 @@ void FieldsInitializer::initializeLevelFromData(
 
    size_t qlen_file = 0;
 #ifdef HAVE_NETCDF3
+   tbox::plog << "Try to read phase values..." << endl;
    NcVar* ncPhase = ncf->get_var( "phase" );
    if ( ncPhase == NULL ) {
       TBOX_ERROR( "Could not read variable 'phase' from input data" << endl );
@@ -174,8 +175,10 @@ void FieldsInitializer::initializeLevelFromData(
    }
 #endif
 #ifdef HAVE_NETCDF4
+   tbox::plog << "Reading fields..." <<endl;
    NcVar ncPhase;
    if ( d_phase_id>=0 ){
+      tbox::plog << "Reading phase values..." << endl;
       ncPhase = ncf->getVar( "phase" );
       if(ncPhase.isNull())
          TBOX_ERROR( "Could not read variable 'phase' from input data" << endl );
@@ -189,13 +192,15 @@ void FieldsInitializer::initializeLevelFromData(
 
    NcVar ncTemp;
    if ( readT() && d_temperature_id>=0 ) {
+      tbox::plog << "Reading temperature values..." << endl;
       ncTemp = ncf->getVar( "temperature" );
       if(ncTemp.isNull())
          TBOX_ERROR( "Could not read variable 'temperature' " <<
                      "from input data" << endl );
    }
    if ( readQ() )
    for ( int ii = 0; ii < d_qlen; ii++ ) {
+      tbox::plog << "Reading quaternion values..." << endl;
       std::ostringstream o;
       o << "quat" << ii+1;
       NcVar ncv=ncf->getVar( o.str() );
@@ -224,6 +229,8 @@ void FieldsInitializer::initializeLevelFromData(
       dims.push_back(ncPhase.getDim(1));
       dims.push_back(ncPhase.getDim(0));
    }else if( readT() ){
+      if(ncTemp.isNull())
+         TBOX_ERROR( "Field 'temperature' was not read" << endl);
       dims.push_back(ncTemp.getDim(2));
       dims.push_back(ncTemp.getDim(1));
       dims.push_back(ncTemp.getDim(0));

source/FreeEnergyFunctions.h

@@ -39,6 +39,7 @@
 #include "Phases.h"
 
 #include <string>
+#include <vector>
 
 class FreeEnergyFunctions
 {
@@ -57,6 +58,13 @@ class FreeEnergyFunctions
    virtual void printEnergyVsComposition(const double temperature, 
                                          const int npts=100 )=0;
 
+   virtual void computeSecondDerivativeFreeEnergy(
+      const double temp,
+      const double* const conc,
+      const PHASE_INDEX pi,
+      std::vector<double>& d2fdc2)=0;
+
+
    virtual bool computeCeqT(
       const double temperature,
       const PHASE_INDEX pi0, const PHASE_INDEX pi1,