EX04 cahn hilliard.edit

import%20marimo%0A%0A__generated_with%20%3D%20%220.23.5%22%0Aapp%20%3D%20marimo.App(width%3D%22medium%22)%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Cahn-Hilliard%20simulation%20with%20regular-solution%20free%20energy%0A%0A%20%20%20%20We%20will%20simulate%20the%20Cahn-Hilliard%20equation%20using%20a%20more%20physical%20potential%20(adapted%20from%20https%3A%2F%2Fweb.tuat.ac.jp%2F~yamanaka%2Fpcoms2019%2FCahn-Hilliard-2d.html).%20The%20governing%20CH%20equation%20is%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cfrac%7B%5Cpartial%20x_B%7D%7B%5Cpartial%20t%7D%20%3D%20%5Cnabla%20%5Ccdot%20%5Cleft(M%20%5Cnabla%20%5Cmu_B%20%5Cright)%0A%20%20%20%20%24%24%0A%0A%20%20%20%20and%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cmu_B%20%3D%20%5Cfrac%7B%5Cpartial%20g%5E%7B%5Ctext%7Bhomo%7D%7D%7D%7B%5Cpartial%20x_B%7D%20-%202%5Ckappa%20%5Cnabla%5E2%20x_B%0A%20%20%20%20%24%24%0A%0A%20%20%20%20-%20Gibbs%20free%20energy%0A%20%20%20%20%24%24%0A%20%20%20%20g%5E%7B%5Ctext%7Bhomo%7D%7D(x_B)%20%3D%20%5Comega%20x_B(1-x_B)%20%2B%20k_B%20T%20%5Cleft%5B(1-x_B)%5Cln(1-x_B)%2Bx_B%5Cln%20x_B%5Cright%5D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20-%20Composition-dependent%20mobility%0A%0A%20%20%20%20%24%24%0A%20%20%20%20M%20%3D%20%5Cfrac%7B%5BD_A%20x_B%20%2B%20D_B%20(1%20-%20x_B)%5D%7D%7B%5Cfrac%7Bk_B%20T%7D%7Bx_B(1-x_B)%7D%20-%202%5Comega%7D%0A%20%20%20%20%24%24%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20setup%20%3D%20mo.md(%22%22%22%0A%20%20%20%20%23%23%20Setup%20the%20potential%20profile%20and%20initial%20concentration%0A%0A%20%20%20%20-%20%24%5C%5Comega%20%2F%20k_B%20T%24%3A%20%7Bomega%7D%0A%20%20%20%20-%20%24T%24%20(K)%3A%20%7BT%7D%0A%20%20%20%20-%20%24x_%7B%7BB0%7D%7D%24%3A%20%7BxB0%7D%0A%20%20%20%20%22%22%22).batch(%0A%20%20%20%20%20%20%20%20omega%20%3D%20mo.ui.slider(%0A%20%20%20%20%20%20%20%20%20%20%20%20start%3D0.5%2C%20stop%3D5.0%2C%20step%3D0.05%2C%20value%3D2.25%2C%20show_value%3DTrue%0A%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%20%20%20%20T%20%3D%20mo.ui.slider(start%3D300%2C%20stop%3D1000%2C%20step%3D10%2C%20value%3D600%2C%20show_value%3DTrue)%2C%0A%20%20%20%20%20%20%20%20xB0%20%3D%20mo.ui.slider(start%3D0.10%2C%20stop%3D0.90%2C%20step%3D0.01%2C%20value%3D0.50%2C%20show_value%3DTrue)%0A%20%20%20%20)%0A%20%20%20%20return%20(setup%2C)%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%208.314%20*%20673%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20%23%23%20definition%20of%20the%20free%20energy%20and%20the%20chemical%20potential%0A%0A%20%20%20%20kB%20%3D%201.380649e-23%0A%20%20%20%20%23%20T%20%3D%20673.0%0A%0A%20%20%20%20nx%20%3D%20128%0A%20%20%20%20ny%20%3D%20128%0A%20%20%20%20dx%20%3D%202.0e-9%0A%20%20%20%20dy%20%3D%202.0e-9%0A%0A%20%20%20%20Da%20%3D%201.0e-22%0A%20%20%20%20Db%20%3D%205.0e-23%0A%0A%20%20%20%20%23%20dt%20%3D%202.0e-2%0A%20%20%20%20dt%20%3D%20(dx*dx%2FDa)*0.1%0A%20%20%20%20print(dt)%0A%20%20%20%20eps%20%3D%201.0e-8%0A%20%20%20%20seed%20%3D%2042%0A%0A%0A%20%20%20%20def%20free_energy(x%2C%20omega_kbT_ratio%3D2.0%2C%20T%3D673.0)%3A%0A%20%20%20%20%20%20%20%20kbT%20%3D%20kB%20*%20T%0A%20%20%20%20%20%20%20%20omega%20%3D%20omega_kbT_ratio%20*%20kbT%0A%20%20%20%20%20%20%20%20x%20%3D%20np.clip(x%2C%20eps%2C%201.0%20-%20eps)%0A%20%20%20%20%20%20%20%20return%20omega%20*%20x%20*%20(1.0%20-%20x)%20%2B%20kbT%20*%20(%0A%20%20%20%20%20%20%20%20%20%20%20%20x%20*%20np.log(x)%20%2B%20(1.0%20-%20x)%20*%20np.log(1.0%20-%20x)%0A%20%20%20%20%20%20%20%20)%0A%0A%0A%20%20%20%20def%20chemical_potential(x%2C%20omega_kbT_ratio%3D2.0%2C%20T%3D673.0)%3A%0A%20%20%20%20%20%20%20%20kbT%20%3D%20kB%20*%20T%0A%20%20%20%20%20%20%20%20omega%20%3D%20omega_kbT_ratio%20*%20kbT%0A%20%20%20%20%20%20%20%20x%20%3D%20np.clip(x%2C%20eps%2C%201.0%20-%20eps)%0A%20%20%20%20%20%20%20%20return%20omega%20*%20(1.0%20-%202.0%20*%20x)%20%2B%20kbT%20*%20np.log(x%20%2F%20(1.0%20-%20x))%0A%0A%0A%20%20%20%20def%20dmu_dc(x%2C%20omega_kbT_ratio%3D2.0%2C%20T%3D673.0)%3A%0A%20%20%20%20%20%20%20%20kbT%20%3D%20kB%20*%20T%0A%20%20%20%20%20%20%20%20omega%20%3D%20omega_kbT_ratio%20*%20kbT%0A%20%20%20%20%20%20%20%20x%20%3D%20np.clip(x%2C%20eps%2C%201.0%20-%20eps)%0A%20%20%20%20%20%20%20%20return%20-2.0%20*%20omega%20%2B%20kbT%20%2F%20(x%20*%20(1.0%20-%20x))%0A%0A%0A%20%20%20%20def%20darken_diffusivity(x%2C%20Da%2C%20Db)%3A%0A%20%20%20%20%20%20%20%20return%20Da%20*%20x%20%2B%20Db%20*%20(1.0%20-%20x)%0A%0A%0A%20%20%20%20def%20mobility(x%2C%20Da%2C%20Db%2C%20omega_kbT_ratio%2C%20T)%3A%0A%20%20%20%20%20%20%20%20return%20darken_diffusivity(x%2C%20Da%2C%20Db)%20%2F%20dmu_dc(x%2C%20omega_kbT_ratio%2C%20T)%0A%0A%0A%20%20%20%20def%20dmdc(x%2C%20Da%2C%20Db%2C%20omega_kbT_ratio%2C%20T%2C%20h%3D1.0e-5)%3A%0A%20%20%20%20%20%20%20%20%23%20Finite%20difference%20of%20dmu%2Fdc%0A%20%20%20%20%20%20%20%20xp%20%3D%20np.clip(x%20%2B%20h%2C%201.0e-8%2C%201.0%20-%201.0e-8)%0A%20%20%20%20%20%20%20%20xm%20%3D%20np.clip(x%20-%20h%2C%201.0e-8%2C%201.0%20-%201.0e-8)%0A%20%20%20%20%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20%20%20%20%20mobility(xp%2C%20Da%2C%20Db%2C%20omega_kbT_ratio%2C%20T)%0A%20%20%20%20%20%20%20%20%20%20%20%20-%20mobility(xm%2C%20Da%2C%20Db%2C%20omega_kbT_ratio%2C%20T)%0A%20%20%20%20%20%20%20%20)%20%2F%20(xp%20-%20xm)%0A%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20chemical_potential%2C%0A%20%20%20%20%20%20%20%20dmdc%2C%0A%20%20%20%20%20%20%20%20dmu_dc%2C%0A%20%20%20%20%20%20%20%20dt%2C%0A%20%20%20%20%20%20%20%20dx%2C%0A%20%20%20%20%20%20%20%20dy%2C%0A%20%20%20%20%20%20%20%20free_energy%2C%0A%20%20%20%20%20%20%20%20kB%2C%0A%20%20%20%20%20%20%20%20mobility%2C%0A%20%20%20%20%20%20%20%20nx%2C%0A%20%20%20%20%20%20%20%20ny%2C%0A%20%20%20%20%20%20%20%20seed%2C%0A%20%20%20%20)%0A%0A%0A%40app.cell%0Adef%20_(chemical_potential%2C%20dmu_dc%2C%20free_energy%2C%20kB%2C%20np%2C%20plt)%3A%0A%20%20%20%20def%20plot_energy(omega_kbT_ratio%3D4.5%2C%20T%3D673.0)%3A%0A%20%20%20%20%20%20%20%20fig%2C%20axes%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(8%2C%204))%0A%20%20%20%20%20%20%20%20ax1%2C%20ax2%20%3D%20axes%5B%3A%5D%0A%20%20%20%20%20%20%20%20xB%20%3D%20np.linspace(0.01%2C%200.99%2C%20200)%0A%20%20%20%20%20%20%20%20g%20%3D%20free_energy(xB%2C%20omega_kbT_ratio%2C%20T)%0A%20%20%20%20%20%20%20%20mu%20%3D%20chemical_potential(xB%2C%20omega_kbT_ratio%2C%20T)%0A%20%20%20%20%20%20%20%20dmudc%20%3D%20dmu_dc(xB%2C%20omega_kbT_ratio%2C%20T)%0A%20%20%20%20%20%20%20%20ax1.plot(xB%2C%20g%20%2F%20kB%20%2F%20T)%0A%20%20%20%20%20%20%20%20ax1.set_xlabel(%22%24x_B%24%22)%0A%20%20%20%20%20%20%20%20ax1.axhline(y%3D0%2C%20ls%3D%22--%22%2C%20color%3D%22grey%22)%0A%0A%20%20%20%20%20%20%20%20ax2.plot(xB%2C%20dmudc%20%2F%20kB%20%2F%20T)%0A%20%20%20%20%20%20%20%20ax1.set_ylabel(%22%24g%5E%7B%7B%5C%5Cmathrm%7B%7Bhomo%7D%7D%7D%7D%24%20(%24k_B%20T%24)%22)%0A%20%20%20%20%20%20%20%20ax2.set_xlabel(%22%24x_B%24%22)%0A%0A%20%20%20%20%20%20%20%20ax2.set_ylabel(%22%24%5C%5Cpartial%5E2%20g%5E%7B%7B%5C%5Cmathrm%7B%7Bhomo%7D%7D%7D%7D%20%2F%20%5C%5Cpartial%20x_B%24%20(%24k_B%20T%24)%22)%0A%20%20%20%20%20%20%20%20ax2.axhline(y%3D0%2C%20ls%3D%22--%22%2C%20color%3D%22grey%22)%0A%20%20%20%20%20%20%20%20ax2.set_ylim(-10%2C%2030)%0A%0A%0A%20%20%20%20%20%20%20%20fig.tight_layout()%0A%20%20%20%20%20%20%20%20return%20fig%0A%0A%20%20%20%20plot_energy()%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(dx%2C%20dy%2C%20np)%3A%0A%20%20%20%20def%20laplacian(c)%3A%0A%20%20%20%20%20%20%20%20%23%20Calculate%20the%20laplacian%20field%20for%20c%0A%20%20%20%20%20%20%20%20ce%20%3D%20np.roll(c%2C%20-1%2C%20axis%3D0)%0A%20%20%20%20%20%20%20%20cw%20%3D%20np.roll(c%2C%20%201%2C%20axis%3D0)%0A%20%20%20%20%20%20%20%20cn%20%3D%20np.roll(c%2C%20-1%2C%20axis%3D1)%0A%20%20%20%20%20%20%20%20cs%20%3D%20np.roll(c%2C%20%201%2C%20axis%3D1)%0A%20%20%20%20%20%20%20%20return%20(ce%20-%202.0%20*%20c%20%2B%20cw)%20%2F%20dx**2%20%2B%20(cn%20-%202.0%20*%20c%20%2B%20cs)%20%2F%20dy**2%0A%0A%20%20%20%20def%20gradient_sq_term(c%2C%20mu)%3A%0A%20%20%20%20%20%20%20%20%23%20Calculate%20the%20shifted%20gradient%20form%0A%20%20%20%20%20%20%20%20ce%20%3D%20np.roll(c%2C%20-1%2C%20axis%3D0)%0A%20%20%20%20%20%20%20%20cw%20%3D%20np.roll(c%2C%20%201%2C%20axis%3D0)%0A%20%20%20%20%20%20%20%20cn%20%3D%20np.roll(c%2C%20-1%2C%20axis%3D1)%0A%20%20%20%20%20%20%20%20cs%20%3D%20np.roll(c%2C%20%201%2C%20axis%3D1)%0A%0A%20%20%20%20%20%20%20%20mue%20%3D%20np.roll(mu%2C%20-1%2C%20axis%3D0)%0A%20%20%20%20%20%20%20%20muw%20%3D%20np.roll(mu%2C%20%201%2C%20axis%3D0)%0A%20%20%20%20%20%20%20%20mun%20%3D%20np.roll(mu%2C%20-1%2C%20axis%3D1)%0A%20%20%20%20%20%20%20%20mus%20%3D%20np.roll(mu%2C%20%201%2C%20axis%3D1)%0A%0A%20%20%20%20%20%20%20%20dc2dx2%20%3D%20((ce%20-%20cw)%20*%20(mue%20-%20muw))%20%2F%20(4.0%20*%20dx**2)%0A%20%20%20%20%20%20%20%20dc2dy2%20%3D%20((cn%20-%20cs)%20*%20(mun%20-%20mus))%20%2F%20(4.0%20*%20dy**2)%0A%20%20%20%20%20%20%20%20return%20dc2dx2%20%2B%20dc2dy2%0A%0A%0A%20%20%20%20return%20gradient_sq_term%2C%20laplacian%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%203.0e-14%20*%201e18%20%2F%205595.322%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20chemical_potential%2C%0A%20%20%20%20dmdc%2C%0A%20%20%20%20dt%2C%0A%20%20%20%20gradient_sq_term%2C%0A%20%20%20%20kB%2C%0A%20%20%20%20laplacian%2C%0A%20%20%20%20mobility%2C%0A%20%20%20%20np%2C%0A%20%20%20%20nx%2C%0A%20%20%20%20ny%2C%0A%20%20%20%20seed%2C%0A)%3A%0A%20%20%20%20rng%20%3D%20np.random.default_rng(seed)%0A%0A%0A%20%20%20%20def%20run_simulation(%0A%20%20%20%20%20%20%20%20nsteps%3D3000%2C%0A%20%20%20%20%20%20%20%20xB0%3D0.5%2C%0A%20%20%20%20%20%20%20%20noise%3D0.01%2C%0A%20%20%20%20%20%20%20%20kappa_kbT_ratio%3D2.5%2C%0A%20%20%20%20%20%20%20%20omega_kbT_ratio%3D2.5%2C%0A%20%20%20%20%20%20%20%20T%3D673.0%2C%0A%20%20%20%20%20%20%20%20Da%3D1e-11%2C%0A%20%20%20%20%20%20%20%20Db%3D1e-13%2C%0A%20%20%20%20%20%20%20%20every%3D600%2C%0A%20%20%20%20)%3A%0A%20%20%20%20%20%20%20%20x%20%3D%20xB0%20%2B%20noise%20*%20rng.standard_normal((nx%2C%20ny))%0A%20%20%20%20%20%20%20%20x%20%3D%20np.clip(x%2C%201.0e-8%2C%201.0%20-%201.0e-8)%0A%0A%20%20%20%20%20%20%20%20kbT%20%3D%20kB%20*%20T%0A%20%20%20%20%20%20%20%20omega%20%3D%20omega_kbT_ratio%20*%20kbT%0A%20%20%20%20%20%20%20%20kappa%20%3D%20kappa_kbT_ratio%20*%20kbT%0A%20%20%20%20%20%20%20%20frames%20%3D%20%5Bx.copy()%5D%0A%0A%20%20%20%20%20%20%20%20for%20step%20in%20range(nsteps)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20mu_chem%20%3D%20chemical_potential(x%2C%20omega_kbT_ratio%2C%20T)%0A%20%20%20%20%20%20%20%20%20%20%20%20mu%20%3D%20mu_chem%20-%20kappa%20*%20laplacian(x)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20M%20%3D%20mobility(x%2C%20Da%2C%20Db%2C%20omega_kbT_ratio%2C%20T)%0A%20%20%20%20%20%20%20%20%20%20%20%20dMdc%20%3D%20dmdc(x%2C%20Da%2C%20Db%2C%20omega_kbT_ratio%2C%20T)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20nabla_mu%20%3D%20laplacian(mu)%0A%20%20%20%20%20%20%20%20%20%20%20%20grad_term%20%3D%20gradient_sq_term(x%2C%20mu)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20dxdt%20%3D%20M%20*%20nabla_mu%20%2B%20dMdc%20*%20grad_term%0A%20%20%20%20%20%20%20%20%20%20%20%20x%20%3D%20x%20%2B%20dt%20*%20dxdt%0A%20%20%20%20%20%20%20%20%20%20%20%20x%20%3D%20np.clip(x%2C%201.0e-8%2C%201.0%20-%201.0e-8)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20(step%20%2B%201)%20%25%20every%20%3D%3D%200%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20frames.append(x.copy())%0A%0A%20%20%20%20%20%20%20%20return%20x%2C%20frames%0A%0A%0A%20%20%20%20final_x%2C%20frames%20%3D%20run_simulation()%0A%20%20%20%20return%20final_x%2C%20frames%0A%0A%0A%40app.cell%0Adef%20_(final_x)%3A%0A%20%20%20%20final_x%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(frames%2C%20plt)%3A%0A%20%20%20%20fig%2C%20ax%20%3D%20plt.subplots(figsize%3D(5%2C%204))%0A%20%20%20%20im%20%3D%20ax.imshow(frames%5B1%5D%2C%20cmap%3D%22bwr%22%2C%20vmin%3D0.0%2C%20vmax%3D1.0)%0A%20%20%20%20ax.set_title(%22final%20composition%20field%22)%0A%20%20%20%20fig.colorbar(im%2C%20ax%3Dax%2C%20label%3Dr%22%24x_B%24%22)%0A%20%20%20%20plt.tight_layout()%0A%20%20%20%20fig%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%0A%20%20%20%20%23%20kappa%20%3D%20mo.ui.slider(%0A%20%20%20%20%23%20%20%20%20%20start%3D0.1e-40%2C%20stop%3D8.0e-40%2C%20step%3D0.1e-40%2C%20value%3D2.0e-40%2C%0A%20%20%20%20%23%20%20%20%20%20label%3Dr%22%24%5Ckappa%24%20%5BJ%20m%24%5E2%24%2Fsite%5D%22%0A%20%20%20%20%23%20)%0A%20%20%20%20%23%20xB0%20%3D%20mo.ui.slider(%0A%20%20%20%20%23%20%20%20%20%20start%3D0.05%2C%20stop%3D0.95%2C%20step%3D0.01%2C%20value%3D0.50%2C%0A%20%20%20%20%23%20%20%20%20%20label%3Dr%22%24x_%7BB0%7D%24%22%0A%20%20%20%20%23%20)%0A%20%20%20%20%23%20noise%20%3D%20mo.ui.slider(%0A%20%20%20%20%23%20%20%20%20%20start%3D0.0%2C%20stop%3D0.05%2C%20step%3D0.001%2C%20value%3D0.01%2C%0A%20%20%20%20%23%20%20%20%20%20label%3D%22initial%20noise%22%0A%20%20%20%20%23%20)%0A%20%20%20%20%23%20nsteps%20%3D%20mo.ui.slider(%0A%20%20%20%20%23%20%20%20%20%20start%3D100%2C%20stop%3D5000%2C%20step%3D100%2C%20value%3D1200%2C%0A%20%20%20%20%23%20%20%20%20%20label%3D%22number%20of%20steps%22%0A%20%20%20%20%23%20)%0A%20%20%20%20%23%20every%20%3D%20mo.ui.slider(%0A%20%20%20%20%23%20%20%20%20%20start%3D5%2C%20stop%3D200%2C%20step%3D5%2C%20value%3D40%2C%0A%20%20%20%20%23%20%20%20%20%20label%3D%22store%20every%20n%20steps%22%0A%20%20%20%20%23%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(setup)%3A%0A%20%20%20%20setup%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_()%3A%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20import%20scipy%0A%20%20%20%20from%20scipy.optimize%20import%20root%0A%0A%20%20%20%20return%20np%2C%20plt%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%0A%20%20%20%20return%20(mo%2C)%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A