Grain orientations for multiorder

CMakeLists.txt

@@ -184,6 +184,7 @@ FortranCInterface_HEADER(
         STORAGE_INDEX_STIFFNESS
         GRADIENT_FLUX
         COMPUTE_FLUX_ISOTROPIC
+        ANISOTROPIC_FLUX
         ANISOTROPIC_GRADIENT_FLUX
         COMPUTERHSPBG
         COMPUTERHSTHREEPHASES

source/CMakeLists.txt

@@ -127,6 +127,7 @@ set(SOURCES ${CMAKE_SOURCE_DIR}/source/KKStools.cc
             ${CMAKE_SOURCE_DIR}/source/PhaseFluxStrategySimple.cc
             ${CMAKE_SOURCE_DIR}/source/PhaseFluxStrategyIsotropic.cc
             ${CMAKE_SOURCE_DIR}/source/PhaseFluxStrategyAnisotropy.cc
+            ${CMAKE_SOURCE_DIR}/source/PhaseFluxStrategyAnisotropyMultiOrder.cc
             ${CMAKE_SOURCE_DIR}/source/TemperatureStrategyFactory.cc
             ${CMAKE_SOURCE_DIR}/source/QuatModel.cc
             ${CMAKE_SOURCE_DIR}/source/QuatRefinePatchStrategy.cc

source/PhaseFluxStrategyAnisotropyMultiOrder.cc

@@ -0,0 +1,73 @@
+// Copyright (c) 2018, Lawrence Livermore National Security, LLC and
+// UT-Battelle, LLC.
+// Produced at the Lawrence Livermore National Laboratory and
+// the Oak Ridge National Laboratory
+// LLNL-CODE-747500
+// All rights reserved.
+// This file is part of AMPE.
+// For details, see https://github.com/LLNL/AMPE
+// Please also read AMPE/LICENSE.
+//
+#include "PhaseFluxStrategyAnisotropyMultiOrder.h"
+#include "QuatFort.h"
+
+#include "SAMRAI/geom/CartesianPatchGeometry.h"
+#include "SAMRAI/pdat/CellData.h"
+#include "SAMRAI/pdat/SideData.h"
+#include "SAMRAI/math/PatchSideDataNormOpsReal.h"
+
+void PhaseFluxStrategyAnisotropyMultiOrder::computeFluxes(
+    const std::shared_ptr<hier::PatchLevel> level, const int phase_id,
+    const int quat_id, const int flux_id)
+{
+   (void)quat_id;  // uses data member d_quat instead
+
+   //  Compute phase "flux" on patches in level.
+   for (hier::PatchLevel::Iterator ip(level->begin()); ip != level->end();
+        ++ip) {
+      std::shared_ptr<hier::Patch> patch = *ip;
+
+      const std::shared_ptr<geom::CartesianPatchGeometry> patch_geom(
+          SAMRAI_SHARED_PTR_CAST<geom::CartesianPatchGeometry,
+                                 hier::PatchGeometry>(
+              patch->getPatchGeometry()));
+      TBOX_ASSERT(patch_geom);
+
+      const double* dx = patch_geom->getDx();
+
+      std::shared_ptr<pdat::CellData<double> > phase(
+          SAMRAI_SHARED_PTR_CAST<pdat::CellData<double>, hier::PatchData>(
+              patch->getPatchData(phase_id)));
+      TBOX_ASSERT(phase);
+
+      std::shared_ptr<pdat::SideData<double> > phase_flux(
+          SAMRAI_SHARED_PTR_CAST<pdat::SideData<double>, hier::PatchData>(
+              patch->getPatchData(flux_id)));
+      TBOX_ASSERT(phase_flux);
+
+      const hier::Box& pbox = patch->getBox();
+      const hier::Index& ifirst = pbox.lower();
+      const hier::Index& ilast = pbox.upper();
+
+      for (int i = 0; i < phase->getDepth(); i++)
+         ANISOTROPIC_FLUX(ifirst(0), ilast(0), ifirst(1), ilast(1),
+#if (NDIM == 3)
+                          ifirst(2), ilast(2),
+#endif
+                          dx, d_epsilon_phase, d_nu, d_knumber,
+                          phase->getPointer(i), phase->getGhostCellWidth()[0],
+                          d_quat[i].data(), d_quat[i].size(),
+                          phase_flux->getPointer(0, i),
+                          phase_flux->getPointer(1, i),
+#if (NDIM == 3)
+                          phase_flux->getPointer(2, i),
+#endif
+                          phase_flux->getGhostCellWidth()[0]);
+
+#ifdef DEBUG_CHECK_ASSERTIONS
+      SAMRAI::math::PatchSideDataNormOpsReal<double> ops;
+      double l2f = ops.L2Norm(phase_flux, pbox);
+      assert(l2f == l2f);
+#endif
+   }
+}

source/PhaseFluxStrategyAnisotropyMultiOrder.h

@@ -0,0 +1,45 @@
+// Copyright (c) 2018, Lawrence Livermore National Security, LLC and
+// UT-Battelle, LLC.
+// Produced at the Lawrence Livermore National Laboratory and
+// the Oak Ridge National Laboratory
+// LLNL-CODE-747500
+// All rights reserved.
+// This file is part of AMPE.
+// For details, see https://github.com/LLNL/AMPE
+// Please also read AMPE/LICENSE.
+//
+//
+#ifndef included_PhaseFluxStrategyAnisotropyMultiOrder
+#define included_PhaseFluxStrategyAnisotropyMultiOrder
+
+#include "PhaseFluxStrategy.h"
+
+#include <vector>
+#include <array>
+
+class PhaseFluxStrategyAnisotropyMultiOrder : public PhaseFluxStrategy
+{
+ public:
+   PhaseFluxStrategyAnisotropyMultiOrder(
+       const double epsilon_phase, const double nu, const int knumber,
+       std::vector<std::array<double, 4>> quat)
+       : d_epsilon_phase(epsilon_phase),
+         d_nu(nu),
+         d_knumber(knumber),
+         d_quat(quat)
+   {
+      tbox::plog << "PhaseFluxStrategyAnisotropyMultiOrder..." << std::endl;
+   }
+
+   void computeFluxes(const std::shared_ptr<hier::PatchLevel> level,
+                      const int phase_id, const int quat_id, const int flux_id);
+
+ private:
+   const double d_epsilon_phase;
+   const double d_nu;
+   const int d_knumber;
+
+   std::vector<std::array<double, 4>> d_quat;
+};
+
+#endif

source/PhaseFluxStrategyFactory.h

@@ -4,6 +4,7 @@
 #include "PhaseFluxStrategyAnisotropy.h"
 #include "PhaseFluxStrategyIsotropic.h"
 #include "PhaseFluxStrategySimple.h"
+#include "PhaseFluxStrategyAnisotropyMultiOrder.h"
 
 class PhaseFluxStrategyFactory
 {
@@ -35,8 +36,15 @@ class PhaseFluxStrategyFactory
       if (epsilon_anisotropy >= 0.) {
          if (!model_parameters.with_orientation())
             TBOX_ERROR("Phase anisotropy requires quaternion orientation");
-         phase_flux_strategy.reset(
-             new PhaseFluxStrategyAnisotropy(epsilon, epsilon_anisotropy, 4));
+         if (model_parameters.norientations() > 0) {
+            phase_flux_strategy.reset(new PhaseFluxStrategyAnisotropyMultiOrder(
+                epsilon, epsilon_anisotropy, 4,
+                model_parameters.orientations()));
+         } else {
+            phase_flux_strategy.reset(
+                new PhaseFluxStrategyAnisotropy(epsilon, epsilon_anisotropy,
+                                                4));
+         }
       } else if (model_parameters.useIsotropicStencil()) {
          phase_flux_strategy.reset(new PhaseFluxStrategyIsotropic(epsilon));
       } else {

source/QuatModelParameters.cc

@@ -866,6 +866,10 @@ void QuatModelParameters::readModelParameters(
    tbox::plog << "norderp_B = " << d_norderp_B << std::endl;
    tbox::plog << "norderp = " << d_norderp << std::endl;
 
+   if (d_norderp > 1) {
+      readOrientations(model_db);
+   }
+
    // Set d_H_parameter to negative value, to turn off orientation terms
    d_H_parameter = model_db->getDoubleWithDefault("H_parameter", -1.);
 
@@ -1256,3 +1260,17 @@ double QuatModelParameters::quatMobilityScaleFactor() const
 
    return tbox::IEEE::getSignalingNaN();
 }
+
+void QuatModelParameters::readOrientations(
+    std::shared_ptr<tbox::Database> model_db)
+{
+   std::shared_ptr<tbox::Database> db(model_db->getDatabase("Orientations"));
+
+   for (int i = 0; i < d_norderp_A + d_norderp_B; i++) {
+      double quat[4];
+      std::string name = "quat" + std::to_string(i);
+      db->getDoubleArray(name, &quat[0], 4);
+      std::array<double, 4> qa{quat[0], quat[1], quat[2], quat[3]};
+      d_orientations.push_back(qa);
+   }
+}

source/QuatModelParameters.h

@@ -124,7 +124,7 @@ class QuatModelParameters
       assert(isp < d_cp.size());
       return d_cp[isp];
    }
-   const std::vector<std::map<short, double> >& cp() const { return d_cp; }
+   const std::vector<std::map<short, double>>& cp() const { return d_cp; }
 
    int ncompositions() const
    {
@@ -546,6 +546,15 @@ class QuatModelParameters
              d_ceq0_solidA[2] * dT * dT;
    }
 
+   int norientations() const { return d_orientations.size(); }
+
+   std::array<double, 4> orientation(const int i) const
+   {
+      return d_orientations[i];
+   }
+
+   std::vector<std::array<double, 4>> orientations() { return d_orientations; }
+
  private:
    void readNumberSpecies(std::shared_ptr<tbox::Database> conc_db);
 
@@ -588,7 +597,7 @@ class QuatModelParameters
    double d_thermal_diffusivity;
    double d_latent_heat;
    // cp for each species
-   std::vector<std::map<short, double> > d_cp;
+   std::vector<std::map<short, double>> d_cp;
    double d_meltingT;
    double d_interface_mobility;
 
@@ -799,8 +808,14 @@ class QuatModelParameters
    double d_ceq0_liquid[3];
    double d_ceq0_solidA[3];
 
+   /*!
+    * Order parameters orientations
+    */
+   std::vector<std::array<double, 4>> d_orientations;
+
    void readMolarVolumes(std::shared_ptr<tbox::Database> db);
    void readEquilibriumCompositions(std::shared_ptr<tbox::Database> conc_db);
+   void readOrientations(std::shared_ptr<tbox::Database> model_db);
 
    void readCahnHilliard(std::shared_ptr<tbox::Database> db);
    void readWangSintering(std::shared_ptr<tbox::Database> db);

source/fortran/2d/quatrhs.m4

@@ -150,6 +150,92 @@ c
       return
       end
 
+c***********************************************************************
+      subroutine anisotropic_flux(
+     &   ifirst0, ilast0, ifirst1, ilast1,
+     &   h, epsilon, nu, knumber,
+     &   phase, ngphase,
+     &   quat, qlen,
+     &   flux0, flux1, ngflux)
+
+      implicit none
+      integer ifirst0, ilast0, ifirst1, ilast1,
+     &     ngphase, ngq, qlen, ngflux, knumber
+      double precision
+     &     phase(CELL2d(ifirst,ilast,ngphase)),
+     &     quat(qlen),
+     &     flux0(SIDE2d0(ifirst,ilast,ngflux)),
+     &     flux1(SIDE2d1(ifirst,ilast,ngflux)),
+     &     nu, h(2)
+
+c     local variables
+      integer i, j
+      double precision dxinv, dyinv, epsilon
+      double precision theta, phi
+      double precision pi, q
+      double precision epsilon2, dphidx, dphidy
+      double precision epstheta, depsdtheta
+c
+      pi = 4.d0*atan(1.d0)
+c
+      epsilon2 = epsilon * epsilon
+
+      dxinv = 1.d0 / h(1)
+      dyinv = 1.d0 / h(2)
+
+      q=quat(1)
+      if( qlen==4 )then
+         phi=2.d0*acos(q)
+      else
+         phi=acos(q)
+      endif
+
+c x faces
+      do j = ifirst1, ilast1
+         do i = ifirst0, ilast0+1
+            dphidx = (phase(i,j) - phase(i-1,j)) * dxinv
+            dphidy = 0.25d0*(phase(i-1,j+1) - phase(i-1,j-1)
+     &                   + phase(i  ,j+1) - phase(i  ,j-1))
+     &                   * dyinv
+            if( abs(dphidx)>1.e-12 )then
+               theta=atan(dphidy/dphidx)
+            else
+               theta=0.5d0*pi
+            endif
+
+            epstheta=epsilon*(1.d0+nu*cos(knumber*(theta-phi)))
+            depsdtheta=-knumber*epsilon*nu*sin(knumber*(theta-phi))
+
+            flux0(i,j) = epstheta*epstheta*dphidx
+     &                 - epstheta*depsdtheta*dphidy
+         enddo
+      enddo
+
+c y faces
+      do j = ifirst1, ilast1+1
+         do i = ifirst0, ilast0
+            dphidx = 0.25*(phase(i+1,j-1) - phase(i-1,j-1)
+     &                   + phase(i+1,j  ) - phase(i-1,j  ))
+     &                   * dxinv
+            dphidy = (phase(i,j) - phase(i,j-1)) * dyinv
+
+            if( abs(dphidx)>1.e-12 )then
+               theta=atan(dphidy/dphidx)
+            else
+               theta=0.5*pi
+            endif
+
+            epstheta=epsilon*(1.+nu*cos(knumber*(theta-phi)))
+            depsdtheta=-knumber*epsilon*nu*sin(knumber*(theta-phi))
+
+            flux1(i,j) = epstheta*epstheta*dphidy
+     &                 + epstheta*depsdtheta*dphidx
+         enddo
+      enddo
+
+      return
+      end
+
 c***********************************************************************
       subroutine anisotropic_gradient_flux(
      &   ifirst0, ilast0, ifirst1, ilast1,

source/fortran/QuatFort.h

@@ -87,6 +87,20 @@ void ANISOTROPIC_GRADIENT_FLUX(const int& ifirst0, const int& ilast0,
 #endif
                                const int& fluxng);
 
+void ANISOTROPIC_FLUX(const int& ifirst0, const int& ilast0, const int& ifirst1,
+                      const int& ilast1,
+#if (NDIM == 3)
+                      const int& ifirst2, const int& ilast2,
+#endif
+                      const double* dx, const double& epsilon, const double& nu,
+                      const int& knumber, double* phase, const int& phaseng,
+                      double* quat, const int& qlen, double* flux0,
+                      double* flux1,
+#if (NDIM == 3)
+                      double* flux2,
+#endif
+                      const int& fluxng);
+
 void COMPUTERHSPBG(const int& ifirst0, const int& ilast0, const int& ifirst1,
                    const int& ilast1,
 #if (NDIM == 3)

tests/AlCuParabolicMultiOrder/2d.input

@@ -27,6 +27,11 @@ ModelParameters {  // required block
 
    epsilon_anisotropy = 0.04 //delta
 
+   Orientations{
+      quat0 = 0,0,0,1
+      quat1 = 1,0,0,0
+   }
+
    Interface {
       sigma = 0.8  // pJ/um^2 = J/m^2
       delta = 0.064 // um