/*********************************************************************************/ /* */ /* Animation of heat equation in a planar domain */ /* */ /* N. Berglund, May 2021 */ /* */ /* 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 heat heat.c */ /* -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -O3 -fopenmp */ /* */ /* To make a video, set MOVIE to 1 and create subfolder tif_heat */ /* 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 */ /* General geometrical parameters */ #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 NX 640 /* number of grid points on x axis */ // #define NY 360 /* number of grid points on y axis */ /* setting NX to WINWIDTH and NY to WINHEIGHT increases resolution */ /* but will multiply run time by 4 */ // #define XMIN -2.0 // #define XMAX 2.0 /* x interval */ #define XMIN -2.5 #define XMAX 1.5 /* x interval */ #define YMIN -1.125 #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ #define JULIA_SCALE 0.5 /* scaling for Julia sets */ /* Choice of the billiard table */ #define B_DOMAIN 26 /* choice of domain shape, see list in global_pdes.c */ #define CIRCLE_PATTERN 0 /* 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 0 /* set to 1 to randomize angle of polygons */ #define LAMBDA -1.0 /* parameter controlling the dimensions of domain */ #define MU 0.1 /* parameter controlling the dimensions of domain */ #define NPOLY 6 /* number of sides of polygon */ #define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */ #define MDEPTH 5 /* depth of computation of Menger gasket */ #define MRATIO 5 /* 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 in sub_wave.c */ /* Physical patameters of wave equation */ // #define DT 0.00001 #define DT 0.000004 // #define DT 0.000002 // #define DT 0.00000002 // #define DT 0.000000005 #define VISCOSITY 10.0 #define T_OUT 2.0 /* outside temperature */ #define T_IN 0.0 /* inside temperature */ // #define T_OUT 0.0 /* outside temperature */ // #define T_IN 2.0 /* inside temperature */ #define SPEED 0.0 /* speed of drift to the right */ /* Boundary conditions, see list in global_pdes.c */ #define B_COND 1 /* Parameters for length and speed of simulation */ #define NSTEPS 1000 /* number of frames of movie */ #define NVID 50 /* number of iterations between images displayed on screen */ // #define NVID 100 /* number of iterations between images displayed on screen */ #define NSEG 100 /* number of segments of boundary */ #define BOUNDARY_WIDTH 1 /* width of billiard boundary */ #define PAUSE 100 /* number of frames after which to pause */ #define PSLEEP 1 /* sleep time during pause */ #define SLEEP1 2 /* initial sleeping time */ #define SLEEP2 1 /* final sleeping time */ /* For debugging purposes only */ #define FLOOR 0 /* set to 1 to limit wave amplitude to VMAX */ #define VMAX 10.0 /* max value of wave amplitude */ /* Field representation */ #define FIELD_REP 1 #define F_INTENSITY 0 /* color represents intensity */ #define F_GRADIENT 1 /* color represents norm of gradient */ #define DRAW_FIELD_LINES 1 /* set to 1 to draw field lines */ #define FIELD_LINE_WIDTH 1 /* width of field lines */ #define N_FIELD_LINES 120 /* number of field lines */ #define FIELD_LINE_FACTOR 120 /* factor controlling precision when computing origin of field lines */ /* Color schemes, see list in global_pdes.c */ #define COLOR_PALETTE 10 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* black background */ #define COLOR_SCHEME 1 /* choice of color scheme */ #define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ // #define SLOPE 0.1 /* sensitivity of color on wave amplitude */ #define SLOPE 0.2 /* sensitivity of color on wave amplitude */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ #define E_SCALE 100.0 /* scaling factor for energy representation */ #define LOG_SCALE 1.0 /* scaling factor for energy log representation */ #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 180.0 /* mean value of hue for color scheme C_HUE */ // #define HUEAMP -180.0 /* amplitude of variation of hue for color scheme C_HUE */ #define HUEMEAN 359.0 /* mean value of hue for color scheme C_HUE */ #define HUEAMP -359.0 /* amplitude of variation of hue for color scheme C_HUE */ // #define HUEMEAN 270.0 /* mean value of hue for color scheme C_HUE */ // #define HUEAMP -130.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 2.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 */ #include "global_pdes.c" #include "sub_wave.c" double courant2; /* Courant parameter squared */ double dx2; /* spatial step size squared */ double intstep; /* integration step */ double intstep1; /* integration step used in absorbing boundary conditions */ void init_gaussian(double x, double y, double mean, double amplitude, double scalex, double *phi[NX], short int * xy_in[NX]) /* initialise field with gaussian at position (x,y) */ { int i, j, in; double xy[2], dist2, module, phase, scale2; scale2 = scalex*scalex; printf("Initialising field\n"); for (i=0; i= 2) phi[i][j] = T_IN*pow(0.75, (double)(in-2)); // else if (in >= 2) phi[i][j] = T_IN*pow(1.0 - 0.5*(double)(in-2), (double)(in-2)); // else if (in >= 2) phi[i][j] = T_IN*(1.0 - (double)(in-2)/((double)MDEPTH))*(1.0 - (double)(in-2)/((double)MDEPTH)); else phi[i][j] = T_OUT; } } void init_julia_set(double *phi[NX], short int * xy_in[NX]) /* change Julia set boundary condition */ { int i, j, in; double xy[2], dist2, module, phase, scale2; // printf("Changing Julia set\n"); for (i=0; i= 2) phi[i][j] = T_IN; } } /*********************/ /* animation part */ /*********************/ void compute_gradient(double *phi[NX], double *nablax[NX], double *nablay[NX]) /* compute the gradient of the field */ { int i, j, iplus, iminus, jplus, jminus; double dx; dx = (XMAX-XMIN)/((double)NX); for (i=0; i NX-1) ij[0] = NX-1; if (ij[1] < 0) ij[1] = 0; if (ij[1] > NY-1) ij[1] = NY-1; nabx = nablax[ij[0]][ij[1]]; naby = nablay[ij[0]][ij[1]]; norm2 = nabx*nabx + naby*naby; if (norm2 > 1.0e-14) { /* avoid too large step size */ if (norm2 < 1.0e-9) norm2 = 1.0e-9; norm = sqrt(norm2); x1 = x1 + delta*nabx/norm; y1 = y1 + delta*naby/norm; } else cont = 0; if (!xy_in[ij[0]][ij[1]]) cont = 0; /* stop if the boundary is hit */ // if (xy_in[ij[0]][ij[1]] != 1) cont = 0; // printf("x1 = %.3lg \t y1 = %.3lg \n", x1, y1); xy_to_pos(x1, y1, pos); glVertex2d(pos[0], pos[1]); i++; } glEnd(); } void draw_wave(double *phi[NX], short int *xy_in[NX], double scale, int time) /* draw the field */ { int i, j, iplus, iminus, jplus, jminus, ij[2], counter = 0; static int first = 1; double rgb[3], xy[2], x1, y1, x2, y2, dx, value, angle, dangle, intens, deltaintens, sum = 0.0; double *nablax[NX], *nablay[NX]; static double linex[N_FIELD_LINES*FIELD_LINE_FACTOR], liney[N_FIELD_LINES*FIELD_LINE_FACTOR], distance[N_FIELD_LINES*FIELD_LINE_FACTOR], integral[N_FIELD_LINES*FIELD_LINE_FACTOR + 1]; for (i=0; i 0) integral[i] = integral[i-1] + intens; else integral[i] = intens; } deltaintens = integral[N_FIELD_LINES*FIELD_LINE_FACTOR-1]/((double)N_FIELD_LINES); // printf("delta = %.5lg\n", deltaintens); i = 0; draw_field_line(linex[0], liney[0], xy_in, nablax, nablay, 0.00002, 100000); for (j = 1; j < N_FIELD_LINES+1; j++) { while ((integral[i] <= j*deltaintens)&&(i < N_FIELD_LINES*FIELD_LINE_FACTOR)) i++; draw_field_line(linex[i], liney[i], xy_in, nablax, nablay, 0.00002, 100000); counter++; } printf("%i lines\n", counter); } for (i=0; i0)&&(i0)&&(j 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 *phi_tmp[NX], short int *xy_in[NX]) /* time step of field evolution */ { evolve_wave_half(phi, phi_tmp, xy_in); evolve_wave_half(phi_tmp, phi, xy_in); } double compute_variance(double *phi[NX], short int * xy_in[NX]) /* compute the variance (total probability) of the field */ { int i, j, n = 0; double variance = 0.0; for (i=1; i= 0.0) sprintf(message, "c = %.5f + %.5f i", julia_x, julia_y); else sprintf(message, "c = %.5f %.5f i", julia_x, julia_y); xy_to_pos(XMIN + 0.1, YMAX - 0.2, pos); write_text(pos[0], pos[1], message); } void set_Julia_parameters(int time, double *phi[NX], short int *xy_in[NX]) { double jangle, cosj, sinj, radius = 0.15; jangle = (double)time*DPI/(double)NSTEPS; // jangle = (double)time*0.001; // jangle = (double)time*0.0001; cosj = cos(jangle); sinj = sin(jangle); julia_x = -0.9 + radius*cosj; julia_y = radius*sinj; init_julia_set(phi, xy_in); printf("Julia set parameters : i = %i, angle = %.5lg, cx = %.5lg, cy = %.5lg \n", time, jangle, julia_x, julia_y); } void set_Julia_parameters_cardioid(int time, double *phi[NX], short int *xy_in[NX]) { double jangle, cosj, sinj, yshift; jangle = pow(1.05 + (double)time*0.00003, 0.333); yshift = 0.02*sin((double)time*PID*0.002); // jangle = pow(1.0 + (double)time*0.00003, 0.333); // jangle = pow(0.05 + (double)time*0.00003, 0.333); // jangle = pow(0.1 + (double)time*0.00001, 0.333); // yshift = 0.04*sin((double)time*PID*0.002); cosj = cos(jangle); sinj = sin(jangle); julia_x = 0.5*(cosj*(1.0 - 0.5*cosj) + 0.5*sinj*sinj); julia_y = 0.5*sinj*(1.0-cosj) + yshift; // julia_x = 0.5*(cosj*(1.0 - 0.5*cosj) + 0.5*sinj*sinj); // julia_y = 0.5*sinj*(1.0-cosj); init_julia_set(phi, xy_in); printf("Julia set parameters : i = %i, angle = %.5lg, cx = %.5lg, cy = %.5lg \n", time, jangle, julia_x, julia_y); } void animation() { double time, scale, dx, var, jangle, cosj, sinj; double *phi[NX], *phi_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