energy density is given by f = f 0 [c( r)]+ κ 2 ∇c r ( )2 . Here, f0 is a bulk free energy density f 0 c r ( ) ! " # $= − A 2 c −c m ( )2 + B 4 c −c m ( )4 + c α 4 c −c α ( )4 + c β 4 c −c β ( )4 , where c m = 1 2 c α +c β ( ), and cα and cβ are the concentrations at which the bulk free energy density has minima (corresponding to the solubilities in the matrix phase and the second phase, respectively). The time evolution of the concentration field c is given by the Cahn-Hilliard equation: ∂c ∂t = ∇⋅ D(c)∇ df 0 dc −κ∇2c % & ' ( ) * + , - . / 0, where D is a diffusivity. Take Cα=0.05, Cβ=0.95, A=2.0, B = A / c α −c m ( )2 , D = D α = D β = 2 / c β −c α ( ) and κ = 2.0. As initial condition, set c=0.45 throughout the domain. Add a noise term with approximately zero mean value by adding δc r ( )=ε cos q ⋅ r ( ), where q = 2, 3 ( ) and ε=0.01. a. The domain is a 2D square of size 200 x 200 units with periodic boundary conditions, x=(0,200), y=(0,200). b. Now the domain is a square, no periodic boundary conditions but zero- flux boundary conditions c. Now the domain is a T-shaped region with zero-flux boundary conditions, a=b=100, c=d=20. Coordinate origin is at the lower left corner near the letter “d” in the figure. d. The domain is the surface of a sphere of radius 100. Now take the noise term to be δc θ,φ ( )=ε cos 233θ ( )sin 239φ ( ), where (θ, ) are the regular polar and azimuthal angles in a spherical coordinate system. Your task: 1. Calculate the time evolution of the concentration – store concentration at time steps to make a movie 2. Plot free energy as a function of time steps until you judge that convergence or a local equilibrium has been reached. 3. Present wall clock time for the calculations, and wall clock time per core used in the calculation. 4. For part a. above, demonstrate that the solution is robust with respect to meshing by refining the mesh (e.g., reduce the mesh size by about a factor of 2 in linear dimensions – use whatever convenient way you have to refine the mesh without exploding the computational time). 2. Ostwald ripening – coupled Cahn-Hilliard and Allen-Cahn equations Expanded problem in that now phase field described by variables ηi is coupled to the concentration field c. The Ginzburg-Landau free energy density is now taken to be Here, f0 is a bulk free energy density, Here, f1 (C) is the free energy density due to the concentration field C with local minima at Cα and Cβ corresponding to the solubilities in the matrix phase and the second phase, respectively; f2 (C,ηi ) is an interaction term between the concentration field and the phase fields, and f3 (ηi ,ηj ) is the free energy density of the phase fields. Simple models for these free energy densities are f = f 0 C r ( ),η1 ,...,ηp ! " # $+ κC 2 ∇C r ( ) ! " # $ 2 + κi 2 i=1 p ∑ ∇ηi r ( ) ! " # $ 2 . f 0 C r ( ),η1 ,...,ηp ! " # $= f 1 (C)+ f 2 i=1 p ∑ C,ηi ( )+ f 3 j≠i p ∑ i=1 p ∑ ηi ,ηj ( ). a b c d where and f 3 ηi ,ηj ( )= εij 2 ηi 2ηj 2,i ≠ j. The time evolution of the system is now given by coupled Cahn-Hilliard and Allen- Cahn (time-dependent GL) equations for the conserved concentration field and the non-conserved phase fields: ∂C ∂t = ∇⋅ D∇ δF δC( r,t) $ % & ' ( ) * + , - . / = ∇⋅ D∇ df 1 dC + df 2 dC −κc ∇2C( r,t) 1 2 3 4 5 6 = D −A+3B(C −C m )2 +3D α (C −C α )2 +3D β (C −C β )2 $ % ' (∇2C −Dγ ηi 2∇2C + 4∇C ⋅∇ηi + 2(C −C α )∇2ηi $ % ' ( i=1 p ∑ − Dκc ∇4C ∂ηi ∂t = −L i δF δηi = df 2 dηi + df 3 dηi −κi ∇2ηi ( r,t) 1 2 3 4 5 6 = L i γ C −C α ( )2 ηi − L i βηi 3 − L i ηi εij j≠i p ∑ ηj 2 + L i κi ∇2ηi . Take Cα=0.05, Cβ=0.95, A=2.0, B = A / C α −C m ( )2 ,γ = 2 / C β −C α ( )2 ,β =1.0, D α = D β =γ /δ2,εij = 3.0, and κi =κj =κC = 2.0. The diffusion constant D is constant and isotropic and the same (unity) for both phases; the mobility-related constants Li are the same (unity) for all phase fields. Here, p=10, and initial concentration c=0.5. Initial values for the phase fields are 0. Add a noise term with approximately zero mean value to the concentration field by adding δc r ( )=ε cos q ⋅ r ( ), where q = 2, 3 ( ) and ε=0.01, and add noise terms to the phase fields by setting ηi = 0.01εi cos2 q i ⋅ r ( ), where εi for each i=1,…10 are the first 10 numbers in the pseudo-random sequence given below. The vectors qi are given by q i = ( 23+i, 149 +i). a. The domain is a 2D square of size 200 x 200 units with periodic boundary conditions. b. Now the domain is a square, no periodic boundary conditions but zero- flux boundary conditions c. Now the domain is a T-shaped region with zero-flux boundary conditions, a=b=100, c=d=20. d. The domain is the surface of a sphere of radius 100. Now take the noise term to be δc θ,φ ( )=ε cos 233θ ( )sin 239φ ( ), and f 1 C ( )= − A 2 (C −C m )2 + B 4 C −C m ( )4 + D α 4 C −C α ( )4 + D β 4 C −C β ( )4 , C m = 1 2 C α +C β ( ), f 2 C,ηi ( )= − γ 2 C −C α ( )2 ηi 2 + β 4 ηi 4, ηi = 0.01εi cos2 θ 23+i ( )sin2 φ 149 +i ( ), where (θ, ) are the regular polar and azimuthal angles in a spherical coordinate system and εi is as before. Your task: 1. Calculate the time evolution of the concentration and phase fields – store concentration and phase fields at time steps to make a movie 2. Plot free energy as function of time steps until you judge convergence or local equilibrium has been reached. 3. Present wall clock time for the calculations, and wall clock time per core used in the calculation. Pseudo-random number sequence 0.979285 0.219812 0.837709 0.695603 0.225115 0.389266 0.585953 0.614471 0.918038 0.518569 Reference Simulation for Microstructured Materials Olle Heinonena • Jonathan Guyerb • James Warrenb • E. Begum Gulsoyc • Peter Voorheesc aArgonne National Laboratory bNational Institute of Standards and Technology cNorthwestern University CHiMaD hosted the first Phase Field Methods Hackathon on October 14, 2015. The motivation for the hackathon came from one of CHiMaD’s missions: to develop community standards for phase field modeling in materials science, and to distribute community codes for phase field modeling. As a first step along this mission, CHiMaD is working on devel- oping and distributing a set of standard problems for phase field modeling. Test the problems beforehand •free energy function was “unusual” •initial condition not good for 24 h run •scale of problem misinterpreted by some Standardize what to report •what are the simulation conditions? •quantitative results for direct comparisons Team Code Custom C/C++ & Matlab Conditions ??? Teams consisting of two students or postdocs were given a set of problems consisting of two kinds of phase field modeling; each set contained sub- problems of increasing difficulty. The problems were (intended to be!) completely defined in terms of initial conditions, geometry, and material pa- rameters. The teams had 24 hours to solve the problems using whatever numerical codes they had at their disposition (with full access to the internet) - there were no codes provided at the hackathon. All attendees were expected to bring their own lap- tops and connect to the servers they regularly use for running their codes. The goal of the hackathon was to see how different codes and different approaches could handle the problems with respect to accuracy and speed. The aim of the hackathon was not to produce winners or losers, but to advance our understanding of phase field modeling: in that context, all results or at- tempts at solving the problems were valuable. In this first iteration, we learned the most about how to define the problems themselves to even be able to have a successful exercise. Number of Iterations 105 0 2 4 6 8 10 Free Energy (arb) -1000 -800 -600 -400 -200 0 200 400 600 Team Code http://www.prisms-center.org Conditions matrix-free finite element static rectilinear grids parallel on 16 cores Code http://fenicsproject.org Conditions finite element static(?) irregular meshes serial -PF Team Code Custom (based on FFTW) Conditions spectral periodic, static, rectilinear grids 1.5% diﬀerence Free energy Wall 3me (min.) Team Code Custom Conditions spectral & finite difference periodic, static, rectilinear grids serial? Problem 1 box! 104 105 106 −0.015 −0.01 −0.005 0 0.005 0.01 t F NoFlux - dx = 1, dt = 0.001 NoFlux - dx = 0.5, dt = 0.001 NoFlux - dx = 1, dt = 0.0005 Periodic - dx = 1, dt = 0.0001 Periodic - dx = 0.5, dt = 0.0001 Periodic - dx = 1, dt = 0.00005 NoFlux! 104 105 106 0.4475 0.448 0.4485 0.449 0.4495 0.45 0.4505 t avg.c NoFlux - dx = 1, dt = 0.001 NoFlux - dx = 0.5, dt = 0.001 NoFlux - dx = 1, dt = 0.0005 Periodic - dx = 1, dt = 0.0001 Periodic - dx = 0.5, dt = 0.0001 Periodic - dx = 1, dt = 0.00005 Periodic! 0.5%%dri)% Question 2! T boundary" C ﬁeld t=150! Random init. cond." C ﬁeld t=150! 1 2 1 2 Finite diﬀerences code! • Uniform square mesh! • Euler time stepping! • 9 point stencil! • Boundary conditions straightforward! • s/iteration! 0.02 Team Code http://www.ctcms.nist.gov/fipy Conditions implicit finite volume static irregular meshes serial FiPy (*This case not done during the 24 h hackathon) Team Code http://mooseframework.org Conditions implicit finite element adaptive irregular meshes parallel This exercise was inspired, in part, by the µMAG standard problems for micromagnetic modeling developed by that community. From their inception in 1995, these standard problems enabled dra- matic improvements in accuracy and consistency of different simulation codes. http://www.ctcms.nist.gov/~rdm/mumag.org.html 1A: Evolution, 1000 k ¢ t Trevor Keller (NIST MSED) MMSP Hackathon Results CHiMaD PFM Workshop 12 / 61 2C: Evolution, 50 k ¢ t concentration magnitude of order Trevor Keller (NIST MSED) MMSP Hackathon Results CHiMaD PFM Workshop 58 / 61 Team Code https://github.com/mesoscale/mmsp Conditions explicit finite difference adaptive rectilinear grids parallel on 64 cores MMSP 1C: Energy Spinodal decomposition on T-square domain with Neumann boundary conditions Trevor Keller (NIST MSED) MMSP Hackathon Results CHiMaD PFM Workshop 26 / 61 We are looking forward to our next round in May 2016.