1 /***********************************************
2 * Sine-Gordon
3 *
4 * Solves the Sine-Gordon Equation.
5 *
6 * Dave Beazley
7 * May 25, 1993
8 *
9 * Modified for SWILL - April 21, 2001
10 ***********************************************/
11
12 #include <stdio.h>
13 #include <math.h>
14 #include <stdlib.h>
15 #include "swill.h"
16
17 #include "gifplot.h"
18
19 #define PI 3.14159265359
20
21 /* Some global variables */
22
23 double *u = 0; /* Positions */
24 double *v = 0; /* Velocities */
25 double *a = 0;
26 double *b = 0;
27 double *c = 0;
28 double *d = 0;
29 double *temp = 0; /* Temporary */
30 int Npoints;
31 double Dt;
32 double Time;
33 double Xmin;
34 double Xmax;
35 int Totalsteps = 0;
36
37 ColorMap *cm = 0; /* Colormap for plotting */
38
39 /* -----------------------------------------------------------------------------
40 * Create a new simulation
41 * ----------------------------------------------------------------------------- */
42
init_sg(int npts)43 void init_sg(int npts) {
44 int i;
45 Npoints = npts;
46 u = (double *) malloc((npts+1)*sizeof(double));
47 v = (double *) malloc((npts+1)*sizeof(double));
48 a = (double *) malloc((npts+1)*sizeof(double));
49 b = (double *) malloc((npts+1)*sizeof(double));
50 c = (double *) malloc((npts+1)*sizeof(double));
51 d = (double *) malloc((npts+1)*sizeof(double));
52
53 for (i = 0; i <= npts; i++) {
54 u[i] = v[i] = a[i] = b[i] = c[i] = d[i] = 0.0;
55 }
56 temp = (double *) malloc((npts+1)*sizeof(double));
57 }
58
59 /************************************************
60 Solves the Sine-Gordon Equation
61 Utt = Uxx - Sin U
62 Ux(0) = 0
63 Ux(L) = 0
64 U(x,0) = f(x)
65 Ut(x,0)= g(x)
66
67 ************************************************/
68
69 /*/ Solve Sine-Gordon Equation for n steps */
70
solve_sg(int nsteps)71 void solve_sg(int nsteps) {
72 double h,k,r,t,p,x,y,h2,k2;
73 void tridiagonal(double *, double *, double *, double *, double *, int);
74 int i,j;
75
76 h = (Xmax - Xmin)/(double) Npoints;
77 k = Dt;
78
79 /* Set up initial conditions */
80
81 h2 = h*h;
82 k2 = k*k;
83
84 /* Calculate tridiagonal Matrix coefficients on first iteration */
85
86 if (!Totalsteps) {
87 a[0] = 0;
88 c[0] = -2*k2;
89 c[Npoints] = 0;
90 a[Npoints] = -2*k2;
91 for (i = 1; i < Npoints; i++) {
92 a[i] = -k2;
93 c[i] = -k2;
94 }
95 for (i = 0; i <= Npoints; i++) {
96 b[i] = 2*(2*h2 + k2);
97 }
98
99 /* Calculate RHS */
100
101 d[0] = 2*k2*u[1] + (4*h2-2*k2)*u[0] + 4*h2*k*v[0] - 2*h2*k2*sin(u[0]);
102 d[Npoints] = 2*k2*u[Npoints-1] + (4*h2-2*k2)*u[Npoints] + 4*h2*k*v[Npoints] - 2*h2*k2*sin(u[Npoints]);
103 for (i = 1; i < Npoints; i++) {
104 d[i] = k2*(u[i+1]+u[i-1])+(4*h2-2*k2)*u[i] +4*h2*k*v[i] - 2*h2*k2*sin(u[i]);
105 }
106
107 tridiagonal(a,b,c,d,v, Npoints+1);
108
109 /* Adjust crank-nicholson coefficients */
110
111 for (i = 0; i <= Npoints; i++) {
112 b[i] = 2*(h2 + k2);
113 }
114 }
115
116 for (i = 1; i <= nsteps; i++) {
117 /* Set up RHS */
118 d[0] = 2*k2*v[1] + (4*h2-2*k2)*v[0] -2*h2*u[0] - 2*h2*k2*sin(v[0]);
119 for (j = 1; j < Npoints; j++) {
120 d[j] = k2*(v[j+1]+v[j-1])+(4*h2-2*k2)*v[j] - 2*h2*u[j] - 2*h2*k2*sin(v[j]);
121 }
122 d[Npoints] = 2*k2*v[Npoints-1] + (4*h2-2*k2)*v[Npoints] - 2*h2*u[Npoints] - 2*h2*k2*sin(v[Npoints]);
123
124 tridiagonal(a,b,c,d,temp, Npoints+1);
125
126 /* Copy solutions */
127
128 for (j = 0; j <= Npoints; j++) {
129 u[j] = v[j];
130 v[j] = temp[j];
131 }
132
133 Time += Dt;
134 Totalsteps++;
135 }
136 }
137
138 /* Equation for a kink with velocity V and center x0 */
eval_kink(double x,double v,double x0,double k)139 double eval_kink(double x, double v, double x0, double k) {
140 double y;
141 y = 4*atan(exp((x-x0)/sqrt(1-v*v))) + 2*PI*k;
142 return(y);
143 }
eval_dkink(double x,double v,double x0)144 double eval_dkink(double x, double v, double x0) {
145 double y;
146 y = -4*v*exp((x-x0)/sqrt(1-v*v))/(sqrt(1-v*v)*(1+exp((x-x0)/sqrt(1-v*v))*exp((x-x0)/sqrt(1-v*v))));
147 return(y);
148 }
149
eval_antikink(double x,double v,double x0,double k)150 double eval_antikink(double x, double v, double x0, double k) {
151 double y;
152 y = 4*atan(exp(-(x-x0)/sqrt(1-v*v))) + 2*PI*k;
153 return(y);
154 }
155
eval_dantikink(double x,double v,double x0)156 double eval_dantikink(double x, double v, double x0) {
157 double y;
158 y = 4*v*exp(-(x-x0)/sqrt(1-v*v))/(sqrt(1-v*v)*(1+exp(-(x-x0)/sqrt(1-v*v))*exp(-(x-x0)/sqrt(1-v*v))));
159 return(y);
160 }
161
162 /* Kink equations */
163
164 /* Add a kink to the initial state */
sg_kink(double vel,double x0,double k)165 void sg_kink(double vel, double x0, double k) {
166
167 double h,x;
168 int i;
169
170 h = (Xmax - Xmin)/ (double) Npoints;
171 for (i = 0; i <= Npoints; i++) {
172 x = Xmin+i*h;
173 u[i] += eval_kink(x,vel,x0,k);
174 v[i] += eval_dkink(x,vel,x0);
175 }
176
177 }
178
sg_antikink(double vel,double x0,double k)179 void sg_antikink(double vel, double x0, double k) {
180
181 double h,x;
182 int i;
183
184 h = (Xmax - Xmin)/ (double) Npoints;
185 for (i = 0; i <= Npoints; i++) {
186 x = Xmin+i*h;
187 u[i] += eval_antikink(x,vel,x0,k);
188 v[i] += eval_dantikink(x,vel,x0);
189 }
190 }
191
192 /* Solve tridiagonal matrix */
193 #define MAXSIZE 128000
194
tridiagonal(double * a,double * b,double * c,double * d,double * u,int N)195 void tridiagonal(double *a, double *b, double *c, double *d, double *u, int N) {
196
197 int i,j,k;
198 static double P[MAXSIZE];
199 static double R[MAXSIZE];
200
201 P[0] = b[0];
202 R[0] = d[0]/b[0];
203
204 for (i = 1; i < N; i++) {
205 P[i] = b[i] - (a[i]*c[i-1])/P[i-1];
206 R[i] = (d[i] - a[i]*R[i-1])/P[i];
207 }
208
209 u[N-1] = R[N-1];
210 for (i = N-2; i >= 0; i--)
211 u[i] = R[i] - (c[i]*u[i+1])/P[i];
212
213
214 }
215
216 /* Print the points out */
print_points()217 void print_points() {
218 int i;
219 double h;
220 h = (Xmax - Xmin)/Npoints;
221 printf("Points: Time = %g\n", Time);
222 for (i = 0; i < Npoints; i++) {
223 printf("%20g %20g\n", i*h, u[i]);
224 }
225 }
226
abort_simulation()227 void abort_simulation() {
228 /* exit(1);*/
229 printf("Really, this would normally abort.\n");
230 }
231
print_info()232 void print_info() {
233 printf("Sine Gordon solver\n\n");
234 printf("Points : %d\n", Npoints);
235 printf("Dt : %g\n", Dt);
236 printf("Xmin : %g\n", Xmin);
237 printf("Xmax : %g\n", Xmax);
238 printf("Current time : %g\n", Time);
239 }
240
241 /* This function makes a plot and writes it to the file f */
make_plot(FILE * f,Plot2D * p)242 void make_plot(FILE *f, Plot2D *p) {
243
244 double x1,x2,y1,y2;
245 double h;
246 char buffer[512000];
247 int len;
248 int i;
249
250 Plot2D_clear(p,BLACK);
251 Plot2D_xaxis(p,0,0,(Xmax-Xmin)/10,4,WHITE);
252 Plot2D_yaxis(p,0,0,5.0,4,WHITE);
253
254 h = (Xmax-Xmin)/Npoints;
255 x1 = Xmin;
256 y1 = u[0];
257 for (i = 1; i < Npoints; i++) {
258 x2 = i*h;
259 y2 = u[i];
260 Plot2D_line(p,x1,y1,x2,y2,YELLOW);
261 x1 = x2;
262 y1 = y2;
263 }
264 sprintf(buffer,"SineGordon : Time = %0.5f", Time);
265
266 FrameBuffer_drawstring(p->frame,10,580, WHITE, BLACK, buffer, HORIZONTAL);
267 len = FrameBuffer_makeGIF(p->frame, cm, buffer, 512000);
268 fwrite(buffer,len,1,f);
269 }
270
main(int argc,char ** argv)271 int main(int argc, char **argv) {
272 int maxsteps;
273 double dt;
274 int npoints;
275 int outf;
276 int i;
277 FrameBuffer *f;
278 FILE *log;
279
280 Plot2D *plot = 0;
281
282 printf("Sine-Gordon Solver\n");
283 printf("Enter npoints : ");
284 scanf("%d", &npoints);
285 printf("Enter timestep : ");
286 scanf("%lf", &dt);
287 printf("Enter nsteps : ");
288 scanf("%d", &maxsteps);
289 printf("Output freq : ");
290 scanf("%d", &outf);
291
292
293 init_sg(npoints);
294 Dt = dt;
295
296
297 Xmax = 100.0;
298 Xmin = 0.0;
299
300 /* Add some kinks and antikinks */
301 sg_kink(0.9,10,0);
302 sg_antikink(0.9,20,-1);
303 sg_kink(-0.9,80,0);
304 sg_antikink(-0.9,90,-1);
305
306 sg_antikink(0.50,45,-1);
307 sg_kink(0.50,55,0);
308
309
310 i = 0;
311 /* print_points(); */
312
313 print_info();
314
315 /* Set up a nice plotting object */
316
317 f = new_FrameBuffer(600,600);
318 cm = new_ColorMap(0);
319 plot = new_Plot2D(f,0,-15,100,15);
320
321
322 /* Initialize the SWILL server */
323 swill_init(8080);
324
325 log = fopen("logfile","w");
326
327 swill_title("Sine-Gordon Solver");
328 swill_log(log);
329
330 swill_handle("stdout:points.txt",print_points,0);
331 swill_handle("stdout:info.txt",print_info,0);
332 swill_handle("stdout:abort",abort_simulation,0);
333 swill_handle("plot.gif",make_plot,plot);
334 swill_file("index.html",0);
335 swill_file("logfile",0);
336 swill_file("sg.c",0);
337 while (i < maxsteps) {
338 solve_sg(outf);
339 /* print_points(); */
340 i += outf;
341 printf("%d\n", Totalsteps);
342 swill_poll();
343 }
344 }
345
346
347
348