cahnhilliard.pyΒΆ
In this script we solve the Cahn-Hiilliard equation, which models the unmixing of two phases under the effect of surface tension.
6from nutils import mesh, function, solver, sample, export, cli, testing
7import numpy, treelog, itertools, enum, typing
8
9class stab(enum.Enum):
10 none = '0' # for educational purposes only
11 linear = '.5 dc^2 (6 - 6 c^2 + 8 c dc - 3 dc^2) / epsilon^2'
12 optimal = '.5 dc^2 (1 - dc^2 / 12) / epsilon^2'
13
14def main(nelems:int, etype:str, btype:str, degree:int, epsilon:typing.Optional[float],
15 contactangle:float, timestep:float, mtol:float, seed:int, circle:bool, stab:stab):
16 '''
17 Cahn-Hilliard equation on a unit square/circle.
18
19 .. arguments::
20
21 nelems [20]
22 Number of elements along domain edge.
23 etype [square]
24 Type of elements (square/triangle/mixed).
25 btype [std]
26 Type of basis function (std/spline), with availability depending on the
27 configured element type.
28 degree [2]
29 Polynomial degree.
30 epsilon []
31 Interface thickness; defaults to an automatic value based on the
32 configured mesh density if left unspecified.
33 contactangle [90]
34 Wall contact angle in degrees.
35 timestep [.01]
36 Time step.
37 mtol [.01]
38 Threshold value for chemical potential peak to peak difference, used as
39 a stop criterion.
40 seed [0]
41 Random seed for the initial condition.
42 circle [no]
43 Select circular domain as opposed to a unit square.
44 stab [linear]
45 Stabilization method (linear/optimal/none).
46 '''
47
48 mineps = 1./nelems
49 if epsilon is None:
50 treelog.info('setting epsilon={}'.format(mineps))
51 epsilon = mineps
52 elif epsilon < mineps:
53 treelog.warning('epsilon under crititical threshold: {} < {}'.format(epsilon, mineps))
54
55 domain, geom = mesh.unitsquare(nelems, etype)
56 bezier = domain.sample('bezier', 5) # sample for plotting
57
58 ns = function.Namespace()
59 if not circle:
60 ns.x = geom
61 else:
62 angle = (geom-.5) * (numpy.pi/2)
63 ns.x = function.sin(angle) * function.cos(angle)[[1,0]] / numpy.sqrt(2)
64 ns.epsilon = epsilon
65 ns.ewall = .5 * numpy.cos(contactangle * numpy.pi / 180)
66 ns.cbasis, ns.mbasis = function.chain([domain.basis('std', degree=degree)] * 2)
67 ns.c = 'cbasis_n ?lhs_n'
68 ns.dc = 'cbasis_n (?lhs_n - ?lhs0_n)'
69 ns.m = 'mbasis_n ?lhs_n'
70 ns.F = '.5 (c^2 - 1)^2 / epsilon^2'
71 ns.dF = stab.value
72 ns.dt = timestep
73
74 nrg_mix = domain.integral('F d:x' @ ns, degree=7)
75 nrg_iface = domain.integral('.5 c_,k c_,k d:x' @ ns, degree=7)
76 nrg_wall = domain.boundary.integral('(abs(ewall) + c ewall) d:x' @ ns, degree=7)
77 nrg = nrg_mix + nrg_iface + nrg_wall + domain.integral('(dF - m dc - .5 dt epsilon^2 m_,k m_,k) d:x' @ ns, degree=7)
78
79 numpy.random.seed(seed)
80 lhs = numpy.random.normal(0, .5, ns.cbasis.shape) # initial condition
81
82 with treelog.iter.plain('timestep', itertools.count()) as steps:
83 for istep in steps:
84
85 E = sample.eval_integrals(nrg_mix, nrg_iface, nrg_wall, lhs=lhs)
86 treelog.user('energy: {0:.3f} ({1[0]:.0f}% mixture, {1[1]:.0f}% interface, {1[2]:.0f}% wall)'.format(sum(E), 100*numpy.array(E)/sum(E)))
87
88 x, c, m = bezier.eval(['x_i', 'c', 'm'] @ ns, lhs=lhs)
89 export.triplot('phase.png', x, c, tri=bezier.tri, clim=(-1,1))
90 export.triplot('chempot.png', x, m, tri=bezier.tri)
91
92 if numpy.ptp(m) < mtol:
93 break
94
95 lhs = solver.optimize('lhs', nrg, arguments=dict(lhs0=lhs), lhs0=lhs, tol=1e-10)
96
97 return lhs
If the script is executed (as opposed to imported), nutils.cli.run()
calls the main function with arguments provided from the command line.
102if __name__ == '__main__':
103 cli.run(main)
Once a simulation is developed and tested, it is good practice to save a few
strategic return values for regression testing. The nutils.testing
module, which builds on the standard unittest
framework, facilitates
this by providing nutils.testing.TestCase.assertAlmostEqual64()
for the
embedding of desired results as compressed base64 data.
111class test(testing.TestCase):
112
113 @testing.requires('matplotlib')
114 def test_initial(self):
115 lhs = main(nelems=3, etype='square', btype='std', degree=2, epsilon=None, contactangle=90, timestep=1, mtol=float('inf'), seed=0, circle=False, stab=stab.linear)
116 self.assertAlmostEqual64(lhs, '''
117 eNoBxAA7/xM3LjTtNYs3MDcUyt41uc14zjo0LzKzNm812jFhNNMzwDYgzbMzV8o0yCM1rzWeypE3Tcnx
118 L07NzTa4NlMyETREyrPIGMxYMl82VDbjy1/M8clZyf3IRjday6XLmMl6NRnJDs1Ayh00WMu1yQHRUDSs
119 MKIz7MoEzM/KCMxwyvjIlzLQyxTJdjQ5yjEwWjX3MTk2n8kwNMbKTsoay1DMWDC8ycM1eTQyyb42NzdK
120 NmLN5skSNs/LXDbnMuw19DNKNREtGTfui1ut''')
121
122 @testing.requires('matplotlib')
123 def test_square(self):
124 lhs = main(nelems=3, etype='square', btype='std', degree=2, epsilon=None, contactangle=90, timestep=1, mtol=.1, seed=0, circle=False, stab=stab.linear)
125 self.assertAlmostEqual64(lhs, '''
126 eNqbZTbHzMHsiGmpCd9V1gszzWaZ2ZjtMQ01eXV+xbk0szSgzAaTDxdNTkue1jbTMpM15TJqP/335PeT
127 100vmyqYaJ3tPNV1svNknmmKqYJR+On3J01Pmp9MMY0y/WIYCOSZn7Q82XCi8UTXiSkn5pxYBISovJYT
128 rSd6T0wD8xae6ATCCSemn5gLlusFwiknZp9YcGIpEE4Ewhkn5p1YfGIFEKLyAN6wcSE=''')
129
130 @testing.requires('matplotlib')
131 def test_contactangle(self):
132 lhs = main(nelems=3, etype='square', btype='std', degree=2, epsilon=None, contactangle=45, timestep=1, mtol=.1, seed=0, circle=False, stab=stab.linear)
133 self.assertAlmostEqual64(lhs, '''
134 eNqzNsszkzZbbfrdOOus6Jlss5lmPmbPTQtNtp6be8bZrNTss6mW6SMDv9OnTokDZRpMbxl7nNE89fTk
135 ItNHpl0mT8+fOzX3ZP7J3yb+ph1G206zn7I+KXWyyOSeibK+1ulzJyVP/joRZhJp0m6yyeSyyXsgDAfy
136 2kw2mlw0eWvyxiTLJNtkgslmk3Mmz4CwzqTeZLbJNpOzJo+AcIrJVJO1JkdMbpi8BsLlJitM9gHNeGLy
137 2eQLkLfSZL/JFZOnJl+BEAAJrlyi''')
138
139 @testing.requires('matplotlib')
140 def test_mixedcircle(self):
141 lhs = main(nelems=3, etype='mixed', btype='std', degree=2, epsilon=None, contactangle=90, timestep=1, mtol=.1, seed=0, circle=True, stab=stab.linear)
142 self.assertAlmostEqual64(lhs, '''
143 eNrTM31uImDqY1puGmwia1prssNY37TERNM01eSOkYuJlck6Q1ED9TP9px+fOmq82FjtfKFJiM6CK70m
144 BsZixmUXgk9XnMo7VX6661zL+cZz58+ln0s6e/PM7DOvjDTOvTz97tS8c6xn9pzYemLHiQMn9p9YDyS3
145 nth4YteJbUCRHUByO5DcfGLDieUnlpyYA2RtP7HpxJ4T64Aih8Bwz4k1QPF5QJ3rgap3ntgCVAHRe+bE
146 biBr5YmDQBMBKJ13Eg==''')