diff --git a/CMakeLists.txt b/CMakeLists.txt index 3bc875af..33f65c69 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -184,6 +184,7 @@ FortranCInterface_HEADER( STORAGE_INDEX_STIFFNESS GRADIENT_FLUX COMPUTE_FLUX_ISOTROPIC + ANISOTROPIC_FLUX ANISOTROPIC_GRADIENT_FLUX COMPUTERHSPBG COMPUTERHSTHREEPHASES diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index 0f88b6e3..1269b976 100644 --- a/source/CMakeLists.txt +++ b/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 diff --git a/source/PhaseFluxStrategyAnisotropyMultiOrder.cc b/source/PhaseFluxStrategyAnisotropyMultiOrder.cc new file mode 100644 index 00000000..d0ef7354 --- /dev/null +++ b/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 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 patch = *ip; + + const std::shared_ptr patch_geom( + SAMRAI_SHARED_PTR_CAST( + patch->getPatchGeometry())); + TBOX_ASSERT(patch_geom); + + const double* dx = patch_geom->getDx(); + + std::shared_ptr > phase( + SAMRAI_SHARED_PTR_CAST, hier::PatchData>( + patch->getPatchData(phase_id))); + TBOX_ASSERT(phase); + + std::shared_ptr > phase_flux( + SAMRAI_SHARED_PTR_CAST, 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 ops; + double l2f = ops.L2Norm(phase_flux, pbox); + assert(l2f == l2f); +#endif + } +} diff --git a/source/PhaseFluxStrategyAnisotropyMultiOrder.h b/source/PhaseFluxStrategyAnisotropyMultiOrder.h new file mode 100644 index 00000000..21cea0dd --- /dev/null +++ b/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 +#include + +class PhaseFluxStrategyAnisotropyMultiOrder : public PhaseFluxStrategy +{ + public: + PhaseFluxStrategyAnisotropyMultiOrder( + const double epsilon_phase, const double nu, const int knumber, + std::vector> 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 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> d_quat; +}; + +#endif diff --git a/source/PhaseFluxStrategyFactory.h b/source/PhaseFluxStrategyFactory.h index 89bad79a..583971fc 100644 --- a/source/PhaseFluxStrategyFactory.h +++ b/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 { diff --git a/source/QuatModelParameters.cc b/source/QuatModelParameters.cc index daef2b73..3189cfef 100644 --- a/source/QuatModelParameters.cc +++ b/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 model_db) +{ + std::shared_ptr 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 qa{quat[0], quat[1], quat[2], quat[3]}; + d_orientations.push_back(qa); + } +} diff --git a/source/QuatModelParameters.h b/source/QuatModelParameters.h index a5df9847..3acb8181 100644 --- a/source/QuatModelParameters.h +++ b/source/QuatModelParameters.h @@ -124,7 +124,7 @@ class QuatModelParameters assert(isp < d_cp.size()); return d_cp[isp]; } - const std::vector >& cp() const { return d_cp; } + const std::vector>& 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 orientation(const int i) const + { + return d_orientations[i]; + } + + std::vector> orientations() { return d_orientations; } + private: void readNumberSpecies(std::shared_ptr conc_db); @@ -588,7 +597,7 @@ class QuatModelParameters double d_thermal_diffusivity; double d_latent_heat; // cp for each species - std::vector > d_cp; + std::vector> 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> d_orientations; + void readMolarVolumes(std::shared_ptr db); void readEquilibriumCompositions(std::shared_ptr conc_db); + void readOrientations(std::shared_ptr model_db); void readCahnHilliard(std::shared_ptr db); void readWangSintering(std::shared_ptr db); diff --git a/source/fortran/2d/quatrhs.m4 b/source/fortran/2d/quatrhs.m4 index 07c31cdc..b5833789 100644 --- a/source/fortran/2d/quatrhs.m4 +++ b/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, diff --git a/source/fortran/QuatFort.h b/source/fortran/QuatFort.h index 94cc0f17..177b56ce 100644 --- a/source/fortran/QuatFort.h +++ b/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) diff --git a/tests/AlCuParabolicMultiOrder/2d.input b/tests/AlCuParabolicMultiOrder/2d.input index d7f33776..62d9486c 100644 --- a/tests/AlCuParabolicMultiOrder/2d.input +++ b/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 diff --git a/tests/AlCuParabolicMultiOrder/test2d.py b/tests/AlCuParabolicMultiOrder/test2d.py index 8d1ef820..0e93fda1 100644 --- a/tests/AlCuParabolicMultiOrder/test2d.py +++ b/tests/AlCuParabolicMultiOrder/test2d.py @@ -5,18 +5,27 @@ print("Test AlCu parabolic multi-order...") +mpicmd = sys.argv[1]+" "+sys.argv[2]+" "+sys.argv[3] +exe = sys.argv[4] +inp = sys.argv[5] +datadir = sys.argv[6] + +data = "2spheres.csv" +src = datadir+'/'+data +try: + os.symlink(src, data) +except FileExistsError: + os.remove(data) + os.symlink(src, data) + #prepare initial conditions file initfilename="160x160.nc" -subprocess.call(["python3", "../../utils/make_nuclei.py", - "--nx", "160", "--ny", "160", "--nz", "1", "-r", "15", - "--center0", "0, 0, 0", - "--concentration-in", "0.03", "--concentration-out", "0.08", +subprocess.call(["python3", "../../utils/make_multi_spheres.py", + "--nx", "160", "--ny", "160", "--nz", "1", + "--spheres", data, + "--concentration-A", "0.03", "--concentration-out", "0.08", initfilename]) -mpicmd = sys.argv[1]+" "+sys.argv[2]+" "+sys.argv[3] -exe = sys.argv[4] -inp = sys.argv[5] - #run AMPE command = "{} {} {}".format(mpicmd,exe,inp) output = subprocess.check_output(command,shell=True) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 95b4baf2..6ef5a038 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -375,6 +375,7 @@ add_test(NAME AlCuParabolicMultiOrder ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${CMAKE_CURRENT_BINARY_DIR}/../source/ampe${NDIM}d ${CMAKE_CURRENT_SOURCE_DIR}/AlCuParabolicMultiOrder/${NDIM}d.input + ${CMAKE_CURRENT_SOURCE_DIR}/AlCuParabolicMultiOrder ) add_test(NAME AlCuCalphadMultiOrder COMMAND ${Python3_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/AlCuCalphadMultiOrder/test${NDIM}d.py