p211
C-program for spinodal decomposition and block-copolymer phase separation CDS model (2D).
The CDS model we have constructed up to this point can be demonstrated easily using the 100 x100 lattice on the laptop.
The following program is annotated (almost surely all the comments are commented out). In this program f(x) is not 1.3 tanh x but is simplified as Ax. That is, exploiting the universality, the map is replaced with a piecewise linear computationally simpler map: f(x) = Ax for (-1/A, 1/A), outside this interval, if x < 0 f(x) = -1, else if x > 0, f(x) = +1. The value of A is requested by the program at the beginning of each run. Choose something like 1.2. If too large, the patterns freeze.
/* spinodal decomposition and block co-polymer system */
#include <stdio.h>
#define L 50 /* lattice size */
#define D 0.35 /* diffusion constant */
#define R 0.1 /* initial amplitude */
#define T 3000 /* time limit number of iterations */
#define DT 50 /* time increment displaying the patterns at every DT steps*/
#define frandom() (float)(((float)random()/(float)( 1 << 30 ))/ 2.0 )
#define max 1. /* display range */
/* A : quench depth the farther away from 1, the deeper the quench
B : block copolymer parameter; this is for the block copplymer cae．B > 0. Smaller B implies longer chains.*/
char *ctab[] = { "e" , " " , "." , "+" , "*" , "#" , "$" , "E" }; /* This is a cheap display converting the field to `gray scale using letters; crowded letters denote larger \phi. */
main()
{
float x[L+ 2 ][L+ 2 ], xt[L+ 2 ][L+ 2 ], fx[L+ 2 ][L+ 2 ];
float fo1, fo2, Map, Diff, Force, A, B, sum;
double atof();
int i, j, k;
char buff[ 128 ];
int seed, atoi();
unsigned int c;
printf( "seed :" );
gets( buff );
seed = atoi(buff);
printf( "A :" ); /* The slope of th map around the origin. A= 1.1 or 1.2 is convenient* /
gets( buff);
A = atof ( buff );
printf( "B :" );/* corresponding to the polymer length^2, something like B = 0.01 is convenient. */
gets( buff);
B = atof (buff );
srandom( seed );
/* initial condition */
for (i= 1 ;i<=L;i++)
{
for (j= 1 ;j<=L;j++)
{
x[i][j] = R*(frandom()- 0.5 );
}
}
/* periodic boundary condition imposed */
for (i= 1 ;i<=L;i++)
{
x[i][L+ 1 ] = x[i][ 1 ];
x[i][ 0 ] = x[i][L];
}
for (i= 0 ;i<=L+ 1 ;i++)
{
x[L+ 1 ][i] = x[ 1 ][i];
x[ 0 ][i] = x[L][i];
}
/* updating */
(k= 0 ;k<=T;k++)
{
/*computation of force: I is computed*/
for (i= 1 ;i<=L;i++)
{
for (j= 1 ;j<=L;j++)
{
Diff = x[i+ 1 ][j] + x[i- 1 ][j] + x[i][j+ 1 ] + x[i][j- 1 ];
Diff = 2.0 *Diff + x[i+ 1 ][j+ 1 ] +
x[i+ 1 ][j- 1 ] + x[i- 1 ][j+ 1 ] + x[i- 1 ][j- 1 ];
Diff = D*(Diff/ 12.0 - x[i][j]);
fo1 = A*x[i][j];
fo2 = (fo1< max)? fo1:max;
Map = (fo2>-max)? fo2:-max;
fx[i][j] = Map + Diff - x[i][j]; /* This is I. */
}
}
for
/* periodic boundary condition for the force */
for (i= 1 ;i<=L;i++)
{
fx[i][L+ 1 ] = fx[i][ 1 ];
fx[i][ 0 ] = fx[i][L];
}
for (j= 0 ;j<=L+ 1 ;j++)
{
fx[L+ 1 ][j] = fx[ 1 ][j];
fx[ 0 ][j] = fx[L][j];
}
/* computation of increment */
for (i= 1 ;i<=L;i++)
{
for (j= 1 ;j<=L;j++)
{
/* Force is negative actual increment */
Force = fx[i+ 1 ][j] + fx[i- 1 ][j] + fx[i][j+ 1 ] + fx[i][j- 1 ];
Force = 2.0 *Force + fx[i+ 1 ][j+ 1 ] +
fx[i+ 1 ][j- 1 ] +fx[i- 1 ][j+ 1 ] + fx[i- 1 ][j- 1 ];
Force = Force/ 12.0 - fx[i][j];
xt[i][j] = (1.- B)*x[i][j] - Force;
}
}
/* periodic boundary condition imposed */
for (i= 1 ;i<=L;i++)
{
xt[i][L+ 1 ] = xt[i][ 1 ];
xt[i][ 0 ] = xt[i][L];
}
for (i= 0 ;i<=L+ 1 ;i++)
{
xt[L+ 1 ][i] = xt[ 1 ][i];
xt[ 0 ][i] = xt[L][i];
}
sum = 0. ;
for (i= 1 ;i<=L;i++)
{
for (j= 1 ;j<=L;j++)
{
sum = sum + xt[i][j];
}
}
/* up dating - data swap */
memcpy( x, xt, (L+ 2 )*(L+ 2 )* sizeof ( float ) );
/* display */
if (k%DT== 0 )
{
printf( "\n\n\ntime= %6d\n" ,k);
printf( "sum = %6.2f \n" , sum);
(i= 1 ;i<=L;i++)
{
if (i% 2 == 0 )
{
for (j= 1 ;j<=L;j++)
{
c = (x[i][j] + max)* 3. /max+ 1. ;
printf( "%s" , ctab[c]);
for
/* if(x[i][j] >0.)
{
printf("$");
}
else
{
printf(" ");
} */
}
printf( "\n" );
}
}
for (i= 20 ;i<= 26 ;i++)/* This part samples actual value; as long as the program works this last part can be commented out. */
{
if (i% 2 == 0 )
{
for (j= 48 ;j<=L;j++)
{
printf( "%f " , x[i][j]);
}
printf( "\n" );
}
}
}
}
}