Custom Iterators
kawin currently comes with two built-in iterators when solving a model: Explicit Euler and 4th order Runga Kutta. However, custom iterators can be made and used in the solve function.
We’ll use the Al-Zr system from the Binary Precipitation example as a use-case.
1 import numpy as np
2 import matplotlib.pyplot as plt
3 from kawin.thermo import BinaryThermodynamics
4 from kawin.precipitation import PrecipitateModel
5 from kawin.precipitation import MatrixParameters, PrecipitateParameters
6
7 #Set up thermodynamics
8 therm = BinaryThermodynamics('AlScZr.tdb', ['AL', 'ZR'], ['FCC_A1', 'AL3ZR'])
9 therm.setGuessComposition(0.24)
10 diff = lambda T: 0.0768 * np.exp(-242000 / (8.314 * T))
11 therm.setDiffusivity(diff, 'FCC_A1')
12
13 a = 0.405e-9
14 matrix = MatrixParameters(['ZR'])
15 matrix.initComposition = 4e-3
16 matrix.volume.setVolume(a**3, 'VA', 4)
17 matrix.nucleationSites.setNucleationDensity(grainSize=1, dislocationDensity=1e15)
18
19 precipitate = PrecipitateParameters('AL3ZR')
20 precipitate.gamma = 0.1
21 precipitate.volume.setVolume(a**3, 'VA', 4)
22 precipitate.nucleation.setNucleationType('dislocations')
23
24 temperature = 450+273.15
Let’s solve the model using the Explicit Euler method for a comparison.
25 from kawin.solver import explicitEulerIterator
26
27 #Set up model with parameters
28 modelEuler = PrecipitateModel(matrix, precipitate, therm, temperature)
29 modelEuler.solve(500*3600, iterator=explicitEulerIterator, verbose=True, vIt=5000)
N Time (s) Sim Time (s) Temperature (K) Matrix Comp
0 0.0e+00 0.0 723 0.4000
Phase Prec Density (#/m3) Volume Frac Avg Radius (m) Driving Force (J/mol)
AL3ZR 0.000e+00 0.0000 0.0000e+00 5.7737e+03
N Time (s) Sim Time (s) Temperature (K) Matrix Comp
3725 1.8e+06 8.2 723 0.0126
Phase Prec Density (#/m3) Volume Frac Avg Radius (m) Driving Force (J/mol)
AL3ZR 1.397e+22 1.5504 6.0868e-09 3.2941e+02
The custom iteration function will need to take in the following parameters:
- f : function that takes in ($t$, $X$) and returns $dX/dt$ (and $\Delta t$)
- t : current time (float)
- X : current $X$ (list - can be a nested list)
- updateX : function that takes in ($X$, $dX/dt$ and $\Delta t$) and returns $X_{new}$
Notes:
- The functions f and updateX are implemented by the Solver object, so these do not need to be implemented
- While $X_{new}$ can be found by $X + dX/dt * \Delta t$, the updateX function will apply any corrections needed to $dX/dt$ defined by the GenericModel to avoid improper values of $X_{new}$
We’ll create a custom function that uses the Midpoint iterative scheme. This goes by: $$ X_{n+1} = X_n + \Delta t * f\left(t + \frac{\Delta t}{2}, X_n + \frac{\Delta t}{2} * \frac{dX_n}{dt}\right) $$
30 def midpointIterator(f, t, Xn, updateX):
31 #Get dXdt at Xn along with dt
32 dXdt, dt = f(t, Xn, True)
33
34 #Calculate X + dXdt * dt/2 and get dXdt at midpoint
35 Xmid = updateX(Xn, dXdt, dt/2)
36 dXdt_mid = f(t, Xmid)
37
38 #Calculate X + dXdt_mid * dt
39 return updateX(Xn, dXdt_mid, dt), dt
We can solve the model with this new iteration by replacing solverType with this function.
40 #Set up new model with same parameters
41 modelMidpoint = PrecipitateModel(matrix, precipitate, therm, temperature)
42 modelMidpoint.solve(500*3600, iterator=midpointIterator, verbose=True, vIt=5000)
N Time (s) Sim Time (s) Temperature (K) Matrix Comp
0 0.0e+00 0.0 723 0.4000
Phase Prec Density (#/m3) Volume Frac Avg Radius (m) Driving Force (J/mol)
AL3ZR 0.000e+00 0.0000 0.0000e+00 5.7737e+03
N Time (s) Sim Time (s) Temperature (K) Matrix Comp
3647 1.8e+06 13.3 723 0.0126
Phase Prec Density (#/m3) Volume Frac Avg Radius (m) Driving Force (J/mol)
AL3ZR 1.386e+22 1.5505 6.0984e-09 3.2681e+02
Then we could compare the results from the Explicit Euler with the Midpoint iterators.
43 from kawin.precipitation.Plot import plotPrecipitateDensity, plotAverageRadius
44
45 fig, ax = plt.subplots(1, 2, figsize=(8,4))
46
47 plotPrecipitateDensity(modelEuler, ax=ax[0], label='Explicit Euler', alpha=0.75)
48 plotPrecipitateDensity(modelMidpoint, ax=ax[0], label='Midpoint', color='k', linestyle=(0,(5,5)), linewidth=2)
49
50 plotAverageRadius(modelEuler, ax=ax[1], label='Explicit Euler', alpha=0.75)
51 plotAverageRadius(modelMidpoint, ax=ax[1], label='Midpoint', color='k', linestyle=(0,(5,5)), linewidth=2)
52
53 ax[1].legend()
