Posts by author wd15

Grain Boundary Motion 4

Bill has corrected the boundary conditions so let's restate the problem. We have (the first equation is missing due to some bug or other in trac).

$ z_t + v z_x = -\left[ \frac{ \kappa_x }{ \left(1 + z_x^2 \right)^{1/2} } \right]_x $

$ \kappa = \left[ \frac{z_x}{\left(1 + z_x^2 \right)^{1/2}} \right]_x $

with a jump at $x=0$, but both $z$ and $\kappa$ are continuous. The jump condition for the first equation is,

$ 2j = \Delta \left( \frac{ \kappa_z }{ \left(1 + z_x^2 \right)^{1/2} } \right) $

This works perfectly as the equation can now be written

$ \int z_t dx + \int \left[ v z_x \right] dx = \int \left[ \frac{ \kappa_x }{ \left(1 + z_x^2 \right)^{1/2} } \right]_x dx + \int \left[\delta \left(x\right) 2 j \right] dx $

in the finite volume formulation. The issue is how to provide a boundary condition for the second equation with the correct form for the jump condition. Define $\psi$ to be the rotation of the grain boundary and $\theta$ the displacement angle of the upper surfaces in the frame of reference where the grain boundary is vertical. The following holds,

$\gamma_{gb} = 2 \gamma_s \sin \theta$

$ z_x^- = \tan \left(\theta + \psi\right) $

$ z_x^+ = \tan \left(\theta - \psi\right) $`

Essentially, we have to eliminate the angles for an expression for only the $z_x$ and the surface energies. Easy enough to do, but what we really want is a jump condition that fits with equation 2. Ideally this would be of the form.

$\beta = \Delta \left( \frac{ z_z }{ \left(1 + z_x^2 \right)^{1/2} } \right)$

where $\beta$ is a constant. This isn't possible with the above expression involving the angles and grain boundaries. It is quite obvious that

$\Delta \left( \frac{ z_z }{ \left(1 + z_x^2 \right)^{1/2} } \right)$

is just the jump in $dz / ds$ where $s$ is an arc length. This condition requires the vertical displacement from the horizontal and hence knowledge of the orientation of the grain boundary. In terms of the angles,

$\Delta \left( \frac{ z_z }{ \left(1 + z_x^2 \right)^{1/2} } \right) = \sin \left(\theta - \psi\right) - \sin \left(\theta + \psi\right) $

From a practical point of view, I suppose we can deal with this boundary conditions simply be iterating at each time step. However, there seems to be something deeply unsatisfying about it.

  • Posted: 2012-06-28 12:09 (Updated: 2012-06-28 12:18)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Grain Boundarly Motion 4

Bill, Has corrected the boundary conditions so let's restate the problem. We have:

$ z_t + v z_x = \left[ \frac{ \kappa_x }{ \left(1 + z_x^2 \right)^{1/2} } \right]_x $

$ \kappa = \[left \frac{z_x}{\left(1 + z_x^2 \right)^{1/2}} \right]_x $

with a jump at $x=0$, but both $z$ and $\kappa$ are continuous. The jump condition for the first equation is,

$ 2j = \Delta \left( \frac{ \kappa_z }{ \left(1 + z_x^2 \right)^{1/2} } \right) $

  • Posted: 2012-06-26 11:50
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Grain Boundary Motion 3

The solution for grain boundary motion with translation is given by,

$ z_- = \frac{2 j}{a^3} - \frac{2 j + a^2 M}{3 a^3} \exp\left(a \xi\right) $

$ z_+ = \exp\left(-a \xi / 2\right) \left[\frac{M}{\sqrt{3} a} \sin\left(a \xi \sqrt{3} / 2\right) + \frac{4 j - a^2 M}{3 a^3} \cos\left(a \xi \sqrt{3} / 2\right)\right] $

with

$ z = \begin{cases} z_- & \text{when} \; x<0 \\ z_+ & \text{when} \; x>0 \end{cases} $

where

$a = (V / B)^{1/3} $

and

$\xi = x - V t$

The comparison with the 1D numerical solution seems to be quite good. The numbers are $L=40$, $N=400$, $B=1$, $M=0.1$, $j=-0.1$ and $V=1$.

source:trunk/boettinger/grainBoundary/translation.png@1561

The code for this is source:trunk/boettinger/grainBoundary/gb.py@1561.

  • Posted: 2012-06-19 15:53 (Updated: 2012-06-28 12:05)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Grain Boundary Motion 2

In blog:GrainBoundaryMotion the calculation of the analytical solution was broken. The complete analytical solution including the power series is

$ z = z_1 + z_2 $

where

$ z_1 = \frac{M}{2} \left( B t \right)^{1 / 4} P \left( \frac{|x|}{\left( B t \right)^{1 / 4} } \right) $

and

$ z_2 = j \left( B t \right)^{3 / 4} W \left( \frac{|x|}{\left( B t \right)^{1 / 4} } \right) $

The power series are given by,

$ P \left( u \right) = \sum_n a_n u^n $

with

$ a_0 = -\frac{1}{2^{1/2} \Gamma \left(5 / 4\right)}$

$a_1 = 1$

$a_2 = -\frac{1}{2^{3/2} \Gamma \left(3 / 4\right)}$

$a_3 = 0$

$a_{n+4} = a_n \frac{n-1}{4 \left(n + 1\right) \left(n + 2\right) \left(n + 3\right) \left(n + 4\right)} $

and

$ W \left( u \right) = \sum_n b_n u^n $

with

$b_0 = \frac{1}{2^{1/2} \Gamma \left(7 / 4\right)}$

$b_1 = 0$

$b_2 = -\frac{1}{2^{3/2} \Gamma \left(5 / 4\right)}$

$b_3 = \frac{1}{6}$

$b_{n+4} = b_n \frac{n-3}{4 \left(n + 1\right) \left(n + 2\right) \left(n + 3\right) \left(n + 4\right)} $

The following image is generated using revision 1559 (source:trunk/boettinger/grainBoundary/gbStationary.py@1559) for $L=40$, $N=400$, $B=1$, $M=0.1$ and $j=-0.1$.

source:trunk/boettinger/grainBoundary/stationary.png@1559

  • Posted: 2012-06-19 14:48
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Grain Boundary Motion

Bill has a simple 1D model for grain boundary motion and hillock formation. The governing equation is given by

$$ z_t + v z_x =- B \left[ \left( 1 + z_x^2 \right)^{-1/2} \kappa_x \right]_x $$

where

$$\kappa = \frac{z_{xx}}{\left(1 + z_x^2\right)^{3/2}}$$

with two jump conditions at $x=0$.

$$z_x^+ - z_x^- = M$$

$$\kappa^+_x - \kappa^-_x = 2 j$$

We can include these conditions in the main equations if we write the equations in the weak form:

$$ \int \left[z_t - v z_x =- B \left[\frac{1}{\sqrt{ 1 + z_x^2}} \kappa_x\right]_x + \frac{B 2 j \delta\left(x\right)}{\sqrt{ 1 + z_x^2}} \right] dx$$

and

$$ \int \left[\kappa \left( 1 + z_x^2 \right)^{3 / 2} = z_{xx} - M \delta\left(x\right) \right] dx$$

The fipy script is in source:trunk/boettinger/grainBoundary/gb.py.

Version 1551 (source:trunk/boettinger/grainBoundary/gb.py@1551) compares against an analytical solution given by:

$$ z_a = \frac{M}{2} \left(Bt\right)^{1/4}P\left(\frac{|x|}{\left(Bt\right)^{1/4}}\right) $$

where

$$P\left(u\right) = \sum_{n=0}^{\infty} a_n u^n$$

where

$$ a_i = (-7.803e-1, 10e-1, -2.866e-1, 0.0, 8.13e-3, 0.0, -2.004e-4, 0.0, 3.625e-6, 0.0, -4.969e-8, 0.0, 5.339e-10, 0.0, -4.655e-12,...$$

The solution is good for small M in the vicinity of the groove and $j=0$.

For values of $j=0$, $M=0.01$ and $B=1$ with a box size of 40 and 200 grid points. The analytical and numerical solutions are in good agreement in the location of the groove at t=10.

source:trunk/boettinger/grainBoundary/plot99.png@1552

Another test is with an analytical solution for $j\ne0$. This looks like

$$z = z_a + j \left(Bt\right)^{3 / 4} W\left(\frac{|x|}{\left(Bt\right)^{1/4}}\right) $$

where

$$ W\left(u\right) = \sum_{n=0}^{\infty} b_n u^n$$

and

$$ b_0 = \frac{1}{\sqrt{2} \Gamma\left(7 / 4\right)} $$
$$ b_1 = 0 $$
$$ b_2 = \frac{-1}{2^{3 / 2} \Gamma\left(5 / 4\right)} $$
$$ b_3 = 1 / 6 $$
$$ b_{n + 4} = b_n \frac{n - 3}{4 \left(n + 1\right) \left(n + 2\right) \left(n + 3\right) \left(n + 4\right)}$$

The code for this comparison is source:trunk/boettinger/grainBoundary/gb.py@1553.

The image below is for $j=0.1$ at t=10 using the first 14 terms in the power series.

source:trunk/boettinger/grainBoundary/plot99.png@1553

The analytical solution does a good job in the region of the groove, but doesn't rise above $z=0$ as it does when $j=0$. Obviously, the power series is not much good if the $u$ argument is much above 1. If we use this as a boundary then we get $x=1.77$ at $t=10$, which seems to correspond with the domain where the solutions are the same.

  • Posted: 2012-05-30 16:43 (Updated: 2012-05-30 18:23)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Notes on Light Propogation using the Eikonal

I want to try and figure out whether it is possible to solve the eikonal equation for a complex function of the wave fronts to simulate absorption of light. Here, we will recap assuming the wave front function is real. The wave equation for the electric field is given by,

$ \nabla^2 \vec{E} + \frac{n^2 \omega^2}{c^2} \vec{E} = 0 $

where we will assume an electric field intensity of

$ \vec{E} = \vec{E}_0 e^{i\left(k R - \omega t \right)} $

where

$k = \frac{\omega}{c} = \frac{2 \pi}{\lambda} $

If we assume a short wavelength than the wave equation reduces to

$ \left[\nabla R\right]\cdot \left[\nabla R\right] = n^2 $

The time averaged Poynting vector is given by,

$\vec{S} = \frac{|\vec{E}|^2}{4 c^2 \mu_0} \nabla R$

and the intensity is $|S|$ given by,

$I = \frac{n^2}{4 c^2 \mu_0} |\vec{E}|^2 $

  • Posted: 2012-03-07 09:34
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Electrical problem including metal diffusion

This post includes the effects of the counter electrode and the metal diffusion. The metal diffusion is included via a solution where the metal gradient doesn't change rapidly at the interface. This is not really very physical, but it will be useful when testing numerical solutions.

At the working electrode $W$:

$ i_T^W = -\kappa \vec{n}^W \cdot \nabla V_{IR}$

$ i_T^W = -c_{DL} \frac{\partial V_{DL}^W}{\partial t} + i_{F}^W $

$ i_F^W = \frac{c_{cu}^0}{ c_{cu}^{\infty}} i_0^W E^W$

$ D_{cu} \vec{n}^W \cdot \nabla c_{cu} = -\frac{i_F^W}{2 F} $

At the counter electrode $C$:

$ i_T^C = -\kappa \vec{n}^C \cdot \nabla V_{IR}$

$ i_T^C = -c_{DL} \frac{\partial V_{DL}^C}{\partial t} + i_{F}^C $

$ i_F^C = i_0^C E^C$

where

$ E = \exp{\left(-\frac{\alpha F V_{DL}}{R T} \right)} - \exp{\left(\frac{\left(2 -\alpha\right) F V_{DL}}{R T} \right)} $

The normals point out of the electrolyte.

In the bulk:

$ \nabla^2 V_{IR} = 0 $

$ \frac{\partial c_{cu}}{\partial t} = D_{cu} \nabla^2 c_{cu} $

and

$V_{CELL} = V_{DL}^W + V_{IR}^W - V_{DL}^C - V_{IR}^C$

$V_{APPLIED} = V_{DL}^W + V_{IR}^W$

$\delta_W V_{IR}^C + \delta_C V_{IR}^W = 0 $

$\delta_C + \delta_W = \delta$

where $\delta_C$ and $\delta_W$ are the distances from the counter and working electrodes to the reference electode, respectively. $\delta$ is the total distance between the electrodes.

The metal diffusion has an analytical solution

$ c_{cu} \left( x, t \right) = \frac{c_x^0}{k} \cos{\left(k x\right)} \exp{\left(-k D_{cu} t\right)} + c_x^0 \left(x - \delta \right) + c_{cu}^{\infty} $

where

$ c_x^0 = \frac{i_F^W}{2 F D_{cu}} $

$ k = \frac{\pi}{2 \delta} $

This solution is valid when $|c_{xt}^0| \ll 1$. This leads to two dimensionless quantities that must be small,

$\Gamma_0 = \frac{\delta i_0^W E^W}{c_{cu}^{\infty} 2 F D_{cu}} \ll 1$

and

$ \Gamma_1 = \Gamma_0 \frac{\alpha F V_{APPLIED}}{R T} \ll 1$

We can now derive an expression for $c_{cu}^0$.

$ c_{cu}^0 = \frac{i_F^W}{2 F D} \left[ \frac{ \exp{\left(-k D_{cu} t\right)} }{k} - \delta \right] + c_{cu}^{\infty} $

Using the above solution we can derive a new expression for $i_F^W$.

$ i_F^W = i_0^W \frac{E^W}{ 1 - \frac{1}{2 F D c_{cu}^{\infty}} \left( \frac{ \exp{\left(-k D_{cu} t\right)} }{k} - \delta \right) i_0^W E^W} $

If we do this the equations can be reduced to two coupled ODEs for $V_{DL}^C$ and $V_{DL}^W$:

$ -\frac{\kappa}{\delta_W} \left(V_{APPLIED} - V_{DL}^W\right) = -c_{DL} \frac{\partial V_{DL}^W}{\partial t} + i_0^W \frac{E^W}{ 1 - \frac{1}{2 F D c_{cu}^{\infty}} \left( \frac{ \exp{\left(-k D_{cu} t\right)} }{k} - \delta \right) i_0^W E^W}$

and

$ \frac{\kappa}{\delta_W} \left(V_{APPLIED} - V_{DL}^W\right) = -c_{DL} \frac{\partial V_{DL}^C}{\partial t} + i_0^C E^C$

The material parameters are:

$i_0^W = i_0^C = 40\;\text{A / m}^2$

$\alpha = 0.4$

$F = 9.6485 \times 10^{-4}\;\text{J / V / mol}$

$V_{APPLIED} = -0.275\;\text{V}$

$R = 8.314\;\text{J / K / mol}$

$T = 298.0\;\text{K}$

$\kappa = 15.26\;\text{A / V / m}$

$\delta = 100 \times10^{-9}\;\text{m}$

$\delta_W = \delta / 2 \;\text{m}$

$C_{DL} = 0.3\;\text{A s / V / m}^2$

$c_{cu}^{\infty} = 250.0\;\text{mol / m}^3$

$\Gamma_1 = 0.012024 $

$\Gamma_2 = -0.0515 $

$D_{cu} = 5 \times 10^{-10} \;\text{m}^2\text{ / s}$

Integrating this gives:

$V_{DL}^W = -0.27499$

$V_{DL}^C = 0.0686$

$i_T^C = 2865.65$

$V_{IR}^W = -9.38943 \times 10^{-6}$

$V_{IR}^C = 9.389434 \times 10^{-6}$

$V_{CELL} = -0.3436$

The code for this is source:trunk/moffat/electrical/misc/1D_metal.py.

The evolution of the $V_{DL}$:

source:trunk/moffat/electrical/misc/voltageVersusTimeCopper.png

and the cupric concentrations at various times:

source:trunk/moffat/electrical/misc/copper.png

  • Posted: 2012-02-21 16:59 (Updated: 2012-03-22 13:26)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Symmetric Electrical Problem

Extending blog:MoffatElectrical1 by making the system totally symmetric

The changes are: 

14c14
< alpha = 0.5
---
> alpha = 0.4
30c30
<     return i0 * (numpy.exp(-alpha * Fbar * V) - numpy.exp((1 - alpha) * Fbar * V))
---
>     return i0 * (numpy.exp(-alpha * Fbar * V) - numpy.exp((2 - alpha) * Fbar * V))
34c34
<                       - (1 - alpha) * Fbar * numpy.exp((1 - alpha) * Fbar * V))
---
>                       - (2 - alpha) * Fbar * numpy.exp((2 - alpha) * Fbar * V))

gives:

$V_{DL}^W = -0.255$

$V_{DL}^C = 0.255$

$i_T^C = 5833$

$V_{IR}^W = -0.0191 $

$V_{IR}^C = 0.0191$

$V_{CELL} = -0.55$

The code for this is source:trunk/moffat/electrical/misc/1D_symmetric.py.

The evolution of the $V_{DL}$:

source:trunk/moffat/electrical/misc/voltageVersusTimeSymmetry.png

  • Posted: 2012-02-14 18:22 (Updated: 2012-03-22 13:28)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Electrical Problem including counter electrode

Including the effects of the counter electrode.

At the working electrode $W$:

$ i_T^W = -\kappa \vec{n}^W \cdot \nabla V_{IR}$

$ i_T^W = -c_{DL} \frac{\partial V_{DL}^W}{\partial t} + i_{F}^W $

$ i_F^W = i_0^W \left[ \exp{\left(-\frac{\alpha F V_{DL}^W}{R T} \right)} - \exp{\left(\frac{\left(2 -\alpha\right) F V_{DL}^W}{R T} \right)} \right]$

At the counter electrode $C$:

$ i_T^C = -\kappa \vec{n}^C \cdot \nabla V_{IR}$

$ i_T^C = -c_{DL} \frac{\partial V_{DL}^C}{\partial t} + i_{F}^C $

$ i_F^C = i_0^C \left[ \exp{\left(-\frac{\alpha F V_{DL}^C}{R T} \right)} - \exp{\left(\frac{\left(2 -\alpha\right) F V_{DL}^C}{R T} \right)} \right]$

The normals point out of the electrolyte.

In the bulk:

$ \nabla^2 V_{IR} = 0 $

and

$V_{CELL} = V_{DL}^W + V_{IR}^W - V_{DL}^C - V_{IR}^C$

$V_{APPLIED} = V_{DL}^W + V_{IR}^W$

$\delta_W V_{IR}^C + \delta_C V_{IR}^W = 0 $

$\delta_C + \delta_W = \delta$

where $\delta_C$ and $\delta_W$ are the distances from the counter and working electrodes to the reference electode, respectively. $\delta$ is the total distance between the electrodes.

If we do this the equations can be reduced to two coupled ODEs for $V_{DL}^C$ and $V_{DL}^W$:

$ -\frac{\kappa}{\delta_W} \left(V_{APPLIED} - V_{DL}^W\right) = -c_{DL} \frac{\partial V_{DL}^W}{\partial t} + i_0^W \left[ \exp{\left(-\frac{\alpha F V_{DL}^W}{R T} \right)} - \exp{\left(\frac{\left(2 -\alpha\right) F V_{DL}^W}{R T} \right)} \right]$

and

$ \frac{\kappa}{\delta_W} \left(V_{APPLIED} - V_{DL}^W\right) = -c_{DL} \frac{\partial V_{DL}^C}{\partial t} + i_0^C \left[ \exp{\left(-\frac{\alpha F V_{DL}^C}{R T} \right)} - \exp{\left(\frac{\left(2 -\alpha\right) F V_{DL}^C}{R T} \right)} \right]$

The material parameters are:

$i_0^W = i_0^C = 40\;\text{A / m}^2$

$\alpha = 0.4$

$F = 9.6485 \times 10^{-4}\;\text{J / V / mol}$

$V_{APPLIED} = -0.275\;\text{V}$

$R = 8.314\;\text{J / K / mol}$

$T = 298.0\;\text{K}$

$\kappa = 15.26\;\text{A / V / m}$

$\delta = 100 \times10^{-6}\;\text{m}$

$\delta_W = 50 \times10^{-6}\;\text{m}$

$C_{DL} = 0.3\;\text{A s / V / m}^2$

Integrating this gives:

$V_{DL}^W = -0.2666$

$V_{DL}^C = 0.0667$

$i_T^C = 2546$

$V_{IR}^W = -0.008345$

$V_{IR}^C = 0.008345$

$V_{CELL} = -0.350$

The code for this is source:trunk/moffat/electrical/misc/1D_counter.py.

The evolution of the $V_{DL}$:

source:trunk/moffat/electrical/misc/voltageVersusTimeCounter.png

  • Posted: 2012-02-14 11:43 (Updated: 2012-03-22 13:29)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

A few changes to Bill's code

Bill wants me to make four changes to his code, see source:trunk/kirkendall/A-V_phase-field.py@1457.

He needs to make four changes.

  • Update $V_m$ to be an expression involving $\sigma_{xx}$. This is done in r1458.
  • Save the data from the simulations for later use and read the saved files and plot the results independently from the simulations. See diff:@1458:1460. This now writes a data file directory of the form data/A-V_phase-field_17_36_53Thu02Feb2012. The data files are stored as numpy zip files .npz and numbered by their time step value. plotdata.py reads and plots from the latest data directory created.

The name of an alternative path can be specified as an argument 

$ python plotdata.py data/A-V_phase-field_17_33_42Thu02Feb2012

This will use the specified directory, while 

$ python plotdata.py

will use the latest created data directory. plotdata.py plots the data for the four main fields using the stored data.

  • Add an improved method for extracting the 0.5 level set. See r1461. This is done using scipy's interp1d. After plotting the four main fields for every stored timestep, plotdata.py will throw up a plot of the interface position.
  • Posted: 2012-02-02 11:45 (Updated: 2012-02-02 18:11)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Light Propogation

Jon is interested in adding a better light calculation to his PV code. At the moment the Beer-Lambert law/equation is being used for light. Firstly we need to understand what is wrong with the current technique. We've been discussing the use of the level set method / fast marching method as an alternative to using FiPy's basic first order upwind scheme, which is quite inefficient and is not working for some problems. Some more thoughts:

  • We don't have a second order accurate implicit scheme in FiPy.
  • Would using a second order explicit scheme be very inefficient. It is a one time cost since it is static, but it could get prohibitively bad for the sort of meshes Jon may need for accurate solutions.
  • FFM will be way more efficient and we have second order schemes available, but not on unstructured meshes.
  • Can we even solve the Beer-Lambert equation using the Eikonal equation?
  • We need to understand what is going wrong with the current solution that Jon showed me in 1D and why it seems to deviate form the analytical. I believe that it is grid spacing.

Given the Beer-Lambert equation,

$ \nabla \cdot \left( \vec{c} \Gamma \right) + \alpha \Gamma = 0 $

let's try and solve a 1D problem with a shock using implicit upwind. See trunk/misc/light.py. With a fixed grid things seem quite good even though this is a 1D scheme. Not sure if this carries over to 1D. Basically the velocity goes from 10 to 1 half way across the domain and we see a shock as expected ($\alpha$ is constant). The results are not so bad.

source:trunk/misc/light.png

The upwind implicit scheme can deal with shocks (at least in 1D). Let's observe how it deals with a large change in the mesh density (see trunk/misc/lightVariableMesh.py). Let's try a 100 fold change. A small jump occurs.

source:trunk/misc/lightVariableMesh.png

A ten fold increase. Seems fine.

source:trunk/misc/lightVariableMesh10.png

A thousand fold increase. Big jump.

source:trunk/misc/lightVariableMesh1000.png

The jumps are not surprising and consistent with a first order scheme. The tail has the right shape even though it has a large displacement. The PV 1D result did not have a tail that modeled the local solution correctly. Maybe that is to do with a gradient of changes rather than 1 jump in the mesh density.

Adding a gradual change example trunk/misc/lightVariableMeshGradual.py shows this

source:trunk/misc/lightVariableMeshGradual.png

This could be the issue with the PV code.

  • Posted: 2012-01-27 12:35 (Updated: 2012-01-27 18:01)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Tom's Electrical Problem

Tom has come up with a new problem that involves solving the electric field. Here are the equations at the electrode interface ($1$):

$ i_T^1 = -\kappa \vec{n} \cdot \nabla V_{IR}$

$ i_T^1 = -c_{DL} \frac{\partial V_{DL}}{\partial t} + i_{F} $

$ i_F = \frac{c_m^i}{c_m^{\infty}} \sum_j \theta_j n i_j^{\infty} F \left( \alpha_j, V_{DL} \right) $

$ V_{CELL} = V_{IR} + V_{DL} $

$ D_m \vec{n} \cdot \nabla c_m = \frac{i_F}{n F} $

and in the bulk:

$ \frac{\partial c_m}{\partial t} = \nabla \cdot D_m \nabla c_m + \nabla \cdot \left[\frac{z F}{R T} D_m c_m \nabla V_{IR} \right] $

$ \nabla^2 V_{IR} = 0$

on the top electrode ($2$):

$ V_{IR} =0 $

$ c_m = c_m^{\infty} $

At first glance, I didn't see that this system is closed. Since the $V_{IR}$ equation is steady state, and surrounded by Neumann boundary conditions it only needs one boundary to fix the value (it is an initial value problem for the slope of $V_{IR}$ rather than a boundary value problem for $V_{IR}$.

Solving for the 1D problem

We have (assuming the normal points out of the electrolyte)

$ i_T = -\kappa V_{IR} / \delta$

where $\delta$ is the distance between the electrodes

$ i_T = i_F - C_{DL} \frac{\partial V_{DL}}{\partial t} $

$ i_F = i_0 \exp{\left(-\frac{\alpha F V_{DL}}{R T} \right)} $

This leads to an ODE for $V_{IR}$

$\frac{\partial V_{IR}}{\partial t} = -\frac{1}{C_{DL}} \left[ i_0 \exp{\left(-\frac{\alpha F \left(V_{CELL} - V_{IR} \right)}{R T} \right)} +\kappa V_{IR} / \delta \right] $

This can be solved with scipy's ODE integrator. To do this I took an example <http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html> and used that. The code is source:trunk/moffat/electrical/1D.py.

The results are

source:trunk/moffat/electrical/voltageVersusTime.png

and

source:trunk/moffat/electrical/currentVersusTime.png

With the following numbers:

$i_0 = 40\;\text{A / m}^2$

$\alpha = 0.4$

$F = 9.6485 \times 10^{-4}\;\text{J / V / mol}$

$V_{CELL} = -0.275\;\text{V}$

$R = 8.314\;\text{J / K / mol}$

$T = 298.0\;\text{K}$

$\kappa = 15.26\;\text{A / V / m}$

$\delta = 100 \times10^{-6}\;\text{m}$

$C_{DL} = 0.3\;\text{A s / V / m}^2$

We get the following steady state values:

$i_F = 2294.7\;\text{A / m}^2$

$V_{IR} = -0.01503\;\text{V}$

$V_{DL} = -0.259\;\text{V}$

  • Posted: 2012-01-12 14:25 (Updated: 2012-02-09 09:54)
  • Author: wd15
  • Categories: (none)
  • Comments (2)

Installing Samrai

In order to use lsmlib in parallel, I need to install SAMRAI. Of course this is easier said than done. I downloaded Samrai 2.4.4 and did the following on bunter: 

$ mkdir build
$ cd build
$ ../configure --prefix=/users/wd15/.virtualenvs/lsmlib --with-F77=/usr/bin/gfortran
checking build system type... x86_64-unknown-linux-gnu
...
config.status: executing default commands
$ make
make NDIM=1 compile || exit 1
...
g++  -g -Wall  -fno-implicit-templates  -I. -I../../../include -I../../../../include  -I/usr/lib/openmpi/include      -c ../../../../source/toolbox/restartdb/HDFDatabaseFactory.C -o HDFDatabaseFactory.o
../../../../source/toolbox/restartdb/HDFDatabaseFactory.C: In member function ‘virtual SAMRAI::tbox::Pointer<SAMRAI::tbox::Database> SAMRAI::tbox::HDFDatabaseFactory::allocate(const std::string&)’:
../../../../source/toolbox/restartdb/HDFDatabaseFactory.C:20: error: ‘HDFDatabase’ was not declared in this scope
../../../../source/toolbox/restartdb/HDFDatabaseFactory.C:20: error: template argument 1 is invalid
../../../../source/toolbox/restartdb/HDFDatabaseFactory.C:20: error: invalid type in declaration before ‘=’ token
../../../../source/toolbox/restartdb/HDFDatabaseFactory.C:20: error: expected type-specifier before ‘HDFDatabase’
../../../../source/toolbox/restartdb/HDFDatabaseFactory.C:20: error: invalid conversion from ‘int*’ to ‘int’
../../../../source/toolbox/restartdb/HDFDatabaseFactory.C:20: error: expected ‘,’ or ‘;’ before ‘HDFDatabase’
../../../../source/toolbox/restartdb/HDFDatabaseFactory.C:21: error: invalid conversion from ‘int’ to ‘SAMRAI::tbox::Database*’
../../../../source/toolbox/restartdb/HDFDatabaseFactory.C:21: error:   initializing argument 1 of ‘SAMRAI::tbox::Pointer<TYPE>::Pointer(TYPE*, bool) [with TYPE = SAMRAI::tbox::Database]’
make[4]: *** [HDFDatabaseFactory.o] Error 1
make[4]: Leaving directory `/users/wd15/packages/SAMRAI/build/source/toolbox/restartdb'
make[3]: *** [libXd] Error 1
make[3]: Leaving directory `/users/wd15/packages/SAMRAI/build/source/toolbox'
make[2]: *** [libXd] Error 1
make[2]: Leaving directory `/users/wd15/packages/SAMRAI/build/source'
make[1]: *** [compile] Error 1
make[1]: Leaving directory `/users/wd15/packages/SAMRAI/build'
make: *** [library] Error 1

So that didn't work. I had the following exchange with someone on the SAMRAI list. I wrote 

Hi, I'm trying to build SAMRAI version 2.4.4. The configure step seems to have no issues, but the compilation step throws the following error.
  $ make
  make NDIM=1 compile || exit 1
  make[1]: Entering directory `/users/wd15/Documents/python/SAMRAI'
  (cd source && make libXd) || exit 1
  ...
  ...
  ...
  g++  -g -Wall  -fno-implicit-templates  -I. -I../../../include -I../../../include  -I/usr/lib/openmpi/include      -c HDFDatabaseFactory.C -o     HDFDatabaseFactory.o
HDFDatabaseFactory.C: In member function ‘virtual SAMRAI::tbox::Pointer<SAMRAI::tbox::Database>  SAMRAI::tbox::HDFDatabaseFactory::allocate(const std::string&)’:
  HDFDatabaseFactory.C:20: error: ‘HDFDatabase’ was not declared in this scope
  HDFDatabaseFactory.C:20: error: template argument 1 is invalid
  HDFDatabaseFactory.C:20: error: invalid type in declaration before ‘=’ token
  HDFDatabaseFactory.C:20: error: expected type-specifier before ‘HDFDatabase’
  HDFDatabaseFactory.C:20: error: invalid conversion from ‘int*’ to ‘int’
  HDFDatabaseFactory.C:20: error: expected ‘,’ or ‘;’ before ‘HDFDatabase’
  HDFDatabaseFactory.C:21: error: invalid conversion from ‘int’ to ‘SAMRAI::tbox::Database*’
  HDFDatabaseFactory.C:21: error:   initializing argument 1 of ‘SAMRAI::tbox::Pointer<TYPE>::Pointer(TYPE*, bool) [with TYPE =   SAMRAI::tbox::Database]’
  make[4]: *** [HDFDatabaseFactory.o] Error 1
  ...
  ...
  ...
  make[1]: *** [compile] Error 1
  make[1]: Leaving directory `/users/wd15/Documents/python/SAMRAI'
  make: *** [library] Error 1
Any ideas what the issue may be?
Other details:
  $ uname -a
  Linux bunter 2.6.32-5-amd64 #1 SMP Mon Oct 3 03:59:20 UTC 2011 x86_64 GNU/Linux
  $ lsb_release -a
  No LSB modules are available.
  Distributor ID: Debian
  Description:    Debian GNU/Linux 6.0.3 (squeeze)
  Release:        6.0.3
  Codename:       squeeze
All other packages are at standard versions for Debian Squeeze including the HDF5 libraries.
Thanks

and I got the reply 

Daniel,
HDFDatabase.h is guarded with a #ifdef HAVE_HDF5.  From the error message it appears that HAVE_HDF5 is not defined.  Grep for HAVE_HDF5 in SAMRAI_config.h to see if it is defined.  I suspect that something went wrong during configuration and you are not getting HDF5.  Another clue that this is likely the problem is that there is no hdf5 include directory on the compile line.
Let us know what you find.
--
Bill Arrighi
Lawrence Livermore National Lab

so I tried a new configure step to include hdf5 

$../configure --prefix=/users/wd15/.virtualenvs/lsmlib --with-F77=/usr/bin/gfortran --with-hdf5
checking build system type... x86_64-unknown-linux-gnu
...
checking for HDF5 installation... /usr
checking whether HDF5 link works... no
configure: error: HDF5 compile/link test failed

Andrew suggested debugging the configure script with the -x flag. The "set -x" was places on line 19422. Any way it showed that the configure script is trying to build and link the following c++ file (conftest.cppa): 

#include "hdf5.h"
#define FILE "file.h5"
int main() {
    hid_t       file_id;   /* file identifier */
    herr_t      status;
    /* Create a new file using default properties. */
    file_id = H5Fcreate(FILE, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
    /* Terminate access to the file. */
    status = H5Fclose(file_id);
}

with 

$ g++  -o conftest -g -O2    conftest.cpp  -lz   -ldl  -lm   -lz -lm
/tmp/ccn31txn.o: In function `main':
/users/wd15/packages/SAMRAI/build/conftest.cpp:10: undefined reference to `H5check_version'
/users/wd15/packages/SAMRAI/build/conftest.cpp:10: undefined reference to `H5Fcreate'
/users/wd15/packages/SAMRAI/build/conftest.cpp:13: undefined reference to `H5Fclose'
collect2: ld returned 1 exit status

if I do 

$ g++  -o conftest -g -O2    conftest.cpp  -lz   -ldl  -lm   -lz -lm -lhdf5
$ ./conftest

it works fine. Looking at the configure code, it looks like it wants a directory for --with-hdf5: 

../configure --prefix=/users/wd15/.virtualenvs/lsmlib --with-F77=/usr/bin/gfortran --with-hdf5=/usr

gets through the config without issue. However make fails with loads of 

g++  -g -Wall  -fno-implicit-templates  -I. -I../../../include -I../../../../include -I/usr/include  -I/usr/lib/openmpi/include      -c ../../../../source/toolbox/restartdb/HDFDatabase.C -o HDFDatabase.o
/usr/include/H5Gpublic.h: In static member function ‘static herr_t SAMRAI::tbox::HDFDatabase::iterateKeys(hid_t, const char*, void*)’:
/usr/include/H5Gpublic.h:148: error: too many arguments to function ‘hid_t H5Gopen1(hid_t, const char*)’
../../../../source/toolbox/restartdb/HDFDatabase.C:160: error: at this point in file
/usr/include/H5Gpublic.h:148: error: too many arguments to function ‘hid_t H5Gopen1(hid_t, const char*)’
../../../../source/toolbox/restartdb/HDFDatabase.C:180: error: at this point in file
../../../../source/toolbox/restartdb/HDFDatabase.C:204: error: invalid conversion from ‘herr_t (**)(hid_t, void*)’ to ‘void**’
/usr/include/H5Epublic.h:212: error: too many arguments to function ‘herr_t H5Eget_auto1(herr_t (**)(void*), void**)’
../../../../source/toolbox/restartdb/HDFDatabase.C:204: error: at this point in file
/usr/include/H5Epublic.h:216: error: too many arguments to function ‘herr_t H5Eset_auto1(herr_t (*)(void*), void*)’
...
make[4]: *** [HDFDatabase.o] Error 1
make[4]: Leaving directory `/users/wd15/packages/SAMRAI/build/source/toolbox/restartdb'
make[3]: *** [libXd] Error 1
make[3]: Leaving directory `/users/wd15/packages/SAMRAI/build/source/toolbox'
make[2]: *** [libXd] Error 1
make[2]: Leaving directory `/users/wd15/packages/SAMRAI/build/source'
make[1]: *** [compile] Error 1
make[1]: Leaving directory `/users/wd15/packages/SAMRAI/build'
make: *** [library] Error 1
  • Posted: 2012-01-11 11:46 (Updated: 2012-01-11 14:48)
  • Author: wd15
  • Categories: (none)
  • Comments (1)

python setup.py develop issues

If python setup.py develop does not seem to be updating easy-install.pth correctly, then try running pip install setuptools --upgrade and rerunning python setup.py develop.

  • Posted: 2011-10-25 17:23
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Installing Ubuntu on the Acer Aspire One Netbook

I purchased the Aspire One at Costco for $300 with the hope that Ubuntu would install smoothly without glitches. Unsurprisingly, it wasn't a straight forward experience. The first problem is that the monitor aspect ratio was wrong. All text was squashed in the vertical plane. Luckily, the first 11.04 update fixed these problems. It also prompted me to install two proprietary drivers, which I did. Before the drivers were installed I couldn't load Unity. After they were installed I could load Unity.

At first this all seemed fine apart from a weird AMD watermark in the bottom right hand corner of the screen. I did the googles and found the following, which fixed the issue,

<http://askubuntu.com/questions/25519/how-to-remove-amd-unsupported-hardware-without-reinstalling-the-driver>

Totally gross, but it works.

This all seemed totally fine, but after awhile the computer just started to freeze and needed to be rebooted. This turned out to be harder to fix in terms of finding answers. I believe that it is somehow connected to the wireless as it always seems to happen when the wireless is starting up, but who knows. Anyway I found this helpful post.

<http://ubuntuforums.org/showpost.php?p=10571078&postcount=16>

This worked (thank god). I had to boot into safe mode in order to have enough time to edit the necessary file. After this things all seem to be fine other than Unity not working. I'm not too bothered about that because I'm perfectly happy with classic mode. In fact classic seems a bit snappier. I'm not really sure about Unity anyway (what does it really add?).

Overall, compared with my now defunct Dell mini, it's a lot faster. Firefox works properly, which it didn't on the mini. I'm a happy bunny once more for $300 and some head scratching.

  • Posted: 2011-09-14 11:55 (Updated: 2011-09-15 11:06)
  • Author: wd15
  • Categories: (none)
  • Comments (1)

Generating Passwords 

   pwgen -Bcny 12 1
  • Posted: 2010-11-23 10:38
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Making lam work across multiple machines

To make lam work across multiple machine take the following steps:

  • Create a bhost file 
      bunter.nist.gov cpu=4
  hobson.nist.gov cpu=4
  alice.nist.gov cpu=4
  • $ export LAMRSH=ssh
  • $ lamboot -ssi boot_rsh_ignore_stderr 1 bhost
  • $ mpirun -np 12 python script.py

This should then work transparently across the three machines. The stdout streams from the other machines didn't appear on the screen, but they were certainly running in parallel.

  • Posted: 2010-10-13 10:14 (Updated: 2010-10-28 14:41)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Rongpei's Core-Shell problem

Rongpei visited from 19th to the 21st of July. He is interested in porting his problem to FiPy. The documents with the required equations are 1 and 2.

  • Posted: 2010-07-22 12:56 (Updated: 2010-07-22 12:56)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

mpi4py, trilinos, lam and i686

We needed both trilinos amd mpi4py working against lam on the i686 machines so we can test effectively and run parallel jobs unopposed. I built mpi4py on bunter by

  • making sure absolutely everything with mpi4py in it was removed from site-packages
  • Adding 
    [i686-lam]
mpi_dir              = /usr/lib/lam
mpicc                = /usr/bin/mpicc.lam
mpicxx               = /usr/bin/mpic++.lam

to mpi.cfg and running Python setup.py build --mpi=i686-lam in the i686-lam virtualenv.

I built trilinos using i686-lam with 

../configure CXX=mpic++.lam CC=mpicc.lam F77=mpif77.lam CXXFLAGS=-O3 CFLAGS=-O3 FFLAGS="-O5 -funroll-all-loops -malign-double" --enable-epetra --enable-aztecoo --enable-pytrilinos --enable-ml --enable-ifpack --enable-amesos --with-gnumake --enable-galeri --enable-mpi --with-mpi=/usr/lib/lam --cache-file=config.cache --prefix=/users/wd15/.virtualenvs/i686-lam
  • Posted: 2010-07-13 16:03 (Updated: 2010-07-13 16:05)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Problem in virtualenvwrapper

I've had some issues with virtualenvwrapper. Whenever I start up a new shell, the shell hangs when virtualenvwrapper is installed. I debugged this by putting a "set -x" at the top of "virtualenvwrapper.sh". "virtualenvwrapper.sh" is sourced from ".bashrc". Anyway, the debug information pointed the finger at line 193 of "virtualenvwrapper.sh" 

    fi
    trap "rm '$file' >/dev/null 2>&1" EXIT
    echo $file

It hangs on the "rm" command as my "rm" command is aliased to "rm -i". I am not sure what the right thing to is. For now I have just set 

   fi
   trap "\rm '$file' >/dev/null 2>&1" EXIT
   echo $file

which deals with the issue. Is the problem with "virtualenvwrapper.sh" because it is assuming that "rm" is the default version or should my aliases come last in the .bashrc file because commands should be the default during the execution of the "rc" type files?

  • Posted: 2010-06-10 11:18 (Updated: 2010-06-10 11:18)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Bill's changes to paper April 6 2010

Stuff that I am not sure about:

  • Think hard about the entropy stuff in the abstract
  • comment about what people believe regarding the viscosity of metals
  • Posted: 2010-04-07 10:32 (Updated: 2010-04-07 10:37)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Convergence

There are three simulations that need to be compared to test convergence 267 (250x360), 271 (500x700) and 272 (125x180). Simulation 272 may not work, I can't remember whether this is an adequate size. I may need to do an even larger simulation, but will need access to more CPUs for that.

  • Posted: 2010-03-03 14:55
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Entropy Production

Given the free energy $f$ per unit volume, the entropy production for our problem should be,

$ S_p = \frac{M}{T^2} \left( \partial_j \left( \mu_1^{NC} - \mu_2^{NC} \right) \right)^2 + M_{\phi} \left( \frac{1} {T} \frac{\partial f}{\partial \phi} - \epsilon_{\phi} \partial_i^2 \phi \right)^2 + \frac{\nu}{2 T} \left( \left(\partial_i u_j \right)^2 + \left( \partial_i u_j \right) \left( \partial_j u_i \right) \right) $

where the $i$ and $j$ are summed. I think the entropy production for the last two terms is correct. Right now, I'm confused by the first term.

  • Posted: 2010-01-20 11:34 (Updated: 2010-01-20 14:29)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Simulations with van Leer convection term

The upwind convection term appeared to help the low concentration simulation work. However, the upwind convection term is inaccurate though. Let's use the van Leer convection term instead.

ID Notes Issues Status $N$ movie
267 base none benson, hobson, alice (stopped) 12
268 $X_1 = 0.1 X_1^{equ}$ none renfield, sancho, kato (stopped) 12
269 $\nu_f = 2 \times 10^{-2}$ crashes luggage (stopped) 12
270 $X_1 = 0.5 X_1^{equ}$ none luggage (stopped) 12
  • Posted: 2009-12-22 13:22 (Updated: 2010-01-14 10:47)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Simulations with a negligible radial displacement

The radial displacement issue has been troubling me. Most of the simulations have had a radial displacement of dr / 10 and some others have dr. Some simulations don't run at all and an instability can clearly be seen to emanate from the r=0 location. I tested the base simulation with dr * 1e-10 just to figure out why I had trouble in the past. Needless to say the simulation seem to be running fine. I have just launched some more simulations with dr * 1e-10 to try and diagnose.

ID Notes Issues Status $N$ movie
264 base None benson, hobson, alice (stopped) 12
265 $X_1 = 0.1 X_1^{equ}$ worked with upwind convection term renfield, sancho, kato (stopped) 12
266 $\nu_f = 2 \times 10^{-2}$ crashed luggage (stopped) 12
  • Posted: 2009-12-17 12:32 (Updated: 2009-12-22 11:31)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Useful Simulations

ID Notes diff Issues Status $N$ movie
222 base diff benson, hobson, alice (stopped) 12
225 $\nu_f = 2 \times 10^{-2}$ diff luggage (stopped) 12
226 $\nu_f = 2 \times 10^{-1}$ diff crashed renfield, sancho, kato (stopped) 12
233 222: $X_1 = 0.1 X_1^{equ}$ $\bar{M}_f = 1 \times 10^{-6}$ diff cluster 7.0 (stopped) 4
240 222: $X_1 = 1.0 X_1^{equ}$ $\bar{M}_f = 1 \times 10^{-6}$ diff crashed luggage (stopped) 12
252 243: $\nu_f = 2 \times 10^{-1}$, shift=dx, $\sigma_{\text{CFL}} = 0.01$ diff cluster 128 (stopped) 4
257 251: $X_1 = 1.0 X_1^{equ}$, $\bar{M}_f = 1 \times 10^{-6}$, shift=2*dx, $\sigma_{\text{CFL}}=0.05$ diff luggage (stopped) 12
  • Posted: 2009-12-15 15:32 (Updated: 2009-12-17 11:59)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Edwin Schwalbach Model

For the last three days Eddie and I have been working on his VLS growth model. We've been trying to implement it in FiPy?. We have successfully implemented a two phase, pure, 1D system. His notes are here

source:trunk/vls/Single_component_perturbation.pdf@1031

The code is here

source:trunk/vls@1031

Essentially, this is simply the van der Waals model with a free energy of the form

$ F_m = w \left( \rho - \rho_l \right)^2 \left( \rho - \rho_v \right)^2 $

We implement the model using the block matrix approach with $u$, $\rho$ and $\mu^{NC}$ as the base variables to "solve for". The $\mu^{NC}$ is linearized with,

$ \mu^{NC} = \mu_0^* + \rho \frac{\partial \mu}{\partial \rho}^* - \kappa \nabla^2 \rho $

where,

$\mu_0^* = \mu^* - \rho^* \frac{\partial \mu}{\partial \rho}^* $

Remember,

$F_m = \frac{m F_v}{\rho}$

and

$P = -F_v + \rho \mu$

and

$\frac{ \partial \mu}{\partial \rho} = \frac{1}{\rho} \frac{ \partial P}{\partial \rho}$

Explicitly,

$P = \frac{2 W \rho^2}{m} \left( \rho - \rho_l \right) \left( \rho - \rho_v \right) \left( 2 \rho - \rho_l - \rho_v \right) $

and

$\frac{\partial P}{\partial \rho} = \frac{2 W \rho}{m} \left( \rho \left( 2 \rho - \rho_l - \rho_v \right)^2 + 2 \left(3 \rho - \rho_l - \rho_v \right) \left( \rho - \rho_l \right) \left( \rho - \rho_v \right) \right) $

Source Terms

Take an equation say,

$ \frac{\partial \phi}{\partial t} = S \left( \phi \right) $

In general $ S \left( \phi \right) $ can be implemented explicitly (i.e. $S \left( \phi^* \right) $ where $*$ is the previous sweep value. Now, to implement this implicitly, we need to linearize,

$ S \left( \phi \right) = S \left( \phi^* \right) - \phi^* \frac{\partial S}{\partial \phi}^* + \phi \frac{\partial S}{\partial \phi}^* $

or

$S = S_c + \phi S_p$

where

$S_p = \frac{\partial S}{\partial \phi}^*$

In general $S_p \le 0$ for stability reasons, though one can often relax this somewhat if the source is not dominant. If the source is not negative, then $S$ should be evaluated entirely explicitly. The reason $S_p$ must be negative is the following. Take,

$ \frac{\partial \phi}{\partial t} = \lambda \phi $

then the discretized implicit form is

$ \phi = \frac{\phi^O}{1 - \lambda \Delta t} $

Say $\phi^O$ is positive, then this equation should increase the value of $\phi$, but if $ \lambda \Delta t > 1$ then we get a negative value of $\phi$. If $\lambda$ is negative, there are no stability criterion. Conversely, if $\lambda<0$ then the explicit scheme,

$ \phi = \left( 1 + \lambda \Delta t \right) \phi^O $

will drive the solution below $0$ for sufficiently negative $\lambda$. For very gentle intro,

clcik here

Pressure velocity coupling

We use the Rhie-Chow interpolation method to deal with this issue since $P$ and $u$ are not on staggered grids but collocated. Essentially, the Rhie-Chow interpolation scheme corrects the face velocities in the continuity equation by the use of a corrected velocity that approximates solving the momentum equation on the faces. Let's take the discretized form of the x-momentum equation:

$ a_P u_P = \sum_{nb} a_{nb} u_{nb} - V_P \rho \nabla_x \mu^{NC} $

where $P$ is a cell value. For small enough grid, we would like

$ a_f u_f - \sum_{nb} a_{nbf} u_{nbf} + V_f \rho_f \nabla_x \mu^{NC} = \overline{a_P u_P} - \overline{\sum_{nb} a_{nb} u_{nb}} + \overline{V_P \rho \nabla_x \mu^{NC}} $

that is to say that we would like the residual of the momentum equation evaluated at the face to be equal to the residual of the averaged residual over the neighboring cells. We assume that the $\sum_{nb}$ terms are unimportant and say that,

$ u_f = \frac{\overline{a_P u_P}}{a_f} - \frac{V_f}{a_f} \rho_f \nabla_x \mu^{NC} + \frac{\overline{V_P}}{a_f} \overline{ \rho \nabla_x \mu^{NC}} $

We approximate $a_f = \overline{a_P} $ and $V_f = \overline{V_P}$ the we get,

$ u_f = \overline{u_P} - \frac{\overline{V_P}} {\overline{a_P}} \left( \rho_f \nabla_x \mu^{NC} - \overline{ \rho \nabla_x \mu^{NC}} \right) $

The continuity equation can be written,

$ \frac{\partial \rho}{\partial t} + \nabla \cdot \left( u_f \rho \right) = 0 $

Substituting in gives a diffusion term for $\mu^{NC}$, which acts to stabilize the solution. See source:trunk/vls/input.py@1031#L78.

Parsitic Currents

Discretize using $\rho \nabla \mu$ rather than $\nabla P$. This induces dissipation assoicated with momentum. See source:trunk/reactiveWetting/paper/jacqmin.pdf@1031.

Movie of solution to equations

  • Posted: 2009-12-11 13:40 (Updated: 2009-12-11 16:04)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Hydrodynamic limit with elevated interface diffusion

Increasing $\bar{M}_i$ certainly makes a big difference.

source:trunk/reactiveWetting/184/dropRadius.png

The image shows the drop radius against time. The inertial time scale scales time and initial radius of the drop scales the spreading distance $r$. There is an order of magnitude change in time between when the spreading occurs for the old hydrodynamic case and the new ones. I will need to plot the spreading rate curves to see how this knocks the solution away from Spelt. Also, there seems to be some difference between $\bar{M}_i=10^{-7}$ and $\bar{M}_i=10^{-5}$. Note that $t_{\text{inertial}}= 7 \times 10 ^{-8}$. The interface diffusion time scale is $t_{\text{ini}}=3\times 10^{-6}$ so we probably need to raise the interface diffusion coefficient by two orders from the fluid value to get it be less than the inertial time scale. If indeed that does matter. It certainly seems to.

  • Posted: 2009-09-25 18:50 (Updated: 2009-12-10 10:58)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Interface Diffusion

The following figure shows melting rates for various $\bar{M}_i$ with the sharp interpolation scheme.

source:trunk/reactiveWetting/146/box.png

Firstly, the melting rates agree well with analytical rates. Secondly, as pointed out in the previous post, $\bar{M}_i$ does not have much influence on the melting. It has some influence at early times, which is evidently related to hydrodynamics. By about $t = 10^{-5}$, the curves for $\bar{M}_f = 10^{-7}$ have almost all collapsed onto each other.

  • Posted: 2009-09-25 18:31 (Updated: 2009-09-25 18:55)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Interpolation Schemes for M

Carrying on from blog:SolidDiffusion1 it is clear that we have to modify the interpolation scheme. It showed that the time scale in the interface has to be fast enough to not influence the bulk diffusion properties. It is clear that for the purposes of diffusion that having the interface diffusion coefficient equal to the bulk diffusion is adequate. However, it may be necessary for hydrodynamic reasons to have faster diffusion in the interface.

source:trunk/reactiveWetting/146/interpolationScheme.png

In this figure $k$ refers to

$\bar{M} = \bar{M}_s^{\phi^k} \bar{M}_f^{1 - \phi^k} $

and $l$ refers to

$\bar{M} = \bar{M}_b + 2^{2l} \left( \bar{M}_i - \bar{M}_b^* / 2 \right) \phi^l \left(1 - \phi \right)^l $

where

$\bar{M}_{b} = \bar{M}_s \phi + \bar{M}_f \left( 1 - \phi \right) $

This red curve demonstrates that allowing the interface to have the fluid diffusion coefficient all the way through gets us back to the bulk limited value. It also demonstrates that the $l$ based interpolation scheme doesn't really work because it over influences the diffusion in the bulk and allows it to have much faster diffusion than it otherwise should do. This is due to the extent that $\phi$ extends out beyond the interface region. As shall be seen, a sharp cut off at $\phi=0.05$ ameliorates this issue.

As shall be seen in the subsequent post, the black line is actually agrees with the analytical solution and is bulk limited. It is also clear that as $l$ goes from 2 to 1 the melting is faster due to this issue.

  • Posted: 2009-09-25 18:05 (Updated: 2009-09-25 18:56)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Analytic Melting Solution

The following os an analytical melting solution from Bill and Geoff

$l = K \sqrt{4 D t}$

where $l$ is the distance the interface has moved. $K$ is given by solving

$K + S \frac{\exp \left( -K^2 \right) }{1 - \text{erf} \left( K \right) } \frac{1}{\sqrt{\pi}} = 0$

and $S$ is given by,

$S = \frac{X_1^l - X_1^{l, \text{equ}}}{X_1^s - X_1^l}$

In this case $X_1^s = X_1^{s, \text{equ}}$.

  • Posted: 2009-09-25 17:40 (Updated: 2009-09-25 17:51)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Influence of Solid Diffusion

The 1D simulations have demonstrated that the solid diffusion coefficient has a major impact on the rate of melting. The following image shows melting with $\bar{M}_s$ being varied.

source:trunk/reactiveWetting/146/solidDiffusion.png

It is clear that $\bar{M}_s$ has a considerable impact on the melting. This is due to the limiting nature of the interface diffusion. The interpolation scheme is

$\bar{M} = \bar{M}_s^{\phi} \bar{M}_f^{1 - \phi}$

which drives the value of the interface diffusion down a lot in the solid. The system is bulk limited as $\bar{M}_s$ approaches $10^{-4}$ as can be seen from the black dashed line. We though that the system was limited by the phase filed motion, but in fact it obviously isn't. The time scale for the phase field equation is given by,

$t_{\phi} = \frac{m T}{M_{\phi} p \left( \phi = 0.5 \right) \left( A_1 - A_2 \right) \rho \Delta X_1} \approx 8 \times 10^{-10}$

which is faster than the bulk diffusion time scale even when $\bar{M}_f = 10^{-2}$,

$t_{bulk} = \frac{l^2}{K^2 4 D} \approx 8 \times 10^{-9}$

where $l$ is the distance moved which I have chosen such that $l=X$. A reasonable choice in this case. Hardly surprising that this system is bulk limited.

So why does the solid diffusion eventually slow things down? Maybe we could say that the interface has to have reached some sort of equilibration before the motion occurs. This would lead to a time scale on the order of,

$t_{int} = \frac{X^2}{D} \approx 3 \times 10^{-11} $

This is probably too fast, it probably also has a prefactor like the $K$ in the bulk. However, it should be faster than the bulk kinetics. This seems obvious now for some reason. Anyway, it makes sense that we need to reduce the solid diffusion by three orders of magnitude before the interface kinetics start to become slower than the bulk kinetics.

  • Posted: 2009-09-25 16:40 (Updated: 2009-09-25 18:56)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Time scales asccoiated with diffusion

Having done a large number of 1D simulations it has become apparent that the comments in blog:SolidDiffusion entry seem to hit the nail on the head and are probably correct. The following image shows the liquid fraction from a number of 1D simulations that start with an system that has a 1/3 of each phase, but with a low concentration is the liquid. The liquid should expand to 1.12 of its original volume.

source:trunk/reactiveWetting/146/solidDiffusion.png

Clearly, the solid diffusion coefficient has no effect on the melting as long as it is above a given threshold. This agrees with our earlier time scale analysis in the previous blog entry. In the following image the interpolation scheme is varied. It is clear that if the fluid diffusion is used in the interface the solid diffusion coefficient does not effect the melting.

source:trunk/reactiveWetting/146/interpolationScheme.png

Given the sharp interpolation scheme, is it still possible to have interface effects. Well, maybe, if $t_{int} > t_{diff}$ we could be in for some problems. To avoid this:

$t_{int} < t_{diff}$

$\frac{X^4 R}{X_1 X_2 \epsilon_i m \rho L^2} < 1$

For the numbers that we are using, the expression on the right hand side is approximately 0.01. That means we are in good shape for the 2D simulations from the perspective of not being interface limited for melting.

However, the absolute value for $t_{int}$ is $6.7 \times 10^{-6}$, which is a lot slower than many of the other time scales in the system. The inertial time scale for example is two orders of magnitude faster. This means that we are in extremely bad shape since interface effects are still way slower than the spreading. Ideally for a physical system, interface effects should be much faster than the spreading. This would indicate that we could be interface limited for spreading, which is bad.

  • Posted: 2009-09-21 11:56 (Updated: 2009-09-21 12:03)
  • Author: wd15
  • Categories: (none)
  • Comments (2)

Melting dependence on solid diffusion coefficient

It is becoming apparent that the rate of melting is not independent of the solid diffusion coefficient. The following link is to a movie where the system has the same viscosity and diffusion coefficient (2E-3 and 2E-2) in all phases.

trunk/reactiveWetting/147/movie/movie.gif

The dynamics is quite fast and the pressure and chemical potential equilibration occur at similar time scales. There is a nice concentration gradient and the phase filed relaxes back as the liquid melts. Note how the chemical potentials equilibrate nicely through the interface and large gradients are quickly removed.

The following is a simulation where the solid diffusion coefficient is 2E-11 (down from 2E-2).

trunk/reactiveWetting/148/movie/movie.gif

From a pure diffusion argument the solid diffusion coefficient should have no influence on the melting rate of the liquid. Clearly in this example it is having an effect. If you look at frames from similar times in both the above simulations, the rate of melting is clearly different and the concentration profile in the liquid has no gradients in the second movie. Also, note that the chemical potential has not equilibrated through the solid liquid interface. This is very interesting.

This does seem somewhat disturbing. However, if we look at the diffusion equation, there are two time scales, the first being the diffusion time scale through the liquid,

$ t_{\text{diff}} = \frac{m \rho L^2}{R \bar{M}_f} $

where $L$ is the length of the liquid phase (it can't be shorter than this because the diffusion equation creates a profile across the entire liquid domain!). The second time scale is and interface time scale associated with the fourth order term,

$ t_{\text{int}} = \frac{X^4}{\bar{M}_s X_1 X_2 \epsilon_i} $

where $X$ is the depth of the interface. $t_{\text{diff}}$ is very fast about 9.27E-9 equivalent to the other time scales in the system such as the phase field and inertial flow scales. Since the chemical potential needs to equilibrate through the entire interface I am going to postulate that I can use \bar{M} in the solid to evaluate $ t_{\text{int}}$ giving a time scale of 2.35E-4. This is very slow. This represents a worse case since $\bar{M}$ in the interface is interpolated between the liquid and solid, but there are still regions of the interface that require a time scale close to this for equilibration.

This argument suggests that $\bar{M}_s$ needs to be large enough or $\bar{M}_l$ needs to be small enough such that $ t_{\text{int}} < t_{\text{diff}}$ and only at this point is melting independent of $\bar{M}_s$. This cross over point for our numbers is roughly when $\bar{M}_s$ is three or four orders of magnitude less than $\bar{M}_f$.

The other thing to note is that for a physical system $X$ is much smaller so for a realistic system $ t_{\text{int}} \ll t_{\text{diff}} $ in general.

Another hand wavy argument may be to say that $ t_{\text{int}} $ represents the time scale for the initial formation of the interface and if this is really slow all bets are off.

I am beginning to realize why it is important to have a more realistic interface depth. The above argument suggests the diffusion coefficient should be large enough such that the time scale associated with the fourth order term in the diffusion equation is always faster than other relevant time scales. All the spreading dynamics will be screwed up otherwise.

  • Posted: 2009-09-17 18:19 (Updated: 2009-09-21 11:22)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

MPI command

mpirun -np 12 -machinefile machinefile ./input.py >> out 2>&1 &

  • Posted: 2009-09-10 18:21 (Updated: 2009-09-10 18:23)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Parallel on the Cluster

Results for parallel simulations on the cluster

  • 1 processor, r005
  • sweep times: 29.3520491123 23.3660550117 23.3256959915 23.1911690235 23.1299791336 22.9121041298 22.8525710106 23.0343871117 22.8907659054 23.039219141 23.1625990868 23.013864994
  • 2 processors, r006 r005
  • sweep times: 15.5924270153 12.1713759899 12.2448968887 12.2202451229 12.2246148586 12.2457139492 12.3183469772 12.3248221874 12.215831995 12.2335240841 12.2320721149
  • 4 processors, r006 r007 r006 r005
  • sweep times: 21.6786100864 6.6802740097 6.63790607452 6.84371995926 6.65304708481 7.14012789726 6.46965503693 6.99537992477 6.83790421486 6.62651705742 6.71980905533 6.77136397362 6.41707897186
  • 8 processors, r007 r008 r006 r007 r006 r008 r007 r005
  • sweep times: 17.380297184 17.7295508385 17.6523389816 17.4138000011 31.6395330429 17.4052929878 24.2372221947 31.4160881042 3.41191792488 3.63104009628 3.52340602875 10.1916849613 17.3036520481 10.5115458965 3.95448207855 3.49768590927
  • 16 processors, r007 r012 r008 r006 r007 r006 r022 r013 r012 r012 r008 r007 r011 r020 r010 r005
  • sweep times: 30.4704580307 15.7536921501 16.4713590145 17.762444973 3.03759098053 3.31603288651 15.6013848782 9.57383298874 43.7733981609 3.11927390099 3.11927390099 15.9374248981 2.89413809776 2.70106816292 2.70635294914 16.7859950066

Results demonstrate:

  • good scaling from 1 to 4 processors
  • something weird happening at 8 processors
  • fastest times for 8 processors does scale correctly, might indicate network traffic
  • similarly, for 16 processors there is a 2.7s sweep
  • Posted: 2009-09-10 14:17
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Melting Simulations

source:trunk/reactiveWetting/128/dropRadiusMelt.png

  • 125 is the base simulation
  • 126 is with 50% reduction of $X_1$.
  • 127 is with $M_{\phi}$ increased by an order of magnitude
  • Not a great deal of variation will take $M_{\phi}$ higher.
  • 126 is has the slowest spreading. makes sense since a reduction in concentration reduces spreading rate (slightly more melting maybe).
  • 127 is faster as $M_{\phi}$ is increased
  • Posted: 2009-09-10 12:14
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Base Simulation with Higher Viscosity

source:trunk/reactiveWetting/128/dropRadiusBase.png

  • 119 is the original base simulation with $\mu_s = 2 \times 10^4$
  • 125 has a higher solid viscosity $\mu_s = 2 \times 10^4$ and a deeper solid region
  • 128 is the same as 125, but with a tolerance of $R = 1 \times 10^{-2}$ instead of $R = 1 \times 10^{-1}$.
  • Figure demonstrates that 128 and 125 map fairly well indicating that $R = 1 \times 10^{-1}$ is adequate
  • 119 has slower spreading dynamics than 125, which is counter intuitive, this may have something to do with the deeper solid region
  • It is clear from the movies for source:trunk/reactiveWetting/119/movie/119.gif and source:trunk/reactiveWetting/125/movie/125.gif, that the solid wave are substantially decreased with the higher solid viscosity.

  • Posted: 2009-09-10 11:45 (Updated: 2009-09-10 11:48)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Dissolution Simulations

The simulations in the table below are being run as I write up the paper, most are very close to production quality.

ID Notes diff Issues Status $N$ movie
119 hydrodynamic limit solid waves, wrong value of X1 close to equilibrium stopped 32 movie
120 removed pytables dependencies to run on rareton diff stopped
121 finer grid density diff better initial conditions required stopped 32 movie
122 $\rho_i = 0.1 \left( \rho_i^s - \rho_i^l \right) \rho_i^l$ and deeper solid diff will induce solidification stopped 6 movie
123 $\rho_i = 0.5 \left( \rho_i^s - \rho_i^l \right) \rho_i^l$ diff drop wets top surface stopped 11
124 $\rho_i = 0.5 \left( \rho_i^s - \rho_i^l \right) \rho_i^l$ and deeper vapor diff will induce solidification stopped 11 movie
125 $\mu_s^{\text{new}} = 100 \mu_s $, hydrodynamic limit diff less waves, but better initial conditions required stopped 12
126 $X_1^{l, \text{new}} = 0.5 X_1^l$ diff no melting stopped 6
127 $M_{\phi}^{\text{new}} = 10 M_{\phi}$ diff increasing mobility will not induce melting stopped 12
128 $R_{\text{tol}}^{\text{new}} = 0.1 R_{\text{tol}}$ diff better initial conditions required stopped 12
129 $\mu_s^{\text{new}} = \mu_l$ and $\bar{M}_s^{\text{new}} = \bar{M}_l$ diff stopped 12
130 $M_{\phi}^{\text{new}} = 10 M_{\phi}$ diff increasing mobility will not induce melting stopped 4
131 $M_{\phi}^{\text{new}} = 10 M_{\phi}$ diff increasing mobility will not induce melting stopped 4
132 $\bar{M}_f = 1 \times 10^{3} \bar{M}_f$ diff X1 has weird profile due to smoothing stopped 12
133 smoothing $\rho$ and $X_1$ diff worked, but fast diffusion not required sropped 4
134 1D to see eventual melting diff not required smaller nx works stopped 4
135 $\bar{M}_f = 1 \times 10^{5} \bar{M}_f$ diff duffusion fast enough, large dissolution impossible stopped 4
136 $\bar{M}_f = 1 \times 10^{5} \bar{M}_f$, 1D, nx=100 diff not required smaller nx works stopped 4
137 1D, nx=50, $X_1^{l, \text{new}} = 0.1 X_1^l$ diff wrong initial conditions stopped 4
138 hydrodynamic limit diff stopped 12
139 1D, nx=10, $X_1^{l, \text{new}} = 0.1 X_1^l$, shift=0.0 diff nx < 29 doesn't seem to work or different shifts stopped 4
140 $X_1^{l, \text{new}} = 0.1 X_1^l$ diff stopped 4
141 $\bar{M}_f = 1 \times 10^{-6} $ diff stopped 12
142 nx=2 diff calculation broke down with nx =2 stopped 4
143 same as 137 but corrected initial conditions diff stopped 4
144 same as 141 but $\bar{M}_f = 1 \times 10^{-5} $ diff stopped 4
145 143: $\bar{M}_s = 1 \times 10^{-2} $ diff stopped 4
146 63: with modifed unsegregatedSolver from 145 $X_1^{l, \text{new}} = 0.1 X_1^l$ and NC1, NC2 = 0 between sweeps diff stopped 1 movie
147 146: $\bar{M}_s = \bar{M}_f = 1 \times 10^{-2}$ diff stopped 1 movie
148 147: $\bar{M}_s = 1 \times 10^{-11}$ diff stopped 1 movie
149 147: $\mu_s = 2 \times 10^{6}$ diff stopped 1 movie
150 148: $\mu_s = 2 \times 10^{6}$ diff stopped 1 movie
151 150: $\bar{M}_s = 1 \times 10^{-1}$ diff stopped 1
152 150: $\bar{M}_s = 1 \times 10^{-3}$ diff stopped 1
153 150: $\bar{M}_s = 1 \times 10^{-4}$ diff stopped 1
154 150: $\bar{M}_s = 1 \times 10^{-5}$ diff stopped 1
155 150: $\bar{M}_s = 1 \times 10^{-6}$ diff stopped 1
156 150: $\bar{M}_s = 1 \times 10^{-7}$ diff stopped 1
157 150: $\bar{M} = \bar{M}_l$ when $\phi < 0.97$ diff stopped 1
158 157: $\bar{M} = \bar{M}_s^{\phi^2} \bar{M}_s^{1 - \phi^2} $ diff stopped 1
159 154: half interface depth and N=1600 diff stopped 1
160 154: $N=1600 $ diff stopped 1
161 159: $\epsilon_i^{\text{new}} = 4 \epsilon_i $ diff stopped 1
162 149: $N=1600 $ diff stopped 1
163 154: $\epsilon_i^{\text{new}} = \epsilon_i / 4 $ diff stopped 1
164 154: $\epsilon_i^{\text{new}} = \epsilon_i / 16 $ diff stopped 1
165 154: $\epsilon_i^{\text{new}} = \epsilon_i / 64 $ diff stopped 1
166 150: $\bar{M}_{\text{int}} = 1 \times 10^{-3}$ diff stopped 1
167 166: $l=1$ diff stopped 1
168 167: box interpolation diff stopped 1
169 167: $\bar{M}_f = 1 \times 10^{-6}$, $l=2$ diff stopped 1
170 167: $\bar{M}_f = 1 \times 10^{-5}$, $l=2$ diff stopped 1
171 170: $\bar{M}_f = 1 \times 10^{-6}$, $l=1$ diff stopped 1
172 170: $\bar{M}_f = 1 \times 10^{-5}$, $l=1$ diff stopped 1
173 167: $\bar{M}_f = 1 \times 10^{-7}$, $l=0.5$ diff stopped 1
174 173: $\bar{M}_f = 1 \times 10^{-6}$, $l=0.5$ diff stopped 1
175 174: $\bar{M}_f = 1 \times 10^{-5}$, $l=0.5$ diff stopped 1
176 175: $\bar{M}_f = 1 \times 10^{-7}$, sharp 0.95 and 0.05 diff stopped 1
177 176: $\bar{M}_f = 1 \times 10^{-6}$ diff stopped 1
178 177: $\bar{M}_f = 1 \times 10^{-5}$ diff stopped 1
179 176: $\bar{M}_i = 1 \times 10^{-2}$ diff stopped 1
180 179: $\bar{M}_i = 1 \times 10^{-4}$ diff stopped 1
181 179: $\bar{M}_i = 1 \times 10^{-5}$ diff stopped 1
182 179: $\bar{M}_i = 1 \times 10^{-6}$ diff stopped 1
184 138: 2D, $\bar{M}_i = 1 \times 10^{-7}$ diff stopped 12
185 184: 2D, $\bar{M}_i = 1 \times 10^{-5}$ diff stopped 12
186 184: 2D, $\bar{M}_i = 1 \times 10^{-3}$ diff stopped 4
187 184: 2D, $\bar{M}_i = 1 \times 10^{-1}$ diff stopped 12
188 186: 2D, $X_1 = 0.5 X_1^{equ}$ diff stopped 12
189 186: 2D, $X_1 = 0.1 X_1^{equ}$ diff calculations broke down stopped 4
190 186: 2D, $X_1 = 0.01 X_1^{equ}$ diff calculations broke down stopped 4
191 189: 2D, $\bar{M}_f = 1 \times 10^{-6}$ diff stopped 4
192 189: 2D, $\bar{M}_f = 1 \times 10^{-5}$ diff stopped 4
193 189: 2D, $\bar{M}_f = 1 \times 10^{-4}$ diff stopped 4
194 189: 2D, $\bar{M}_f = 1 \times 10^{-3}$ diff stopped 4
195 186: 2D, $X_1 = 0.3 X_1^{equ}$, $R_{tol} = 1 \times 10^{-2}$ diff stopped 12
196 195: 2D, $X_1 = 1.0 X_1^{equ}$ diff stopped 12
197 196: 2D, factor=2 diff stopped 32
198 185: 2D, $\mu = \mu_s^{\phi^4} \mu_l^{\left(1 - \phi^4\right)}$ diff stopped 12
199 185: 2D, $\mu = \mu_s^{\phi^6} \mu_l^{\left(1 - \phi^6\right)}$ diff stopped 12
200 185: 2D, $\mu = \mu_s^{\phi^{10}} \mu_l^{\left(1 - \phi^{10}\right)}$ diff stopped 12
201 185: 2D, $\mu = \mu_s^{\phi^8} \mu_l^{\left(1 - \phi^8\right)}$ diff stopped 12
202 199: 2D, $\mu = \mu_s^{\phi^2} \mu_l^{\left(1 - \phi^2\right)}$ diff stopped 4
203 199: 2D, $\mu = \mu_s^{\phi^3} \mu_l^{\left(1 - \phi^3\right)}$ diff stopped 4
204 198: 2D, $X_1 = 0.5 X_1^{equ}$ $\bar{M}_f = 1 \times 10^{-7}$ diff stopped 12
205 198: 2D, $X_1 = 0.1 X_1^{equ}$ $\bar{M}_f = 1 \times 10^{-7}$ $\Delta t < 1 \times 10^{-11} $ diff crashed benson, hobson, alice (stopped) 12
206 198: 2D, $X_1 = 1.0 X_1^{equ}$ $\bar{M}_f = 1 \times 10^{-6}$ diff crashed stopped 12
207 198: 2D, $X_1 = 0.5 X_1^{equ}$ $\bar{M}_f = 1 \times 10^{-6}$ diff crashed stopped 4
208 198: 2D, $X_1 = 0.1 X_1^{equ}$ $\bar{M}_f = 1 \times 10^{-6}$ diff stopped 4
209 198: 2D, $X_1 = 0.1 X_1^{equ}$ $\bar{M}_f = 1 \times 10^{-5}$ diff crashed cluster 411 (stopped) 4
210 198: 2D, $X_1 = 0.1 X_1^{equ}$ $\bar{M}_f = 1 \times 10^{-4}$ diff cluster 412 (stopped) 4
211 198: 2D, $X_1 = 0.1 X_1^{equ}$ $\bar{M}_f = 1 \times 10^{-3}$ diff cluster 413 (stopped) 4
212 198: 2D, $X_1 = 1.0 X_1^{equ}$ $\bar{M}_f = 1 \times 10^{-7}$ $\nu_f = 2 \times 10^{-2}$ $\Delta t < 1 \times 10^{-11} $ diff crashed renfield, sancho, kato (stopped) 12
213 212: 2D, $X_1 = 1.0 X_1^{equ}$ $\bar{M}_f = 1 \times 10^{-7}$ $\nu_f = 2 \times 10^{-1}$ 5000 iterations diff crashed cluster 418 (stopped) 4
214 205: 2D, $X_1 = 0.1 X_1^{equ}$ $\bar{M}_f = 1 \times 10^{-7}$ $\nu_f = 2 \times 10^{-3}$ $R_{\text{tol}} = 1 \times 10^{-2} $ diff crashed poole (stopped) 4
215 212: 2D, $X_1 = 1.0 X_1^{equ}$ $\bar{M}_f = 1 \times 10^{-7}$ $\nu_f = 2 \times 10^{-2}$ $R_{\text{tol}} = 1 \times 10^{-2} $ diff crashed luggage (stopped) 12
216 205: 2D, $X_1 = 0.1 X_1^{equ}$ $\bar{M}_f = 1 \times 10^{-7}$ $\nu_f = 2 \times 10^{-3}$ $R_{\text{tol}} = 1 \times 10^{-1} $ $\Delta_{\text{CFL}} = 5 \times 10^{-2} $ diff crashed benson, hobson, alice (stopped) 12
217 212: 2D, $X_1 = 1.0 X_1^{equ}$ $\bar{M}_f = 1 \times 10^{-7}$ $\nu_f = 2 \times 10^{-2}$ $R_{\text{tol}} = 1 \times 10^{-1} $ $\Delta_{\text{CFL}} = 5 \times 10^{-2} $ diff crashed renfield, sancho, kato (stopped) 12
218 198: $ \bar{M} = \bar{M}_f^{1 - \phi^m} \bar{M}_s^{\phi^m} $ diff crashed luggage (stopped) 12
219 218: $ m = 2 $ diff cluster 424 (stopped) 4
220 218: $ m = 3 $ diff cluster 425 (stopped) 4
221 218: $ m = 4 $ diff cluster 426 (stopped) 4

New base simulation

ID Notes diff Issues Status $N$ movie
222 base diff benson, hobson, alice 12
223 $X_1 = 0.5 X_1^{equ}$ diff None luggage (stopped) 12
224 $X_1 = 0.1 X_1^{equ}$ diff crashed luggage (stopped) 12
225 $\nu_f = 2 \times 10^{-2}$ diff luggage 12
226 $\nu_f = 2 \times 10^{-1}$ diff crashed renfield, sancho, kato (stopped) 12
227 224: $X_1 = 0.1 X_1^{equ}$, $R_{\text{solver}}=1 \times 10^{-9}$, $\sigma_{\text{CFL}} = 0.05$ diff crashed luggage (stopped) 12
228 225: $\nu_f = 2 \times 10^{-2}$, $R_{\text{solver}}=1 \times 10^{-9}$, $\sigma_{\text{CFL}} = 0.05$ diff crashed renfield, sancho, kato (stopped) 12
229 226: $\nu_f = 2 \times 10^{-1}$, $R_{\text{tol}}=1 \times 10^{-3}$, $\alpha_{\text{relax}} = 1.0$ diff crashed cluster (stopped) 4
230 222: $X_1 = 2.0 X_1^{equ}$ diff crashed cluster 12.0 (stopped) 4
231 222: $\nu_f = 1 \times 10^{-2}$ diff crashed oddjob (stopped) 4
232 222: $\nu_f = 2 \times 10^{-2}$ $\bar{M}_f = 1 \times 10^{-6}$ diff crashed oddjob (stopped) 4
233 222: $X_1 = 0.1 X_1^{equ}$ $\bar{M}_f = 1 \times 10^{-6}$ diff cluster 7.0 4
234 225: $\sigma_{\text{CFL}} = 0.02$ diff crashed oddjob (stopped) 3
235 225: $\nu_f = 4 \times 10^{-3}$ diff None cluster 2.0 (stopped) 4
236 225: $\alpha = 0.2$ diff crashed cluster 31.0 (stopped) 4
237 224: $\sigma_{\text{CFL}} = 0.02$ diff crashed renfield, sancho, kato (stopped) 12
238 224: $\alpha = 0.2$ diff crashed cluster 3.0 (stopped) 4
239 230: $X_1 = 1.5 X_1^{equ}$ diff crashed cluster 4.0 (stopped) 4
240 222: $X_1 = 1.0 X_1^{equ}$ $\bar{M}_f = 1 \times 10^{-6}$ diff crashed luggage (stopped) 12
241 223: $X_1 = 0.5 X_1^{equ}$ $\bar{M}_f = 1 \times 10^{-6}$ diff crashed cluster 5.0 (stopped) 4
242 235: $X_1 = 1.0 X_1^{equ}$ $\nu_f = 8 \times 10^{-3}$ diff crashed cluster 6.0 (stopped) 4
243 226: $\nu_f = 2 \times 10^{-1}$, shift=dx diff None renfield, sancho, kato (stopped) 12
244 226: $\nu_f = 2 \times 10^{-1}$, shift=10 * dx diff None luggage (stopped) 12
245 226: $\nu_f = 2 \times 10^{-1}$, shift=100 * dx diff None cluster 28 (stopped) 4
246 226: $\nu_f = 2 \times 10^{-1}$ diff None cluster 29 (stopped) 1
247 225: $\nu_f = 2 \times 10^{-2}$, shift=dx diff crashed luggage (stopped) 12
248 224: $X_1 = 0.1 X_1^{equ}$, shift=dx diff crashed cluster 30 (stopped) 4
249 247: $\nu_f = 2 \times 10^{-2}$, shift=2 * dx diff crashed renfield, sancho, kato (stopped) 12
250 247: $\nu_f = 2 \times 10^{-2}$, $\sigma_{\text{CFL}} = 0.05$, shift=dx diff crashed luggage (stopped) 12
251 240: $X_1 = 1.0 X_1^{equ}$, $\bar{M}_f = 1 \times 10^{-6}$, shift=dx diff crashed cluster 31 (stopped) 4
252 243: $\nu_f = 2 \times 10^{-1}$, shift=dx, $\sigma_{\text{CFL}} = 0.01$ diff cluster 32 4
253 243: $\nu_f = 2 \times 10^{-4}$, shift=dx diff crashed cluster 33 (stopped) 4
254 253: $\nu_f = 2 \times 10^{-4}$, shift=dx, $\nu_s = 2 \times 10^{5}$ diff crashed renfield, sancho, kato (stopped) 12
255 254: $\nu_f = 2 \times 10^{-4}$, shift=dx, $\nu_s = 2 \times 10^{4}$ diff crashed renfield, sancho, kato (stopped) 12
256 255: $\nu_f = 2 \times 10^{-4}$, shift=dx, $\nu_s = 2 \times 10^{7}$ diff crashed cluster 34 (stopped) 4
257 251: $X_1 = 1.0 X_1^{equ}$, $\bar{M}_f = 1 \times 10^{-6}$, shift=2*dx, $\sigma_{\text{CFL}}=0.05$ diff renfield, sancho, kato 12
258 256: $\nu_f = 2 \times 10^{-4}$, shift=dx, $\nu_s = 2 \times 10^{3}$ diff crashed cluster 35 (stopped) 4
259 248: $X_1 = 0.1 X_1^{equ}$, shift=2 * dx diff crashed cluster 36 (stopped) 4
260 250: $\nu_f = 2 \times 10^{-2}$, $\sigma_{\text{CFL}} = 0.05$, shift=2 * dx diff crashed renfield, sancho, kato (stopped) 12
261 259: $X_1 = 0.1 X_1^{equ}$, shift=10 * dx diff crashed luggage (stopped) 12
262 257: $X_1 = 1.0 X_1^{equ}$, $\bar{M}_f = 1 \times 10^{-6}$, shift=10*dx, $\sigma_{\text{CFL}}=0.05$ diff crashed bunter (stopped) 4
263 260: $\nu_f = 2 \times 10^{-2}$, $\sigma_{\text{CFL}} = 0.05$, shift=10 * dx diff crashed renfield, sancho, kato (stopped) 12
  • Posted: 2009-08-27 13:02 (Updated: 2009-12-15 15:24)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Fixing LAMMPI bug

Had a weird problem with the following error when using LAMMPI. 

-----------------------------------------------------------------------------
The selected RPI failed to initialize during MPI_INIT.  This is a
fatal error; I must abort.
This occurred on host poole (n0).
The PID of failed process was 29080 (MPI_COMM_WORLD rank: 0)
-----------------------------------------------------------------------------
-----------------------------------------------------------------------------
One of the processes started by mpirun has exited with a nonzero exit
code.  This typically indicates that the process finished in error.
If your process did not finish in error, be sure to include a "return
0" or "exit(0)" in your C code before exiting the application.
PID 29081 failed on node n0 (127.0.0.1) with exit status 1.
-----------------------------------------------------------------------------

The solution involved using the command "ipcs" and "ipcrm" and is explained here:

<http://www.riothamus.net/drupal/?q=node/12>

The above post has a perl script to run it use 

exec `./ipcrmAll.pl -s`
  • Posted: 2009-08-13 11:43 (Updated: 2009-08-13 11:55)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Parallel on the cluster

The table below shows sweep times for the reactive wetting simulation on the cluster (Simulation 122). The results are awful as was probably to be expected due to the bandwidth issues with the cluster.

Nproc sweep 1 sweep 2 sweep 3
1 17.43 13.65 13.72
2 9.65 7.34 7.40
4 19.23 19.76 7.26
8 23.11 20.72 21.44
16 29.33 100.40 126.25
  • Posted: 2009-07-28 16:49 (Updated: 2009-07-28 16:50)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Simulation 121

Simulation 121 is a fine grid version of the base case. Copied from 119 with the grid density doubled.

  • Posted: 2009-07-28 14:56 (Updated: 2009-07-28 14:57)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Simulations on Raritan

  • We copied simulation 119 to 120 to remove the HDF5 calls. At this stage we don't have HDF5 installed due to build issues with szip

Commands for raritan

This submits a job with 4 nodes and 1 processor per node

qsub -q infiniband -lnodes=4:ppn=1 -W x=NACCESSPOLICY:SINGLEJOB ./myscript

Simulations

  • 8-128 processors
  • 8-32 processor per node basis
  • dt from 1e-15 to something
  • Problem size 256x256, 512x512, 1024x1024

Simulation set 1, $\Delta t = 1 \times 10^{-15}$, 256x256

$N$: number of nodes (each node has X processors) $P$: number of processors being used on each node $t_i$ : wall clock time at sweep $i$

$N$ $P$ $t_{1}$ $t_{2}$ $t_{3}$ $t_{4}$
4 1

  • Posted: 2009-07-21 10:27 (Updated: 2009-07-21 10:51)
  • Author: wd15
  • Categories: (none)
  • Comments (22)

Equilibrium calculations

Currently running 1D equilibrium calculations to determine thermodynamics which has an equilibrium contact angle of $\Pi / 18$. This is done by varying $B$. Simulations 116, 117 and 118 are being used for this calculation. B_0 = 25418.6582832

$B$ $\gamma_{lv}$ $\gamma_{sl}$ $\gamma_{sv}$ $\theta$
$4B_0$ 18.98 31.74 44.16 $0.32 \Pi$
$8B_0$ 18.53 31.95 43.74 $0.28 \Pi$
$16B_0$ 18.50 32.15 45.45 $0.24 \Pi$
$50B_0$
  • Posted: 2009-07-20 17:58 (Updated: 2009-07-24 16:05)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Spatial convergence with a smaller drop

The radius of the drop in the following simulations is 4.3 $\times 10^{-7}$ m and $\Delta x=1 \times 10^{-8}$ m. The previous simulations had drop with radius 1.1 $\times 10^{-6}$ m and $\Delta x=1.1 \times 10^{-8}$ m.

ID Machine $N$ scheme
104 alice (4) stopped 150 central difference
105 renfield (4) 10589 300 central difference
106 luggage (32) 40986 600 central difference
107 benson (4) 15976 150 power law
108 poole (4) 15934 300 power law
109 600 power law
110 sancho (4) stopped 150 van Leer
111 alice (4) 11750 300 van Leer
112 benson (4) 19216 150 QUICK
  • Posted: 2009-05-12 15:58 (Updated: 2009-05-18 17:45)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Spatial Convergence

The Figure below shows triple junction position against time for thesimulations 93, 99, 100, 101, 102, 103. What is clear is that they don't all agree. In and off itself this is not a huge problem. As long as we can show that for some grid resolution their is convergence. If under resolved calculations spread in a similar manner and have the same trends we can get away with running parameter search-type simulations at under resolved conditions. Simulation 102 is running with a 900 x 900 grid. I think the drop size can be reduced for these convergence calculations. It is fairly clear that the weird spurious velocities at the triple junction go away when the van Leer or power law schemes are used.

source:trunk/reactiveWetting/93/radiusVtime.png

Central Difference N=225

source:trunk/reactiveWetting/93/velocity.png

Power Law N=225

source:trunk/reactiveWetting/99/velocity.png

Power Law N=450

source:trunk/reactiveWetting/100/velocity.png

van Leer N=225

source:trunk/reactiveWetting/101/velocity.png

  • Posted: 2009-05-12 10:39 (Updated: 2009-05-12 11:01)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Simulations using the Van Leer convection term

ID Machine $N$ movie
101 poole stopped (4) 225 movie
103 benson stopped (4) 450 movie
  • Posted: 2009-05-08 18:08 (Updated: 2009-05-12 15:59)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Spatial Convergence

ID Machine $N$ movie
99 poole stopped (4) 225 movie
100 luggage stopped (32) 450 movie
102 luggage stopped (32) 900 movie
  • Posted: 2009-05-08 17:40 (Updated: 2009-05-12 15:59)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Simulation using the power law convection term

ID Machine movie
99 poole stopped (4) movie
  • Posted: 2009-05-08 15:59 (Updated: 2009-05-12 11:13)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Meeting With Frank

Why has the reactive wetting problem taken so damn long?

It's really difficult. No one has modeled multicomponent van der Waals liquid-vapor (and we have a solid too).

  • source:trunk/reactiveWetting/paper/sandia1.png
  • source:trunk/reactiveWetting/paper/sandia2.png
  • flow with difficult source terms in the Navier-Stokes, standard codes don't deal with this
    • parasitic currents
    • pressure-correction algorithm needed to be invented
    • compressible due to nature of pressure density coupling (most work has perfect gas relationship for pressure based approach)
    • low Mach number compressible density coupled (rather than pressure) flow
    • multicomponent decoupling
  • didn't have parallel capability
  • thermo
  • optimal choice of primitive variables ($\rho$, $\rho_1$, $\rho_2$, $X_1$)
  • choice of diffusion coefficient
  • code complexity and difficult to organize, it has to solve a variety of different problems
  • did not have state of the art preconditioners and solvers
  • viscosity and diffusion ratios (interpolation schemes) (requires good preconditioners and solvers)
  • density trapping

Some of these issues are avoided with the three phase field approach, density is defined by the phase field variable

$\phi_s + \phi_l + \phi_v = 1$

but the van der Waals model ($P = P_0 \left( \rho \right) + \frac{\lambda}{2}\left( \nabla \rho \right)^2$) has fundamental advantages

  • a more fundamental physical description
  • handles the critical temperature
  • compression waves, which are interesting
  • no arbitrary phase field equations

What is the current situation?

What did we do to fix it?

  • full matrix method
  • Trilinos preconditioners and solvers
  • Exponential interpolation scheme

$\mu = \mu_s^{\phi} \mu_f^{1 - \phi}$

  • Full implementation of Trilinos as the "backend" for fipy, state of the art preconditioners, solvers and parallel capability
  • Using the full matrix method
    • allows us to leverage full capability of Trilinos preconditioners
  • Parallel version
  • reason to use FiPy? is to have easy ways to add new terms and components
  • Talking to Aaron Lott (preconditioner expert in maths and computer div.)

Will I write a paper soon?

Yes. Here is the current draft source:trunk/reactiveWetting/paper/paper.pdf

What did we do wrong?

  • in essence, this is a very deep and difficult mathematical model
  • ignore pressure for "real" physical model and write intermediate papers
  • two computer guys should be working together on each problem, never alone as we currently do
    • better code organization

Will all this work carry over into other problems?

Yes. Certainly the parallel, full matrix and Trilinos capabilities are immediately applicable to

  • Photo-Voltaics
  • Superconformal Ni Deposition

Other problems that we could possibly tackle include

  • Zeta-potential
  • VLS nano-wire growth
  • Vehicle Light-weighting
  • Ink Jet Process
  • mem-resistive technologies
  • VLS nano-wire growth
  • Crystal Plasticity with Lyle
  • Martensites

What other techniques do we still need to implement?

  • get all stuff back into fipy that I have developed for the Reactive Wetting code.
  • second order accurate unstructured capability
  • adapativity
  • Newton iterations
  • block preconditioners (working with Aaron Lott)
  • Posted: 2009-04-29 11:47 (Updated: 2009-04-30 09:20)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Simulations to Test Tolerance

Rerunning some of the viscosity simulations

ID Machine $|R_2|$ movie
91 luggage (14) stopped $1\times 10^{-1}$ movie
96 luggage (15) finished $1\times 10^{-2}$ movie
97 luggage finished $5\times 10^{-1}$ movie
98 benson stopped $1\times 10^{-3}$ movie
  • Posted: 2009-04-20 12:24 (Updated: 2009-05-12 11:15)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Testing Advection Schemes

As a start to getting a higher order advection scheme in FiPy I've been looking at what we already have. The following test is for the advection at 45 degrees to a Cartesian grid. The analytical solution is this on a 110 x 60 grid: source:trunk/hofv/analytical.png

The Van Leer scheme solotion on a 110 x 60 grid:

source:trunk/hofv/vanleer.png

The power law scheme on a 110 x 60 grid:

source:trunk/hofv/powerlaw.png

The actual convergence rates based on the L1 norm (on 5 grids with 11 * N x 6 * N cells with N=(1,2,4,8,16)):

source:trunk/hofv/convergence.png

I implemented the QUICK scheme without limiting, is essentially the Van Leer scheme without the limiting.

source:trunk/hofv/QUICK.png

  • Posted: 2009-04-16 15:33 (Updated: 2009-04-16 17:23)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Timings on Luggage

  • Luggage had one additional job running.

$\Delta t = 1 \times 10^{-11}$

900 x 900 $\Delta t = 1 \times 10^{-11}$

sweep 1 2 4 8 16 32 60 62
1 677.69 382.98 230.94 107.78 56.23 30.74 16.25 16.53
2 x 303.69 192.95 90.93 51.63 28.6 14.85 14.43

1800 x 1800 $\Delta t = 1 \times 10^{-11}$

sweep 16 32 62 63
1 231.93 141.68 93.18 87.18
2 187.85 111.25 69.97 63.80
  • Posted: 2009-04-13 15:26
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Full Parallel Simulations

Rerunning some of the viscosity simulations

ID Machine $\mu_s$ movie
91 luggage (14) stopped $2\times 10^3$ movie
92 finished $2\times 10^2$ movie
93 luggage (15) stopped $2\times 10^4$ movie
94 luggage (15) finished $2\times 10^1$ movie
95 luggage (15) finished $2\times 10^5$ movie
  • Posted: 2009-04-07 10:58 (Updated: 2009-05-08 17:30)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Chat with Jim --- Reactive Wetting Direction

Three main directions

  • compare with Cox theory

  • compare with Warren, Boettin, Roosen for spreading
  • compare with Pismen & Pomeau
  • Posted: 2009-04-06 11:14
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Full Parallel Implementation of the Reactive Wetting

  • The numbers are wall clock times in seconds, cpu clock time is hard to get with mpi jobs
  • The number in parentheses is the matrix solution time.
  • Luggage was heavily loaded when these calculations were done
  • Poole has two additional processes running.

900 x 900 $\Delta t = 1 \times 10^{-12}$

sweep 1 (poole) 2 (poole) 4 (poole) 8 (luggage) 16 (luggage) 32 (luggage) 64 (luggage)
1 x x x 189.89 (1.53) 90.41 (1.71) 46.48 (0.43) 34.40 (0.94)
2 x x x 210.94 (5.00) 75.16 (5.29) 43.60 (2.55) 34.48 (4.67)

450 x 450 $\Delta t = 1 \times 10^{-12}$

sweep 1 (poole) 2 (poole) 4 (poole) 8 (poole) 8 (luggage) 16 (luggage) 32 (luggage) 64 (luggage)
1 73.86 (0.64) 39.23 (0.332) 20.62 (0.189) 11.18 (0.102) 26.59 (0.220) 32.05 (0.264) 29.24 (1.08) 13.37 (0.24)
2 61.97 (3.13) 30.82 (1.198) 16.07 (0.742) 8.62 (0.461) 23.74 (0.535) 19.855 (0.640) 10.22 (0.596) 7.51 (0.795)

225 x 225 $\Delta t = 1 \times 10^{-12}$

sweep 1 (poole) 2 (poole) 4 (poole) 8 (poole) 8 (luggage) 16 (luggage) 32 (luggage) 64 (luggage)
1 19.16 (0.16) 10.16 (0.086) 5.57 (0.04) 3.44 (0.026) 6.01 (0.076) 13.82 (0.12) 9.66 (0.14) 16.87 (0.12)
2 14.51 (0.55) 7.56 (0.31) 3.92 (0.18) 2.03 (0.13) 5.01 (0.19) 8.36 (0.14) 5.91 (0.22) x

$\Delta t = 1 \times 10^{-11}$

900 x 900 $\Delta t = 1 \times 10^{-11}$

sweep 1 (poole) 2 (poole) 4 (poole) 8 (luggage) 16 (luggage) 32 (luggage) 64 (luggage)
1 x x x 136.75 (1.80) 130.58 (1.69) 47.51 (1.87) 36.69 (0.976)
2 x x x 253.05 (133.37) 194.08 (74.12) 142.07 (84.92) 72.17 (47.37)

450 x 450 $\Delta t = 1 \times 10^{-11}$

sweep 1 (poole) 2 (poole) 4 (poole) 8 (poole) 8 (luggage) 16 (luggage) 32 (luggage) 64 (luggage)
1 81.49 (1.07) 40.69 (0.51) 21.46 (0.30) 11.19 (0.16) 55.53 (1.82) 14.31 (0.39) 16.32 (0.24) 19.46 (0.52)
2 158.42 (97.25) 80.21 (49.89) 41.32 (25.99) 23.89 (15.94) 83.85 (20.27) 32.60 (20.04) 18.87 (11.13) 35.52 (28.13)

225 x 225 $\Delta t = 1 \times 10^{-11}$

sweep 1 (poole) 2 (poole) 4 (poole) 8 (poole
1 20.03 (0.27) 10.21 (0.12) 5.57 (0.075) 3.95 (0.041)
2 38.44 (23.78 ) 19.25 (11.97) 10.08 (6.39) 6.10 (4.16)
  • Posted: 2009-04-02 15:15 (Updated: 2009-04-03 15:44)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Parallel Efficiency Data

The following table shows run times and memory usage for the reactive wetting problem on a number of nodes. The time step is $\delta t = 1 \times 10^{-10}$. The machine is poole, which has 8 nodes and 32GB.

category no mpi mpi 1 mpi 2 mpi 4 mpi 8
sweep 1 (s) 85 85 74 53 39
sweep 2 (s) 186 180 156 105 81
sweep 3 (s) 289 284 243 161 136
memory high (GB) 0.5 1.1 3.4 6.4 10.2
memory stable (GB) 0.5 1.1 3.4 3.9 2.8

The memory usage results are very strange. Initially the memory spikes at a high value somewhere in the initialization phase (memory high) and then for the rest of the simulation (all subsequent sweeps) the memory stays at or around the stable value (memory stable). The stable value seems to be independent of the number of processors, but is higher when there is more than one processor. The memory usage is taken from top in the section with Mem:XXXX total, YYYY used. I'm using YYYY minus the initial YYYY to get the number.

  • Posted: 2009-03-17 13:16
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Numpy on the Altix

Thus far, I've done the following

I've added

[mkl] library_dirs = /opt/intel/Compiler/11.0/074/mkl/lib/64 lapack_libs = mkl_lapack mkl_libs = mkl_intel_lp64, mkl_intel_thread, mkl_lapack, mkl_core

to the site.cfg file in numpy and the built numpy with the following command

python setup.py config --compiler=intel build_clib --fcompiler=intel build_ext --compiler=intel --verbose build install --prefix=${USR}

Currently I'm getting an import error:

>>> import numpy

Traceback (most recent call last):

File "<stdin>", line 1, in ? File "/users/wd15/usr/ia64/2.6.16.60-0.23.PTF.403865.0-default//lib/python2.4/site-packages/numpy/init.py", line 143, in ?

import linalg

File "/users/wd15/usr/ia64/2.6.16.60-0.23.PTF.403865.0-default//lib/python2.4/site-packages/numpy/linalg/init.py", line 47, in ?

from linalg import *

File "/users/wd15/usr/ia64/2.6.16.60-0.23.PTF.403865.0-default//lib/python2.4/site-packages/numpy/linalg/linalg.py", line 29, in ?

from numpy.linalg import lapack_lite

ImportError?: /opt/intel/Compiler/11.0/081/mkl/lib/64/libmkl_intel_thread.so: undefined symbol: kmpc_end_critical

  • Posted: 2009-03-11 17:18 (Updated: 2009-03-11 17:19)
  • Author: wd15
  • Categories: (none)
  • Comments (3)

Compiling Trilinos on the Altix with the Intel compilers

This looks promising (from a trilinos mailing list entry):

../configure CC=icc CXX=icc F77=ifort CFLAGS="-O3 -fPIC" CXXFLAGS="-O3 -fPIC -LANG:std -LANG:ansi -DMPI_NO_CPPBIND" FFLAGS="-O3 -fPIC" --prefix=${USR} --with-install="/usr/bin/install -p" --with-blas="-L/opt/intel/Compiler/11.0/081/mkl/lib/64 -lmkl -lpthread" --with-lapack="-L/opt/intel/Compiler/11.0/081/mkl/lib/64 -lmkl_lapack -lpthreads" --enable-mpi --enable-amesos --enable-ifpack --enable-shared --enable-aztecoo --enable-epetra --enable-epetraext --enable-external --enable-ml --enable-threadpool --enable-thyra --enable-stratimikos --enable-triutils --enable-galeri --enable-zoltan --cache-file=config.cache

This is the "normal" thing we do on 64 bit:

../configure CXXFLAGS="-O3 -fPIC" CFLAGS="-O3 -fPIC" FFLAGS="-O5 -funroll-all-loops -fPIC" F77="g77" --enable-epetra --enable-aztecoo --enable-pytrilinos --enable-ml --enable-ifpack --enable-amesos --with-gnumake --enable-galeri --cache-file=config.cache --with-swig=$PATH_TO_SWIG_EXECUTABLE --prefix=$LOCAL_INSTALLATION_DIR

  • Posted: 2009-03-11 11:31 (Updated: 2009-03-18 15:09)
  • Author: wd15
  • Categories: (none)
  • Comments (12)

Preconditioning

Chat with Aaron. We need to do the following:

  • figure out max time step for trilinos jacobi, block jacobi, mL and block ML (where I change the A that is passed, only the diagonal blocks)
  • look at overall solution time as well
  • number of iterations
  • Posted: 2009-02-27 13:12
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Parallel Simulations

More simulations with one in parallel running on poole with 6 processors.

ID Machine $\mu_s$ $k$ $N$ $L$ (m) $\Delta t_{\text{max}}$ movie
84 cluster 1525 $2\times 10^2$ 1 225 $2.5\times 10^{-6}$ ? movie
85 cluster 1524 $2\times 10^2$ 1 450 $5\times 10^{-6}$ ? movie
86 cluster 1523 $2\times 10^1$ 1 225 $2.5\times 10^{-6}$ ? movie
87 cluster 1522 $2\times 10^3$ 1 225 $2.5\times 10^{-6}$ ? movie
88 cluster 1521 $2\times 10^{12}$ 1 225 $2.5\times 10^{-6}$ ? movie
89 poole 20791 $2\times 10^{12}$ 2 225 $2.5\times 10^{-6}$ ? movie

  • Posted: 2009-02-18 17:02 (Updated: 2009-03-18 15:55)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Working 2D simulations with whole matrix method (unsegregated solver)

Finally have working 2D simulations with the whole matrix method. Simulation source:trunk/reactiveWetting/65 and source:trunk/reactiveWetting/66 both have $\mu=2\times10^{-3}$, but $\bar{M}=1\times10^{-6}$ and $1\times10^{-7}$, respectively. Movies are in the source directories for velocity, density, concentration and phase. The sharp interface simple analytic solution is plotted in the images. It is questionable whether the results have reached a final equilibrium, there is a sift from the analytic curve, especially for simulation 65. The diffusion is an order of magnitude slower in this case. It is probable that this simulation has to run for longer in order to see convergence to the equilibrium result. The final equilibrium shapes for simulation 65 and 66 are in source:trunk/reactiveWetting/65/conc0005400.png and source:trunk/reactiveWetting/66/conc0005400.png, respectively.

Convergence

Figures source:trunk/reactiveWetting/65/velocity.png and source:trunk/reactiveWetting/66/velocity.png show the maximum velocity against time. Neither simulation has completely converged, it may be that the drops still need to adjust to equilibrium or the velocities could be parasitic in nature.

Next steps

  • Run a simulation with $\bar{M}=1\times10^{-4}$. This may converge a little better. (67)
  • Run a simulation with $\bar{M}=1\times10^{-5}$ but for a grid of 200 * 200 (68)
  • Profile for speed and adjust.
  • axi-symmetric drops.

  • Posted: 2009-01-16 17:17 (Updated: 2009-04-29 16:04)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

1D solid liquid vapor simulation

source:trunk/reactiveWetting/63/movie.gif

This is a 1D solid-liquid-vapor simulation with the unsegregated solver. The parameter values are $\bar{M}= 1 \times 10^{-4}$ and $\nu = 1 \times 10^{-3}$ everywhere.

  • Posted: 2008-12-16 10:52 (Updated: 2009-04-29 16:10)
  • Author: wd15
  • Categories: (none)
  • Comments (2)

Unsegregated Reactive Wetting

Equilibrium calculations

ID Changes Notes
47 with 48 Used the trilinos solver and newer simulations use superlu
52 with 47 Used the trilinos solver and newer simulations use superlu
53 with 47 Used the trilinos solver and newer simulations use superlu
60 with 47 using superlu

Convergence calculations

$\mu=2 \times 10^{-2}$

ID Changes Notes
48 with 47 $N=100$
49 with 48 $N=200$
50 with 48 $N=400$
51 with 48 $N=800$

Experiment with new solver

Trying to optimize the solver so it doesn't crash. Using superlu really helps.

ID Changes Notes
54 with 51 $N=100$

Convergence check with new method

ID Changes Notes
55 with 54 $N=100$
56 with 55 $N=200$
57 with 55 $N=400$
58 with 55 $N=800$
59 with 55 $N=50$

Cleaning up input files

ID Changes Notes
61 with 55 $N=100$

Unsegregated solver class and term multiplication

Simulation 61 represents a binary liquid-vapor using the Unsegregated solver class and the term multiplication branch. It seems to work.

  • Posted: 2008-10-06 15:50 (Updated: 2008-11-18 13:11)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Reactive Wetting Equations in Cylindrical Coordinates

$\partial_j = \left( \partial_r, \partial_z \right) $

Note

We need to remove the radial distance $r$ and embed it in the term's coefficient. For a general equation in cylindrical coordinates that are axisymmetric the general equation can be written,

$ \frac{\partial \phi}{\partial t} + \frac{1}{r} \partial_r \left( r u_r \phi \right) + \partial_z \left( u_z \phi \right) = \frac{1}{r} \partial_r \left( r \Gamma \partial_r \phi \right) + \partial_z \left( \Gamma \partial_z \phi \right)$

Multiply through by $r$ to get,

$ \frac{\partial \left( r \phi \right)}{\partial t} + \partial_r \left( r u_r \phi \right) + \partial_z \left( r u_z \phi \right) = \partial_r \left( r \Gamma \partial_r \phi \right) + \partial_z \left( r \Gamma \partial_z \phi \right)$

which is,

$ \frac{\partial \left( r \phi \right) }{\partial t} + \partial_j \left( r u_j \phi \right) = \partial_j \left( r \Gamma \partial_j \phi \right) $

Thus, in general, it is very easy to recast equations in cylindrical axi-symmetric coordinates with fipy. The fourth order term in the concentration equations may be an issue and the correction equation also has a fourth order term.

Fourth Order Term

$\frac{1}{r} \partial_r \left( r \Gamma_1 \partial_r \left[ \frac{1}{r} \partial_r \left( r \Gamma_2 \partial_r \phi \right) + \partial_z \left( \Gamma_1 \partial_z \phi \right) \right] \right) + \partial_z \left( \Gamma_1 \partial_z \left[ \frac{1}{r} \partial_r \left( r \Gamma_2 \partial_r \phi \right) + \partial_z \left( \Gamma_1 \partial_z \phi \right) \right] \right) $

This can be rewritten as

$\frac{1}{r} \partial_r \left( r \Gamma_1 \partial_r \left[ \frac{1}{r} \partial_r \left( r \Gamma_2 \partial_r \phi \right) + \frac{1}{r} \partial_z \left( r \Gamma_1 \partial_z \phi \right) \right] \right) + \frac{1}{r} \partial_z \left( r \Gamma_1 \partial_z \left[ \frac{1}{r} \partial_r \left( r \Gamma_2 \partial_r \phi \right) + \frac{1}{r} \partial_z \left( r \Gamma_1 \partial_z \phi \right) \right] \right) $

There is an awkward $\frac{1}{r}$ that is not next to a $\Gamma$ and thus can not just multiply a regular coefficient. The outer $\frac{1}{r}$ gets dealt with because we multiply all the terms in the equation by $r$. The above can be rewritten as,

$\frac{1}{r} \partial_j \left( \Gamma_1 \left[\partial_j - \frac{\hat{r}_j}{r} \right] \partial_k \left[ r \Gamma_2 \partial_k \phi \right] \right)$

Unfortunately, The last term needs to be added as a source term. If we had a third order convection term then it could be added implicitly.

Continuity

$ \frac{\partial \left(r \rho \right)}{\partial t} + \partial_j \left( r u_j \rho \right) = 0 $

Velocity

$ \frac{\partial \left(r \rho u_i \right)}{\partial t} + \partial_j \left( r u_j u_i \rho \right) = \partial_j \left( r \mu \left[ \partial_j u_i + \partial_i u_j \right] \right) - r \rho_k \partial_i \mu_k + \epsilon_k T r \rho_k \partial_i \left( \frac{1}{r} \partial_j \left(r \partial_j \rho_k \right) \right) $

Concentration Equation

$ \frac{\partial \left(r \rho X_2 \right)}{\partial t} + \partial_j \left( r u_j X_2 \rho \right) = \partial_j \left( \frac{ r \bar{M}}{T} X_1 X_2 \left[ \left(\frac{ \left( A_2 - A_1 \right) }{m} - \frac{ \left( e_2 - e_1 \right) }{m^2} \rho \right) \partial_j p + \left( 1 - p \right) \frac{ \left( e_2 - e_1 \right)}{m^2} \partial_j \rho \right] \right) $

$ + \partial_j \left( \frac{ r \bar{M} R }{m} \partial_j X_2 \right) - \partial_j \left( \bar{M} X_1 X_2 \left[ \epsilon_1 + \epsilon_2 \right] \left[\partial_j - \frac{\hat{r}_j}{r} \right] \partial_k \left[ r \rho \partial_k X_2 \right] \right) $

$+ \partial_j \left( r \bar{M} X_1 X_2 \partial_j \left[ \frac{1}{r} \partial_k \left( r \left[ \epsilon_1 X_1 - \epsilon_2 X_2 \right] \partial_k \rho \right) \right] \right)$

The last term can remain in its original form as it is a convection term.

Alternate Method

The formulation above could lead to some coding issues and maybe other difficulties concerning the fourth order source terms. A more expedient and relatively trivial method that is probably no less accurate is to define a CylindricalGrid2D object that has its mesh volumes and face areas modified accordingly. I just realized that the discretiztion is exactly equivalent in the finite volume method when using a cylindrically modified grid and rewriting the equations.

Reactive Wetting Paper

Outline

  1. Introduction
    1. Boiler plate - motivation
    2. Summary of physical parameters coinciding with typical applications (water-ice-salt, solder alloys)
    3. Thermodynamics
  2. Model
    1. Transport
    2. Thermodynamics
    3. Surface Tensions
  3. Parameters Choices - Numerics
    1. Governing equations summary.
    2. Numerical Approach
      1. Solution algorithm
      2. Parasitic Currents
    3. One dimensional liquid vapor analysis
      1. Dimensionless form of the pure-liquid vapor
      2. Convergence of a one-dimensional pure liquid-vapor system
        1. $\delta$,
        2. $Re$
        3. $M$,
        4. $Ca$
        5. $\kappa$ liquid
      3. Convergence of a one-dimensional binary liquid vapor system
        1. $Pe$, $Pe \rightarrow \infty$, What is the meaning?
        2. $Re > 1$ issues.
    4. Two dimensional Binary Solid-Liquid-Vapor
      1. Obtaining Plausible contact angles
      2. $\kappa$ solid
      3. Approximation of a solid with $\mu \rightarrow \infty$, density trapping.
      4. Comparison with analytical solution in equilibrium
  4. Results
    1. Extracting interface locations and junctions with phase field methods
      1. Ridge detection
      2. Wheeler method
      3. Villanueva method
    2. $\theta_{\text{contact}}$ versus $t$.
    3. Cox comparisons - look at later cox work with surfactants.
    4. $R\left(t\right)$ for different initial conditions.
    5. Pure Hydrodynamics with and without reactive wetting ($\theta_{\text{liquid}}$ vs. $\vec{v}$)
    6. Pretty pictures
    7. Compare with Gustav
      1. Steady state versus transient
      2. Compare streamlines
    8. $\theta$ versus diffusion
      1. $\bar{M}\ll 1$ versus $\bar{M}=0$
      2. $\bar{M}_{\phi}=0$
    9. Compare with S. Troian
    10. Spherical cap deviation
    11. Best measures of angle and shape
    12. Compare with Jim and Bill's right angle assumption

Copies of notes

Warren's notes

Wheeler's notes

Contact Angle Metrics

Given the simulations from here, we can predict the triple point motion and contact angles.

The following figure show the triple point position against time normalized by the initial radius of the drop.

source:/trunk/reactiveWetting/contact.png@latest

The triple point position is given by

$\vec{r}_T$

where

$P_T\left(\vec{r}_T \right) = \underset{\vec{r}}{\min}\,P_T\left(\vec{r}\right)$

and

$P_T\left( \vec{r}; \phi_T, X_{2T}, \rho_T \right) = 1 + \frac{ \left(\phi \left( \vec{r} \right) - \phi_T\right)^2 }{ \underset{\vec{r}}{\max}\,\left(\phi \left( \vec{r} \right) - \phi_T \right)^2} + \frac{ \left(X_2 \left( \vec{r} \right) - X_{2T}\right)^2 }{ \underset{\vec{r}}{\max}\,\left(X_2 \left( \vec{r} \right) - X_{2T} \right)^2} + \frac{ \left(\rho \left( \vec{r} \right) - \rho_T\right)^2 }{ \underset{\vec{r}}{\max}\,\left(\rho \left( \vec{r} \right) - \rho_T \right)^2}$

and

$\phi_T = \frac{\phi_l + \phi_v + \phi_s}{3}$

$X_{2T} = \frac{X_{2l} + X_{2v} + X2_{2s}}{3}$

$\rho_T = \frac{\rho_l + \rho_v + \rho_s}{3}$

The roughness of the curves in the figure is due to the zero order interpolation of $P_T$. Currently, the center of the cell with the minimum value is deemed to be the triple point. Higher order interpolation should improve the jagged nature of the curves. Need at least second order interpolation.

This movie shows the four points that are chosen to determine the triple point angles. The position of the four points seems reasonable.

The solid vapor interface point, $\vec{r}_{sv}$ such that,

$P_{sv}\left(\vec{r}_{sv} \right) = \underset{\vec{r}}{\min}\,P_{sv}\left(\vec{r}; \phi_{sv}, X_{2sv}, \rho_{sv} \right) C\left (\vec{r} ; \vec{r}_T \right)$

where

$C \left( \vec{r} ; \vec{r}_T \right) = 1 + 10 \left(\frac{|\vec{r} - \vec{r}_T| - R_0}{R_0}\right)^2$

with

$\phi_{sv} = \frac{\phi_v + \phi_s}{2}$

$X_{2sv} = \frac{X_{2v} + X2_{2s}}{2}$

$\rho_{sv} = \frac{\rho_v + \rho_s}{2}$

The radius of the circle $R_0$ is chosen to be $20 \Delta x$. We'll need to see how contact angle varies with this parameter.

Contact Angles

The links in the table below show contact angles versus time using the method described above.

ID $R$ $\mu_L$ Figure
2 $5\times10^{-7}$ $10^4$ figure
4 $1\times10^{-6}$ $10^4$ figure
6 $1\times10^{-6}$ $10^3$ figure
7 $5\times10^{-7}$ $10^2$ figure
8 $1\times10^{-6}$ $10^2$ figure
9 $5\times10^{-7}$ $10^3$ figure
10 $2\times10^{-6}$ $10^4$ figure
12 $2\times10^{-6}$ $10^3$ figure

Thoughts:

  • The initial contact angles in these figures seem to be somewhat wrong. The solid contact angle should be $\pi$ while the other angles should be $\pi / 2$.
  • Initial contact angle should improve with bigger droplets and this isn't happening.
  • The liquid and vapor angles both seem to drift away from the equilibrium values.
  • The curves are extremely jagged.
  • The surface tensions used to calculate the equilibrium contact angle could be for a different thermo. Need to double check.
  • Even though the system is not in equilibrium the contact angles should probably be in equilibrium.

What next?

  • Figure out why the initial angles are wrong.
  • Plot the equilibrium solution against that the last frame in the simulation to figure out how close to equilibrium.
  • Use higher order interpolation for the contact points to improve the smoothness.
  • Making $R_0$ smaller may improve the results.
  • Test the equilibrium contact angle in the case where the solid has the same viscosity as the liquid, so that the system reaches equilibrium. Probably need to run that case again.
  • run a full sphere on a solid surface that doesn't melt (correct equilibrium concentrations). Use this to compare with cox stuff.
  • vary $R_0$ for one frame of a simulation to see how the angle changes.

Jobs

Two new jobs

ID PID machine status $R$ ($\mu$m) movie Angles
25 13638.0 cluster stopped 0.4 movie
26 13639.0 cluster R 0.8 movie figure
27 13740.0 cluster R 0.35 movie figure
28 12539 poole R 1.8 movie figure
29 14230.0 cluster R 0.35 movie figure

Various plots for the simulations above,

Issues:

  • The $Ca$ number seems to be larger than Cox's theory predicts.
  • The drop off is slower than Cox's theory predicts.
  • Cox's theory is for high $Pe$ number. Need to run some simulations in this regime.
  • Need to improve data sampling. This is relatively easy to do by refining the grid and interpolating the base fields.
  • axi-symmetric may also be required

Peclet Number

A rough estimate for the Peclet Number $Pe$ is given by

$Pe = \frac{X m V_{\text{tp}} \rho}{\bar{M} R}$

If we take an interface length for $X$ of $~10^{-7}$. At early times the Peclet number is 100 in the solid and 0.01 in the liquid. This indicates that if we give the solid value of $\bar{M}$ to the fluid then the Peclet number will be large everywhere and will hopefully allow a better fit to Cox's theory. Simulation 29 in the table above is running with $\bar{M}_{f}=10^{-11}$

Jobs with lower Peclet number

ID PID machine status $R$ ($\mu$m) $\bar{M}_f$ $\mu_f$ $\Delta t_{\text{max}}$ movie Angles
29 14230.0 cluster unstable 0.35 $1 \times 10^{-11}$ $1 \times 10^{4}$ $1 \times 10^{-7}$ movie figure
30 14256.0 cluster finished 0.35 $1 \times 10^{-8}$ $1 \times 10^{4}$ $1 \times 10^{-4}$ movie figure
31 14257.0 cluster finished 0.35 $1 \times 10^{-9}$ $1 \times 10^{4}$ $1 \times 10^{-4}$ movie figure
32 14258.0 cluster unstable 0.35 $1 \times 10^{-10}$ $1 \times 10^{4}$ $1 \times 10^{-4}$ movie figure
33 14280.0 cluster R 0.35 0 $1 \times 10^{4}$ $1 \times 10^{-4}$ movie figure
34 14296.0 cluster R 0.35 0 $1 \times 10^{-3}$ $1 \times 10^{-4}$ movie figure

Higher order interpolation

It works!

High Peclet number

Can't seem to run the code with a Peclet number of about 10 and above. This is certainly connected with other issues. The figure below shows the Cox theory fit for decreasing $\bar{M}$ corresponding to increasing Peclet number. As the Peclet number is increasing, the curvature of the simulation points is increasing. This seems to show it is getting closer to the theory, though it is hard to tell. It is probable that as the Peclet number passes through 1 there will be a shift from one spreading regime to another.

Mbar=0

Far out. By solving for $\rho_1$ and $\rho_2$ separately and with corrections the system can be solved when $\bar{M}=0$. It's interesting that the system remains stable. The movie for this simulations is movie. I can't use this for Cox's theory as the system does not seem to movie towards the equilibrium contact angles. What are the equilibrium contact angles in the case when $\bar{M}=0$. It's also interesting that the pressure remains higher in the fluid state.

Low viscosity

The next step is try the $\bar{M}=0$ system with a much lower viscosity in the fluid. Running a job with $\mu_f=1 \times 10^{-3}$. Working on this issue made me think about modeling the system by solving for $\rho_1$ and $\rho_2$ instead. Simulation 35 is my attempt at this issue. I added RHS flux back in to each of these equations and did the variational derivative. Unfortunately this wasn't really successful at low viscosity. Same sort of issues as before.

  • Posted: 2008-06-06 16:54 (Updated: 2009-04-29 17:46)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Simulations with varying solver tolerance

ID status condor ID $R$ Tolerance Step Time (s) Steps Distance (R) Sim Time (s) $\left|\left|\frac{\phi - \phi^*}{\phi^*}\right|\right|_{\infty}$ $\left|\left|\frac{\rho - \rho^*}{\rho^*}\right|\right|_{\infty}$ $\left|\left|\frac{\rho_2 - \rho_2^*}{\rho_2^*}\right|\right|_{\infty}$ Bill's movie
18 stopped 12087.0 $5\times10^{-7}$ $1\times10^{-7}$ 120 460 0.13 0.00256 $7.0\times10^{-9}$ $1.4\times10^{-9}$ $1.5\times10^{-9}$ movie
19 stopped 12088.0 $5\times10^{-7}$ $1\times10^{-4}$ 96 460 0.13 0.00256 $1.2\times10^{-6}$ $9.9\times10^{-7}$ $9.7\times10^{-7}$ movie
20 stopped 12089.0 $5\times10^{-7}$ $1\times10^{-1}$ 66 460 0.13 0.00256 $6.1\times10^{-3}$ $5.9\times10^{-3}$ $5.9\times10^{-3}$ movie
21 stopped 12090.0 $5\times10^{-7}$ $1\times10^{-10}$ 156 460 0.13 0.00256 $0$ $0$ $0$ movie
22 running 12093.0 $5\times10^{-7}$ $1$ -- -- -- -- -- -- -- movie

Notes

  • These simulations need to run for longer to ascertain the long term effects of larger solver tolerances.
  • According to the Trilinos manuals, the default stopping criterion is $\frac{\left|\left|\vec{r}\right|\right|_2}{\left|\left|\vec{r}^{(0)}\right|\right|_2}$.
  • I'm still confused about how the tolerance relates to the error. Certainly I know that the error is linear with residual (Ae=r). This would say that an order of magnitude decrease in the tolerance leads to an order of magnitude decrease in error. However if the tolerance is 1 then no iterations should be done, which would be pathological. I am running a tolerance = 1 case, so we'll see.
  • Posted: 2008-04-17 22:42 (Updated: 2008-04-21 16:13)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Altix and !FiPy. What can we expect?

What are the expected performance benefits for FiPy using the Altix machine?

The following questions are pertinent.

  • Which examples?
  • What proportion of time is spent in routines that could benefit from a parallel capability?
  • How long do the jobs currently take on one processor?
  • How much memory do the jobs currently use?
job ID PID Status Radius Drop Size N Proportion Solver Max Memory (Gb) Steps Time Per Step (s) Time (days) Simulation time (s) Altix (k=1.6) (days) Poole (k=1.6)
14 stopped $5\times10^{-7}$ 34000 0.774 (@) 0.542 0.41 450 185 0.96 0.0016 0.21 0.28
15 stopped $1\times10^{-6}$ 136000 0.806 0.579 0.99 680 863 6.8 0.00226 1.5 2.0
16 poole 4430 running $2\times10^{-6}$ 544000 0.822 (+) 0.562 (+) 3.3 929 (%) 3336.0 35.9 (%) 0.00275 (%) 8.0 11
17 poole 4714 running $4\times10^{-6}$ 2176000 0.798 (&) 0.496 (&) 13 50667 (#) 11034.0 190 (647.1(#)) 0.0293 (#) 42 56
$6\times10^{-6}$ 4896000 28 502 112 148
$2.5\times10^{-5}$ $8.5\times10^{7}$ 483 15420 3440 4552
  • (+) still need many more data points
  • (%) reasonable extrapolation, waiting for more data points
  • (#) very bad extrapolation, waiting for more data points
  • (&) only at 200 steps, this number rises until 600 steps
  • (@) steps are too quick to get a good number

Notes

The simulations wall clock time requires the triple point to cover a distance of R / 10.

The memory usage seems to go like $\alpha + \beta R^2$ with $\alpha=0.217$ and $\beta=7.73\times10^{11}$.

The simulation times are estimated with a power law, $\alpha N^{\beta}$ with $\alpha=4.69e-6$ and $\beta=1.2$.

The following figure shows the time spent in routines that could benefit from a parallel capability .

A quick study of solver iterations and tolerance may also yield gains.

Other functions that are costly, but are not readily reimplemented in parallel.

Radius Drop Size N FillComplete? InsertGlobalValues? array
$5\times10^{-7}$ 34000 0.068 0.039 0.015
$1\times10^{-6}$ 136000 0.052 0.024 0.017
$2\times10^{-6}$ 544000 0.052 0.026 0.017
$4\times10^{-6}$ 2176000 0.057 0.027 0.020
  • Posted: 2008-04-14 20:39 (Updated: 2008-04-22 19:10)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Poor weave results on the altix

The times below are very bad. I have no idea why this is the case. Poole had great speed ups.

N Threads speedup speedup factor raw times (s)
500 0 1.00 1.00 0.50
500 1 0.33 0.00 1.50
500 2 0.65 0.65 0.76
500 4 1.29 1.13 0.38
500 8 2.49 1.36 0.20
500 16 4.43 1.45 0.11
500 32 6.08 1.43 0.08
500 64 6.35 1.36 0.08
1000 0 1.00 1.00 12.37
1000 1 0.62 0.00 20.06
1000 2 1.05 1.05 11.76
1000 4 1.88 1.37 6.57
1000 8 2.75 1.40 4.50
1000 16 2.91 1.31 4.25
1000 32 3.11 1.25 3.98
1000 64 3.38 1.23 3.66
2000 0 1.00 1.00 179.58
2000 1 0.74 0.00 244.13
2000 2 1.43 1.43 125.62
2000 4 2.91 1.71 61.63
2000 8 4.20 1.61 42.77
2000 16 4.97 1.49 36.16
2000 32 4.89 1.37 36.74
2000 64 5.87 1.34 30.60
3000 0 1.00 1.00 510.73
3000 1 0.71 0.00 719.87
3000 2 1.26 1.26 404.31
3000 4 2.30 1.52 221.79
3000 8 3.82 1.56 133.80
3000 16 3.98 1.41 128.48
3000 32 4.07 1.32 125.33
3000 64 4.80 1.30 106.35
  • Posted: 2008-04-11 16:56
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Posting code

Posting some code because it can't be shared using svn. 

1       import numpy
2 	from scipy import weave
3
4 	def multiply(N, nthreads):
5
6 	    i, j = numpy.indices((N, N), 'd')
7 	    c = numpy.zeros((N, N), 'd')
8
9 	    a = i + j
10 	    b = i * j
11
12 	    code_pragma = """
13
14 	    int i, j, k;
15
16 	    #pragma omp parallel shared(a,b,c,nthreads) private(i,j,k) num_threads(nthreads)
17 	      {
18
19 	      #pragma omp for schedule (dynamic)
20 	      for (i=0; i<N; i++)
21 	        for(j=0; j<N; j++)
22 	          for (k=0; k<N; k++)
23 	            c[i * N + j] += a[i * N + k] * b[k * N + j];
24
25 	      }
26
27 	    """
28
29 	    code_no_pragma = """
30
31 	    int i, j, k;
32
33 	    for (i=0; i<N; i++)
34 	      for(j=0; j<N; j++)
35 	          for (k=0; k<N; k++)
36 	            c[i * N + j] += a[i * N + k] * b[k * N + j];
37
38 	    """
39
40 	    args = {'N' : N, 'a' : a, 'b' : b, 'c' : c, 'nthreads' : nthreads}
41
42 	    if nthreads == 0:
43 	        code = code_no_pragma
44 	    else:
45 	        code = code_pragma
46
47 	    import time
48
49 	    t0 = time.time()
50
51 	    weave.inline(code,
52 	                 args.keys(),
53 	                 local_dict=args,
54 	                 type_converters=None,
55 	                 compiler = 'gcc',
56 	                 force=0,
57 	                 verbose=0,
58 	                 extra_compile_args =['-O3 -fopenmp'],
59 	                 extra_link_args=['-O3 -fopenmp'],
60 	                 headers=['<omp.h>'])
61
62 	    return time.time() - t0
63
64 	from fipy.tools import dump
65
66 	results = {}
67
68 	dump.write(results, 'datafile.gz')
69
70 	for N in (500, 1000, 2000, 3000, 4000, 5000, 6000):
71
72 	    results = dump.read('datafile.gz')
73 	    results[N] = {}
74 	    dump.write(results, 'datafile.gz')
75
76 	    for nthreads in range(9):
77
78 	        results = dump.read('datafile.gz')
79 	        results[N][nthreads] = []
80 	        dump.write(results, 'datafile.gz')
81
82 	        for sim in range(5):
83
84 	            results = dump.read('datafile.gz')
85 	            results[N][nthreads].append(multiply(N, nthreads))
86 	            dump.write(results, 'datafile.gz')
87
88
  • Posted: 2008-04-07 18:27 (Updated: 2008-04-07 18:30)
  • Author: wd15
  • Categories: (none)
  • Comments (0)

openmp and weave on poole

Some results for openmp on poole. The test is for an N by N matrix multiplication. The code is here. The speed up factor is certainly improving as the threads increase. The result for N=4000 are nonsense. Probably Ed running stuff. I should really have recorded resources available during each simulation by keeping a record with time.clock() and time.time().

N speedup factor for 8 threads
500 1.45
1000 1.78
2000 1.80
3000 1.82
4000 1.30
5000 1.76

More of the data.

500 0 1.00 1.00 0.67
500 1 0.52 0.00 1.30
500 2 0.97 0.97 0.69
500 3 1.41 1.24 0.48
500 4 1.82 1.35 0.37
500 5 2.22 1.41 0.30
500 6 2.65 1.46 0.25
500 7 2.86 1.45 0.23
500 8 3.03 1.45 0.22
1000 0 1.00 1.00 11.52
1000 1 0.85 0.00 13.55
1000 2 1.65 1.65 6.98
1000 3 2.35 1.72 4.89
1000 4 3.05 1.75 3.78
1000 5 3.84 1.79 3.00
1000 6 4.43 1.78 2.60
1000 7 5.12 1.79 2.25
1000 8 5.69 1.78 2.03
2000 0 1.00 1.00 113.21
2000 1 0.89 0.00 127.71
2000 2 1.69 1.69 67.07
2000 3 2.46 1.77 45.99
2000 4 3.21 1.79 35.27
2000 5 3.96 1.81 28.56
2000 6 4.64 1.81 24.38
2000 7 5.22 1.80 21.68
2000 8 5.81 1.80 19.48
3000 0 1.00 1.00 385.15
3000 1 0.89 0.00 431.11
3000 2 1.70 1.70 226.09
3000 3 2.49 1.78 154.88
3000 4 3.24 1.80 118.80
3000 5 4.02 1.82 95.81
3000 6 4.79 1.83 80.46
3000 7 5.43 1.83 70.93
3000 8 6.03 1.82 63.85
4000 0 1.00 1.00 2706.37
4000 1 0.78 0.00 3488.43
4000 2 1.31 1.31 2066.20
4000 3 1.77 1.44 1526.49
4000 4 2.07 1.44 1308.90
4000 5 2.29 1.43 1182.51
4000 6 2.28 1.38 1185.99
4000 7 2.16 1.32 1252.05
4000 8 2.21 1.30 1224.98
5000 0 1.00 1.00 1909.29
5000 1 0.78 0.00 2457.51
5000 2 1.61 1.61 1187.01
5000 3 2.45 1.76 780.55
5000 4 3.24 1.80 589.95
5000 5 3.96 1.81 481.70
5000 6 4.47 1.78 427.20
5000 7 5.01 1.78 380.78
5000 8 5.47 1.76 349.23

Machine details 

   $uname -a
   Linux poole 2.6.18-6-amd64 #1 SMP Sun Feb 10 17:50:19 UTC 2008 x86_64 GNU/Linux

8 of these 

processor       : 0
vendor_id       : AuthenticAMD
cpu family      : 15
model           : 65
model name      : Dual-Core AMD Opteron(tm) Processor 8220
stepping        : 3
cpu MHz         : 2812.978
cache size      : 1024 KB
physical id     : 0
siblings        : 2
core id         : 0
cpu cores       : 2
fpu             : yes
fpu_exception   : yes
cpuid level     : 1
wp              : yes
flags           : fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush mmx fxsr sse sse2 ht syscall nx mmxext fxsr_opt rdtscp lm 3dnowext 3dnow pni cx16 lahf_lm cmp_legacy svm cr8_legacy
bogomips        : 5630.17
TLB size        : 1024 4K pages
clflush size    : 64
cache_alignment : 64
address sizes   : 40 bits physical, 48 bits virtual
power management: ts fid vid ttp tm stc

and memory 

MemTotal:     32964928 kB
  • Posted: 2008-04-04 22:06
  • Author: wd15
  • Categories: (none)
  • Comments (0)

Curvature expressions in !FiPy

Curvature Expression FiPy Expression
iso $ \nabla \cdot \left( \frac{ \nabla \phi }{|\nabla \phi| } \right) $ (phi.getFaceGrad() / phi.getFaceGrad().getMag()).getDivergence()
mean $ \nabla \cdot \left( \frac{ \nabla \phi }{ \left(1 + |\nabla \phi|^2 \right)^{1/2} } \right) $ (phi.getFaceGrad() / numerix.sqrt(1 + phi.getFaceGrad().getMag()**2)).getDivergence()
Gauss $ \frac{ \phi_{xx} \phi_{yy} - \phi_{xy}^2 }{\left(1 + |\nabla \phi|^2 \right)^2} $ ((phi.getFaceGrad() * (1,0)).getDivergence() * (phi.getFaceGrad() * (0,1)).getDivergence() - phi.getFaceGrad()[::-1].getDivergence() / 2) / (1 + phi.getFaceGrad().getMag()**2)**2

Is there a vector form for the Gaussian curvature? It may help the FiPy expression.

Very good differential geometry book http://www.archive.org/details/differentialgeom003681mbp.

  • Posted: 2008-04-04 20:25
  • Author: wd15
  • Categories: (none)
  • Comments (0)

openmp on the altix

The following are some timing results on the Altix. The test case is a matrix multiplication. The time goes like N**3 where N is a side of the matrix. The numbers are averaged over 5 runs.

N Threads speedups scale factor raw times (s)
1000 1 1.00 1.00 9.99
1000 2 1.69 1.69 5.93
1000 4 2.86 1.69 3.23
1000 8 4.66 1.67 2.09
1000 16 8.05 1.68 1.17
1000 32 8.08 1.52 1.05
1000 64 9.53 1.46 0.95
2000 1 1.00 1.00 134.38
2000 2 1.79 1.79 75.21
2000 4 3.26 1.81 38.96
2000 8 5.85 1.80 22.89
2000 16 10.95 1.82 12.28
2000 32 13.35 1.68 8.03
2000 64 18.32 1.62 6.81
4000 1 1.00 1.00 1139.32
4000 2 1.77 1.77 640.28
4000 4 3.24 1.80 348.68
4000 8 5.90 1.81 190.78
4000 16 11.03 1.82 102.81
4000 32 16.39 1.75 63.89
4000 64 22.61 1.68 44.12
500 1 1.00 1.00 0.32
500 2 1.49 1.49 0.21
500 4 2.91 1.71 0.11
500 8 2.44 1.35 0.11
500 16 2.91 1.31 0.11
500 32 1.46 1.08 0.21
500 64 1.46 1.07 0.21

The important consideration is the scale factor and how that holds up as the threads increase. The C code used on the Altix is below. Still having problems with fipy as distutils is missing. As a point of reference I found this quote in an article, "More typical code will have a lower limit; 1.7x-1.8x are generally considered very good speedup numbers for code run on two threads" (http://cache-www.intel.com/cd/00/00/31/64/316421_316421.pdf). There is a start up time associated with threading. It is recommended that threads are maintained during the duration of a programming running. This may be difficult to achieve with weave.

Note:

The results need to be tested against unthreaded code. 

/******************************************************************************
* FILE: omp_mm.c
* DESCRIPTION:
*   OpenMp Example - Matrix Multiply - C Version
*   Demonstrates a matrix multiply using OpenMP. Threads share row iterations
*   according to a predefined chunk size.
* AUTHOR: Blaise Barney
* LAST REVISED: 06/28/05
******************************************************************************/
#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
int main (int argc, char *argv[])
{
  int   tid, nthreads, i, j, k, chunk;
int N=4000;
              //  int loop;
  //  int loops=1;
  double *a;
  a = malloc(N * N * sizeof(double *));
  double *b;
  b = malloc(N * N * sizeof(double *));
  double *c;
  c = malloc(N * N * sizeof(double *));
  printf("finished allocating\n");
  chunk = 10;                    /* set loop iteration chunk size */
  /*** Spawn a parallel region explicitly scoping all variables ***/
#pragma omp parallel shared(a,b,c,nthreads,chunk) private(tid,i,j,k)
  {
  tid = omp_get_thread_num();
  if (tid == 0)
    {
    nthreads = omp_get_num_threads();
    printf("Starting matrix multiple example with %d threads\n",nthreads);
    //printf("Initializing matrices...\n");
    }
  /*** Initialize matrices ***/
#pragma omp for schedule (runtime)
  //  #pragma omp for schedule (static, chunk)
  for (i=0; i<N; i++)
    for (j=0; j<N; j++)
      a[i * N + j]= i+j;
  //#pragma omp for schedule (static, chunk)
#pragma omp for schedule (runtime)
  for (i=0; i<N; i++)
    for (j=0; j<N; j++)
      b[i * N + j]= i*j;
#pragma omp for schedule (runtime)
  //#pragma omp for schedule (static, chunk)
  for (i=0; i<N; i++)
    for (j=0; j<N; j++)
      c[i * N + j]= 0;
  /*** Do matrix multiply sharing iterations on outer loop ***/
  /*** Display who does which iterations for demonstration purposes ***/
  //  printf("Thread %d starting matrix multiply...\n",tid);
#pragma omp for schedule (runtime)
  //  #pragma omp for schedule (static, chunk)
  //  for(loop=0; loop<loops; loop++)
    for(i=0; i<N; i++)
      for(j=0; j<N; j++)
        for (k=0; k<N; k++)
          c[i * N + j] += a[i * N + k] * b[k * N + j];
  }   /*** End of parallel region ***/
/*** Print results ***/
//
//printf("******************************************************\n");
//printf("Result Matrix:\n");
//for (i=0; i<NRA; i++)
//  {
//  for (j=0; j<NCB; j++)
//    printf("%6.2f   ", c[i][j]);
//  printf("\n");
//  }
//printf("******************************************************\n");
//printf ("Done.\n");
}

openmp/weave timings.

A matrix multiplication in weave really scales well with openmp. The code is here. The observed speedup is almost perfect with two threads.

This code, involving large array multiplications of size N, has the following speedups with two threads.

N Speedup
1E7 1.51
1E6 1.37
1E5 1.39
1E4 1.0

It should be noted that the number of loops in python increased inversely with the size of the array.

The question remains whether we can get speed ups for smaller arrays typically used in FiPy.

openmp and weave

The following steps are required to build openmp to work with weave. I've tried this on poole and rosie. I followed the follwoing steps

1) Install mpfr version 2.3.1 

        $ ../configure --prefix=${USR}
        $ make
        $ make install

2) Get gcc version 4.3 or 4.2.4 (when it is released). This version is needed otherwise you will get the "ImportError: libgomp.so.1: shared object cannot be dlopen()ed" error. gomp was not set up for dynamic loading in early gcc openmp compatible versions.

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=28482

3) Create a new directory and configure with 

    $ ../configure --prefix=${USR} --disable-multilib

where USR is in your local directories somewhere.

4) make; make install

The code to test weave and openmp.

A C code that actually works in parallel with openmp.

  • Posted: 2008-03-20 20:05 (Updated: 2008-03-20 21:20)
  • Author: wd15
  • Categories: fipy
  • Comments (0)

Streamlines

How to make a contour plot of the streamlines for a flow field in 2D?

Given a velocity field $\vec{u}$ we need a field $\psi$ that is constant in directions that are tangent to the velocity field. This leads to the equation,

$$ \vec{u} \cdot \nabla \psi = 0$$
The relationships
$$ u = \frac{ \partial \psi }{ \partial y} $$
and
$$v = -\frac{ \partial \psi }{ \partial x} $$
satisfy the equation. From the above equations we can say,
$$ \nabla^2 \psi = \frac{\partial u}{\partial y} - \frac{\partial v}{\partial x} $$
Solving this gives the stream-function. If we have no inflow conditions on regular domain walls, then we can choose $\psi=0$ on these walls. Since these walls have no circulation, we have made the choice that $\psi=0$ represents no circulation, which is sensible. The circulation orientation is given by the sign of $\psi$.

See an example of streamlines

FiPy code to compute streamlines

2D simulations integrated with subversion

ID status condor ID $\Delta t_{\text{max}}$ $R$ $\sigma_{\text{CFL}}$ $\bar{M}_S$ $\bar{M}_L$ $\mu_S$ $\mu_L$ Jim's movie Bill's movie
1 unstable $5\times10^{-7}$ $2\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ movie movie
2 running 12095.0 $1\times 10^{-4}$ $5\times10^{-7}$ $5\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ movie movie
3 stopped $1\times10^{-6}$ $5\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ movie movie
4 running 12096.0 $1\times10^{-6}$ $5\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ movie movie
5 unstable $5\times10^{-7}$ $2\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^3$ movie movie
6 running 12097.0 $1\times10^{-6}$ $2\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^3$ movie movie
7 stopped $1\times 10^{-4}$ $5\times10^{-7}$ $2\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^2$ movie movie
8 running 12100.0 $1\times 10^{-4}$ $1\times10^{-6}$ $2\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^2$ movie movie
9 unstable $1\times 10^{-4}$ $5\times10^{-7}$ $2\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^3$ movie movie
10 running poole 5545 $1\times 10^{-4}$ $2\times10^{-6}$ $5\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ movie movie
11 running poole 9411 $1\times 10^{-4}$ $4\times10^{-6}$ $5\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ movie movie
12 stopped $1\times 10^{-4}$ $2\times10^{-6}$ $2\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^3$ movie movie
13 stopped $1\times 10^{-4}$ $4\times10^{-6}$ $5\times10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^3$ movie movie

2D simulations

Table of simulations

ID status condor ID Initial T $\Delta t_{\text{max}}$ $R$ $\sigma_{\text{CFL}}$ $\bar{M}_S$ $\bar{M}_L$ $\mu_S$ $\mu_L$ data directory movie link source files
H2D-restart H2D 550 $1\times 10^{-7}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ movie
I2D-restart I2D 550 $1\times 10^{-8}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ movie
J2D-restart running 9099.0 J2D 600 $1\times 10^{-7}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataFri-Feb-29-12:25:35-2008 movie files
K2D-restart running 9098.0 K2D 600 $1\times 10^{-6}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataFri-Feb-29-12:25:35-2008 movie files
L2D unstable 9077.0 650 $5 \times 10^{-7}$ $10^{-1}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataWed-Feb-27-16:59:24-2008 movie files
L2D1 unstable 9087.0 650 $5 \times 10^{-7}$ $5 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataThu-Feb-28-11:27:18-2008 movie files
M2D running 9088.0 650 $1 \times 10^{-6}$ $5 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataThu-Feb-28-11:28:02-2008 movie files
N2D swapping 9089.0 650 $2 \times 10^{-6}$ $5 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataThu-Feb-28-11:31:03-2008 movie files
O2D unstable 9090.0 650 $5 \times 10^{-7}$ $5 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^3$ dataFri-Feb-29-01:27:24-2008 movie files
O2D1 running 9108.0 650 $5 \times 10^{-7}$ $2 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^3$ dataFri-Feb-29-14:35:04-2008 movie files
P2D running 9111.0 650 $1 \times 10^{-6}$ $2 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^3$ dataSat-Mar--1-07:01:06-2008 movie files
Q2D swapping 9112.0 650 $2 \times 10^{-6}$ $2 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^3$ dataSun-Mar--2-22:30:18-2008 movie files
R2D unstable 9083.0 650 $5 \times 10^{-7}$ $5 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^2$ dataThu-Feb-28-07:43:56-2008 movie files
R2D1 running 9093.0 650 $5 \times 10^{-7}$ $2 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^2$ dataThu-Feb-28-11:37:49-2008 movie files
S2D running 9094.0 650 $1 \times 10^{-6}$ $2 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^2$ dataThu-Feb-28-11:38:37-2008 movie files
T2D swapping 9095.0 650 $2 \times 10^{-6}$ $2 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^2$ dataThu-Feb-28-11:41:06-2008 movie files

Upper bound on CFL condition

The simulations below do not seem to have a bounded $\Delta t$. They do, however, have an upper bound for $\sigma_{\text{CFL}}$ that decreases dramatically as $\mu_L$ is reduced. Of course, at lower $\mu_L$ the upper bound for $\Delta t$ is probably due to some other numerical factor unrelated to $\sigma_{\text{CFL}}$.

Matrix of simulations

ID status $\Delta t_{\text{max}}$ $L$ $f_V$ $f_L$ $f_S$ $\sigma_{\text{CFL}}$ $\bar{M}_S$ $\bar{M}_L$ $\mu_S$ $\mu_L$ data directory movie link
N finished None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-1}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataMon-Feb-25-14:55:59-2008 movie
N3 unstable None $5 \times 10^{-6}$ 1/4 1/2 1/4 $2\times 10^{-1}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataWed-Feb-27-12:27:04-2008 movie
O finished None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-1}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^3$ dataMon-Feb-25-17:58:15-2008 movie
P unstable None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-1}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^2$ dataMon-Feb-25-18:30:26-2008 movie
Q finished None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^2$ dataMon-Feb-25-18:57:53-2008 movie
Q1 finished None $5 \times 10^{-6}$ 1/4 1/2 1/4 $5 \times 10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^2$ dataWed-Feb-27-14:51:24-2008 movie
R unstable None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^1$ dataTue-Feb-26-10:13:19-2008 movie
S unstable None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-3}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^1$ dataTue-Feb-26-11:45:36-2008 movie

2D simulations with reduced temperature

Matrix of 2D simulations with reduced temperature

ID status condor ID Initial T $\Delta t_{\text{max}}$ $R$ $\sigma_{\text{CFL}}$ $\bar{M}_s$ $\bar{M}_l$ $\mu_s$ $\mu_l$ data directory movie link
A2D finished 6952.0 650 $1\times 10^{-7}$ $3.33 \times10^{-7}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataWed-Jan-16-00:48:07-2008 movie
B2D unstable 7672.0 A2D 550 $1\times 10^{-7}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataFri-Jan-25-18:17:34-2008 movie
C2D finished 7828.0 650 $1\times 10^{-7}$ $6.66 \times10^{-7}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataWed-Jan-30-11:53:40-2008 movie
D2D finished 7836.0 650 $1\times 10^{-7}$ $1.33 \times10^{-6}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataWed-Jan-30-12:19:30-2008 movie
E2D finished 7848.0 650 $1\times 10^{-7}$ $1 \times10^{-6}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataWed-Jan-30-13:52:02-2008 movie
F2D finished 7859.0 650 $1\times 10^{-7}$ $5 \times10^{-7}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataWed-Jan-30-14:09:54-2008 movie
G2D finished 7869.0 650 $1\times 10^{-7}$ $2 \times10^{-6}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataWed-Jan-30-14:15:37-2008 movie
H2D finished 8541.0 F2D 550 $1\times 10^{-7}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataThu-Feb-14-12:09:33-2008 movie
I2D finished 8582.0 F2D 550 $1\times 10^{-8}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataFri-Feb-15-14:55:08-2008 movie
J2D finished 8803.0 F2D 600 $1\times 10^{-7}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataTue-Feb-19-10:20:01-2008 movie
K2D finished 8805.0 F2D 600 $1\times 10^{-6}$ $10^{-1}$ $10^{-9}$ $10^{-9}$ $10^{4}$ $10^4$ dataTue-Feb-19-10:20:56-2008 movie

Most of these jobs ran for a good period of time, but were brought down with server problems on 02/14/08 and 02/22/08. When they restarted the condor system began them from scratch, I killed them and will restart any useful ones from the latest point in the data directory above. However, it is probably best to keep processors for the new set of 2D jobs as these will probably have unbounded time steps.

Efficiency

It is interesting to compare time versus time steps for some of the 2D jobs and the new 1D jobs. $N$ represents number of time steps.

Simulations E2D, F2D and G2D all lie on the same curve as expected. Elapsed time in simulations N and O is orders of magnitude greater than in the 2D simulations for an equivalent number of time steps. This is not the case in simulation O as it has a more restrictive CFL condition (see Table).

1D simulation with vapor

Previously, we've managed to get the 1D liquid-solid system working with all the desired numbers (apart from the liquid viscosity) with only a CFL time step stability criterion. We now need to include the vapor and see how the time step is affected. I'm assuming from previous experience and intuition that adding the vapor will require a far stricter stability limit.

Simulation With Vapor

As per usual my intuition proves to be false. The time step size is completely unrestricted. I'm extremely excited. We may be able to do the 2D simulations in orders of magnitude less real time. There are also some extra simulations to look at initial liquid fraction and domain size. It is clear from simulation N that diffusion based interface motion is occurring at these parameter values.

Diffusion Driven Interface Motion

Examining simulations N, O, Q, it is evident that interface motion is driven by diffusion. Notice also that the pressure equilibrates sooner as the viscosity is reduced and it seems that the concentration gradients are more pronounced. In simulation Q, the pressure is equilibrated much earlier, before significant interface motion occurs and concentration gradients form.

Stability

To maintain stability for lower $\mu_l$, $\sigma_{\text{CFL}}$ is being reduced. Simulations P, R, S were all unstable. The lowest value of $\mu_l$ that was successfully computed is $10^2$, simulation Q. The following figure shows simulation time against number of required steps ($N$). It is clear that as $\sigma_{\text{CFL}}$ the run times are significantly increased.

Interface Motion

Interface motion plot for simulation N, O and Q. At least between $\mu_l=10^4$ and $\mu_l=10^2$ the solid liquid interface motion is mostly unaffected by $\mu_l$.

Matrix of simulations

ID status $\Delta t_{\text{max}}$ $L$ $f_V$ $f_L$ $f_S$ $\sigma_{\text{CFL}}$ $\bar{M}_s$ $\bar{M}_l$ $\mu_s$ $\mu_l$ data directory movie link
H finished None $10^{-6}$ 0 1/2 1/2 $10^{-1}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataThu-Feb-21-17:46:23-2008 movie
L finished None $10^{-6}$ 1/3 1/3 1/3 $10^{-1}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataMon-Feb-25-11:44:28-2008 movie
M finished None $10^{-6}$ 1/4 1/2 1/4 $10^{-1}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataMon-Feb-25-12:08:13-2008 movie
N finished None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-1}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataMon-Feb-25-14:55:59-2008 movie
O finished None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-1}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^3$ dataMon-Feb-25-17:58:15-2008 movie
P unstable None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-1}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^2$ dataMon-Feb-25-18:30:26-2008 movie
Q running None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^2$ dataMon-Feb-25-18:57:53-2008 movie
R unstable None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-2}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^1$ dataTue-Feb-26-10:13:19-2008 movie
S unstable None $5 \times 10^{-6}$ 1/4 1/2 1/4 $10^{-3}$ $10^{-11}$ $10^{-7}$ $10^{12}$ $10^1$ dataTue-Feb-26-11:45:36-2008 movie

Realistic Solid Behavior

The following plot includes the curves for each of the 1D simulations (see Table).

It is clear that $\bar{M}_s$ has more influence on the interface motion than $\mu_s$ during melting.

The next step is to try simulation H with the vapor added in 1D to figure out a satisfactory time step.

Also, it would be interesting to know how $\mu_l$ influences the interface motion.

Finite Size Effects

Having examined simulations I and K and compared them with B and D (see Table), it seems that there may be finite size effects in these simulations. The interface motions are way off. For example simulation B and I use the same parameter set yet the interface motion is over two orders of magnitude slower for case I. Some simulations need to be performed to ascertain whether there is an eventual convergence of the interface motion as the domain is expanded.

Table of Simulations

ID status $\Delta t_{\text{max}}$ $L$ $\sigma_{\text{CFL}}$ $\bar{M}_s$ $\bar{M}_l$ $\mu_s$ $\mu_l$ data directory movie link
B finished None $10^{-6}$ 0.1 $10^{-7}$ $10^{-7}$ $10^4$ $10^4$ dataWed-Feb-20-14:09:53-2008 movie
I finished None $10^{-5}$ 0.1 $10^{-7}$ $10^{-7}$ $10^{4}$ $10^4$ dataThu-Feb-21-18:08:24-2008 movie
J finished None $2\times10^{-6}$ 0.1 $10^{-7}$ $10^{-7}$ $10^{4}$ $10^4$ dataFri-Feb-22-15:14:20-2008 movie

Chat with Jim

After a chat with Jim we both came to the realization that curves on a liquid fraction against time plot for varying domain sizes do not have to correspond at all. That raises the question about how to determine whether we have finite size effects. There does not seem to be an easy way to do this. Certainly plotting interface position as a function of time may be more enlightening. The velocity of the interface at early times should be independent, but at long times I don't think that is true.

Absolute Interface Position

It is clear from the figure that all our fears are unfounded and this whole episode was a waste of time. Anyhow, it is interesting to notice that there is more "ringing" in the larger domain.

High solid viscosity

Jim wants to understand why we cannot model a rigid solid interface with a high viscosity. He needs the story for his talk. Suggestion is to run a simulation with a very high viscosity in the solid and a liquid that is not at equilibrium to induce melting of the solid

Initial Condition

The initial condition is a phase fraction of 1/2 liquid and 1/2 solid. The solid is at the three phase equilibrium values, $\rho_1^s$ and $\rho_2^s$. The initial densities in the liquid region are given by,

$ \rho_1 = \rho_1^l + \lambda \left( \rho_1^s - \rho_1^l \right) $

$ \rho_2 = \rho_2^l + \lambda \left( \rho_2^s - \rho_2^l \right) $

In this case $\lambda=-0.1$. This induces melting of the solid. For the numbers that we use in the example the final liquid fraction is given by 0.55. Using this approach there is a lower bound for $\lambda$ and thus an upper bound on the final liquid fraction.

Regular solid viscosity

When $\mu_s=1000$ the solid behaves much like a liquid. The movie is given here for $\Deltat=1\times10^{-6}$ (dir: dataWed-Feb-20-10:37:23-2008). In this simulation there is no clear distinction between interface motion due to pressure equilibration and motion due to the concentration gradient. In simulations where the temperature is reduced there is a clear distinction between these regimes of interface motion. One issue is that the total interface motion is equivalent to the interface width. It is definitely a good idea to run simulations with a small interface width compared with the total interface motion.

Interestingly enough, we can run this simulation with no upper bound on time step other than CFL. The parasitic velocities eventually give an upper bound less than 1. Here is the movie (dir: dataWed-Feb-20-14:09:53-2008)

High solid viscosity

When $\mu_s=1\times10^{12}$ the simulation shows similar interface motion. Here is the movie (dir: dataWed-Feb-20-14:36:18-2008). A simulation has also been run with a bounded $\Delta t$ (dir: dataWed-Feb-20-10:47:36-2008)

Liquid Fraction plot

The following plot is the liquid fraction against time for the whole matrix of simulations. The simulation with $\Delta t < \Delta t_{\text{CFL}}$, $\bar{M}_s=10^{12}$ and $\Delta t<\Delta t_{\text{CFL}}$ shows rather strange long time behavior. $\Delta t\approx5$ at long times, which could be too large.

Low diffusion in the solid

To get closer to real solid behavior, the value of $\bar{M}_s$ can be reduced. In this movie, $\bar{M}_s=10^{-15}$ (dir: dataWed-Feb-20-18:22:20-2008).

Matrix of simulations

ID status $\Delta t_{\text{max}}$ $L$ $\sigma_{\text{CFL}}$ $\bar{M}_s$ $\bar{M}_l$ $\mu_s$ $\mu_l$ data directory movie link
A finished $10^{-6}$ $10^{-6}$ 0.1 $10^{-7}$ $10^{-7}$ $10^4$ $10^4$ dataWed-Feb-20-10:37:23-2008 movie
B finished None $10^{-6}$ 0.1 $10^{-7}$ $10^{-7}$ $10^4$ $10^4$ dataWed-Feb-20-14:09:53-2008 movie
C finished $10^{-6}$ $10^{-6}$ 0.1 $10^{-7}$ $10^{-7}$ $10^{12}$ $10^4$ dataWed-Feb-20-10:47:36-2008 movie
D finished None $10^{-6}$ 0.1 $10^{-7}$ $10^{-7}$ $10^{12}$ $10^4$ dataWed-Feb-20-14:36:18-2008 movie
E finished $10^{-6}$ $10^{-6}$ 0.1 $10^{-15}$ $10^{-7}$ $10^{12}$ $10^4$ dataWed-Feb-20-18:22:20-2008 movie
F finished None $10^{-6}$ 0.1 $10^{-15}$ $10^{-7}$ $10^{12}$ $10^4$ dataThu-Feb-21-17:16:51-2008 long short
G finished $10^{-6}$ $10^{-6}$ 0.1 $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataThu-Feb-21-17:40:48-2008 movie
H finished None $10^{-6}$ 0.1 $10^{-11}$ $10^{-7}$ $10^{12}$ $10^4$ dataThu-Feb-21-17:46:23-2008 movie
I finished None $10^{-5}$ 0.1 $10^{-7}$ $10^{-7}$ $10^{4}$ $10^4$ dataThu-Feb-21-18:08:24-2008 movie
K finished None $10^{-5}$ 0.1 $10^{-7}$ $10^{-7}$ $10^{12}$ $10^4$ dataThu-Feb-21-18:11:13-2008 movie

2008/02/20/16.43

Test