Upgrade to Pro — share decks privately, control downloads, hide ads and more …

SymPhas: Symbolic Algebra for Phase-Field Simulations

SymPhas: Symbolic Algebra for Phase-Field Simulations

97d945680ed363e4cce48666d41c586e?s=128

Daniel Wheeler

July 21, 2022
Tweet

More Decks by Daniel Wheeler

Other Decks in Science

Transcript

  1. SymPhas: Symbolic Algebra for Phase-Field Simulations Steven Silber Dr. Mikko

    Karttunen Western University, London, ON, Canada 1
  2. Agenda 1. SymPhas 2. Problem Representation 3. Symbolic Algebra 4.

    Simulation Examples 2 S. A. Silber and M. Karttunen, “SymPhas —General Purpose Software for Phase- Field, Phase-Field Crystal, and Reaction- Diffusion Simulations,” Adv. Theory Simul., vol. 5, Art. no. 1, Dec. 2021.
  3. What is SymPhas? • C++ API & software • Numerical

    solutions to phase-field problems • P.F. Model → system of PDEs (1st order in time) • Symbolic algebra manipulation of phase-field models • User-defined P.F. Model → solver-enabled representation 3
  4. What is SymPhas? • Phase-field problems: • 1, 2, or

    3 dimensions • Any number of order parameters • Real-, complex- or vector- valued order parameters • Phase-field, phase-field crystal or reaction-diffusion • I.e. User-defined P.F. Models 4
  5. Objectives of SymPhas • Accessible • …to non-developers & those

    unfamiliar with C++ E.g. easy-to-use API + extensive documentation (Doxygen) • Efficient/high performance • Parallelization: Single-CPU • Efficient code design 5
  6. Objectives of SymPhas 1. Unconstrained definition of phase-field models 2.

    Comprehensive symbolic algebra functionality 3. Modular numerical solvers 6
  7. 4 basic modules: 2 additional modules: → lib → sym

    → io → datatypes → sol → conf • Module selection and compilation process: → CMake build process allows the user to: → Select modules & compilation parameters for SymPhas → Install SymPhas library to OS → Import SymPhas to other CMake projects SymPhas API & Implementation 7
  8. SymPhas API & Implementation • User defined models and workflow:

    User-defined driver file using SymPhas API (as few as 22 lines) Phase-field models implemented via macro-based grammar User-implemented solvers simplified with C++ macros • Included in software: • Driver files (including a “smallest working” implementation) • Model definitions (phase-field & phase-field crystal) • Solvers (forward Euler + semi-implicit Fourier spectral) 8
  9. SymPhas API & Implementation • Design with best practises and

    maintenance-friendly techniques → Reducing coupling/increasing cohesion → Single responsibility principle → Common namespaces → Informative function names → Data groups → API specific type aliases + Extensive documentation (using Doxygen) Object design User facing functionality 9
  10. SymPhas API & Implementation • Design implements high performance features:

    → Parallelization with OpenMP & C++ execution header → Layered object design to streamline & minimize data passing → Minimal branching for all P.F. model simulations → Template metaprogramming 10
  11. 2 | Problem Encapsulation 11

  12. Problem Representation • C++ → Object Oriented Programming (OOP) •

    Generic programming AKA template metaprogramming → Everything is a type …As opposed to inheritance and virtual dispatch → C++ templates and generics → “Template argument deduction” 12
  13. Problem Representation Uniform Grid Uniform Grid Uniform Grid 13

  14. Problem Representation System Uniform Grid System Uniform Grid System Uniform

    Grid 14
  15. Group of Systems Problem Representation System Uniform Grid System Uniform

    Grid System Uniform Grid 15
  16. Group of Systems Problem Representation System Uniform Grid System Uniform

    Grid System Uniform Grid Equation of Motion Equation of Motion Equation of Motion 16
  17. Group of Systems Problem Representation System Uniform Grid System Uniform

    Grid System Uniform Grid Equation of Motion Equation of Motion Equation of Motion 𝐸(𝜓1 , 𝜓2 , 𝜓3 ) 17
  18. Group of Systems Problem Representation System Uniform Grid System Uniform

    Grid System Uniform Grid Equation of Motion Equation of Motion Equation of Motion 𝐸(𝜓1 , 𝜓2 , 𝜓3 ) 𝜓1 𝜓2 𝜓3 18
  19. Working Data Group of Systems Problem Representation Working Data Working

    Data Working Data System System System Uniform Grid Uniform Grid Uniform Grid Equation of Motion Equation of Motion Equation of Motion 𝐸(𝜓1 , 𝜓2 , 𝜓3 ) 𝜓1 𝜓2 𝜓3 Solver 19
  20. Working Data Group of Systems Problem Encapsulation Working Data Working

    Data Working Data System Uniform Grid System Uniform Grid System Uniform Grid Equation of Motion Equation of Motion Equation of Motion 𝐸(𝜓1 , 𝜓2 , 𝜓3 ) 𝜓1 𝜓2 𝜓3 Solver 20
  21. Model Group of Systems Problem Representation Working Data Working Data

    Working Data System Uniform Grid System Uniform Grid System Uniform Grid Equation of Motion Equation of Motion Working Data Equation of Motion 𝐸(𝜓1 , 𝜓2 , 𝜓3 ) 𝜓1 𝜓2 𝜓3 Solver 21
  22. Problem Representation • OOP provides flexibility: • Different Solver implementations

    (incl. user-defined solvers) • Different System specializations based on Solver type • Different Grid specializations based on System type • Different Model types for phase-field and phase-field crystal problems 22
  23. Problem Representation • User-defined phase-field models: 23 MODEL(GRAYSCOTT, (SCALAR, SCALAR),

    MODEL_DEF_PREAMBLE((auto prr = psi * rho * rho;), dpsi = c1 * lap(psi) - prr + c3 * (one - psi), drho = c2 * lap(rho) + prr - (c3 + c4) * rho) ) LINK_WITH_NAME(GRAYSCOTT, GRAYSCOTT) MODEL(MB, (SCALAR), MODEL_DEF( dpsi = -bilap(psi) - lap((c1 - c2 * psi * psi) * psi)) ) LINK_WITH_NAME(MB, MODELB) MODEL(MC, (SCALAR, SCALAR), MODEL_DEF_PREAMBLE( ( auto psi3 = psi * psi * psi; auto rho3 = rho * rho * rho; ), dpsi = lap(psi) + c1 * psi - c2 * psi3 + lit(2.) * c5 * psi * rho, drho = -bilap(rho) - lap(c3 * rho - c4 * rho3 + c5 * psi * psi)) ) LINK_WITH_NAME(MC, MODELC) DEFINE_MODEL_FIELD_NAMES(MC, ("psi", "m")) MODEL(MA, (SCALAR), MODEL_DEF( dpsi = lap(psi) + (c1 - c2 * psi * psi) * psi) ) LINK_WITH_NAME(MA, MODELA)
  24. 3 | Symbolic Algebra 24

  25. Symbolic Algebra • Representation of equation of motion? → Of

    the form… • How the equation will be used: → Evaluated multiple times (e.g., each solver loop) → Modified to fit the design of a given numerical solver 25
  26. Symbolic Algebra Expression tree • Expression trees fulfill these needs:

    evaluable & modifiable → Recursive evaluation beginning at root → Each node & subtree is also an expression tree = 1 = 2 + 1 = 3 = 6 ÷ 3 = 2 Recursive evaluation + + 𝜓 8 × 4 1 2 ÷ 6 + 26
  27. Symbolic Algebra • Expression tree which is nonbranching at runtime

    → virtual dispatch deterministic control flow → Deterministic control ⇨ object types are fixed • Deterministic control flow via expression templates → ⇓ Runtime ⇑ Compile Time Since: C++ is a compiled language, Then: Types must be known at point of invocation … notion of “compile-time constant” expression trees 27
  28. Symbolic Algebra • Implementation: → Curiously Recurring Template Pattern (CRTP)

    → Defines an interface that every specialization must implement • Basic terms: Node term → OpExpression Variable term → Data instance (i.e. order parameter) Literal term → A number Binary term → An operation 28
  29. Symbolic Algebra • Binary Operators → OpBinaryAdd → OpBinarySub →

    OpBinaryMul → OpBinaryDiv • Implements operator overloads • Order of operations of overloaded operator is always followed 29
  30. Symbolic Algebra • Common terms for an equation of motion:

    → OpLiteral 𝑎 𝑛𝑢𝑚𝑏𝑒𝑟 → OpLVariable 𝜓 → OpNLVariable 𝜓 ∗ 𝜓 ∗ ⋯ = 𝜓𝑛 → OpFuncDerivative ∇2(𝑎𝑛 𝑒𝑥𝑝𝑟𝑒𝑠𝑠𝑖𝑜𝑛) • More control over representing derivatives (PFC): → OpOperator (1 + ∇2) ⋅ (𝑎𝑛 𝑒𝑥𝑝𝑟𝑒𝑠𝑠𝑖𝑜𝑛) 30
  31. Cast to “recurring” or inherited type, E template<typename E> struct

    OpExpression { auto eval(int n) const { return cast().eval(n); } void print(char* out) const { return cast().print(out); } auto& cast() const { return *static_cast<E const*>(this); } ... Symbolic Algebra CRTP Base Class 31
  32. Symbolic Algebra CRTP Specialization template<typename T, typename G> struct OpLVariable

    : OpExpression<OpLVariable<T, G>> { T value; G data; OpLVariable(T value, G data) : value{ value }, data{ data } {} auto eval(iter_type n) const { return value * BaseData<G>::evaluate(data, n); } ... 32
  33. Symbolic Algebra Return type selection template<typename T, typename G> struct

    OpLVariable : OpExpression<OpLVariable<T, G>> { T value; G data; OpLVariable(T value, G data) : value{ value }, data{ data } {} auto eval(iter_type n) const { return value * BaseData<G>::evaluate(data, n); } ... 33
  34. Symbolic Algebra • The CRTP is used with template argument

    deduction template<typename E1, typename E2> auto add(OpExpression<E1> const& a, OpExpression<E2> const& b) { ... } Template argument deduction 34
  35. Symbolic Algebra • Library defines: → Distribution of multiplied terms

    → Collection of like terms → Multiplication by 0 or 1 → Addition/Subtraction by 0 → … 35
  36. Symbolic Algebra • Library defines: → … → Defining a

    function with arguments → Factoring to cancel terms in a division → Transmuting terms → “Pretty printing” of expressions → Custom symbolic algebra terms 36
  37. Symbolic Algebra • Custom symbolic algebra terms 37 DEFINE_BASE_TYPE_CONDITION((typename K_t),

    (K_t), expr::is_K_type<K_t>::value, K<order_K_type<K_t>::value>) DEFINE_SYMBOL_ID_CONDITION((typename K_t), (K_t), (expr::is_K_type<K_t>::value), { static char* name = new_k_name(order_K_type<K_t>::value); return name; }); ALLOW_COMBINATION_CONDITION((typename K_t), (K_t), expr::is_K_type<K_t>::value); RESTRICT_MULTIPLICATION_CONDITION((typename K_t), (K_t), expr::is_K_type<K_t>::value); ALLOW_DIVISION_CONDITION((size_t O1, size_t O2), (K<O1>, K<O2>), O2 > O1); template<size_t O> template<typename T, size_t D> struct K<O>::WaveVectorData : Grid<T, D>, K<O> { ... }; New object
  38. Symbolic Algebra • Custom symbolic algebra terms 38 DEFINE_BASE_TYPE_CONDITION((typename K_t),

    (K_t), expr::is_K_type<K_t>::value, K<order_K_type<K_t>::value>) template<size_t O> template<typename T, size_t D> struct K<O>::WaveVectorData : Grid<T, D>, K<O> { ... }; “Type identity” of the object: → i.e. the “class” or “identity” of the object in the symbolic algebra
  39. Symbolic Algebra • Custom symbolic algebra terms 39 DEFINE_SYMBOL_ID_CONDITION((typename K_t),

    (K_t), (expr::is_K_type<K_t>::value), { static char* name = new_k_name(order_K_type<K_t>::value); return name; }); template<size_t O> template<typename T, size_t D> struct K<O>::WaveVectorData : Grid<T, D>, K<O> { ... }; Character name/representation of the object when printed
  40. Symbolic Algebra • Custom symbolic algebra terms 40 ALLOW_COMBINATION_CONDITION((typename K_t),

    (K_t), expr::is_K_type<K_t>::value); RESTRICT_MULTIPLICATION_CONDITION((typename K_t), (K_t), expr::is_K_type<K_t>::value); ALLOW_DIVISION_CONDITION((size_t O1, size_t O2), (K<O1>, K<O2>), O2 > O1); template<size_t O> template<typename T, size_t D> struct K<O>::WaveVectorData : Grid<T, D>, K<O> { ... }; Rules of symbolic algebra
  41. 4 | Simulation Examples 41

  42. Simulation Examples • Model A (Allen-Cahn) • Model B (Cahn-Hilliard)

    • Model C (coupled A & B) • PFC (conserved dynamics) • Validate dynamics of A & B using Porod’s Law • Structure Factor ⇨ measures the amount of “solidification” 𝑆 𝑘 ~ 𝑘𝑑+1 −1 Solved using spectral solver 42
  43. Simulation Examples • Runtime tests show: • Good scaling between

    2D & 3D 43 S. A. Silber and M. Karttunen, “SymPhas —General Purpose Software for Phase-Field, Phase-Field Crystal, and Reaction-Diffusion Simulations,” Adv. Theory Simul., vol. 5, Art. no. 1, Dec. 2021.
  44. Model C MODEL(MC, (SCALAR, SCALAR), MODEL_PREAMBLE_DEF( ( auto psi3 =

    c2 * psi * psi * psi, auto m3 = c4 * m * m * m) ), dpsi = lap(psi) + c1 * psi - psi3 + lit(2.) * c5 * psi * m, dm = -bilap(m) - lap(c3 * m - m3 + c5 * psi * psi))) coupled conserved and non-conserved dynamics; eutectic growth 44
  45. Model C 2D coupled conserved and non-conserved dynamics; eutectic growth

    45
  46. Model C 3D coupled conserved and non-conserved dynamics; eutectic growth

    4,000 20,000 4,000 20,000 46
  47. Phase-Field Crystal Model (PFC) PFC_TYPE(PC, DEFAULTS( DEFAULT_DYNAMIC(PFC_CONSERVED) ), (SCALAR)) Introduces

    elasticity; atomistic model 47
  48. PFC 2D 700,000 iterations ≈ 700 diffusion time lengths 48

  49. PFC 2D 500 1,000 2,000 4,000 8,000 16,000 32,000 64,000

    49
  50. PFC 3D 70,000 iterations 50

  51. Model A non-conserved dynamics; chemical potential driven phase transition MODEL(MA,

    (SCALAR), MODEL_DEF( dpsi = lap(psi) + (c1 - c2 * psi * psi) * psi)) 51
  52. Model A 2D non-conserved dynamics; chemical potential driven phase transition

    52
  53. Model A 3D non-conserved dynamics; chemical potential driven phase transition

    400 2000 800 53
  54. Model A Dynamics 54

  55. Model B MODEL(MB, (SCALAR), MODEL_DEF( dm = -bilap(m) - lap(c1

    - c2 * m * m) * m)) conserved dynamics; chemical potential driven phase transition 55
  56. Model B 2D conserved dynamics; chemical potential driven phase transition

    56
  57. Model B 3D conserved dynamics; chemical potential driven phase transition

    4,000 20,000 8,000 57
  58. Model B Dynamics 58

  59. Acknowledgements • Western University • NSERC CGSM Scholarship • Mitacs

    GlobalLink Travel Grant • Aalto University 59
  60. Appendix Slides 60

  61. • Expression trees fulfill these needs … a construct that

    is both evaluable & modifiable Symbolic Algebra 1 2 + + ÷ 6 𝜓 8 × + 4 61
  62. SymPhas API & Implementation • A SymPhas implementation in 22

    lines: #include "symphas.h" MODEL(EX, (SCALAR), MODEL_DEF(dop(1) = lap(op(1)) + (c1 - lit (4.) * c2 * op(1) * op(1)) * op(1))) int main(int argc , char* argv []) { double dt = 0.1; symphas::problem_parameters_type pp{ 1 }; symphas::b_data_type bdata; symphas::init_data_type tdata{ Inside::UNIFORM, { -1, 1 } }; symphas::interval_data_type vdata; symphas::interval_element_type interval; interval.set_interval_count(0, 128, 128); bdata[Side::LEFT] = BoundaryType::PERIODIC; bdata[Side::RIGHT] = BoundaryType::PERIODIC; bdata[Side::TOP] = BoundaryType::PERIODIC; bdata[Side::BOTTOM] = BoundaryType::PERIODIC; vdata[Axis::X] = interval; vdata[Axis::Y] = interval; pp.set_boundary_data(&bdata); pp.set_initial_data(&tdata); pp.set_interval_data(&vdata); pp.set_problem_time_step(dt); model_EX_t<2, SolverSP<Stencil2d2h<5, 9, 6>>> model{ pp }; symphas::find_solution(model, dt, 50);} Define new Model . ⊳ Solve phase-field model ⊳ Boundary conditions ⊳ Initial conditions … ⊳ System size … . ⊳ Initialize problem … ⊳ 62
  63. Problem Representation Uniform Grid Uniform Grid Uniform Grid 63

  64. Problem Representation System Uniform Grid System Uniform Grid System Uniform

    Grid 64 *Changes based on use of io module
  65. 4 | Numerical Solvers 65

  66. Numerical Solvers • Phase-field data, 𝜓, represented by a regular

    grid → Cartesian coordinate system,Δ𝑥 (× Δ𝑦 (× Δ𝑧)) → Each element represents the average of the phase-field Numerical solution in 1D denoted: 𝜓𝑗 𝑛 ≔ 𝜓 𝑥 = 𝑗Δ𝑥,𝑡 = 𝑛Δ𝑡 • SymPhas implements two numerical solvers: → Forward-Time Central-Space (FTCS) (Forward Euler) Solver → Semi-implicit Fourier spectral solver Δ𝑦 Δ𝑥 Δz 𝑥 𝑦 𝑧 66
  67. Numerical Solvers: Forward Euler • Derivatives approximated with a finite

    difference stencil → A linear combination of grid elements around a given point → Stencil is characterized by: Dimension of the system Order of accuracy of the approximation Derivative order that is approximated (Also number of grid points in the stencil) ∇2𝜓 2D uniform grid 𝒪(ℎ4) 𝒪(ℎ2) 𝚫𝒙 𝚫𝒚 𝚫𝐱 = 𝚫𝐲 = 𝒉 67
  68. Numerical Solvers: Forward Euler & periodic BC • Approximation of

    𝝏 𝝏𝒕 𝝍 ≈ 𝜓𝑗 𝑛+1 − 𝜓𝑗 𝑛 Δ𝑡 • Approximation of 𝛁𝟐𝝍 ≈ 𝜓𝑗+1 𝑛 − 2𝜓𝑗 𝑛 + 𝜓𝑗−1 𝑛 ℎ2 Resulting Scheme: ∇2𝜓 1D grid 𝒉 ⊲ 𝒪 Δ𝑡 accurate ⊲ 𝒪 ℎ2 accurate 68
  69. Numerical Solvers: Forward Euler • Derivatives approximated with a finite

    difference stencil → A stencil: template object can be readily extended specializations interchangeable with generalizations The whole point: template<size_t DD, size_t OA> struct SelfSelectingStencil 69 len_type dims[]{ 10, 10 }; auto stencil = SelfSelectingStencil<2, 2>::Points<9, 13, 6>{dims, 1.0};
  70. Numerical Solvers: Forward Euler • Derivatives approximated with a finite

    difference stencil → A stencil combines an interface and general implementation 70 template<typename Sp> struct Stencil { template <typename T> auto laplacian(T * const v) const { return (*static_cast<const Sp*>(this)).laplacian(v); } template <typename T> auto bilaplacian(T * const v) const ... template <typename T> auto gradlaplacian(T * const v) const ... template<typename T> auto gradient(T * const v) const }; template<size_t DD, size_t OA = DEFAULT_STENCIL_ACCURACY> struct GeneralizedStencil { len_type dims[DD]; double divh; GeneralizedStencil(const len_type* dims, double h) : divh{ 1. / h } { std::copy(dims, dims + DD, this->dims); } template<size_t OD, typename T> auto apply(T* const v) const { return *v; } }; +
  71. Numerical Solvers: Spectral Solver • Given phase-field problem: • Solution

    is determined in Fourier space → Forms a natural basis set for computational domain (FFT & grids) • Apply the Fourier transform and solve as a linear ODE 71
  72. Numerical Solvers: Spectral Solver • Given phase-field problem: Nonlinear terms

    72
  73. Numerical Solvers: Spectral Solver • Given phase-field problem: Apply Fourier

    Transform… linear operators ⇢ linear functions ∇𝑗→ −𝑘𝑗 ∇𝑚𝑖→ −𝑘𝑚𝑖 73
  74. Numerical Solvers: Spectral Solver • Solve linear ODE: 74

  75. Numerical Solvers: Spectral Solver • Solve linear ODE: ෠ 𝜓

    𝑘 𝑡 + Δ𝑡 = 𝑒−𝐿 𝑘𝑗 Δ𝑡 ෠ 𝜓 𝑘 𝑡 75
  76. Numerical Solvers: Spectral Solver • Solve linear ODE: ෠ 𝜓

    𝑘 𝑡 + Δ𝑡 = 𝑒−𝐿 𝑘𝑗 Δ𝑡 ෠ 𝜓 𝑘 𝑡 Implicit relation for the next time index 76
  77. Numerical Solvers: Spectral Solver • Scheme: Semi-implicit approximation First order

    accurate in time; 𝒪(Δ𝑡) In space, “exponential convergence” 77
  78. Numerical Solvers: Spectral Solver • Determine scheme: Semi-implicit spectral Fourier

    scheme First order accurate in time; 𝒪(Δ𝑡) In space, “exponential convergence” 78
  79. Numerical Solvers: Spectral Solver • Solve linear ODE: “Semi-implicit operators”

    79
  80. Numerical Solvers: Spectral Solver • Semi-implicit Fourier spectral scheme: •

    SymPhas automatically constructs the operators for any problem! 80