/*********************************************************************************/ /* */ /* Animation of wave equation in a planar domain */ /* */ /* N. Berglund, december 2012, may 2021 */ /* */ /* UPDATE 24/04: distinction between damping and "elasticity" parameters */ /* UPDATE 27/04: new billiard shapes, bug in color scheme fixed */ /* UPDATE 28/04: code made more efficient, with help of Marco Mancini */ /* */ /* Feel free to reuse, but if doing so it would be nice to drop a */ /* line to nils.berglund@univ-orleans.fr - Thanks! */ /* */ /* compile with */ /* gcc -o wave_billiard wave_billiard.c */ /* -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -O3 -fopenmp */ /* */ /* OMP acceleration may be more effective after executing */ /* export OMP_NUM_THREADS=2 in the shell before running the program */ /* */ /* To make a video, set MOVIE to 1 and create subfolder tif_wave */ /* It may be possible to increase parameter PAUSE */ /* */ /* create movie using */ /* ffmpeg -i wave.%05d.tif -vcodec libx264 wave.mp4 */ /* */ /*********************************************************************************/ /*********************************************************************************/ /* */ /* NB: The algorithm used to simulate the wave equation is highly paralellizable */ /* One could make it much faster by using a GPU */ /* */ /*********************************************************************************/ #include #include #include #include #include #include #include /* Sam Leffler's libtiff library. */ #include #define MOVIE 0 /* set to 1 to generate movie */ #define WINWIDTH 1280 /* window width */ #define WINHEIGHT 720 /* window height */ #define NX 1280 /* number of grid points on x axis */ #define NY 720 /* number of grid points on y axis */ #define YMID 360 /* mid point of display */ // #define XMIN -1.777777778 // #define XMAX 1.777777778 /* x interval */ // #define YMIN -1.0 // #define YMAX 1.0 /* y interval for 9/16 aspect ratio */ #define XMIN -2.0 #define XMAX 2.0 /* x interval */ #define YMIN -1.125 #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ #define JULIA_SCALE 1.0 /* scaling for Julia sets */ /* Choice of the billiard table */ #define B_DOMAIN 20 /* choice of domain shape, see list in global_pdes.c */ #define B_DOMAIN_B 20 /* choice of domain shape, see list in global_pdes.c */ #define CIRCLE_PATTERN 2 /* pattern of circles, see list in global_pdes.c */ #define CIRCLE_PATTERN_B 11 /* pattern of circles, see list in global_pdes.c */ #define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */ #define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */ #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ #define RANDOM_POLY_ANGLE_B 0 /* set to 1 to randomize angle of polygons */ #define XDEP_POLY_ANGLE 0 /* set to 1 to rotate polygons depending on x coordinate */ #define XDEP_POLY_ANGLE_B 1 /* set to 1 to rotate polygons depending on x coordinate */ #define POLY_ROTATION_ANGLE -0.645 /* rotation angle for |x|=1 in units of Pi/2 */ #define HEX_NONUNIF_COMPRESSSION 0.15 /* compression factor for HEX_NONUNIF pattern */ #define HEX_NONUNIF_COMPRESSSION_B -0.15 /* compression factor for HEX_NONUNIF pattern */ #define LAMBDA 0.75 /* parameter controlling the dimensions of domain */ #define MU 0.03 /* parameter controlling the dimensions of domain */ #define MUB 0.03 /* parameter controlling the dimensions of domain */ #define NPOLY 3 /* number of sides of polygon */ #define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */ #define APOLY_B 0.335 /* angle by which to turn polygon, in units of Pi/2 */ #define MDEPTH 4 /* depth of computation of Menger gasket */ #define MRATIO 3 /* ratio defining Menger gasket */ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ #define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */ #define FOCI 1 /* set to 1 to draw focal points of ellipse */ #define NGRIDX 15 /* number of grid point for grid of disks */ #define NGRIDY 20 /* number of grid point for grid of disks */ #define X_SHOOTER -0.2 #define Y_SHOOTER -0.6 #define X_TARGET 0.4 #define Y_TARGET 0.7 /* shooter and target positions in laser fight */ #define ISO_XSHIFT_LEFT -1.65 #define ISO_XSHIFT_RIGHT 0.4 #define ISO_YSHIFT_LEFT -0.05 #define ISO_YSHIFT_RIGHT -0.05 #define ISO_SCALE 0.85 /* coordinates for isospectral billiards */ /* You can add more billiard tables by adapting the functions */ /* xy_in_billiard and draw_billiard below */ /* Physical parameters of wave equation */ #define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */ #define OSCILLATE_LEFT 0 /* set to 1 to add oscilating boundary condition on the left */ #define OSCILLATE_TOPBOT 0 /* set to 1 to enforce a planar wave on top and bottom boundary */ #define OMEGA 0.0 /* frequency of periodic excitation */ #define AMPLITUDE 0.025 /* amplitude of periodic excitation */ #define COURANT 0.02 /* Courant number */ #define COURANTB 0.004 /* Courant number in medium B */ #define GAMMA 0.0 /* damping factor in wave equation */ #define GAMMAB 1.0e-8 /* damping factor in wave equation */ #define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */ #define GAMMA_TOPBOT 1.0e-6 /* damping factor on boundary */ #define KAPPA 0.0 /* "elasticity" term enforcing oscillations */ #define KAPPA_SIDES 5.0e-4 /* "elasticity" term on absorbing boundary */ #define KAPPA_TOPBOT 0.0 /* "elasticity" term on absorbing boundary */ /* The Courant number is given by c*DT/DX, where DT is the time step and DX the lattice spacing */ /* The physical damping coefficient is given by GAMMA/(DT)^2 */ /* Increasing COURANT speeds up the simulation, but decreases accuracy */ /* For similar wave forms, COURANT^2*GAMMA should be kept constant */ /* Boundary conditions, see list in global_pdes.c */ // #define B_COND 2 #define B_COND 3 /* Parameters for length and speed of simulation */ #define NSTEPS 3750 /* number of frames of movie */ #define NVID 25 /* number of iterations between images displayed on screen */ #define NSEG 100 /* number of segments of boundary */ #define INITIAL_TIME 200 /* time after which to start saving frames */ #define COMPUTE_ENERGIES 1 /* set to 1 to compute and print energies */ #define BOUNDARY_WIDTH 2 /* width of billiard boundary */ #define PAUSE 1000 /* number of frames after which to pause */ #define PSLEEP 1 /* sleep time during pause */ #define SLEEP1 1 /* initial sleeping time */ #define SLEEP2 1 /* final sleeping time */ #define END_FRAMES 100 /* number of still frames at end of movie */ /* Parameters of initial condition */ #define INITIAL_AMP 0.2 /* amplitude of initial condition */ #define INITIAL_VARIANCE 0.002 /* variance of initial condition */ #define INITIAL_WAVELENGTH 0.1 /* wavelength of initial condition */ /* Plot type, see list in global_pdes.c */ #define PLOT 1 /* Color schemes */ #define COLOR_PALETTE 14 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* background */ #define BLACK_TEXT 0 /* set to 1 to write text in black instead of white */ #define COLOR_SCHEME 3 /* choice of color scheme, see list in global_pdes.c */ #define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ #define SLOPE 10.0 /* sensitivity of color on wave amplitude */ #define PHASE_FACTOR 1.0 /* factor in computation of phase in color scheme P_3D_PHASE */ #define PHASE_SHIFT 0.0 /* shift of phase in color scheme P_3D_PHASE */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ #define E_SCALE 500.0 /* scaling factor for energy representation */ #define LOG_SCALE 1.5 /* scaling factor for energy log representation */ #define LOG_SHIFT 1.0 /* shift of colors on log scale */ #define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */ #define COLORHUE 260 /* initial hue of water color for scheme C_LUM */ #define COLORDRIFT 0.0 /* how much the color hue drifts during the whole simulation */ #define LUMMEAN 0.5 /* amplitude of luminosity variation for scheme C_LUM */ #define LUMAMP 0.3 /* amplitude of luminosity variation for scheme C_LUM */ #define HUEMEAN 220.0 /* mean value of hue for color scheme C_HUE */ #define HUEAMP -220.0 /* amplitude of variation of hue for color scheme C_HUE */ #define DRAW_COLOR_SCHEME 0 /* set to 1 to plot the color scheme */ #define COLORBAR_RANGE 4.0 /* scale of color scheme bar */ #define COLORBAR_RANGE_B 12.0 /* scale of color scheme bar for 2nd part */ #define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ /* For debugging purposes only */ #define FLOOR 0 /* set to 1 to limit wave amplitude to VMAX */ #define VMAX 5.0 /* max value of wave amplitude */ /* the following constants are only used by wave_billiard and wave_3d so far */ #define COMPARISON 0 /* set to 1 to compare two different patterns */ #define OSCILLATION_SCHEDULE 3 /* oscillation schedule, see list in global_pdes.c */ #define ACHIRP 0.2 /* acceleration coefficient in chirp */ #define DAMPING 0.0 /* damping of periodic excitation */ /* end of constants only used by wave_billiard and wave_3d */ /* for compatibility with sub_wave and sub_maze */ #define NXMAZE 7 /* width of maze */ #define NYMAZE 7 /* height of maze */ #define MAZE_MAX_NGBH 4 /* max number of neighbours of maze cell */ #define RAND_SHIFT 24 /* seed of random number generator */ #define MAZE_XSHIFT 0.0 /* horizontal shift of maze */ #define ADD_POTENTIAL 0 #define POT_MAZE 7 #define POTENTIAL 0 /* end of constants only used by sub_wave and sub_maze */ #include "global_pdes.c" /* constants and global variables */ #include "sub_maze.c" /* support for generating mazes */ #include "sub_wave.c" /* common functions for wave_billiard, heat and schrodinger */ #include "wave_common.c" /* common functions for wave_billiard, wave_comparison, etc */ #include "sub_wave_comp.c" /* some functions specific to wave_comparison */ double courant2, courantb2; /* Courant parameters squared */ double compute_energy_x(int i, double *phi[NX], double *psi[NX], short int *xy_in[NX]) /* compute energy in column i */ { double energy = 0.0; int j; for (j=0; j= COL_TURBO) color_scheme_asym(COLOR_SCHEME, energy, scale, time, rgb); else color_scheme(COLOR_SCHEME, energy, scale, time, rgb); break; } case (P_MIXED): { if (j > NY/2) color_scheme(COLOR_SCHEME, phi[i][j], scale, time, rgb); else color_scheme(COLOR_SCHEME, compute_energy(phi, psi, xy_in, i, j), scale, time, rgb); break; } } glColor3f(rgb[0], rgb[1], rgb[2]); glVertex2i(i, j); glVertex2i(i+1, j); glVertex2i(i+1, j+1); glVertex2i(i, j+1); } } glEnd (); /* compute and plot energies */ for (i=0; i= jmid) jplus -= jmid; jminus = (j-1); if (jminus < 0) jminus += jmid; } /* imposing linear wave on top and bottom by making Laplacian 1d */ if (OSCILLATE_TOPBOT) { if (j == NY-1) jminus = NY-1; else if (j == 0) jplus = 0; } delta = phi_in[iplus][j] + phi_in[iminus][j] + phi_in[i][jplus] + phi_in[i][jminus] - 4.0*phi_in[i][j]; x = phi_in[i][j]; y = psi_in[i][j]; /* evolve phi */ if ((B_COND == BC_PERIODIC)||(B_COND == BC_DIRICHLET)) phi_out[i][j] = -y + 2*x + cc*delta - KAPPA*x - gamma*(x-y); else if ((B_COND == BC_ABSORBING)||(B_COND == BC_ABS_REFLECT)) { if ((i>0)&&(i0)&&(j0)&&(i VMAX) phi_out[i][j] = VMAX; if (phi_out[i][j] < -VMAX) phi_out[i][j] = -VMAX; if (psi_out[i][j] > VMAX) psi_out[i][j] = VMAX; if (psi_out[i][j] < -VMAX) psi_out[i][j] = -VMAX; } } } } // printf("phi(0,0) = %.3lg, psi(0,0) = %.3lg\n", phi[NX/2][NY/2], psi[NX/2][NY/2]); } void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX], short int *xy_in[NX]) /* time step of field evolution */ /* phi is value of field at time t, psi at time t-1 */ { int i, j, iplus, iminus, jplus, jminus, jmid = NY/2; double delta, x, y, c, cc, gamma; static long time = 0; static double tc[NX][NY/2], tcc[NX][NY/2], tgamma[NX][NY/2]; static short int first = 1; time++; /* initialize tables with wave speeds and dissipation */ if (first) { for (i=0; i VMAX) phi_out[i][j] = VMAX; if (phi_out[i][j] < -VMAX) phi_out[i][j] = -VMAX; } } } // printf("phi(0,0) = %.3lg, psi(0,0) = %.3lg\n", phi[NX/2][NY/2], psi[NX/2][NY/2]); } void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *xy_in[NX]) /* time step of field evolution */ /* phi is value of field at time t, psi at time t-1 */ { // evolve_wave_half_old(phi, psi, phi_tmp, psi_tmp, xy_in); // evolve_wave_half_old(phi_tmp, psi_tmp, phi, psi, xy_in); evolve_wave_half(phi, psi, tmp, xy_in); evolve_wave_half(tmp, phi, psi, xy_in); evolve_wave_half(psi, tmp, phi, xy_in); } void animation() { double time, scale, energies[6], top_energy, bottom_energy; double *phi[NX], *psi[NX], *tmp[NX]; short int *xy_in[NX]; int i, j, s; /* Since NX and NY are big, it seemed wiser to use some memory allocation here */ for (i=0; i= INITIAL_TIME) save_frame(); else printf("Initial phase time %i of %i\n", i, INITIAL_TIME); /* it seems that saving too many files too fast can cause trouble with the file system */ /* so this is to make a pause from time to time - parameter PAUSE may need adjusting */ if (i % PAUSE == PAUSE - 1) { printf("Making a short pause\n"); sleep(PSLEEP); s = system("mv wave*.tif tif_wave/"); } } } if (MOVIE) { for (i=0; i