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

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. 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
  12. Group of Systems Problem Representation System Uniform Grid System Uniform

    Grid System Uniform Grid Equation of Motion Equation of Motion Equation of Motion 16
  13. 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
  14. 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
  15. 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
  16. 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
  17. 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
  18. 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
  19. 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)
  20. 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
  21. 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
  22. 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
  23. 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
  24. Symbolic Algebra • Binary Operators → OpBinaryAdd → OpBinarySub →

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

    → OpLiteral 𝑎 𝑛𝑢𝑚𝑏𝑒𝑟 → OpLVariable 𝜓 → OpNLVariable 𝜓 ∗ 𝜓 ∗ ⋯ = 𝜓𝑛 → OpFuncDerivative ∇2(𝑎𝑛 𝑒𝑥𝑝𝑟𝑒𝑠𝑠𝑖𝑜𝑛) • More control over representing derivatives (PFC): → OpOperator (1 + ∇2) ⋅ (𝑎𝑛 𝑒𝑥𝑝𝑟𝑒𝑠𝑠𝑖𝑜𝑛) 30
  26. 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
  27. 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
  28. 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
  29. 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
  30. Symbolic Algebra • Library defines: → Distribution of multiplied terms

    → Collection of like terms → Multiplication by 0 or 1 → Addition/Subtraction by 0 → … 35
  31. 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
  32. 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
  33. 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
  34. 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
  35. 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
  36. 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
  37. 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.
  38. 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
  39. Model A non-conserved dynamics; chemical potential driven phase transition MODEL(MA,

    (SCALAR), MODEL_DEF( dpsi = lap(psi) + (c1 - c2 * psi * psi) * psi)) 51
  40. Model B MODEL(MB, (SCALAR), MODEL_DEF( dm = -bilap(m) - lap(c1

    - c2 * m * m) * m)) conserved dynamics; chemical potential driven phase transition 55
  41. • Expression trees fulfill these needs … a construct that

    is both evaluable & modifiable Symbolic Algebra 1 2 + + ÷ 6 𝜓 8 × + 4 61
  42. 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
  43. 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
  44. 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
  45. Numerical Solvers: Forward Euler & periodic BC • Approximation of

    𝝏 𝝏𝒕 𝝍 ≈ 𝜓𝑗 𝑛+1 − 𝜓𝑗 𝑛 Δ𝑡 • Approximation of 𝛁𝟐𝝍 ≈ 𝜓𝑗+1 𝑛 − 2𝜓𝑗 𝑛 + 𝜓𝑗−1 𝑛 ℎ2 Resulting Scheme: ∇2𝜓 1D grid 𝒉 ⊲ 𝒪 Δ𝑡 accurate ⊲ 𝒪 ℎ2 accurate 68
  46. 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};
  47. 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; } }; +
  48. 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
  49. Numerical Solvers: Spectral Solver • Given phase-field problem: Apply Fourier

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

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

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

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

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

    SymPhas automatically constructs the operators for any problem! 80