From de54999d9883f124ca338f116f9df3f68d11df02 Mon Sep 17 00:00:00 2001 From: Nils Berglund <83530463+nilsberglund-orleans@users.noreply.github.com> Date: Sat, 8 Jul 2023 23:28:23 +0200 Subject: [PATCH] Add files via upload --- global_3d.c | 20 +- global_ljones.c | 11 +- global_particles.c | 2 +- global_pdes.c | 12 +- lennardjones.c | 246 +++++--- particle_billiard.c | 33 +- rde.c | 621 +++++++++++++++----- sub_hashgrid.c | 1 + sub_lj.c | 472 ++++++++++++++- sub_rde.c | 1350 ++++++++++++++++++++++++++++++++++++++----- sub_wave.c | 323 ++++++++++- sub_wave_3d.c | 78 ++- sub_wave_3d_rde.c | 176 ++++++ sub_wave_comp.c | 132 ++++- wave_3d.c | 163 +++--- wave_billiard.c | 150 ++--- wave_common.c | 3 +- wave_comparison.c | 272 ++++++--- wave_energy.c | 6 - 19 files changed, 3393 insertions(+), 678 deletions(-) diff --git a/global_3d.c b/global_3d.c index 203ef5c..14a3de8 100644 --- a/global_3d.c +++ b/global_3d.c @@ -27,6 +27,7 @@ #define E_SCHRODINGER 5 /* Schrodinger equation */ #define E_EULER_INCOMP 6 /* incompressible Euler equation */ #define E_EULER_COMP 7 /* compressible Euler equation */ +#define E_SHALLOW_WATER 8 /* shallow water equation */ /* Choice of potential */ @@ -51,7 +52,14 @@ #define GF_ELLIPSE 2 /* repelling ellipse */ #define GF_AIRFOIL 3 /* curved repelling ellipse */ #define GF_WING 4 /* wing shape */ -#define GF_COMPUTE_FROM_BC 5 /* compute force field as gradient of bc_field */ +#define GF_COMPUTE_FROM_BC 5 /* compute force field as gradient of bc_field2 */ + +/* Choice of water depth for shallow water equation */ + +#define SH_CIRCLE 1 /* circular shallower obstacle */ +#define SH_CIRCLES 2 /* shallow obstacle specified by CIRCLE_PATTERN */ +#define SH_COAST 3 /* depth varying with x-coordinate */ +#define SH_COAST_MONOTONE 4 /* depth decreasing with x-coordinate */ /* macros to avoid unnecessary computations in 3D plots */ @@ -74,6 +82,7 @@ int global_time = 0; +double max_depth = 1.0; /* structure used for color and height representations */ /* possible extra fields: zfield, cfield, interpolated coordinates */ @@ -111,11 +120,20 @@ typedef struct double cos_angle; /* cos of angle between normal vector and direction of light */ double log_vorticity; /* logarithm of vorticity (for Euler equation) */ double Lpressure; /* Laplacian of pressure (for Euler equation) */ + double height; /* wave height (for shallow wave equation) */ double dxu, dyu, dxv, dyv; /* gradient of velocity field (for compressible Euler equation) */ double rgb[3]; /* RGB color code */ double *p_zfield[2]; /* pointers to z field (second pointer for option DOUBLE_MOVIE) */ double *p_cfield[2]; /* pointers to color field (second pointer for option DOUBLE_MOVIE) */ + double depth; /* water depth */ + double cos_depth_angle; /* cos of angle of depth profile */ + double gradx, grady; /* gradient of water depth */ } t_rde; +typedef struct +{ + double depth; /* water depth */ + double gradx, grady; /* gradient of water depth */ +} t_swater_depth; diff --git a/global_ljones.c b/global_ljones.c index 83b6a2e..7bb9adc 100644 --- a/global_ljones.c +++ b/global_ljones.c @@ -44,6 +44,7 @@ #define O_GENUS_TWO 2 /* obstacles in corners of L-shape domeain (for genus 2 b.c.) */ #define O_POOL_TABLE 3 /* obstacles around pockets of pool table */ #define O_HLINE_HOLE_SPOKES 181 /* tips of spokes for S_HLINE_HOLE_SPOKES segment pattern */ +#define O_CIRCLE 4 /* one circle at the origin */ /* pattern of additional repelling segments */ #define S_RECTANGLE 0 /* segments forming a rectangle */ @@ -64,6 +65,7 @@ #define S_DAM_WITH_HOLE 13 /* segments forming a dam in which a hole can open */ #define S_DAM_WITH_HOLE_AND_RAMP 14 /* segments forming a dam in which a hole can open */ #define S_MAZE 15 /* segments forming a maze */ +#define S_MAZE_DIAG 151 /* segments forming a maze with diagonally opposed exits */ #define S_EXT_RECTANGLE 16 /* particles outside a rectangle */ #define S_DAM_BRICKS 17 /* dam made of several bricks */ #define S_HLINE_HOLE 18 /* horizontal line with a hole in the bottom */ @@ -84,6 +86,7 @@ #define I_VICSEK 8 /* Vicsek-type interaction */ #define I_VICSEK_REPULSIVE 9 /* Vicsek-type interaction with harmonic repulsion */ #define I_VICSEK_SPEED 10 /* Vicsek-type interaction with speed adjustment */ +#define I_VICSEK_SHARK 11 /* Vicsek-type interaction with speed adjustment, and one shark */ /* Boundary conditions */ @@ -103,6 +106,7 @@ #define BC_GENUS_TWO 14 /* surface of genus 2, obtained by identifying opposite sides of an L shape */ #define BC_ABSORBING 20 /* "no-return" boundary conditions outside BC area */ #define BC_REFLECT_ABS 21 /* reflecting on lower boundary, and "no-return" boundary conditions outside BC area */ +#define BC_REFLECT_ABS_BOTTOM 22 /* absorbing on lower boundary, and reflecting elsewhere */ /* Regions for partial thermostat couplings */ @@ -122,6 +126,7 @@ #define RCK_DISC 0 /* disc-shaped rocket */ #define RCK_RECT 1 /* rectangular rocket */ #define RCK_RECT_HAT 2 /* rectangular rocket with a hat */ +#define RCK_RECT_BAR 3 /* rectangular rocket with a hat and a separating bar */ /* Nozzle shapes */ @@ -131,6 +136,7 @@ #define NZ_CONE 3 /* cone-shaped nozzle */ #define NZ_TRUMPET 4 /* trumpet-shaped nozzle */ #define NZ_BROAD 5 /* broad straight nozzle */ +#define NZ_DELAVAL 6 /* a type of de Laval nozzle */ #define NZ_NONE 99 /* no nozzle */ /* Types of chemical reactions */ @@ -151,6 +157,7 @@ #define CHEM_ABCDABC 13 /* reactions A + B <-> C, A + C <-> D */ #define CHEM_BZ 14 /* simplified Belousov-Zhabotinski reaction with 6 types (Oregonator) */ #define CHEM_BRUSSELATOR 15 /* Brusselator oscillating reaction */ +#define CHEM_ABDACBE 16 /* A + B -> D, A + C -> B + E */ /* Initial conditions for chemical reactions */ @@ -164,6 +171,7 @@ #define IC_BZ 6 /* initial state for BZ reaction */ #define IC_SIGNX 7 /* type 1 or 2 depending on sign of x */ #define IC_TWOROCKETS 8 /* type 1 or 2 depending on rocket position */ +#define IC_TWOROCKETS_TWOFUELS 9 /* type 1 and 2 or 1 and 3 depending on rocket */ /* Plot types */ @@ -179,6 +187,7 @@ #define P_DIFF_NEIGHB 9 /* colors represent number of neighbours of different type */ #define P_THERMOSTAT 10 /* colors show which particles are coupled to the thermostat */ #define P_INITIAL_POS 11 /* colors depend on initial position of particle */ +#define P_NUMBER 12 /* colors depend on particle number */ /* Color schemes */ @@ -206,7 +215,7 @@ #define COL_TURBO_CYCLIC 101 /* TURBO color palette (by Anton Mikhailov) corrected to be cyclic, beta */ -#define VICSEK_INT ((INTERACTION == I_VICSEK)||(INTERACTION == I_VICSEK_REPULSIVE)||(INTERACTION == I_VICSEK_SPEED)) +#define VICSEK_INT ((INTERACTION == I_VICSEK)||(INTERACTION == I_VICSEK_REPULSIVE)||(INTERACTION == I_VICSEK_SPEED)||(INTERACTION == I_VICSEK_SHARK)) typedef struct { diff --git a/global_particles.c b/global_particles.c index 41b1499..319be34 100644 --- a/global_particles.c +++ b/global_particles.c @@ -58,7 +58,7 @@ double x_shooter = -0.2, y_shooter = -0.6, x_target = 0.4, y_target = 0.7; #define D_ALT_REU 11 /* alternating between star and Reuleaux */ #define D_ANGLE 12 /* angular sector */ #define D_LSHAPE 13 /* L-shaped billiard for conical singularity */ -#define D_GENUSN 14 /* polygon with identifies opposite sides */ +#define D_GENUSN 14 /* polygon with identified opposite sides */ #define D_PARABOLAS 15 /* polygon with parabolic sides */ #define D_PENROSE 16 /* Penrose solution to illumination problem */ diff --git a/global_pdes.c b/global_pdes.c index e4ad7d0..6674c62 100644 --- a/global_pdes.c +++ b/global_pdes.c @@ -69,6 +69,7 @@ #define D_MAZE 53 /* maze */ #define D_MAZE_CLOSED 54 /* closed maze */ #define D_MAZE_CHANNELS 541 /* maze with two channels attached */ +#define D_MAZE_CHANNELS_INT 542 /* maze with two channels attached, distance defined from interior of cells */ #define D_CHESSBOARD 55 /* chess board configuration */ #define D_TRIANGLE_TILES 56 /* triangular tiling */ #define D_HEX_TILES 57 /* honeycomb tiling */ @@ -128,6 +129,8 @@ #define IOR_RANDOM_WELLS 5 /* random superposition of "wells" */ #define IOR_PERIODIC_WELLS_ROTATING 6 /* periodic superposition rotating in time */ #define IOR_PERIODIC_WELLS_ROTATING_LARGE 7 /* periodic superposition rotating in time, larger area */ +#define IOR_POISSON_WELLS 8 /* wells located on a random Poisson disc process */ +#define IOR_PPP_WELLS 9 /* wells located on a Poisson point process */ /* Boundary conditions */ @@ -229,12 +232,17 @@ #define Z_EULER_LPRESSURE 53 /* Laplacian of pressure */ #define Z_EULER_PRESSURE 54 /* pressure */ -/* for Euler compressible Euler equation */ +/* for compressible Euler equation */ #define Z_EULER_DENSITY 60 /* density */ #define Z_EULER_SPEED 61 /* norm of velocity */ #define Z_EULERC_VORTICITY 62 /* vorticity of velocity */ #define Z_EULER_DIRECTION 63 /* direction of velocity */ -#define Z_EULER_DIRECTION_SPEED 64 /* hut for direction of velocity, luminosity for speed */ +#define Z_EULER_DIRECTION_SPEED 64 /* hue for direction of velocity, luminosity for speed */ + +/* for shallow water equation */ +#define Z_SWATER_HEIGHT 70 /* height */ +#define Z_SWATER_SPEED 71 /* speed */ +#define Z_SWATER_DIRECTION_SPEED 74 /* hue for direction of velocity, luminosity for speed */ /* special boundary conditions for Euler equation */ #define BCE_TOPBOTTOM 1 /* special flow at top and bottom */ diff --git a/lennardjones.c b/lennardjones.c index a854d13..5848c39 100644 --- a/lennardjones.c +++ b/lennardjones.c @@ -36,15 +36,15 @@ #include #include -#define MOVIE 0 /* set to 1 to generate movie */ -#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */ +#define MOVIE 0 /* set to 1 to generate movie */ +#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */ #define SAVE_MEMORY 1 /* set to 1 to save memory while saving frames */ #define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */ #define TIME_LAPSE 0 /* set to 1 to add a time-lapse movie at the end */ /* so far incompatible with double movie */ #define TIME_LAPSE_FACTOR 3 /* factor of time-lapse movie */ -#define TIME_LAPSE_FIRST 0 /* set to 1 to show time-lapse version first */ +#define TIME_LAPSE_FIRST 1 /* set to 1 to show time-lapse version first */ #define SAVE_TIME_SERIES 0 /* set to 1 to save time series of particle positions */ @@ -58,10 +58,18 @@ #define YMIN -1.125 #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ -#define INITXMIN -1.97 -#define INITXMAX 1.97 /* x interval for initial condition */ +#define INITXMIN -2.0 +#define INITXMAX 2.0 /* x interval for initial condition */ #define INITYMIN -1.1 #define INITYMAX 1.1 /* y interval for initial condition */ +// #define INITYMAX 7.0 /* y interval for initial condition */ +// #define INITYMAX 9.0 /* y interval for initial condition */ + +#define ADDXMIN -1.97 +#define ADDXMAX -0.8 /* x interval for adding particles */ +#define ADDYMIN -1.125 +#define ADDYMAX 1.125 /* y interval for adding particles */ + #define BCXMIN -2.0 #define BCXMAX 2.0 /* x interval for boundary condition */ @@ -71,39 +79,42 @@ #define OBSXMIN -2.0 #define OBSXMAX 2.0 /* x interval for motion of obstacle */ -#define CIRCLE_PATTERN 8 /* pattern of circles, see list in global_ljones.c */ +#define CIRCLE_PATTERN 1 /* pattern of circles, see list in global_ljones.c */ -#define ADD_FIXED_OBSTACLES 0 /* set to 1 do add fixed circular obstacles */ -#define OBSTACLE_PATTERN 181 /* pattern of obstacles, see list in global_ljones.c */ +#define ADD_INITIAL_PARTICLES 0 /* set to 1 to add a second type of particles */ +#define CIRCLE_PATTERN_B 1 /* pattern of circles for additional particles */ + +#define ADD_FIXED_OBSTACLES 1 /* set to 1 do add fixed circular obstacles */ +#define OBSTACLE_PATTERN 4 /* pattern of obstacles, see list in global_ljones.c */ #define ADD_FIXED_SEGMENTS 0 /* set to 1 to add fixed segments as obstacles */ -#define SEGMENT_PATTERN 181 /* pattern of repelling segments, see list in global_ljones.c */ -#define ROCKET_SHAPE 2 /* shape of rocket combustion chamber, see list in global_ljones.c */ -#define ROCKET_SHAPE_B 2 /* shape of second rocket */ -#define NOZZLE_SHAPE 2 /* shape of nozzle, see list in global_ljones.c */ -#define NOZZLE_SHAPE_B 4 /* shape of nozzle for second rocket, see list in global_ljones.c */ +#define SEGMENT_PATTERN 151 /* pattern of repelling segments, see list in global_ljones.c */ +#define ROCKET_SHAPE 3 /* shape of rocket combustion chamber, see list in global_ljones.c */ +#define ROCKET_SHAPE_B 3 /* shape of second rocket */ +#define NOZZLE_SHAPE 6 /* shape of nozzle, see list in global_ljones.c */ +#define NOZZLE_SHAPE_B 6 /* shape of nozzle for second rocket, see list in global_ljones.c */ #define TWO_TYPES 0 /* set to 1 to have two types of particles */ -#define TYPE_PROPORTION 0.3 /* proportion of particles of first type */ +#define TYPE_PROPORTION 0.5 /* proportion of particles of first type */ #define SYMMETRIZE_FORCE 1 /* set to 1 to symmetrize two-particle interaction, only needed if particles are not all the same */ #define CENTER_PX 0 /* set to 1 to center horizontal momentum */ #define CENTER_PY 0 /* set to 1 to center vertical momentum */ #define CENTER_PANGLE 0 /* set to 1 to center angular momentum */ -#define INTERACTION 10 /* particle interaction, see list in global_ljones.c */ -#define INTERACTION_B 1 /* particle interaction for second type of particle, see list in global_ljones.c */ -#define SPIN_INTER_FREQUENCY 1.0 /* angular frequency of spin-spin interaction */ -#define SPIN_INTER_FREQUENCY_B 2.0 /* angular frequency of spin-spin interaction for second particle type */ +#define INTERACTION 2 /* particle interaction, see list in global_ljones.c */ +#define INTERACTION_B 2 /* particle interaction for second type of particle, see list in global_ljones.c */ +#define SPIN_INTER_FREQUENCY 4.0 /* angular frequency of spin-spin interaction */ +#define SPIN_INTER_FREQUENCY_B 4.0 /* angular frequency of spin-spin interaction for second particle type */ #define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */ #define NPOISSON 100 /* number of points for Poisson C_RAND_POISSON arrangement */ -#define PDISC_DISTANCE 1.3 /* minimal distance in Poisson disc process, controls density of particles */ +#define PDISC_DISTANCE 2.8 /* minimal distance in Poisson disc process, controls density of particles */ #define PDISC_CANDIDATES 100 /* number of candidates in construction of Poisson disc process */ #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */ -#define LAMBDA 0.5 /* parameter controlling the dimensions of domain */ -#define MU 0.035 /* parameter controlling radius of particles */ -#define MU_B 0.012 /* parameter controlling radius of particles of second type */ +#define LAMBDA 0.8 /* parameter controlling the dimensions of domain */ +#define MU 0.015 /* parameter controlling radius of particles */ +#define MU_B 0.015 /* parameter controlling radius of particles of second type */ #define NPOLY 25 /* number of sides of polygon */ #define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */ #define MDEPTH 4 /* depth of computation of Menger gasket */ @@ -111,8 +122,8 @@ #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 120 /* number of grid point for grid of disks */ -#define NGRIDY 51 /* number of grid point for grid of disks */ +#define NGRIDX 140 /* number of grid point for grid of disks */ +#define NGRIDY 70 /* number of grid point for grid of disks */ #define EHRENFEST_RADIUS 0.9 /* radius of container for Ehrenfest urn configuration */ #define EHRENFEST_WIDTH 0.035 /* width of tube for Ehrenfest urn configuration */ #define TWO_CIRCLES_RADIUS_RATIO 0.8 /* ratio of radii for S_TWO_CIRCLES_EXT segment configuration */ @@ -125,11 +136,12 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 3800 /* number of frames of movie */ -#define NVID 50 /* number of iterations between images displayed on screen */ +#define NSTEPS 2800 /* number of frames of movie */ +// #define NSTEPS 200 /* number of frames of movie */ +#define NVID 60 /* number of iterations between images displayed on screen */ #define NSEG 250 /* number of segments of boundary */ -#define INITIAL_TIME 10 /* time after which to start saving frames */ -#define OBSTACLE_INITIAL_TIME 200 /* time after which to start moving obstacle */ +#define INITIAL_TIME 0 /* time after which to start saving frames */ +#define OBSTACLE_INITIAL_TIME 150 /* time after which to start moving obstacle */ #define BOUNDARY_WIDTH 1 /* width of particle boundary */ #define LINK_WIDTH 2 /* width of links between particles */ #define CONTAINER_WIDTH 4 /* width of container boundary */ @@ -148,17 +160,20 @@ /* Plot type, see list in global_ljones.c */ #define PLOT 4 -#define PLOT_B 6 /* plot type for second movie */ +#define PLOT_B 11 /* plot type for second movie */ #define DRAW_BONDS 1 /* set to 1 to draw bonds between neighbours */ #define COLOR_BONDS 1 /* set to 1 to color bonds according to length */ #define FILL_TRIANGLES 0 /* set to 1 to fill triangles between neighbours */ #define ALTITUDE_LINES 0 /* set to 1 to add horizontal lines to show altitude */ -#define COLOR_SEG_GROUPS 0 /* set to 1 to collor segment groups differently */ +#define COLOR_SEG_GROUPS 1 /* set to 1 to collor segment groups differently */ +#define N_PARTICLE_COLORS 200 /* number of colors for P_NUMBER color scheme */ /* Color schemes */ -#define COLOR_PALETTE 0 /* Color palette, see list in global_ljones.c */ +#define COLOR_PALETTE 18 /* Color palette, see list in global_ljones.c */ +#define COLOR_PALETTE_ANGLE 10 /* Color palette for angle representation */ +#define COLOR_PALETTE_INITIAL_POS 10 /* Color palette for initial position representation */ #define BLACK 1 /* background */ @@ -175,7 +190,7 @@ #define HUEMEAN 220.0 /* mean value of hue for color scheme C_HUE */ #define HUEAMP -50.0 /* amplitude of variation of hue for color scheme C_HUE */ -#define PRINT_PARAMETERS 1 /* set to 1 to print certain parameters */ +#define PRINT_PARAMETERS 0 /* set to 1 to print certain parameters */ #define PRINT_TEMPERATURE 0 /* set to 1 to print current temperature */ /* particle properties */ @@ -184,54 +199,57 @@ #define ENERGY_HUE_MAX 50.0 /* color of saturated particle */ #define PARTICLE_HUE_MIN 359.0 /* color of original particle */ #define PARTICLE_HUE_MAX 0.0 /* color of saturated particle */ -#define PARTICLE_EMAX 1.2e3 /* energy of particle with hottest color */ -#define HUE_TYPE0 260.0 /* hue of particles of type 0 */ -#define HUE_TYPE1 210.0 /* hue of particles of type 1 */ -#define HUE_TYPE2 160.0 /* hue of particles of type 2 */ -#define HUE_TYPE3 290.0 /* hue of particles of type 3 */ +#define PARTICLE_EMAX 6000.0 /* energy of particle with hottest color */ +#define HUE_TYPE0 60.0 /* hue of particles of type 0 */ +#define HUE_TYPE1 280.0 /* hue of particles of type 1 */ +#define HUE_TYPE2 140.0 /* hue of particles of type 2 */ +#define HUE_TYPE3 200.0 /* hue of particles of type 3 */ #define RANDOM_RADIUS 0 /* set to 1 for random circle radius */ #define DT_PARTICLE 3.0e-6 /* time step for particle displacement */ -#define KREPEL 0.1 /* constant in repelling force between particles */ -#define EQUILIBRIUM_DIST 2.0 /* Lennard-Jones equilibrium distance */ -#define EQUILIBRIUM_DIST_B 2.0 /* Lennard-Jones equilibrium distance for second type of particle */ +#define KREPEL 1.0 /* constant in repelling force between particles */ +#define EQUILIBRIUM_DIST 5.0 /* Lennard-Jones equilibrium distance */ +#define EQUILIBRIUM_DIST_B 5.0 /* Lennard-Jones equilibrium distance for second type of particle */ #define REPEL_RADIUS 15.0 /* radius in which repelling force acts (in units of particle radius) */ -#define DAMPING 0.0 /* damping coefficient of particles */ -#define INITIAL_DAMPING 50.0 /* damping coefficient of particles during initial phase */ +#define DAMPING 100.0 /* damping coefficient of particles */ +#define OMEGA_INITIAL 10.0 /* initial angular velocity range */ +#define INITIAL_DAMPING 5.0 /* damping coefficient of particles during initial phase */ #define DAMPING_ROT 100.0 /* dampint coefficient for rotation of particles */ -#define PARTICLE_MASS 0.25 /* mass of particle of radius MU */ -#define PARTICLE_MASS_B 0.5 /* mass of particle of radius MU */ -#define PARTICLE_INERTIA_MOMENT 0.5 /* moment of inertia of particle */ +#define PARTICLE_MASS 2.5 /* mass of particle of radius MU */ +#define PARTICLE_MASS_B 1.0 /* mass of particle of radius MU */ +#define PARTICLE_INERTIA_MOMENT 0.2 /* moment of inertia of particle */ #define PARTICLE_INERTIA_MOMENT_B 0.02 /* moment of inertia of second type of particle */ -#define V_INITIAL 20.0 /* initial velocity range */ -#define OMEGA_INITIAL 5.0 /* initial angular velocity range */ +#define V_INITIAL 0.0 /* initial velocity range */ +#define OMEGA_INITIAL 10.0 /* initial angular velocity range */ #define VICSEK_VMIN 1.0 /* minimal speed of particles in Vicsek model */ #define VICSEK_VMAX 40.0 /* minimal speed of particles in Vicsek model */ #define THERMOSTAT 0 /* set to 1 to switch on thermostat */ #define VARY_THERMOSTAT 0 /* set to 1 for time-dependent thermostat schedule */ #define SIGMA 5.0 /* noise intensity in thermostat */ -#define BETA 0.002 /* initial inverse temperature */ +#define BETA 0.02 /* initial inverse temperature */ #define MU_XI 0.01 /* friction constant in thermostat */ #define KSPRING_BOUNDARY 1.0e7 /* confining harmonic potential outside simulation region */ #define KSPRING_OBSTACLE 1.0e11 /* harmonic potential of obstacles */ -#define NBH_DIST_FACTOR 3.0 /* radius in which to count neighbours */ +#define NBH_DIST_FACTOR 2.5 /* radius in which to count neighbours */ +// #define NBH_DIST_FACTOR 4.0 /* radius in which to count neighbours */ #define GRAVITY 0.0 /* gravity acting on all particles */ -#define GRAVITY_X 0.0 /* horizontal gravity acting on all particles */ +#define GRAVITY_X 5000.0 /* horizontal gravity acting on all particles */ #define INCREASE_GRAVITY 0 /* set to 1 to increase gravity during the simulation */ #define GRAVITY_SCHEDULE 2 /* type of gravity schedule, see list in global_ljones.c */ #define GRAVITY_FACTOR 100.0 /* factor by which to increase gravity */ #define GRAVITY_INITIAL_TIME 200 /* time at start of simulation with constant gravity */ #define GRAVITY_RESTORE_TIME 700 /* time at end of simulation with gravity restored to initial value */ #define KSPRING_VICSEK 0.2 /* spring constant for I_VICSEK_SPEED interaction */ +#define VICSEK_REPULSION 10.0 /* repulsion between particles in Vicsek model */ #define ROTATION 1 /* set to 1 to include rotation of particles */ -#define COUPLE_ANGLE_TO_THERMOSTAT 0 /* set to 1 to couple angular degrees of freedom to thermostat */ +#define COUPLE_ANGLE_TO_THERMOSTAT 1 /* set to 1 to couple angular degrees of freedom to thermostat */ #define DIMENSION_FACTOR 1.0 /* scaling factor taking into account number of degrees of freedom */ -#define KTORQUE 2.0e3 /* force constant in angular dynamics */ +#define KTORQUE 50.0 /* force constant in angular dynamics */ +#define KTORQUE_BOUNDARY 1.0e6 /* constant in torque from the boundary */ #define KTORQUE_B 10.0 /* force constant in angular dynamics */ #define KTORQUE_DIFF 150.0 /* force constant in angular dynamics for different particles */ -#define KTORQUE_BOUNDARY 1.0e6 /* constant in torque from the boundary */ #define DRAW_SPIN 0 /* set to 1 to draw spin vectors of particles */ #define DRAW_SPIN_B 0 /* set to 1 to draw spin vectors of particles */ #define DRAW_CROSS 1 /* set to 1 to draw cross on particles of second type */ @@ -240,7 +258,7 @@ #define QUADRUPOLE_RATIO 0.6 /* anisotropy in quadrupole potential */ #define INCREASE_BETA 0 /* set to 1 to increase BETA during simulation */ -#define BETA_FACTOR 0.5 /* factor by which to change BETA during simulation */ +#define BETA_FACTOR 0.02 /* factor by which to change BETA during simulation */ #define N_TOSCILLATIONS 1.5 /* number of temperature oscillations in BETA schedule */ #define NO_OSCILLATION 1 /* set to 1 to have exponential BETA change only */ #define MIDDLE_CONSTANT_PHASE 2000 /* final phase in which temperature is constant */ @@ -257,7 +275,7 @@ #define CENTER_VIEW_ON_OBSTACLE 0 /* set to 1 to center display on moving obstacle */ #define RESAMPLE_Y 0 /* set to 1 to resample y coordinate of moved particles (for shock waves) */ #define NTRIALS 2000 /* number of trials when resampling */ -#define OBSTACLE_RADIUS 0.12 /* radius of obstacle for circle boundary conditions */ +#define OBSTACLE_RADIUS 0.35 /* radius of obstacle for circle boundary conditions */ #define FUNNEL_WIDTH 0.25 /* funnel width for funnel boundary conditions */ #define OBSTACLE_XMIN 0.0 /* initial position of obstacle */ #define OBSTACLE_XMAX 3.0 /* final position of obstacle */ @@ -267,7 +285,7 @@ #define N_T_AVERAGE 1 /* size of temperature averaging window */ #define MAX_PRESSURE 3.0e10 /* pressure shown in "hottest" color */ #define PARTIAL_THERMO_COUPLING 0 /* set to 1 to couple only some particles to thermostat */ -#define PARTIAL_THERMO_REGION 4 /* region for partial thermostat coupling (see list in global_ljones.c) */ +#define PARTIAL_THERMO_REGION 1 /* region for partial thermostat coupling (see list in global_ljones.c) */ #define PARTIAL_THERMO_SHIFT 0.2 /* distance from obstacle at the right of which particles are coupled to thermostat */ #define PARTIAL_THERMO_WIDTH 0.5 /* vertical size of partial thermostat coupling */ #define PARTIAL_THERMO_HEIGHT 0.25 /* vertical size of partial thermostat coupling */ @@ -279,11 +297,11 @@ #define MASS_PART_BOTTOM 10000.0 /* mass of particles at bottom */ #define NPART_BOTTOM 100 /* number of particles at the bottom */ -#define ADD_PARTICLES 0 /* set to 1 to add particles */ +#define ADD_PARTICLES 0 /* set to 1 to add particles */ #define ADD_TIME 0 /* time at which to add first particle */ -#define ADD_PERIOD 10000 /* time interval between adding further particles */ -#define N_ADD_PARTICLES 20 /* number of particles to add */ -#define FINAL_NOADD_PERIOD 200 /* final period where no particles are added */ +#define ADD_PERIOD 10 /* time interval between adding further particles */ +#define N_ADD_PARTICLES 2 /* number of particles to add */ +#define FINAL_NOADD_PERIOD 0 /* final period where no particles are added */ #define SAFETY_FACTOR 2.0 /* no particles are added at distance less than MU*SAFETY_FACTOR of other particles */ #define TRACER_PARTICLE 0 /* set to 1 to have a tracer particle */ @@ -306,49 +324,51 @@ #define MOVE_BOUNDARY 0 /* set to 1 to move repelling segments, due to force from particles */ #define SEGMENTS_MASS 40.0 /* mass of collection of segments */ #define DEACTIVATE_SEGMENT 1 /* set to 1 to deactivate last segment after a certain time */ -#define SEGMENT_DEACTIVATION_TIME 200 /* time at which to deactivate last segment */ -#define RELEASE_ROCKET_AT_DEACTIVATION 1 /* set to 1 to limit segments velocity before segment release */ +#define SEGMENT_DEACTIVATION_TIME 20 /* time at which to deactivate last segment */ +#define RELEASE_ROCKET_AT_DEACTIVATION 0 /* set to 1 to limit segments velocity before segment release */ #define SEGMENTS_X0 1.5 /* initial position of segments */ #define SEGMENTS_Y0 0.0 /* initial position of segments */ #define SEGMENTS_VX0 0.0 /* initial velocity of segments */ #define SEGMENTS_VY0 0.0 /* initial velocity of segments */ #define DAMP_SEGS_AT_NEGATIVE_Y 0 /* set to 1 to dampen segments when y coordinate is negative */ -#define MOVE_SEGMENT_GROUPS 1 /* set to 1 to group segments into moving units */ -#define SEGMENT_GROUP_MASS 1000.0 /* mass of segment group */ +#define MOVE_SEGMENT_GROUPS 0 /* set to 1 to group segments into moving units */ +#define SEGMENT_GROUP_MASS 500.0 /* mass of segment group */ #define SEGMENT_GROUP_I 1000.0 /* moment of inertia of segment group */ #define SEGMENT_GROUP_DAMPING 0.0 /* damping of segment groups */ -#define GROUP_REPULSION 1 /* set to 1 for groups of segments to repel each other */ -#define KSPRING_GROUPS 1.0e11 /* harmonic potential between segment groups */ +#define GROUP_REPULSION 0 /* set to 1 for groups of segments to repel each other */ +#define KSPRING_GROUPS 5.0e11 /* harmonic potential between segment groups */ #define GROUP_WIDTH 0.05 /* interaction width of groups */ -#define GROUP_G_REPEL 1 /* set to 1 to add repulsion between centers of mass of groups */ +#define GROUP_G_REPEL 0 /* set to 1 to add repulsion between centers of mass of groups */ #define GROUP_G_REPEL_RADIUS 1.2 /* radius within which centers of mass of groups repel each other */ -#define TRACK_SEGMENT_GROUPS 1 /* set to 1 for view to track group of segments */ +#define TRACK_SEGMENT_GROUPS 0 /* set to 1 for view to track group of segments */ #define TRACK_X_PADDING 2.0 /* distance from x boundary where tracking starts */ #define POSITION_DEPENDENT_TYPE 0 /* set to 1 to make particle type depend on initial position */ #define POSITION_Y_DEPENDENCE 0 /* set to 1 for the separation between particles to be horizontal */ +#define POSITION_DEP_SIGN -1.0 /* sign in position dependence condition */ +#define POSITION_DEP_X -0.625 /* threshold value for position-dependent type */ #define PRINT_ENTROPY 0 /* set to 1 to compute entropy */ #define REACTION_DIFFUSION 0 /* set to 1 to simulate a chemical reaction (particles may change type) */ -#define RD_REACTION 15 /* type of reaction, see list in global_ljones.c */ +#define RD_REACTION 16 /* type of reaction, see list in global_ljones.c */ #define RD_TYPES 5 /* number of types in reaction-diffusion equation */ -#define RD_INITIAL_COND 2 /* initial condition of particles */ -#define REACTION_DIST 3.5 /* maximal distance for reaction to occur */ -#define REACTION_PROB 0.5 /* probability controlling reaction term */ +#define RD_INITIAL_COND 9 /* initial condition of particles */ +#define REACTION_DIST 4.0 /* maximal distance for reaction to occur */ +#define REACTION_PROB 1.0 /* probability controlling reaction term */ #define DISSOCIATION_PROB 0.002 /* probability controlling dissociation reaction */ #define CENTER_COLLIDED_PARTICLES 0 /* set to 1 to recenter particles upon reaction (may interfere with thermostat) */ -#define EXOTHERMIC 0 /* set to 1 to make reaction exo/endothermic */ -#define DELTA_EKIN 500.0 /* change of kinetic energy in reaction */ +#define EXOTHERMIC 1 /* set to 1 to make reaction exo/endothermic */ +#define DELTA_EKIN 2000.0 /* change of kinetic energy in reaction */ #define COLLISION_TIME 15 /* time during which collisions are shown */ #define PRINT_PARTICLE_NUMBER 0 /* set to 1 to print total number of particles */ #define PLOT_PARTICLE_NUMBER 0 /* set to 1 to make of plot of particle number over time */ #define PARTICLE_NB_PLOT_FACTOR 0.5 /* expected final number of particles over initial number */ -#define PRINT_LEFT 1 /* set to 1 to print certain parameters at the top left instead of right */ +#define PRINT_LEFT 0 /* set to 1 to print certain parameters at the top left instead of right */ #define PLOT_SPEEDS 0 /* set to 1 to add a plot of obstacle speeds (e.g. for rockets) */ #define PLOT_TRAJECTORIES 0 /* set to 1 to add a plot of obstacle trajectories (e.g. for rockets) */ -#define VMAX_PLOT_SPEEDS 0.6 /* vertical scale of plot of obstacle speeds */ +#define VMAX_PLOT_SPEEDS 0.25 /* vertical scale of plot of obstacle speeds */ #define EHRENFEST_COPY 0 /* set to 1 to add equal number of larger particles (for Ehrenfest model) */ @@ -360,20 +380,20 @@ #define WALL_VMAX 100.0 /* max speed of wall */ #define WALL_TIME 0 /* time during which to keep wall */ -#define NXMAZE 10 /* width of maze */ -#define NYMAZE 10 /* height of maze */ +#define NXMAZE 12 /* width of maze */ +#define NYMAZE 12 /* height of maze */ #define MAZE_MAX_NGBH 4 /* max number of neighbours of maze cell */ -#define RAND_SHIFT 200 /* seed of random number generator */ -#define MAZE_XSHIFT 0.0 /* horizontal shift of maze */ +#define RAND_SHIFT 4 /* seed of random number generator */ +#define MAZE_XSHIFT 0.5 /* horizontal shift of maze */ #define MAZE_WIDTH 0.01 /* width of maze walls */ #define FLOOR_FORCE 1 /* set to 1 to limit force on particle to FMAX */ -#define FMAX 1.0e10 /* maximal force */ +#define FMAX 1.0e9 /* maximal force */ #define FLOOR_OMEGA 0 /* set to 1 to limit particle momentum to PMAX */ #define PMAX 1000.0 /* maximal force */ -#define HASHX 60 /* size of hashgrid in x direction */ -#define HASHY 30 /* size of hashgrid in y direction */ +#define HASHX 100 /* size of hashgrid in x direction */ +#define HASHY 50 /* size of hashgrid in y direction */ #define HASHMAX 100 /* maximal number of particles per hashgrid cell */ #define HASHGRID_PADDING 0.1 /* padding of hashgrid outside simulation window */ @@ -382,6 +402,8 @@ #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 */ +#define LIMIT_ENERGY 0 /* set to 1 to limit energy, when there is no thermostat */ + #define NO_WRAP_BC ((BOUNDARY_COND != BC_PERIODIC)&&(BOUNDARY_COND != BC_PERIODIC_CIRCLE)&&(BOUNDARY_COND != BC_PERIODIC_TRIANGLE)&&(BOUNDARY_COND != BC_KLEIN)&&(BOUNDARY_COND != BC_PERIODIC_FUNNEL)&&(BOUNDARY_COND != BC_BOY)&&(BOUNDARY_COND != BC_GENUS_TWO)) #define PERIODIC_BC ((BOUNDARY_COND == BC_PERIODIC)||(BOUNDARY_COND == BC_PERIODIC_CIRCLE)||(BOUNDARY_COND == BC_PERIODIC_FUNNEL)||(BOUNDARY_COND == BC_PERIODIC_TRIANGLE)) #define TWO_OBSTACLES ((SEGMENT_PATTERN == S_TWO_CIRCLES_EXT)||(SEGMENT_PATTERN == S_TWO_ROCKETS)) @@ -826,7 +848,8 @@ void evolve_segments(t_segment segment[NMAXSEGMENTS], int time) if (y < YMIN + padding) fy[group] += KSPRING_BOUNDARY*(YMIN + padding - y); else if (y > YMAX - padding) fy[group] -= KSPRING_BOUNDARY*(y - YMAX + padding); } - else if (BOUNDARY_COND == BC_REFLECT_ABS) /* add force from simulation boundary */ + else if ((BOUNDARY_COND == BC_REFLECT_ABS)||(BOUNDARY_COND == BC_REFLECT_ABS_BOTTOM)) + /* add force from simulation boundary */ { y = 0.5*(segment[i].y1 + segment[i].y2); if (y < YMIN) fy[group] += KSPRING_BOUNDARY*(YMIN - y); @@ -899,7 +922,7 @@ void evolve_segment_groups(t_segment segment[NMAXSEGMENTS], int time, t_group_se /* new version of evolve_segments that takes the group structure into account */ { double fx[NMAXGROUPS], fy[NMAXGROUPS], torque[NMAXGROUPS], dx[NMAXGROUPS], dy[NMAXGROUPS], dalpha[NMAXGROUPS]; - double x, y, dx0, dy0, padding, proj, distance, f, xx[2], yy[2], xmean = 0.0, ymean = 0.0; + double x, y, dx0, dy0, padding, proj, distance, f, xx[2], yy[2], xmean = 0.0, ymean = 0.0, ymax = 0.0; int i, j, k, group = 0; static double maxdepth, saturation_depth, xmax; static int first = 1; @@ -944,7 +967,8 @@ void evolve_segment_groups(t_segment segment[NMAXSEGMENTS], int time, t_group_se if (y < YMIN + padding) fy[group] += KSPRING_BOUNDARY*(YMIN + padding - y); else if (y > YMAX - padding) fy[group] -= KSPRING_BOUNDARY*(y - YMAX + padding); } - else if (BOUNDARY_COND == BC_REFLECT_ABS) /* add force from simulation boundary */ + else if ((BOUNDARY_COND == BC_REFLECT_ABS)||(BOUNDARY_COND == BC_REFLECT_ABS_BOTTOM)) + /* add force from simulation boundary */ { y = 0.5*(segment[i].y1 + segment[i].y2); if (y < YMIN) fy[group] += KSPRING_BOUNDARY*(YMIN - y); @@ -1047,10 +1071,14 @@ void evolve_segment_groups(t_segment segment[NMAXSEGMENTS], int time, t_group_se { xmean += segment_group[group].xc; ymean += segment_group[group].yc; + if (segment_group[group].yc > ymax) ymax = segment_group[group].yc; } xmean = xmean/((double)(ngroups-1)); ymean = ymean/((double)(ngroups-1)); + /* bias towards ymax */ + ymean = 0.75*ymax + 0.25*ymean; + if (ymean > ytrack) ytrack = ymean; if (xmean > xmax) xtrack = xmean - xmax; @@ -1063,7 +1091,7 @@ void evolve_segment_groups(t_segment segment[NMAXSEGMENTS], int time, t_group_se void animation() { double time, scale, diss, rgb[3], dissip, gradient[2], x, y, dx, dy, dt, xleft, xright, - a, b, length, fx, fy, force[2], totalenergy = 0.0, krepel = KREPEL, pos[2], prop, vx, beta = BETA, xi = 0.0, xmincontainer = BCXMIN, xmaxcontainer = BCXMAX, torque, torque_ij, fboundary = 0.0, pleft = 0.0, pright = 0.0, entropy[2], mean_energy, gravity = GRAVITY, speed_ratio, ymin, ymax, delta_energy, speed, ratio = 1.0, ratioc; + a, b, length, fx, fy, force[2], totalenergy = 0.0, krepel = KREPEL, pos[2], prop, vx, beta = BETA, xi = 0.0, xmincontainer = BCXMIN, xmaxcontainer = BCXMAX, torque, torque_ij, fboundary = 0.0, pleft = 0.0, pright = 0.0, entropy[2], mean_energy, gravity = GRAVITY, speed_ratio, ymin, ymax, delta_energy, speed, ratio = 1.0, ratioc, cum_etot = 0.0, emean = 0.0; double *qx, *qy, *px, *py, *qangle, *pangle, *pressure, *obstacle_speeds; int i, j, k, n, m, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q, p1, q1, p2, q2, total_neighbours = 0, min_nb, max_nb, close, wrapx = 0, wrapy = 0, nactive = 0, nadd_particle = 0, nmove = 0, nsuccess = 0, @@ -1117,13 +1145,6 @@ void animation() lj_final_position = fopen("lj_final_position.dat", "w"); } - /* initialise positions and radii of circles */ - init_particle_config(particle); - init_hashgrid(hashgrid); - - xshift = OBSTACLE_XMIN; - speed_ratio = (double)(25*NVID)*DT_PARTICLE; - if (ADD_FIXED_OBSTACLES) init_obstacle_config(obstacle); if (ADD_FIXED_SEGMENTS) init_segment_config(segment); @@ -1133,6 +1154,17 @@ void animation() group_speeds = (t_group_data *)malloc(ngroups*(INITIAL_TIME + NSTEPS)*sizeof(t_group_data)); } + /* initialise positions and radii of circles */ + init_particle_config(particle); + + /* add some particles, beta */ + if (ADD_INITIAL_PARTICLES) add_particle_config(particle, -0.6, 1.6, -1.0, 1.0, MU_B); + + init_hashgrid(hashgrid); + + xshift = OBSTACLE_XMIN; + speed_ratio = (double)(25*NVID)*DT_PARTICLE; + if (RECORD_PRESSURES) for (i=0; i 0.0)&&(speed < VICSEK_VMIN)) speed = VICSEK_VMIN; + if ((INTERACTION == I_VICSEK_SHARK)&&(particle[j].type == 2)) speed *= 1.75; if (speed > VICSEK_VMAX) speed = 0.5*(speed + VICSEK_VMAX); particle[j].vx = speed*cos(particle[j].angle); particle[j].vy = speed*sin(particle[j].angle); speed = module2(px[j],py[j]); if ((VICSEK_VMIN > 0.0)&&(speed < VICSEK_VMIN)) speed = VICSEK_VMIN; + if ((INTERACTION == I_VICSEK_SHARK)&&(particle[j].type == 2)) speed *= 1.75; if (speed > VICSEK_VMAX) speed = 0.5*(speed + VICSEK_VMAX); px[j] = speed*cos(particle[j].angle); py[j] = speed*sin(particle[j].angle); + + } /* add gravity */ @@ -1284,6 +1320,7 @@ void animation() /* timestep of thermostat algorithm */ totalenergy = evolve_particles(particle, hashgrid, qx, qy, qangle, px, py, pangle, beta, &nactive, &nsuccess, &nmove, i < INITIAL_TIME); + /* evolution of lid coordinate */ if (BOUNDARY_COND == BC_RECTANGLE_LID) evolve_lid(fboundary); if (BOUNDARY_COND == BC_RECTANGLE_WALL) @@ -1378,6 +1415,25 @@ void animation() } printf("Mean kinetic energy: %.3f\n", totalenergy/(double)ncircles); + + if ((!THERMOSTAT)&&(LIMIT_ENERGY)) + { + if (cum_etot > 0.0) + { + emean = cum_etot/(double)(i+1); + if (totalenergy > 10.0*emean) + { + reset_energy(particle, px, py, totalenergy, emean); + totalenergy = 0.0; + for (j=0; j= INITIAL_TIME) { - if (TIME_LAPSE_FIRST) + if ((TIME_LAPSE)&&(TIME_LAPSE_FIRST)) { if ((TIME_LAPSE)&&((i - INITIAL_TIME)%TIME_LAPSE_FACTOR == 0)&&(!DOUBLE_MOVIE)) { diff --git a/particle_billiard.c b/particle_billiard.c index dd53478..7f109d0 100644 --- a/particle_billiard.c +++ b/particle_billiard.c @@ -35,15 +35,18 @@ #include #include -#define MOVIE 0 /* set to 1 to generate movie */ +#define MOVIE 1 /* set to 1 to generate movie */ #define SAVE_MEMORY 1 /* set to 1 to save memory when writing tiff images */ #define INVERT_COUNTER 0 /* set to 1 to save frames in inverse order */ -#define WINWIDTH 1280 /* window width */ +// #define WINWIDTH 1280 /* window width */ +#define WINWIDTH 720 /* window width */ #define WINHEIGHT 720 /* window height */ -#define XMIN -1.5 -#define XMAX 2.5 /* x interval */ +// #define XMIN -1.5 +// #define XMAX 2.5 /* x interval */ +#define XMIN -1.125 +#define XMAX 1.125 /* x interval */ #define YMIN -1.125 #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ @@ -51,7 +54,7 @@ /* Choice of the billiard table, see global_particles.c */ -#define B_DOMAIN 31 /* choice of domain shape */ +#define B_DOMAIN 30 /* choice of domain shape */ #define CIRCLE_PATTERN 1 /* pattern of circles */ #define POLYLINE_PATTERN 10 /* pattern of polyline */ @@ -83,7 +86,7 @@ /* Simulation parameters */ // #define NPART 10 /* number of particles */ -#define NPART 5000 /* number of particles */ +#define NPART 50000 /* number of particles */ #define NPARTMAX 100000 /* maximal number of particles after resampling */ #define LMAX 0.01 /* minimal segment length triggering resampling */ #define DMIN 0.02 /* minimal distance to boundary for triggering resampling */ @@ -91,10 +94,10 @@ #define SHOWTRAILS 0 /* set to 1 to keep trails of the particles */ #define HEATMAP 1 /* set to 1 to show heat map of particles */ #define DRAW_FINAL_HEATMAP 1 /* set to 1 to show final heat map of particles */ -#define DRAW_HEATMAP_HISTOGRAM 1 /* set to 1 to draw a histogram of particle distribution in heat map */ -#define NBIN_FACTOR 8.0 /* constant controlling number of bins in histogram */ +#define DRAW_HEATMAP_HISTOGRAM 0 /* set to 1 to draw a histogram of particle distribution in heat map */ +#define NBIN_FACTOR 6.0 /* constant controlling number of bins in histogram */ #define DRAW_HEATMAP_PARTICLES 1 /* set to 1 to draw particles in heat map */ -#define HEATMAP_MAX_PART_BY_CELL 5 /* set to positive value to draw only limited number of particles in cell */ +#define HEATMAP_MAX_PART_BY_CELL 50 /* set to positive value to draw only limited number of particles in cell */ #define PLOT_HEATMAP_AVERAGE 1 /* set to 1 to plot average number of particles in heat map */ #define SHOWZOOM 0 /* set to 1 to show zoom on specific area */ #define PRINT_PARTICLE_NUMBER 0 /* set to 1 to print number of particles */ @@ -105,8 +108,8 @@ #define TEST_INITIAL_COND 0 /* set to 1 to allow only initial conditions that pass a test */ -#define NSTEPS 22000 /* number of frames of movie */ -#define TIME 2000 /* time between movie frames, for fluidity of real-time simulation */ +#define NSTEPS 1300 /* number of frames of movie */ +#define TIME 3000 /* time between movie frames, for fluidity of real-time simulation */ // #define DPHI 0.000002 /* integration step */ #define DPHI 0.00002 /* integration step */ #define NVID 25 /* number of iterations between images displayed on screen */ @@ -120,12 +123,12 @@ /* Colors and other graphical parameters */ -#define COLOR_PALETTE 17 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */ #define NCOLORS 500 /* number of colors */ #define COLORSHIFT 0 /* hue of initial color */ #define COLOR_HUEMIN 0 /* minimal color hue */ -#define COLOR_HUEMAX 160 /* maximal color hue */ +#define COLOR_HUEMAX 300 /* maximal color hue */ #define RAINBOW_COLOR 1 /* set to 1 to use different colors for all particles */ #define FLOWER_COLOR 0 /* set to 1 to adapt initial colors to flower billiard (tracks vs core) */ #define NSEG 100 /* number of segments of boundary */ @@ -145,8 +148,8 @@ #define SLEEP1 1 /* initial sleeping time */ #define SLEEP2 1 /* final sleeping time */ -#define NXMAZE 36 /* width of maze */ -#define NYMAZE 36 /* height of maze */ +#define NXMAZE 18 /* width of maze */ +#define NYMAZE 18 /* height of maze */ #define MAZE_MAX_NGBH 8 /* max number of neighbours of maze cell */ #define RAND_SHIFT 15 /* seed of random number generator */ #define MAZE_XSHIFT 0.0 /* horizontal shift of maze */ diff --git a/rde.c b/rde.c index 1054bb9..91c05ca 100644 --- a/rde.c +++ b/rde.c @@ -48,14 +48,21 @@ #define WINWIDTH 1920 /* window width */ #define WINHEIGHT 1150 /* window height */ -#define NX 960 /* number of grid points on x axis */ -#define NY 575 /* number of grid points on y axis */ +// #define NY 575 /* number of grid points on y axis */ +// #define NX 960 /* number of grid points on x axis */ +// #define NY 575 /* number of grid points on y axis */ // #define WINWIDTH 1280 /* window width */ // #define WINHEIGHT 720 /* window height */ -// #define NX 640 /* number of grid points on x axis */ -// #define NY 360 /* 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 */ +// #define NX 320 /* number of grid points on x axis */ +// #define NY 180 /* number of grid points on y axis */ +// #define XMIN -1.2 +// #define XMAX 1.2 /* x interval */ +// #define YMIN -1.2 +// #define YMAX 1.2 /* y interval for 9/16 aspect ratio */ #define XMIN -2.0 #define XMAX 2.0 /* x interval */ #define YMIN -1.197916667 @@ -63,7 +70,7 @@ /* Choice of simulated equation */ -#define RDE_EQUATION 7 /* choice of reaction term, see list in global_3d.c */ +#define RDE_EQUATION 8 /* choice of reaction term, see list in global_3d.c */ #define NFIELDS 3 /* number of fields in reaction-diffusion equation */ #define NLAPLACIANS 0 /* number of fields for which to compute Laplacian */ @@ -72,13 +79,16 @@ #define ADD_FORCE_FIELD 0 /* set to 1 to add a foce field (for compressible Euler equation) */ #define POTENTIAL 7 /* type of potential or vector potential, see list in global_3d.c */ #define FORCE_FIELD 5 /* type of force field, see list in global_3d.c */ -#define ADD_CORIOLIS_FORCE 1 /* set to 1 to add Coriolis force (quasigeostrophic Euler equations) */ +#define ADD_CORIOLIS_FORCE 0 /* set to 1 to add Coriolis force (quasigeostrophic Euler equations) */ +#define VARIABLE_DEPTH 1 /* set to 1 for variable depth in shallow water equation */ +#define SWATER_DEPTH 4 /* variable depth in shallow water equation */ #define ANTISYMMETRIZE_WAVE_FCT 0 /* set tot 1 to make wave function antisymmetric */ #define ADAPT_STATE_TO_BC 0 /* to smoothly adapt initial state to obstacles */ -#define OBSTACLE_GEOMETRY 72 /* geometry of obstacles, as in B_DOMAIN */ -#define BC_STIFFNESS 50.0 /* controls region of boundary condition control */ +#define OBSTACLE_GEOMETRY 7 /* geometry of obstacles, as in B_DOMAIN */ // #define BC_STIFFNESS 100.0 /* controls region of boundary condition control */ +#define BC_STIFFNESS 50.0 /* controls region of boundary condition control */ +#define CHECK_INTEGRAL 1 /* set to 1 to check integral of first field */ #define JULIA_SCALE 0.5 /* scaling for Julia sets */ @@ -86,14 +96,14 @@ #define B_DOMAIN 999 /* choice of domain shape, see list in global_pdes.c */ -#define CIRCLE_PATTERN 99 /* pattern of circles, see list in global_pdes.c */ +#define CIRCLE_PATTERN 8 /* 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 0.6 /* parameter controlling the dimensions of domain */ -#define MU 0.08 /* parameter controlling the dimensions of domain */ +#define LAMBDA 0.5 /* parameter controlling the dimensions of domain */ +#define MU 0.06 /* parameter controlling the dimensions of domain */ #define NPOLY 5 /* number of sides of polygon */ #define APOLY 2.0 /* angle by which to turn polygon, in units of Pi/2 */ #define MDEPTH 7 /* depth of computation of Menger gasket */ @@ -101,8 +111,8 @@ #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 NGRIDX 6 /* number of grid point for grid of disks */ +#define NGRIDY 8 /* number of grid point for grid of disks */ #define REVERSE_TESLA_VALVE 1 /* set to 1 to orient Tesla valve in blocking configuration */ #define X_SHOOTER -0.2 @@ -123,7 +133,12 @@ #define DT 0.00000025 -#define VISCOSITY 2.0 +// #define VISCOSITY 0.0001 +#define VISCOSITY 1.5e-5 +// #define VISCOSITY 5.0e-4 +// #define VISCOSITY 1.0e-3 +#define DISSIPATION 1.0e-8 +// #define DISSIPATION 1.0e-7 #define RPSA 0.75 /* parameter in Rock-Paper-Scissors-type interaction */ #define RPSLZB 0.75 /* second parameter in Rock-Paper-Scissors-Lizard-Spock type interaction */ @@ -138,9 +153,8 @@ #define BZQ 0.0008 /* parameter in BZ equation */ #define BZF 1.2 /* parameter in BZ equation */ #define B_FIELD 10.0 /* magnetic field */ -#define G_FIELD 0.004 /* gravity/constant in repulsive field from obstacles */ -// #define G_FIELD 1.0e-4 /* gravity/constant in repulsive field from obstacles */ -// #define G_FIELD 2.0e-5 /* gravity/constant in repulsive field from obstacles */ +#define G_FIELD 0.75 /* gravity */ +#define BC_FIELD 1.0e-5 /* constant in repulsive field from obstacles */ #define AB_RADIUS 0.2 /* radius of region with magnetic field for Aharonov-Bohm effect */ #define K_EULER 50.0 /* constant in stream function integration of Euler equation */ #define K_EULER_INC 0.5 /* constant in incompressible Euler equation */ @@ -181,11 +195,22 @@ #define IN_OUT_BC_FACTOR 0.001 /* factor of convex combination between old and new flow */ #define BC_FLOW_TYPE 1 /* type of initial condition */ /* see list in global_pdes.c */ -#define IN_OUT_FLOW_MIN_AMP 0.45 /* amplitude of in-flow/out-flow boundary conditions (for Euler equation) - min value */ -#define IN_OUT_FLOW_AMP 0.45 /* amplitude of in-flow/out-flow boundary conditions (for Euler equation) - max value */ +#define IN_OUT_FLOW_MIN_AMP 0.25 /* amplitude of in-flow/out-flow boundary conditions (for Euler equation) - min value */ +#define IN_OUT_FLOW_AMP 0.25 /* amplitude of in-flow/out-flow boundary conditions (for Euler equation) - max value */ #define LAMINAR_FLOW_MODULATION 0.01 /* asymmetry of laminar flow */ #define LAMINAR_FLOW_YPERIOD 1.0 /* period of laminar flow in y direction */ -#define PRESSURE_GRADIENT 0.3 /* amplitude of pressure gradient for Euler equation */ +#define PRESSURE_GRADIENT 0.2 /* amplitude of pressure gradient for Euler equation */ + +#define SWATER_MIN_HEIGHT 2.0 /* min height of initial condition for shallow water equation */ +// #define DEPTH_FACTOR 0.0125 /* proportion of min height in variable depth */ +// #define DEPTH_FACTOR 0.001 /* proportion of min height in variable depth */ +// #define DEPTH_FACTOR 0.0012 /* proportion of min height in variable depth */ +// #define DEPTH_FACTOR 0.0015 /* proportion of min height in variable depth */ +// #define DEPTH_FACTOR 0.002 /* proportion of min height in variable depth */ +// #define DEPTH_FACTOR 0.005 /* proportion of min height in variable depth */ +// #define DEPTH_FACTOR 0.01 /* proportion of min height in variable depth */ +#define DEPTH_FACTOR 0.015 /* proportion of min height in variable depth */ +#define TANH_FACTOR 1.0 /* steepness of variable depth */ #define EULER_GRADIENT_YSHIFT 0.0 /* y-shift in computation of gradient in Euler equation */ @@ -193,15 +218,22 @@ #define B_COND 1 +#define B_COND_LEFT 0 +#define B_COND_RIGHT 0 +#define B_COND_TOP 1 +#define B_COND_BOTTOM 1 + + /* Parameters for length and speed of simulation */ -#define NSTEPS 2250 /* number of frames of movie */ -// #define NSTEPS 500 /* number of frames of movie */ -#define NVID 100 /* number of iterations between images displayed on screen */ +#define NSTEPS 2100 /* number of frames of movie */ +// #define NSTEPS 500 /* number of frames of movie */ +// #define NVID 100 /* number of iterations between images displayed on screen */ +#define NVID 30 /* number of iterations between images displayed on screen */ #define ACCELERATION_FACTOR 1.0 /* factor by which to increase NVID in course of simulation */ #define DT_ACCELERATION_FACTOR 1.0 /* factor by which to increase time step in course of simulation */ #define MAX_DT 0.024 /* maximal value of integration step */ -#define NSEG 100 /* number of segments of boundary */ +#define NSEG 999 /* number of segments of boundary */ #define BOUNDARY_WIDTH 2 /* width of billiard boundary */ #define PAUSE 100 /* number of frames after which to pause */ @@ -215,7 +247,7 @@ /* Visualisation */ -#define PLOT_3D 0 /* controls whether plot is 2D or 3D */ +#define PLOT_3D 1 /* controls whether plot is 2D or 3D */ #define ROTATE_VIEW 0 /* set to 1 to rotate position of observer */ #define ROTATE_ANGLE 360.0 /* total angle of rotation during simulation */ @@ -224,22 +256,27 @@ /* Plot type - color scheme */ -#define CPLOT 64 -#define CPLOT_B 62 +#define CPLOT 70 +#define CPLOT_B 74 /* Plot type - height of 3D plot */ -#define ZPLOT 61 /* z coordinate in 3D plot */ -#define ZPLOT_B 64 /* z coordinate in second 3D plot */ +#define ZPLOT 70 /* z coordinate in 3D plot */ +#define ZPLOT_B 71 /* z coordinate in second 3D plot */ #define AMPLITUDE_HIGH_RES 1 /* set to 1 to increase resolution of P_3D_AMPLITUDE plot */ #define SHADE_3D 1 /* set to 1 to change luminosity according to normal vector */ #define NON_DIRICHLET_BC 0 /* set to 1 to draw only facets in domain, if field is not zero on boundary */ #define WRAP_ANGLE 1 /* experimental: wrap angle to [0, 2Pi) for interpolation in angle schemes */ #define FADE_IN_OBSTACLE 0 /* set to 1 to fade color inside obstacles */ +#define FADE_WATER_DEPTH 0 /* set to 1 to make wave color depth-dependent */ #define DRAW_OUTSIDE_GRAY 0 /* experimental - draw outside of billiard in gray */ -#define ADD_POTENTIAL_TO_Z 1 /* set to 1 to add the external potential to z-coordinate of plot */ -#define ADD_POT_CONSTANT 0.35 /* constant in front of added potential */ +#define ADD_POTENTIAL_TO_Z 0 /* set to 1 to add the external potential to z-coordinate of plot */ +#define ADD_POT_CONSTANT 0.35 /* constant added to wave height */ +#define ADD_HEIGHT_CONSTANT -2.0 /* constant added to wave height */ +#define DRAW_DEPTH 1 /* set to 1 to draw water depth */ +#define DEPTH_SCALE 0.75 /* vertical scaling of depth plot */ +#define DEPTH_SHIFT -0.015 /* vertical shift of depth plot */ #define PLOT_SCALE_ENERGY 0.05 /* vertical scaling in energy plot */ @@ -269,26 +306,28 @@ /* Color schemes, see list in global_pdes.c */ -#define COLOR_PALETTE 17 /* Color palette, see list in global_pdes.c */ -#define COLOR_PALETTE_B 10 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE_B 10 /* Color palette, see list in global_pdes.c */ +// #define COLOR_PALETTE_B 0 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* black background */ #define COLOR_SCHEME 3 /* choice of color scheme */ -#define COLOR_PHASE_SHIFT 0.0 /* phase shift of color scheme, in units of Pi */ +#define COLOR_PHASE_SHIFT 0.5 /* phase shift of color scheme, in units of Pi */ #define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ #define SLOPE 1.0 /* sensitivity of color on wave amplitude */ -#define VSCALE_AMPLITUDE 15.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ +#define VSCALE_AMPLITUDE 10.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ #define CURL_SCALE 0.000015 /* scaling factor for curl representation */ #define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */ -#define SLOPE_SCHROD_LUM 50.0 /* sensitivity of luminosity on module, for color scheme Z_ARGUMENT */ +#define SLOPE_SCHROD_LUM 200.0 /* sensitivity of luminosity on module, for color scheme Z_ARGUMENT */ #define MIN_SCHROD_LUM 0.2 /* minimal luminosity in color scheme Z_ARGUMENT*/ #define VSCALE_PRESSURE 2.0 /* additional scaling factor for color scheme Z_EULER_PRESSURE */ #define PRESSURE_SHIFT 10.0 /* shift for color scheme Z_EULER_PRESSURE */ #define PRESSURE_LOG_SHIFT -2.5 /* shift for color scheme Z_EULER_PRESSURE */ +#define VSCALE_WATER_HEIGHT 0.4 /* vertical scaling of water height */ #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 */ @@ -301,25 +340,28 @@ #define LOG_SCALE 0.5 /* scaling factor for energy log representation */ #define LOG_SHIFT 1.0 #define LOG_MIN 1.0e-3 /* floor value for log vorticity plot */ -#define VSCALE_SPEED 10.0 /* additional scaling factor for color scheme Z_EULER_SPEED */ +#define VSCALE_SPEED 200.0 /* additional scaling factor for color scheme Z_EULER_SPEED */ #define VMEAN_SPEED 0.0 /* mean value around which to scale for color scheme Z_EULER_SPEED */ -#define SHIFT_DENSITY 2.6 /* shift for color scheme Z_EULER_DENSITY */ -#define VSCALE_DENSITY 7.5 /* additional scaling factor for color scheme Z_EULER_DENSITY */ +#define SHIFT_DENSITY 8.5 /* shift for color scheme Z_EULER_DENSITY */ +#define VSCALE_DENSITY 3.0 /* additional scaling factor for color scheme Z_EULER_DENSITY */ #define VSCALE_VORTICITY 20.0 /* additional scaling factor for color scheme Z_EULERC_VORTICITY */ #define VORTICITY_SHIFT 0.0 /* vertical shift of vorticity */ -#define ZSCALE_SPEED 1.0 /* additional scaling factor for z-coord Z_EULER_SPEED */ +#define ZSCALE_SPEED 0.3 /* additional scaling factor for z-coord Z_EULER_SPEED and Z_SWATER_SPEED */ +#define VSCALE_SWATER 300.0 /* additional scaling factor for color scheme Z_EULER_DENSITY */ -#define NXMAZE 13 /* width of maze */ -#define NYMAZE 13 /* height of 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 3 /* seed of random number generator */ #define MAZE_XSHIFT 0.0 /* horizontal shift of maze */ -#define MAZE_WIDTH 0.03 /* half width of maze walls */ +#define MAZE_WIDTH 0.04 /* half width of maze walls */ #define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */ #define COLORBAR_RANGE 2.5 /* scale of color scheme bar */ #define COLORBAR_RANGE_B 2.5 /* scale of color scheme bar for 2nd part */ #define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ +#define CIRC_COLORBAR 0 /* set to 1 to draw circular color scheme */ +#define CIRC_COLORBAR_B 1 /* set to 1 to draw circular color scheme */ /* only for compatibility with wave_common.c */ #define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */ @@ -355,23 +397,23 @@ double u_3d[2] = {0.75, -0.45}; /* projections of basis vectors for REP_AXO_ double v_3d[2] = {-0.75, -0.45}; double w_3d[2] = {0.0, 0.015}; double light[3] = {0.816496581, 0.40824829, 0.40824829}; /* vector of "light" direction for P_3D_ANGLE color scheme */ -double observer[3] = {8.0, 8.0, 10.0}; /* location of observer for REP_PROJ_3D representation */ +double observer[3] = {8.0, 8.0, 7.0}; /* location of observer for REP_PROJ_3D representation */ int reset_view = 0; /* switch to reset 3D view parameters (for option ROTATE_VIEW) */ -#define Z_SCALING_FACTOR 1.25 /* overall scaling factor of z axis for REP_PROJ_3D representation */ +#define Z_SCALING_FACTOR 35.0 /* overall scaling factor of z axis for REP_PROJ_3D representation */ #define XY_SCALING_FACTOR 1.7 /* overall scaling factor for on-screen (x,y) coordinates after projection */ #define ZMAX_FACTOR 1.0 /* max value of z coordinate for REP_PROJ_3D representation */ #define XSHIFT_3D 0.0 /* overall x shift for REP_PROJ_3D representation */ -#define YSHIFT_3D -0.1 /* overall y shift for REP_PROJ_3D representation */ +#define YSHIFT_3D 0.1 /* overall y shift for REP_PROJ_3D representation */ #define BORDER_PADDING 0 /* distance from boundary at which to plot points, to avoid boundary effects due to gradient */ /* For debugging purposes only */ #define FLOOR 1 /* set to 1 to limit wave amplitude to VMAX */ -#define VMAX 1000.0 /* max value of wave amplitude */ +#define VMAX 100.0 /* max value of wave amplitude */ #define TEST_GRADIENT 0 /* print norm squared of gradient */ #define REFRESH_B (ZPLOT_B != ZPLOT)||(CPLOT_B != CPLOT) /* to save computing time, to be improved */ -#define COMPUTE_WRAP_ANGLE ((WRAP_ANGLE)&&((cplot == Z_ANGLE_GRADIENT)||(cplot == Z_ANGLE_GRADIENTX)||(cplot == Z_ARGUMENT)||(cplot == Z_ANGLE_GRADIENTX))) +#define COMPUTE_WRAP_ANGLE ((WRAP_ANGLE)&&((cplot == Z_ANGLE_GRADIENT)||(cplot == Z_ANGLE_GRADIENTX)||(cplot == Z_ARGUMENT)||(cplot == Z_ANGLE_GRADIENTX)||(cplot == Z_EULER_DIRECTION_SPEED)||(cplot == Z_SWATER_DIRECTION_SPEED))) #define PRINT_PARAMETERS ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)||(PRINT_PROBABILITIES)||(PRINT_NOISE)||(PRINT_FLOW_SPEED)||(PRINT_AVERAGE_SPEED)) #define COMPUTE_PRESSURE ((ZPLOT == Z_EULER_PRESSURE)||(CPLOT == Z_EULER_PRESSURE)||(ZPLOT_B == Z_EULER_PRESSURE)||(CPLOT_B == Z_EULER_PRESSURE)) @@ -519,16 +561,16 @@ void compute_gfield(int i, int j, double *gx, double *gy) { r = module2(x,y) + 1.0e-2; f = 0.5*(1.0 - tanh(BC_STIFFNESS*(r - LAMBDA))); - *gx = G_FIELD*f*x/r; - *gy = G_FIELD*f*y/r; + *gx = BC_FIELD*f*x/r; + *gy = BC_FIELD*f*y/r; break; } case (GF_ELLIPSE): { r = module2(x/LAMBDA,y/MU) + 1.0e-2; f = 0.5*(1.0 - tanh(BC_STIFFNESS*(r - 1.0))); - *gx = G_FIELD*f*x/(LAMBDA*LAMBDA); - *gy = G_FIELD*f*y/(MU*MU); + *gx = BC_FIELD*f*x/(LAMBDA*LAMBDA); + *gy = BC_FIELD*f*y/(MU*MU); break; } case (GF_AIRFOIL): @@ -536,8 +578,8 @@ void compute_gfield(int i, int j, double *gx, double *gy) y1 = y + a*x*x; r = module2(x/LAMBDA,y1/MU) + 1.0e-2; f = 0.5*(1.0 - tanh(BC_STIFFNESS*(r - 1.0))); - *gx = G_FIELD*f*(x/(LAMBDA*LAMBDA) + a*y1/(MU*MU)); - *gy = G_FIELD*f*y1/(MU*MU); + *gx = BC_FIELD*f*(x/(LAMBDA*LAMBDA) + a*y1/(MU*MU)); + *gy = BC_FIELD*f*y1/(MU*MU); break; } case (GF_WING): @@ -554,8 +596,8 @@ void compute_gfield(int i, int j, double *gx, double *gy) y1 = y + a*x*x; r = module2(x/LAMBDA,y1/(MU*x1)) + 1.0e-2; f = 0.5*(1.0 - tanh(BC_STIFFNESS*(r - 1.0))); - *gx = G_FIELD*f*(x/(LAMBDA*LAMBDA) + 2.0*a*x*y1/(MU*MU*x1*x1) - y1*y1/(MU*MU*x1*x1*x1)); - *gy = G_FIELD*f*y1/(MU*MU*x1*x1); + *gx = BC_FIELD*f*(x/(LAMBDA*LAMBDA) + 2.0*a*x*y1/(MU*MU*x1*x1) - y1*y1/(MU*MU*x1*x1*x1)); + *gy = BC_FIELD*f*y1/(MU*MU*x1*x1); // *gx = 0.1*G_FIELD*f*(x/(LAMBDA*LAMBDA) + 2.0*a*x*y1/(MU*MU*x1*x1) - y1*y1/(MU*MU*x1*x1*x1)); // *gy = 0.1*G_FIELD*f*y1/(MU*MU*x1*x1); // hx = x/(LAMBDA*LAMBDA) + 2.0*a*x*y1/(MU*MU*x1*x1) - y1*y1/(MU*MU*x1*x1*x1); @@ -613,8 +655,8 @@ void initialize_gfield(double gfield[2*NX*NY], double bc_field[NX*NY], double bc #pragma omp parallel for private(i,j) for (i=1; idepth = SWATER_MIN_HEIGHT*DEPTH_FACTOR*h; + rde->gradx = x*z; + rde->grady = y*z; + break; + } + case (SH_CIRCLES): + { + h = 0.0; + z = 0.0; + r = 10.0; + /* compute minimal distance to circles */ + for (n = 0; n < ncircles; n++) + for (nx = -1; nx < 2; nx++) + for (ny = -1; ny < 2; ny++) + { + x1 = circles[n].xc + (double)nx*(XMAX-XMIN); + y1 = circles[n].yc + (double)ny*(YMAX-YMIN); + r0 = module2(x - x1,y - y1)/circles[n].radius; + if (r0 < r) r = r0; + } + h = tanh(TANH_FACTOR*(r-1.0)); + h = 0.5*(h + 1.0); + rde->depth = SWATER_MIN_HEIGHT*DEPTH_FACTOR*h; + break; + } + case (SH_COAST): + { + x1 = PI*(x + 0.1 - 0.1*cos(PI*y/YMAX))/XMAX; + d = 3.0*sin(x1); + h = tanh(-TANH_FACTOR*d); + h = 0.5*(h + 1.0); + rde->depth = SWATER_MIN_HEIGHT*DEPTH_FACTOR*h; + break; + } + case (SH_COAST_MONOTONE): + { + x1 = PI*(x + 0.1 - 0.2*cos(PI*y/YMAX))/XMAX; + h = tanh(-TANH_FACTOR*x1); + h = 0.5*(h + 1.0); + rde->depth = SWATER_MIN_HEIGHT*DEPTH_FACTOR*h; + break; + } + } +} + +double initialize_water_depth(t_rde rde[NX*NY]) +{ + int i, j, ncircles; + double dx, dy, min, max, pscal, norm, vz = 0.01; + + if (SWATER_DEPTH == SH_CIRCLES) ncircles = init_circle_config_pattern(circles, CIRCLE_PATTERN); + + #pragma omp parallel for private(i,j) + for (i=0; i max) max = rde[i*NY + j].gradx; + if (rde[i*NY + j].gradx < min) min = rde[i*NY + j].gradx; + } + } + printf("gradx min = %.3lg, max = %.3lg\n", min, max); + + min = 0.0; + max = 0.0; + #pragma omp parallel for private(i,j) + for (i=0; i max) max = rde[i*NY + j].grady; + if (rde[i*NY + j].grady < min) min = rde[i*NY + j].grady; + } + } + printf("grady min = %.3lg, max = %.3lg\n", min, max); + + min = 0.0; + max = 0.0; + #pragma omp parallel for private(i,j) + for (i=0; i max) max = rde[i*NY + j].depth; + if (rde[i*NY + j].depth < min) min = rde[i*NY + j].depth; + } + } + printf("Depth min = %.3lg, max = %.3lg\n", min, max); + + +// for (j=0; j= SMOOTHEN_PERIOD) smooth = 0; } - smooth++; - if (smooth >= SMOOTHEN_PERIOD) smooth = 0; + break; } - } - /* compute gradients of fields for compressible Euler equation */ - else if (RDE_EQUATION == E_EULER_COMP) - { - nabla_rho = (double *)malloc(2*NX*NY*sizeof(double)); -// nabla_u = (double *)malloc(2*NX*NY*sizeof(double)); -// nabla_v = (double *)malloc(2*NX*NY*sizeof(double)); - compute_gradient_euler_test(phi_in[0], nabla_rho, xy_in); - compute_velocity_gradients(phi_in, rde); -// compute_gradient_euler_test(phi_in[1], nabla_u, xy_in); -// compute_gradient_euler_test(phi_in[2], nabla_v, xy_in); - - if (SMOOTHEN_VELOCITY) /* beta: try to reduce formation of ripples */ + /* compute gradients of fields for compressible Euler equation */ + case (E_EULER_COMP): { - if (smooth == 0) + nabla_rho = (double *)malloc(2*NX*NY*sizeof(double)); + compute_gradient_euler_test(phi_in[0], nabla_rho, xy_in); + compute_velocity_gradients(phi_in, rde, xy_in); + + if (SMOOTHEN_VELOCITY) /* beta: try to reduce formation of ripples */ + { + if (smooth == 0) + { + delta_u = (double *)malloc(NX*NY*sizeof(double)); + delta_v = (double *)malloc(NX*NY*sizeof(double)); + compute_laplacian_rde(phi_in[1], delta_u, xy_in); + compute_laplacian_rde(phi_in[2], delta_v, xy_in); + #pragma omp parallel for private(i) + for (i=0; i= SMOOTHEN_PERIOD) smooth = 0; + } + break; + } + + case (E_SHALLOW_WATER): + { + nabla_eta = (double *)malloc(2*NX*NY*sizeof(double)); + compute_gradient_euler_test(phi_in[0], nabla_eta, xy_in); + compute_velocity_gradients(phi_in, rde, xy_in); + + if (VISCOSITY > 0.0) { delta_u = (double *)malloc(NX*NY*sizeof(double)); delta_v = (double *)malloc(NX*NY*sizeof(double)); compute_laplacian_rde(phi_in[1], delta_u, xy_in); compute_laplacian_rde(phi_in[2], delta_v, xy_in); - #pragma omp parallel for private(i) - for (i=0; i= SMOOTHEN_PERIOD) smooth = 0; + break; } - } + } if (TEST_GRADIENT) { test = 0.0; @@ -899,6 +1175,48 @@ void evolve_wave_half(double *phi_in[NFIELDS], double *phi_out[NFIELDS], short i } break; } + case (E_SHALLOW_WATER): + { + eta = phi_in[0][i*NY+j]; + if (eta == 0.0) eta = 1.0e-6; + u = phi_in[1][i*NY+j]; + v = phi_in[2][i*NY+j]; + etax = nabla_eta[i*NY+j]; + etay = nabla_eta[NX*NY+i*NY+j]; + + ux = rde[i*NY+j].dxu; + uy = rde[i*NY+j].dyu; + vx = rde[i*NY+j].dxv; + vy = rde[i*NY+j].dyv; + + phi_out[0][i*NY+j] = eta - intstep*(u*etax + v*etay + eta*(ux + vy)); + phi_out[1][i*NY+j] = u - intstep*(u*ux + v*uy + G_FIELD*etax); + phi_out[2][i*NY+j] = v - intstep*(u*vx + v*vy + G_FIELD*etay); + + if (VISCOSITY > 0.0) + { + phi_out[1][i*NY+j] += intstep*VISCOSITY*delta_u[i*NY+j]; + phi_out[2][i*NY+j] += intstep*VISCOSITY*delta_v[i*NY+j]; + } + if (DISSIPATION > 0.0) + { + phi_out[1][i*NY+j] -= intstep*DISSIPATION*u; + phi_out[2][i*NY+j] -= intstep*DISSIPATION*v; + } + if (ADD_FORCE_FIELD) + { + phi_out[1][i*NY+j] += intstep*gfield[i*NY+j]; + phi_out[2][i*NY+j] += intstep*gfield[NX*NY+i*NY+j]; + } + if (VARIABLE_DEPTH) + { + phi_out[0][i*NY+j] -= intstep*rde[i*NY+j].depth*(ux + vy); + phi_out[0][i*NY+j] -= intstep*rde[i*NY+j].gradx*u; + phi_out[0][i*NY+j] -= intstep*rde[i*NY+j].grady*v; + } + + break; + } } } } @@ -942,6 +1260,15 @@ void evolve_wave_half(double *phi_in[NFIELDS], double *phi_out[NFIELDS], short i { free(nabla_rho); } + else if (RDE_EQUATION == E_SHALLOW_WATER) + { + free(nabla_eta); + if (VISCOSITY > 0.0) + { + free(delta_u); + free(delta_v); + } + } if (COMPUTE_PRESSURE) { @@ -1142,7 +1469,7 @@ void print_parameters(double *phi[NFIELDS], t_rde rde[NX*NY], short int xy_in[NX } } -void draw_color_bar_palette(int plot, double range, int palette, int fade, double fade_value) +void draw_color_bar_palette(int plot, double range, int palette, int circular, int fade, double fade_value) { double width = 0.14; // double width = 0.2; @@ -1151,7 +1478,8 @@ void draw_color_bar_palette(int plot, double range, int palette, int fade, doubl { if (ROTATE_COLOR_SCHEME) draw_color_scheme_palette_3d(XMIN + 0.3, YMIN + 0.1, XMAX - 0.3, YMIN + 0.1 + width, plot, -range, range, palette, fade, fade_value); -// draw_color_scheme_palette_3d(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette, fade, fade_value); + else if (circular) + draw_circular_color_scheme_palette_3d(XMAX - 2.0*width, YMAX - 2.0*width, 1.5*width, plot, -range, range, palette, fade, fade_value); else draw_color_scheme_palette_3d(XMAX - 1.5*width, YMIN + 0.1, XMAX - 0.5*width, YMAX - 0.1, plot, -range, range, palette, fade, fade_value); } @@ -1159,7 +1487,6 @@ void draw_color_bar_palette(int plot, double range, int palette, int fade, doubl { if (ROTATE_COLOR_SCHEME) draw_color_scheme_palette_fade(XMIN + 0.8, YMIN + 0.1, XMAX - 0.8, YMIN + 0.1 + width, plot, -range, range, palette, fade, fade_value); -// draw_color_scheme_palette_fade(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette, fade, fade_value); else draw_color_scheme_palette_fade(XMAX - 1.5*width, YMIN + 0.1, XMAX - 0.5*width, YMAX - 0.1, plot, -range, range, palette, fade, fade_value); } @@ -1238,12 +1565,13 @@ void viewpoint_schedule(int i) void animation() { double time = 0.0, scale, dx, var, jangle, cosj, sinj, sqrintstep, - intstep0, viscosity_printed, fade_value, noise = NOISE_INTENSITY; + intstep0, viscosity_printed, fade_value, noise = NOISE_INTENSITY, x, y, sign; double *phi[NFIELDS], *phi_tmp[NFIELDS], *potential_field, *vector_potential_field, *tracers, *gfield, *bc_field, *bc_field2; short int *xy_in; int i, j, k, s, nvid, field; static int counter = 0; t_rde *rde; +// t_swater_depth *water_depth; /* Since NX and NY are big, it seemed wiser to use some memory allocation here */ for (i=0; i ymin - radius)&&(particles[n].xc < xmax + radius)&&(particles[n].xc > xmin - radius)) particles[n].active = 1; + else particles[n].active = 0; + } + ncircles += NGRIDX*(NGRIDY+1); + break; + } + case (C_POISSON_DISC): + { + ncirc0 = ncircles; + printf("Generating new Poisson disc sample\n"); + /* generate first circle */ + particles[ncirc0].xc = (xmax - xmin)*(double)rand()/RAND_MAX + xmin; + particles[ncirc0].yc = (ymax - ymin)*(double)rand()/RAND_MAX + ymin; + active_poisson[0] = 1; +// // particles[0].active = 1; + n_p_active = 1; + newcircles = 1; + + while ((n_p_active > 0)&&(ncircles < NMAXCIRCLES)) + { + /* randomly select an active circle */ + i = rand()%(newcircles); + while (!active_poisson[i]) i = rand()%(ncircles); +// printf("Starting from circle %i at (%.3f,%.3f)\n", i, particles[i].xc, particles[i].yc); + /* generate new candidates */ + naccepted = 0; + for (j=0; j= dpoisson*dpoisson); + /* new circle is in domain */ + far = far*(x < xmax)*(x > xmin)*(y < ymax)*(y > ymin); +// far = far*(vabs(x) < LAMBDA)*(y < INITYMAX)*(y > INITYMIN); + } + if (far) /* accept new circle */ + { + printf("New circle at (%.3f,%.3f) accepted\n", x, y); + particles[ncircles].xc = x; + particles[ncircles].yc = y; + particles[ncircles].radius = radius; + particles[ncircles].active = 1; + active_poisson[ncircles] = 1; + ncircles++; + ncirc0++; + n_p_active++; + naccepted++; + } +// else printf("Rejected\n"); + } + if (naccepted == 0) /* inactivate circle i */ + { +// printf("No candidates work, inactivate circle %i\n", i); + active_poisson[i] = 0; + n_p_active--; + } + printf("%i active circles\n", n_p_active); + } + + printf("Generated %i circles\n", ncircles); + break; + } + default: + { + printf("Function init_circle_config not defined for this pattern \n"); + } + } +} + void init_people_config(t_person people[NMAXCIRCLES]) /* initialise particle configuration */ { @@ -1021,6 +1140,17 @@ void init_obstacle_config(t_obstacle obstacle[NMAXOBSTACLES]) add_obstacle(0.5*(double)i, YMIN + 0.3, radius, obstacle); break; } + case (O_CIRCLE): + { + n = 0; + obstacle[n].xc = 0.0; + obstacle[n].yc = 0.0; + obstacle[n].radius = OBSTACLE_RADIUS; + obstacle[n].active = 1; + nobstacles = 1; + break; + } + default: { printf("Function init_obstacle_config not defined for this pattern \n"); @@ -1162,7 +1292,7 @@ void add_circle_to_segments(double x, double y, double r, int nsegs, double angl double nozzle_width(double x, double width, int nozzle_shape) /* width of bell-shaped nozzle */ { - double lam = 0.5*LAMBDA; + double lam = 0.5*LAMBDA, a, b; if (x >= 0.0) return(width); else switch (nozzle_shape) { @@ -1176,6 +1306,16 @@ double nozzle_width(double x, double width, int nozzle_shape) if (-x < 0.1) return(width - (0.5 - width)*x/0.1); else return(0.5); } + case (NZ_DELAVAL): + { + a = 1.5; + b = 0.05; + return(sqrt(width*width - 0.5*x) + a*b*x*(1.0 + x)/(b + x*x)); +// a = (sqrt(width*width+0.5) - width)/sqrt(0.5); +// c = (a*sqrt(0.5*h) - width)/(h*h); +// if (-x < h) return(width + c*x*x); +// else return(width + a * sqrt(-0.5*x)); + } default: return(0.0); } } @@ -1185,7 +1325,7 @@ void add_rocket_to_segments(t_segment segment[NMAXSEGMENTS], double x0, double y /* add one or several rocket_shaped set of segments */ { int i, j, cycle = 0, nsegments0; - double angle, dx, x1, y1, x2, y2, nozx, nozy, a, b; + double angle, dx, x1, y1, x2, y2, nozx, nozy, a, b, c, w; nsegments0 = nsegments; @@ -1234,6 +1374,25 @@ void add_rocket_to_segments(t_segment segment[NMAXSEGMENTS], double x0, double y add_rotated_angle_to_segments(x0-a, nozy, x0-nozx, nozy, 0.02, segment, 0); break; } + case (RCK_RECT_BAR): /* rectangular chamber with a hat and separating bar */ + { + a = 0.5*LAMBDA; + b = (0.49*PI-0.25)*LAMBDA; + c = 0.5*a; + w = 0.025; + add_rotated_angle_to_segments(x0+nozx, nozy, x0+nozx, nozy+0.5*c, 0.02, segment, 0); + add_rotated_angle_to_segments(x0+nozx, nozy+0.5*c, x0+nozx+0.5*c, nozy+c, 0.02, segment, 0); + add_rotated_angle_to_segments(x0+nozx+0.5*c, nozy+c, x0+a, nozy+c, 0.02, segment, 0); + add_rotated_angle_to_segments(x0+a, nozy+c, x0+a, nozy+b+c, 0.02, segment, 0); + add_rotated_angle_to_segments(x0+a, nozy+b+c, x0, nozy+b+a+c, 0.02, segment, 0); + add_rotated_angle_to_segments(x0-w, nozy+a+c, x0-w, nozy+b+a+c, 2.0*w, segment, 0); + add_rotated_angle_to_segments(x0, nozy+b+a+c, x0-a, nozy+b+c, 0.02, segment, 0); + add_rotated_angle_to_segments(x0-a, nozy+b+c, x0-a, nozy+c, 0.02, segment, 0); + add_rotated_angle_to_segments(x0-a, nozy+c, x0-nozx-0.5*c, nozy+c, 0.02, segment, 0); + add_rotated_angle_to_segments(x0-nozx-0.5*c, nozy+c, x0-nozx+w, nozy+0.5*c, 0.02, segment, 0); + add_rotated_angle_to_segments(x0-nozx+w, nozy+0.5*c, x0-nozx+w, nozy, 0.02, segment, 0); + break; + } } dx = LAMBDA/(double)(nsides); @@ -1275,7 +1434,7 @@ void add_rocket_to_segments(t_segment segment[NMAXSEGMENTS], double x0, double y for (i=nsegments0; i0)||(j!=NYMAZE/2))&&(maze[n].west)) add_rectangle_to_segments(x1, y1, x1 - width, y1 + dy, segment, 0); + if (diag) + { + if (((i>0)||(j0)||(j!=NYMAZE/2))&&(maze[n].west)) add_rectangle_to_segments(x1, y1, x1 - width, y1 + dy, segment, 0); if (maze[n].south) add_rectangle_to_segments(x1, y1, x1 + dx, y1 - width, segment, 0); } @@ -1305,15 +1468,29 @@ int init_maze_segments(t_segment segment[NMAXSEGMENTS]) add_rectangle_to_segments(YMIN + padding + MAZE_XSHIFT, YMAX - padding, YMAX - padding + MAZE_XSHIFT, YMAX - padding - width, segment, 0); /* right side of maze */ - y1 = YMIN + padding + dy*((double)NYMAZE/2); x1 = YMAX - padding + MAZE_XSHIFT; - add_rectangle_to_segments(x1, YMIN - 1.0, x1 - width, y1 - dy, segment, 0); - add_rectangle_to_segments(x1, y1, x1 - width, YMAX + 1.0, segment, 0); + if (diag) + { + y1 = YMIN + padding + dy; + add_rectangle_to_segments(x1, y1, x1 - width, YMAX + 10.0, segment, 0); + } + else + { + y1 = YMIN + padding + dy*((double)NYMAZE/2); + add_rectangle_to_segments(x1, YMIN - 1.0, x1 - width, y1 - dy, segment, 0); + add_rectangle_to_segments(x1, y1, x1 - width, YMAX + 1.0, segment, 0); + } /* left side of maze */ x1 = YMIN + padding + MAZE_XSHIFT; add_rectangle_to_segments(x1, YMIN - 1.0, x1 - width, YMIN + padding, segment, 0); - add_rectangle_to_segments(x1, YMAX - padding, x1 - width, YMAX + 1.0, segment, 0); + add_rectangle_to_segments(x1, YMAX - padding, x1 - width, YMAX + 10.0, segment, 0); + + if (diag) + { + add_rotated_angle_to_segments(XMIN, YMAX - 0.5*dy, x1, YMAX - dy - 2.0*width, width, segment, 0); + add_rectangle_to_segments(XMIN, YMAX - 0.5*dy, XMIN - width, YMAX + 10.0, segment, 0); + } free(maze); } @@ -1822,7 +1999,20 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS]) } case (S_MAZE): { - init_maze_segments(segment); + init_maze_segments(segment, 0); + + cycle = 0; + for (i=0; i 0.95*a) return(0); + if (vabs(x) < 0.1*a) return(0); + if (y1 < 0.05 + c) return(0); + if (y1 < b - 0.05 + c) return(1); + return(y1 < a + b - 0.05 + c - vabs(x)); + } } } @@ -2293,6 +2495,14 @@ int in_segment_region(double x, double y, t_segment segment[NMAXSEGMENTS]) } return(1); } + case (S_MAZE_DIAG): + { + for (i=0; i 0.0) + { +// f = VICSEK_REPULSION*harmonic_force(distance, particle[k]); + f = VICSEK_REPULSION*coulomb_force(distance, particle[k]); + force[0] = f*ca; + force[1] = f*sa; + } + else + { + force[0] = 0.0; + force[1] = 0.0; + } + } + break; + } } if (ROTATION) @@ -3053,6 +3294,27 @@ int compute_particle_interaction(int i, int k, double force[2], double *torque, else *torque = sin(particle[k].angle - particle[i].angle); break; } + case (I_VICSEK_SHARK): + { + if (dist_scaled > 10.0) *torque = 0.0; + else if (particle[k].type == particle[i].type) /* fish adjusting direction */ + { + if (twrapx||twrapy) *torque = sin(-particle[k].angle - particle[i].angle); + else *torque = sin(particle[k].angle - particle[i].angle); + } + else if (particle[k].type == 2) /* fish fleeing a shark */ + { + alpha = argument(ca,sa); + particle[i].angle = alpha + PI; + *torque = 0.0; + } + else /* shark tracking fish */ + { + *torque = cos(particle[k].angle)*sa - sin(particle[k].angle)*ca; + } + + break; + } default: { spin_f = particle[i].spin_freq; @@ -3266,6 +3528,9 @@ int add_particle(double x, double y, double vx, double vy, double mass, short in // particle[i].spin_freq = SPIN_INTER_FREQUENCY_B; // } + if ((PLOT == P_NUMBER)||(PLOT_B == P_NUMBER)) + particle[i].color_hue = 360.0*(double)(i%N_PARTICLE_COLORS)/(double)N_PARTICLE_COLORS; + ncircles++; printf("Added particle at (%.3lg, %.3lg)\n", x, y); @@ -3442,6 +3707,13 @@ void compute_particle_colors(t_particle particle, int plot, double rgb[3], doubl *width = BOUNDARY_WIDTH; break; } + case (P_NUMBER): + { + hue = particle.color_hue; + *radius = particle.radius; + *width = BOUNDARY_WIDTH; + break; + } } switch (plot) { @@ -3466,6 +3738,13 @@ void compute_particle_colors(t_particle particle, int plot, double rgb[3], doubl hsl_to_rgb_twilight(hue, 0.9, 0.5, rgby); break; } + case (P_ANGLE): + { + hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_ANGLE); + hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_ANGLE); + hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_ANGLE); + break; + } case (P_DIRECT_ENERGY): { hsl_to_rgb_twilight(hue, 0.9, 0.5, rgb); @@ -3486,6 +3765,13 @@ void compute_particle_colors(t_particle particle, int plot, double rgb[3], doubl hsl_to_rgb_twilight(hue, 0.9, 0.5, rgby); break; } + case (P_INITIAL_POS): + { + hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_INITIAL_POS); + hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_INITIAL_POS); + hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_INITIAL_POS); + break; + } default: { hsl_to_rgb(hue, 0.9, 0.5, rgb); @@ -3988,6 +4274,29 @@ int draw_special_particle(t_particle particle, double xc1, double yc1, double ra } break; } + case (CHEM_ABDACBE): + { + if (particle.type == 4) + { + for (wsign = -1; wsign <= 1; wsign+=2) + { + x1 = xc1 + (double)wsign*0.7*radius*cos(particle.angle); + y1 = yc1 + (double)wsign*0.7*radius*sin(particle.angle); + if (wsign == 1) + { + if (fill) draw_colored_polygon(x1, y1, 1.2*MU_B, nsides, angle + APOLY*PID, rgb); + else draw_polygon(x1, y1, 1.2*MU_B, nsides, angle + APOLY*PID); + } + else + { + if (fill) draw_colored_polygon(x1, y1, 1.2*MU, nsides, angle + APOLY*PID, rgb); + else draw_polygon(x1, y1, 1.2*MU, nsides, angle + APOLY*PID); + } + } + return(0); + } + break; + } } return(1); } @@ -4011,7 +4320,7 @@ void draw_one_particle(t_particle particle, double xc, double yc, double radius, /* specific shapes for chemical reactions */ if (REACTION_DIFFUSION) cont = draw_special_particle(particle, xc1, yc1, radius, angle, nsides, rgb, 1); - if ((particle.interaction == I_LJ_QUADRUPOLE)||(particle.interaction == I_LJ_DIPOLE)||(particle.interaction == I_VICSEK)||(particle.interaction == I_VICSEK_REPULSIVE)||(particle.interaction == I_VICSEK_SPEED)) + if ((particle.interaction == I_LJ_QUADRUPOLE)||(particle.interaction == I_LJ_DIPOLE)||(particle.interaction == I_VICSEK)||(particle.interaction == I_VICSEK_REPULSIVE)||(particle.interaction == I_VICSEK_SPEED)||(particle.interaction == I_VICSEK_SHARK)) draw_colored_rhombus(xc1, yc1, radius, angle + APOLY*PID, rgb); else if (cont) draw_colored_polygon(xc1, yc1, radius, nsides, angle + APOLY*PID, rgb); @@ -4041,7 +4350,7 @@ void draw_one_particle(t_particle particle, double xc, double yc, double radius, glColor3f(1.0, 1.0, 1.0); if (REACTION_DIFFUSION) cont = draw_special_particle(particle, xc1, yc1, radius, angle, nsides, rgb, 0); - if ((particle.interaction == I_LJ_QUADRUPOLE)||(particle.interaction == I_LJ_DIPOLE)||(particle.interaction == I_VICSEK)||(particle.interaction == I_VICSEK_REPULSIVE)||(particle.interaction == I_VICSEK_SPEED)) + if ((particle.interaction == I_LJ_QUADRUPOLE)||(particle.interaction == I_LJ_DIPOLE)||(particle.interaction == I_VICSEK)||(particle.interaction == I_VICSEK_REPULSIVE)||(particle.interaction == I_VICSEK_SPEED)||(particle.interaction == I_VICSEK_SHARK)) draw_rhombus(xc1, yc1, radius, angle + APOLY*PID); else if (cont) draw_polygon(xc1, yc1, radius, nsides, angle + APOLY*PID); @@ -5646,6 +5955,30 @@ double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacl return(fperp); } + case (BC_REFLECT_ABS_BOTTOM): + { + /* add harmonic force outside screen */ + padding = MU; + x = particle[j].xc; + y = particle[j].yc; + + if (y < BCYMIN) + { + particle[j].active = 0; + particle[j].vx = 0.0; + particle[j].vy = 0.0; + particle[j].xc = BCXMAX + 2.0*padding; + particle[j].yc = BCYMIN - 2.0*padding; + } + else + { + if (y > BCYMAX + padding) particle[j].fy -= KSPRING_BOUNDARY*(y - BCYMAX - padding); + if (x > BCXMAX - padding) particle[j].fx -= KSPRING_BOUNDARY*(x - BCXMAX + padding); + else if (x < BCXMIN + padding) particle[j].fx += KSPRING_BOUNDARY*(BCXMIN - x + padding); + } + + return(0.0); + } } } @@ -5730,9 +6063,8 @@ int reorder_particles(t_particle particle[NMAXCIRCLES], double py[NMAXCIRCLES], } -int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY], - t_obstacle obstacle[NMAXOBSTACLES], double px[NMAXCIRCLES], double py[NMAXCIRCLES], double pangle[NMAXCIRCLES], int tracer_n[N_TRACER_PARTICLES], - t_segment segment[NMAXSEGMENTS]) +int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY], t_obstacle obstacle[NMAXOBSTACLES], double px[NMAXCIRCLES], double py[NMAXCIRCLES], double pangle[NMAXCIRCLES], int tracer_n[N_TRACER_PARTICLES], +t_segment segment[NMAXSEGMENTS]) /* initialize all particles, obstacles, and the hashgrid */ { int i, j, k, n, nactive = 0; @@ -5747,6 +6079,11 @@ int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgr particle[i].type = 2; particle[i].radius = MU_B; } + if ((INTERACTION == I_VICSEK_SHARK)&&(i==1)) + { + particle[i].type = 2; + particle[i].radius = MU_B; + } particle[i].neighb = 0; particle[i].diff_neighb = 0; @@ -5781,6 +6118,13 @@ int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgr particle[i].vx = V_INITIAL*gaussian(); particle[i].vy = V_INITIAL*gaussian(); + + if ((INTERACTION == I_VICSEK_SHARK)&&(i==1)) + { + particle[i].vx *= 1000.0; + particle[i].vy *= 1000.0; + } + particle[i].energy = (particle[i].vx*particle[i].vx + particle[i].vy*particle[i].vy)*particle[i].mass_inv; px[i] = particle[i].vx; @@ -5802,6 +6146,8 @@ int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgr if ((PLOT == P_INITIAL_POS)||(PLOT_B == P_INITIAL_POS)) particle[i].color_hue = 360.0*(particle[i].yc - INITYMIN)/(INITYMAX - INITYMIN); + else if ((PLOT == P_NUMBER)||(PLOT_B == P_NUMBER)) + particle[i].color_hue = 360.0*(double)(i%N_PARTICLE_COLORS)/(double)N_PARTICLE_COLORS; } /* initialize dummy values in case particles are added */ for (i=ncircles; i < NMAXCIRCLES; i++) @@ -5892,7 +6238,7 @@ int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgr /* position-dependent particle type */ if (POSITION_DEPENDENT_TYPE) for (i=0; i 1) delta_n = 1; if (EXOTHERMIC) *delta_e = (double)(delta_n)*DELTA_EKIN; return(ncollisions); } @@ -7240,11 +7620,30 @@ int update_types(t_particle particle[NMAXCIRCLES], t_collision *collisions, int } } return(ncollisions); - } + } + case (CHEM_ABDACBE): + { + for (i=0; i emax) + { + printf("Particle %i at (%.3lg, %.3lg) has energy %.5lg, resetting to 0\n", i, particle[i].xc, particle[i].yc, particle[i].energy); + particle[i].vx = 0.0; + particle[i].vy = 0.0; + px[i] = 0.0; + py[i] = 0.0; + particle[i].energy = 0.0; + } + printf("\n\n"); + +} diff --git a/sub_rde.c b/sub_rde.c index ce8cf68..7b78421 100644 --- a/sub_rde.c +++ b/sub_rde.c @@ -699,6 +699,127 @@ void init_pressure_gradient_flow(double vx, double pmax, double pmin, double *ph } } +double init_gaussian_wave(double x, double y, double amp, double radius, double min, double *phi[NFIELDS], short int xy_in[NX*NY]) +/* initialise gaussian wave height for shallow water equation */ +{ + int i, j; + double xy[2], var, a; + + var = radius*radius; + a = amp/(sqrt(DPI)*radius); + + for (i=0; i 0.5*(LAMBDA + 1.0)) f = 0.5*(1.0 + tanh(BC_STIFFNESS*(1.0 - r))); + else f = 0.5*(1.0 + tanh(BC_STIFFNESS*(r - LAMBDA))); + bc_field[i*NY+j] = f; + f = 0.5*(1.0 + tanh(0.5*BC_STIFFNESS*(0.99 - r))); + bc_field2[i*NY+j] = f; + } + break; + } case (D_WING): { a = 0.4; @@ -850,11 +1000,21 @@ void initialize_bcfield(double bc_field[NX*NY], double bc_field2[NX*NY], t_recta nsides = init_polyrect_euler(polyrect, D_MAZE); break; } + case (D_MAZE_CLOSED): + { + nsides = init_polyrect_euler(polyrect, D_MAZE_CLOSED); + break; + } case (D_MAZE_CHANNELS): { nsides = init_polyrect_euler(polyrect, D_MAZE_CHANNELS); break; } + case (D_MAZE_CHANNELS_INT): + { + nsides = init_polyrect_euler(polyrect, D_MAZE_CHANNELS_INT); + break; + } case (D_TESLA): { a = 0.16; @@ -947,7 +1107,7 @@ void initialize_bcfield(double bc_field[NX*NY], double bc_field2[NX*NY], t_recta } } - if ((OBSTACLE_GEOMETRY == D_MAZE)||(OBSTACLE_GEOMETRY == D_MAZE_CHANNELS)) + if ((OBSTACLE_GEOMETRY == D_MAZE)||(OBSTACLE_GEOMETRY == D_MAZE_CLOSED)||(OBSTACLE_GEOMETRY == D_MAZE_CHANNELS)) { d0 = 0.25*MAZE_WIDTH; fmin = 0.5*(1.0 - tanh(d0*BC_STIFFNESS)); @@ -1012,7 +1172,78 @@ void initialize_bcfield(double bc_field[NX*NY], double bc_field2[NX*NY], t_recta else f = 0.5*(1.0 + tanh(BC_STIFFNESS*(distance - 1.25*MAZE_WIDTH))); bc_field[i*NY+j] = f; - if (distance >= d0) f = 0.5*(1.0 + tanh(BC_STIFFNESS*(distance - 1.5*MAZE_WIDTH))); + if (distance >= d0) f = 0.5*(1.0 + tanh(0/BC_STIFFNESS*(distance - 1.5*MAZE_WIDTH))); + bc_field2[i*NY+j] = f; +// printf("distance = %.5lg, bcfield = %.5lg\n", distance, f); + } + } + } + else if (OBSTACLE_GEOMETRY == D_MAZE_CHANNELS_INT) + { + d0 = 2.0*MAZE_WIDTH; + fmin = 0.5*(1.0 - tanh(d0*BC_STIFFNESS)); + for (i=0; i 4.0*MAZE_WIDTH)&&(height > 4.0*MAZE_WIDTH)) + { + if (x < x1) + { + if (y < y1) d = module2(x - x1, y - y1); + else if (y > y2) d = module2(x - x1, y - y2); + else d = x1 - x; + } + else if (x > x2) + { + if (y < y1) d = module2(x - x2, y - y1); + else if (y > y2) d = module2(x - x2, y - y2); + else d = x - x2; + } + else + { + if (y < y1) d = y1 - y; + else if (y > y2) d = y - y2; + else d = 0.0; + } + } + else if (length > height) + { + mid = 0.5*(polyrect[s].y1 + polyrect[s].y2); + if ((x > x1)&&(x < x2)) d = vabs(y - mid); + else if (x <= x1) d = module2(x - x1, y - mid); + else d = module2(x - x2, y - mid); + } + else + { + mid = 0.5*(polyrect[s].x1 + polyrect[s].x2); + if ((y > y1)&&(y < y2)) d = vabs(x - mid); + else if (y <= y1) d = module2(x - mid, y - y1); + else d = module2(x - mid, y - y2); + } + if (d < distance) distance = d; + } + if (distance < d0) f = 1.0; + else f = 0.5*(1.0 - tanh(BC_STIFFNESS*(distance - 1.25*MAZE_WIDTH))); + bc_field[i*NY+j] = f; + + if (distance >= d0) f = 0.5*(1.0 - tanh(0.5*BC_STIFFNESS*(distance - 1.5*MAZE_WIDTH))); bc_field2[i*NY+j] = f; // printf("distance = %.5lg, bcfield = %.5lg\n", distance, f); } @@ -1025,15 +1256,22 @@ void adapt_state_to_bc(double *phi[NFIELDS], double bc_field[NX*NY], short int x /* apply smooth modulation to adapt initial state to obstacles */ { int i, j, field; - double xy[2], r, f; -// double ratio = 1.0e-1; + double xy[2], r, f, integral = 0.0, factor; + + if (CHECK_INTEGRAL) + { + integral = 0.0; + +// #pragma omp parallel for private(i) + for (i=0; i 360.0) value -= 360.0; + hsl_to_rgb_palette(value, 0.9, 0.5, rgb, palette); + break; + } + } } @@ -2169,10 +3085,7 @@ double adjust_field(double z, double pot) } -double compute_interpolated_colors_rde(int i, int j, t_rde rde[NX*NY], double potential[NX*NY], double palette, int cplot, - double rgb_e[3], double rgb_w[3], double rgb_n[3], double rgb_s[3], - double *z_sw, double *z_se, double *z_nw, double *z_ne, - int fade, double fade_value, int movie) +double compute_interpolated_colors_rde(int i, int j, t_rde rde[NX*NY], double potential[NX*NY], double palette, int cplot, double rgb_e[3], double rgb_w[3], double rgb_n[3], double rgb_s[3], double *z_sw, double *z_se, double *z_nw, double *z_ne, int fade, double fade_value, int movie) { int k; double cw, ce, cn, cs, c_sw, c_se, c_nw, c_ne, c_mid, ca, z_mid; @@ -2226,7 +3139,8 @@ double compute_interpolated_colors_rde(int i, int j, t_rde rde[NX*NY], double po compute_field_color_rde(cs, cplot, palette, rgb_s); // if ((cplot == Z_ARGUMENT)||(cplot == Z_REALPART)) - if (cplot == Z_ARGUMENT) +// if (cplot == Z_ARGUMENT) + if ((cplot == Z_ARGUMENT)||(cplot == Z_EULER_DIRECTION_SPEED)||(cplot == Z_SWATER_DIRECTION_SPEED)) { lum = tanh(SLOPE_SCHROD_LUM*rde[i*NY+j].field_norm) + MIN_SCHROD_LUM; for (k=0; k<3; k++) @@ -2262,6 +3176,45 @@ double compute_interpolated_colors_rde(int i, int j, t_rde rde[NX*NY], double po } +double compute_depth_colors_rde(int i, int j, t_rde rde[NX*NY], double potential[NX*NY], double palette, int cplot, double rgb_e[3], double rgb_w[3], double rgb_n[3], double rgb_s[3], double *z_sw, double *z_se, double *z_nw, double *z_ne, int fade, double fade_value, int movie) +{ + int k; + double cw, ce, cn, cs, c_sw, c_se, c_nw, c_ne, c_mid, ca, z_mid; + double lum; + + *z_sw = -DEPTH_SCALE*rde[i*NY+j].depth + DEPTH_SHIFT; + *z_se = -DEPTH_SCALE*rde[(i+1)*NY+j].depth + DEPTH_SHIFT; + *z_nw = -DEPTH_SCALE*rde[i*NY+j+1].depth + DEPTH_SHIFT; + *z_ne = -DEPTH_SCALE*rde[(i+1)*NY+j+1].depth + DEPTH_SHIFT; + + z_mid = 0.25*(*z_sw + *z_se + *z_nw + *z_ne); + + if (SHADE_3D) + { +// printf("cos angle = %.3lg\n", rde[i*NY+j].cos_depth_angle); + ca = rde[i*NY+j].cos_depth_angle; + ca = (ca + 1.0)*0.2 + 0.1; + for (k=0; k<3; k++) + { + rgb_e[k] = ca; + rgb_w[k] = ca; + rgb_n[k] = ca; + rgb_s[k] = ca; + } + } + if (fade) + for (k=0; k<3; k++) + { + rgb_e[k] *= fade_value; + rgb_w[k] *= fade_value; + rgb_n[k] *= fade_value; + rgb_s[k] *= fade_value; + } + + return(z_mid); +} + + void compute_rde_fields(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot, int cplot, t_rde rde[NX*NY]) /* compute the necessary auxiliary fields */ { @@ -2330,6 +3283,14 @@ void compute_rde_fields(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot, compute_vorticity(rde); break; } + case (E_SHALLOW_WATER): + { + if (zplot == Z_SWATER_HEIGHT) adjust_height(phi, rde); + if ((zplot == Z_SWATER_SPEED)||(cplot == Z_SWATER_SPEED)||(zplot == Z_SWATER_DIRECTION_SPEED)||(cplot == Z_SWATER_DIRECTION_SPEED)) + compute_speed(phi, rde); + if ((zplot == Z_SWATER_DIRECTION_SPEED)||(cplot == Z_SWATER_DIRECTION_SPEED)) + compute_direction(phi, rde); + } default : break; } } @@ -2490,6 +3451,24 @@ void init_zfield_rde(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot, t_ for (i=0; i= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_MEAN_ENERGY): + { + value = dy_e*(double)(j)*100.0/E_SCALE; + if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_LOG_ENERGY): + { + value = LOG_SCALE*log(dy_e*(double)(j)*100.0/E_SCALE); + color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_LOG_MEAN_ENERGY): + { + value = LOG_SHIFT + LOG_SCALE*log(dy_e*(double)(j)*100.0/E_SCALE); + color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_ENERGY_FLUX): + { + value = dy_e*(double)(j); + if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_TOTAL_ENERGY_FLUX): + { + value = dy_phase*(double)(j); + color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb); + break; + } + case (P_PHASE): + { + value = min + 1.0*dy*(double)(j); + amp = 0.5*color_amplitude_linear(value, 1.0, 1); + while (amp > 1.0) amp -= 2.0; + while (amp < -1.0) amp += 2.0; + amp_to_rgb(0.5*(1.0 + amp), rgb); + break; + } + case (Z_EULER_VORTICITY): + { + value = min + 1.0*dy*(double)(j); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (Z_EULER_LOG_VORTICITY): + { + value = min + 1.0*dy*(double)(j); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (Z_EULER_VORTICITY_ASYM): + { + value = min + 1.0*dy*(double)(j); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (Z_EULER_LPRESSURE): + { + value = min + 1.0*dy*(double)(j); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (Z_EULER_PRESSURE): + { + value = min + 1.0*dy*(double)(j); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (Z_EULER_DENSITY): + { + value = min + 1.0*dy*(double)(j); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (Z_EULER_SPEED): + { + value = min + 1.0*dy*(double)(j); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (Z_EULERC_VORTICITY): + { + value = min + 1.0*dy*(double)(j); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + default: + { + value = min + 1.0*dy*(double)(j); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + } + if (fade) for (k=0; k<3; k++) rgb[k] *= fade_value; + glColor3f(rgb[0], rgb[1], rgb[2]); + + xy_to_pos(x1 + radius*cos(dphi*(double)j), y1 + radius*sin(dphi*(double)j), xy); + glVertex2d(xy[0], xy[1]); + } + xy_to_pos(x1 + radius*cos(dphi), y1 + radius*sin(dphi), xy); + glVertex2d(xy[0], xy[1]); + glEnd (); + + if (fade) glColor3f(fade_value, fade_value, fade_value); + else glColor3f(1.0, 1.0, 1.0); + glLineWidth(BOUNDARY_WIDTH*3/2); + glEnable(GL_LINE_SMOOTH); + + dphi = DPI/(double)NSEG; + glBegin(GL_LINE_LOOP); + for (j = 0; j < NSEG; j++) + { + xy_to_pos(x1 + radius*cos(dphi*(double)j), y1 + radius*sin(dphi*(double)j), xy); + glVertex2d(xy[0], xy[1]); + glVertex2d(xy[0], xy[1]); + } + glEnd (); +} + + void print_speed(double speed, int fade, double fade_value) { char message[100]; @@ -5550,12 +5793,13 @@ void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double ior_angle) /* compute variable index of refraction */ /* should be at some point merged with 3D version in suv_wave_3d.c */ { - int i, j, k, n, inlens; + int i, j, k, n, inlens, ncircles; double courant2 = COURANT*COURANT, courantb2 = COURANTB*COURANTB, lambda1, mu1; double u, v, u1, x, y, xy[2], norm2, speed, r2, c, salpha, h, ll, ca, sa, x1, y1, dx, dy, sum, sigma, x0, y0, rgb[3]; double xc[NGRIDX*NGRIDY], yc[NGRIDX*NGRIDY], height[NGRIDX*NGRIDY]; static double xc_stat[NGRIDX*NGRIDY], yc_stat[NGRIDX*NGRIDY], sigma_stat; static int first = 1; + t_circle circles[NMAXCIRCLES]; rgb[0] = 1.0; rgb[1] = 1.0; @@ -5749,7 +5993,9 @@ void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double ior_angle) } // #pragma omp parallel for private(i,j) - for (i=0; i 1.0) height[n] = 1.0; + if (height[n] < 0.0) height[n] = 0.0; + } + + for (n = 0; n 1.0) height[n] = 1.0; + if (height[n] < 0.0) height[n] = 0.0; + } + + for (n = 0; n 1.0) height[n] = 1.0; + if (height[n] < 0.0) height[n] = 0.0; + } + + for (n = 0; n 1.0) height[n] = 1.0; + if (height[n] < 0.0) height[n] = 0.0; + } + + for (n = 0; n= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_3D_LOG_ENERGY): + { + value = LOG_SCALE*dy_e*(double)(j)*100.0/E_SCALE; + color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_3D_TOTAL_ENERGY): + { + value = dy_e*(double)(j)*100.0/E_SCALE; + if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_3D_LOG_TOTAL_ENERGY): + { + value = LOG_SCALE*dy_e*(double)(j)*100.0/E_SCALE; + color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_3D_MEAN_ENERGY): + { + value = dy_e*(double)(j)*100.0/E_SCALE; + if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_3D_LOG_MEAN_ENERGY): + { + value = LOG_SCALE*dy_e*(double)(j)*100.0/E_SCALE; + color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_3D_PHASE): + { + value = dy_phase*(double)(j); + color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb); + break; + } + case (P_3D_FLUX_INTENSITY): + { + value = dy_e*(double)(j)*100.0/E_SCALE; + if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_3D_FLUX_DIRECTION): + { + value = dy_phase*(double)(j); + color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb); + break; + } + case (Z_EULER_VORTICITY): + { + value = min + 1.0*dy*(double)(j); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (Z_EULER_LOG_VORTICITY): + { + value = min + 1.0*dy*(double)(j); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (Z_EULER_VORTICITY_ASYM): + { + value = min + 1.0*dy*(double)(j); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (Z_EULER_LPRESSURE): + { + value = min + 1.0*dy*(double)(j); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (Z_EULERC_VORTICITY): + { + value = min + 1.0*dy*(double)(j); + printf("Palette value %.3lg\n", value); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (Z_SWATER_DIRECTION_SPEED): + { + value = dy_phase*(double)(j); + color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb); + break; + } + + } + if (fade) for (k=0; k<3; k++) rgb[k] *= fade_value; + glColor3f(rgb[0], rgb[1], rgb[2]); + x = cos(dphi*(double)j)*(double)ir; + y = sin(dphi*(double)j)*(double)ir; + ij_to_xy(ic + x, jc + y, xy); + glVertex2d(xy[0], xy[1]); + } + draw_vertex_ij(ic + ir, jc); + glEnd (); + + if (fade) glColor3f(fade_value, fade_value, fade_value); + else glColor3f(1.0, 1.0, 1.0); + glLineWidth(BOUNDARY_WIDTH*3/2); + glEnable(GL_LINE_SMOOTH); + + dphi = DPI/(double)NSEG; + glBegin(GL_LINE_LOOP); + for (j = 0; j < NSEG; j++) + { + x = cos(dphi*(double)j)*(double)ir; + y = sin(dphi*(double)j)*(double)ir; + xy[0] = XMIN + ((double)ic + x)*(XMAX-XMIN)/((double)NX); + xy[1] = YMIN + ((double)jc + y)*(YMAX-YMIN)/((double)NY); + glVertex2d(xy[0], xy[1]); + } + glEnd (); + +} + diff --git a/sub_wave_comp.c b/sub_wave_comp.c index 22d46df..c0f0932 100644 --- a/sub_wave_comp.c +++ b/sub_wave_comp.c @@ -483,6 +483,11 @@ int xy_in_billiard_half(double x, double y, int domain, int pattern, int top) int i, j, k, k1, k2, condition, type; switch (domain) { + case (D_NOTHING): + { + return(1); + break; + } case (D_MENGER): { x1 = 0.5*(x+1.0); @@ -744,7 +749,12 @@ void draw_billiard_half(int domain, int pattern, int top) /* draws the bill /* Do nothing */ break; } - default: + case (D_NOTHING): + { + /* Do nothing */ + break; + } + default: { printf("Function draw_billiard not defined for this billiard \n"); } @@ -799,6 +809,31 @@ void init_wave_flat_comp( double *phi[NX], double *psi[NX], short int * xy_in[NX } +void add_circular_wave_comp(double factor, double x, double y, double *phi[NX], double *psi[NX], short int * xy_in[NX], int half) +/* add drop at (x,y) to the field with given prefactor */ +{ + int i, j; + double xy[2], dist2; + + for (i=0; i= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, energy, scale, time, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, energy, scale, time, rgb); + break; + } + case (P_MIXED): + { + if (j > NY/2) color_scheme_palette(COLOR_SCHEME, palette, phi[i][j], scale, time, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, compute_energy(phi, psi, xy_in, i, j), scale, time, rgb); + break; + } + case (P_MEAN_ENERGY): + { + energy = compute_energy(phi, psi, xy_in, i, j); + total_energy[i][j] += energy; + if (COLOR_PALETTE >= COL_TURBO) + color_scheme_asym_palette(COLOR_SCHEME, palette, total_energy[i][j]/(double)(time+1), scale, time, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, total_energy[i][j]/(double)(time+1), scale, time, rgb); + break; + } + case (P_LOG_ENERGY): + { + energy = compute_energy(phi, psi, xy_in, i, j); + color_scheme(COLOR_SCHEME, LOG_SHIFT + LOG_SCALE*log(energy), scale, time, rgb); + break; + } +// case (P_LOG_MEAN_ENERGY): +// { +// energy = compute_energy(phi, psi, xy_in, i, j); +// if (energy == 0.0) energy = 1.0e-20; +// total_energy[i][j] += energy; +// color_scheme(COLOR_SCHEME, LOG_SCALE*log(total_energy[i][j]/(double)(time+1)), scale, time, rgb); +// break; +// } + } + +// if (PLOT == P_AMPLITUDE) +// color_scheme(COLOR_SCHEME, phi[i][j], scale, time, rgb); +// else if (PLOT == P_ENERGY) +// { +// energy = compute_energy(phi, psi, xy_in, i, j); +// if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym(COLOR_SCHEME, energy, scale, time, rgb); +// else color_scheme(COLOR_SCHEME, energy, scale, time, rgb); +// } +// else if (PLOT == 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); +// } + if (fade) for (k=0; k<3; k++) rgb[k] *= fade_value; + glColor3f(rgb[0], rgb[1], rgb[2]); + + glVertex2i(i, j); + glVertex2i(i+size, j); + glVertex2i(i+size, j+size); + glVertex2i(i, j+size); + } + } + + glEnd (); + + /* draw horizontal mid line */ + glColor3f(1.0, 1.0, 1.0); + glBegin(GL_LINE_STRIP); + xy_to_pos(XMIN, 0.5*(YMIN+YMAX), pos); + glVertex2d(pos[0], pos[1]); + xy_to_pos(XMAX, 0.5*(YMIN+YMAX), pos); + glVertex2d(pos[0], pos[1]); + glEnd(); +} + void compute_energy_tblr(double *phi[NX], double *psi[NX], short int *xy_in[NX], double energies[6]) /* compute total energy in top/bottom left/right boxes */ { diff --git a/wave_3d.c b/wave_3d.c index 3dfa63d..a79ba10 100644 --- a/wave_3d.c +++ b/wave_3d.c @@ -54,21 +54,19 @@ #define WINHEIGHT 1150 /* window height */ // // // #define NX 2500 /* number of grid points on x axis */ // // // #define NY 1250 /* number of grid points on y axis */ +#define NX 1920 /* number of grid points on x axis */ +#define NY 1150 /* number of grid points on y axis */ // #define NX 2840 /* number of grid points on x axis */ // #define NY 2300 /* number of grid points on y axis */ // #define NX 2500 /* number of grid points on x axis */ // #define NY 1250 /* number of grid points on y axis */ -#define NX 1700 /* number of grid points on x axis */ -#define NY 1700 /* number of grid points on y axis */ +// #define NX 1700 /* number of grid points on x axis */ +// #define NY 1700 /* number of grid points on y axis */ -// #define XMIN -2.0 -// #define XMAX 2.0 /* x interval */ -// #define YMIN -1.041666667 -// #define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */ -#define XMIN -1.5 -#define XMAX 1.5 /* x interval */ -#define YMIN -1.5 -#define YMAX 1.5 /* y interval for 9/16 aspect ratio */ +#define XMIN -1.0 +#define XMAX 3.0 /* x interval */ +#define YMIN -1.041666667 +#define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */ #define HIGHRES 0 /* set to 1 if resolution of grid is double that of displayed image */ @@ -94,7 +92,7 @@ /* Choice of the billiard table */ -#define B_DOMAIN 999 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 10 /* choice of domain shape, see list in global_pdes.c */ #define CIRCLE_PATTERN 2 /* pattern of circles or polygons, see list in global_pdes.c */ @@ -103,16 +101,16 @@ #define CIRCLE_PATTERN_B 0 /* second pattern of circles or polygons */ #define VARIABLE_IOR 0 /* set to 1 for a variable index of refraction */ -#define IOR 101 /* choice of index of refraction, see list in global_pdes.c */ +#define IOR 9 /* choice of index of refraction, see list in global_pdes.c */ #define IOR_TOTAL_TURNS 1.0 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */ #define MANDEL_IOR_SCALE -0.05 /* parameter controlling dependence of IoR on Mandelbrot escape speed */ #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 NPOISSON 1000 /* number of points for Poisson C_RAND_POISSON arrangement */ #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ -#define LAMBDA 0.5 /* parameter controlling the dimensions of domain */ -#define MU 0.5 /* parameter controlling the dimensions of domain */ +#define LAMBDA 0.1197916667 /* parameter controlling the dimensions of domain */ +#define MU 0.035 /* parameter controlling the dimensions of domain */ #define NPOLY 6 /* number of sides of polygon */ #define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */ #define MDEPTH 7 /* depth of computation of Menger gasket */ @@ -120,8 +118,8 @@ #define MANDELLEVEL 2000 /* iteration level for Mandelbrot set */ #define MANDELLIMIT 20.0 /* limit value for approximation of Mandelbrot set */ #define FOCI 1 /* set to 1 to draw focal points of ellipse */ -#define NGRIDX 8 /* number of grid point for grid of disks */ -#define NGRIDY 8 /* number of grid point for grid of disks */ +#define NGRIDX 30 /* number of grid point for grid of disks */ +#define NGRIDY 18 /* number of grid point for grid of disks */ #define X_SHOOTER -0.2 #define Y_SHOOTER -0.6 @@ -149,8 +147,8 @@ #define AMPLITUDE 0.8 /* amplitude of periodic excitation */ #define ACHIRP 0.2 /* acceleration coefficient in chirp */ #define DAMPING 0.0 /* damping of periodic excitation */ -#define COURANT 0.06 /* Courant number */ -#define COURANTB 0.06 /* Courant number in medium B */ +#define COURANT 0.1 /* Courant number */ +#define COURANTB 0.01 /* Courant number in medium B */ #define GAMMA 0.0 /* damping factor in wave equation */ #define GAMMAB 1.0e-6 /* damping factor in wave equation */ #define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */ @@ -164,7 +162,7 @@ /* For similar wave forms, COURANT^2*GAMMA should be kept constant */ #define ADD_OSCILLATING_SOURCE 1 /* set to 1 to add an oscillating wave source */ -#define OSCILLATING_SOURCE_PERIOD 2 /* period of oscillating source */ +#define OSCILLATING_SOURCE_PERIOD 30 /* period of oscillating source */ #define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */ #define ADD_WAVE_PACKET_SOURCES 0 /* set to 1 to add several sources emitting wave packets */ @@ -181,9 +179,9 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 2450 /* number of frames of movie */ -// #define NSTEPS 500 /* number of frames of movie */ -#define NVID 5 /* number of iterations between images displayed on screen */ +#define NSTEPS 2300 /* number of frames of movie */ +// #define NSTEPS 300 /* number of frames of movie */ +#define NVID 4 /* number of iterations between images displayed on screen */ // #define NVID 10 /* number of iterations between images displayed on screen */ #define NSEG 1000 /* number of segments of boundary */ #define INITIAL_TIME 0 /* time after which to start saving frames */ @@ -201,18 +199,19 @@ /* Parameters of initial condition */ -// #define INITIAL_AMP 0.02 /* amplitude of initial condition */ -#define INITIAL_AMP 0.007 /* amplitude of initial condition */ -#define INITIAL_VARIANCE 0.0003 /* variance of initial condition */ -#define INITIAL_WAVELENGTH 0.015 /* wavelength of initial condition */ +#define INITIAL_AMP 2.0 /* amplitude of initial condition */ +// #define INITIAL_VARIANCE 0.000025 /* variance of initial condition */ +#define INITIAL_VARIANCE 0.00005 /* variance of initial condition */ +#define INITIAL_WAVELENGTH 0.05 /* wavelength of initial condition */ + /* Plot type, see list in global_pdes.c */ #define ZPLOT 103 /* wave height */ #define CPLOT 103 /* color scheme */ -#define ZPLOT_B 112 -#define CPLOT_B 113 /* plot type for second movie */ +#define ZPLOT_B 107 +#define CPLOT_B 107 /* plot type for second movie */ #define CHANGE_LUMINOSITY 1 /* set to 1 to let luminosity depend on energy flux intensity */ #define FLUX_WINDOW 30 /* size of averaging window of flux intensity */ @@ -242,8 +241,8 @@ /* Color schemes */ -#define COLOR_PALETTE 15 /* Color palette, see list in global_pdes.c */ -#define COLOR_PALETTE_B 10 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE_B 13 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* background */ @@ -252,11 +251,11 @@ #define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ #define SLOPE 1.0 /* sensitivity of color on wave amplitude */ #define VSCALE_AMPLITUDE 2.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ -#define VSCALE_ENERGY 15.0 /* additional scaling factor for color scheme P_3D_ENERGY */ +#define VSCALE_ENERGY 10.0 /* additional scaling factor for color scheme P_3D_ENERGY */ #define PHASE_FACTOR 20.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 240.0 /* scaling factor for energy representation */ +#define E_SCALE 20.0 /* scaling factor for energy representation */ #define LOG_SCALE 0.75 /* scaling factor for energy log representation */ #define LOG_SHIFT 0.5 /* shift of colors on log scale */ #define LOG_ENERGY_FLOOR -10.0 /* floor value for log of (total) energy */ @@ -281,7 +280,7 @@ #define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */ #define COLORBAR_RANGE 2.5 /* scale of color scheme bar */ -#define COLORBAR_RANGE_B 5.0 /* scale of color scheme bar for 2nd part */ +#define COLORBAR_RANGE_B 0.5 /* scale of color scheme bar for 2nd part */ #define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ #define SAVE_TIME_SERIES 0 /* set to 1 to save wave time series at a point */ @@ -289,7 +288,7 @@ #define ADD_POTENTIAL 1 /* set to 1 to add potential to z coordinate */ // #define POT_MAZE 7 #define POTENTIAL 10 -#define POT_FACT 400.0 +#define POT_FACT 20.0 /* end of constants only used by sub_wave and sub_maze */ @@ -303,14 +302,14 @@ double u_3d[2] = {0.75, -0.45}; /* projections of basis vectors for REP_AXO_ double v_3d[2] = {-0.75, -0.45}; double w_3d[2] = {0.0, 0.015}; double light[3] = {0.816496581, -0.40824829, 0.40824829}; /* vector of "light" direction for P_3D_ANGLE color scheme */ -double observer[3] = {8.0, 8.0, 8.0}; /* location of observer for REP_PROJ_3D representation */ +double observer[3] = {8.0, 6.0, 5.0}; /* location of observer for REP_PROJ_3D representation */ int reset_view = 0; /* switch to reset 3D view parameters (for option ROTATE_VIEW) */ -#define Z_SCALING_FACTOR 0.5 /* overall scaling factor of z axis for REP_PROJ_3D representation */ -#define XY_SCALING_FACTOR 1.8 /* overall scaling factor for on-screen (x,y) coordinates after projection */ +#define Z_SCALING_FACTOR 0.1 /* overall scaling factor of z axis for REP_PROJ_3D representation */ +#define XY_SCALING_FACTOR 1.7 /* overall scaling factor for on-screen (x,y) coordinates after projection */ #define ZMAX_FACTOR 1.0 /* max value of z coordinate for REP_PROJ_3D representation */ -#define XSHIFT_3D -0.1 /* overall x shift for REP_PROJ_3D representation */ -#define YSHIFT_3D 0.2 /* overall y shift for REP_PROJ_3D representation */ +#define XSHIFT_3D 0.2 /* overall x shift for REP_PROJ_3D representation */ +#define YSHIFT_3D 0.5 /* overall y shift for REP_PROJ_3D representation */ #include "global_pdes.c" /* constants and global variables */ @@ -1067,8 +1066,8 @@ void animation() // for (j=1; j 0.0) wave_source[k].sign = 1; -// else wave_source[k].sign = -1; - wave_source[k].sign = 1; - } - first_source = 0; - - for (k=0; k<25; k++) - wave_source[k].phase += 1.4*sin(0.7*(1.0 + wave_source[k].xc*cos(angle)/0.05 + wave_source[k].yc*sin(angle)/0.05)); - - angle = DPI*(double)i/(double)NSTEPS; -// phase_shift = 0.02 + 0.08*(double)i/(double)NSTEPS; -// printf("Phase shift = %.3lg\n", phase_shift); - - for (k=0; k<25; k++) - { -// wave_source[k].phase += 0.07; - wave_source[k].phase += phase_shift; - wave_source[k].phase -= 1.4*sin(0.7*(1.0 + wave_source[k].xc*cos(angle)/0.05 + wave_source[k].yc*sin(angle)/0.05)); - - if (wave_source[k].phase > 1.0) - { - add_circular_wave_mod((double)wave_source[k].sign*wave_source[k].amp, wave_source[k].xc, wave_source[k].yc, phi, psi, xy_in); - printf("Adding wave at (%.2lg, %.2lg)\n", wave_source[k].xc, wave_source[k].yc); - wave_source[k].phase -= 1.0; - wave_source[k].sign *= -1; - } - } +// phase_shift = 1/30.0; +// // phase_shift = 1/60.0; +// +// if (first_source) for (k=0; k<25; k++) +// { +// angle = 0.0; +// omega = DPI/25.0; +// wave_source[k].xc = 0.05*cos((double)k*omega);; +// wave_source[k].yc = 0.05*sin((double)k*omega); +// wave_source[k].phase = 0.99 - 1.4*sin(0.7*(1.0 + wave_source[k].xc/0.05)); +// wave_source[k].amp = 1.0; +// // if (wave_source[k].phase > 0.0) wave_source[k].sign = 1; +// // else wave_source[k].sign = -1; +// wave_source[k].sign = 1; +// } +// first_source = 0; +// +// for (k=0; k<25; k++) +// wave_source[k].phase += 1.4*sin(0.7*(1.0 + wave_source[k].xc*cos(angle)/0.05 + wave_source[k].yc*sin(angle)/0.05)); +// +// angle = DPI*(double)i/(double)NSTEPS; +// // phase_shift = 0.02 + 0.08*(double)i/(double)NSTEPS; +// // printf("Phase shift = %.3lg\n", phase_shift); +// +// for (k=0; k<25; k++) +// { +// // wave_source[k].phase += 0.07; +// wave_source[k].phase += phase_shift; +// wave_source[k].phase -= 1.4*sin(0.7*(1.0 + wave_source[k].xc*cos(angle)/0.05 + wave_source[k].yc*sin(angle)/0.05)); +// +// if (wave_source[k].phase > 1.0) +// { +// add_circular_wave_mod((double)wave_source[k].sign*wave_source[k].amp, wave_source[k].xc, wave_source[k].yc, phi, psi, xy_in); +// printf("Adding wave at (%.2lg, %.2lg)\n", wave_source[k].xc, wave_source[k].yc); +// wave_source[k].phase -= 1.0; +// wave_source[k].sign *= -1; +// } +// } // for (j=0; j 1.0) - { - add_circular_wave((double)wave_source[k].sign*wave_source[k].amp, wave_source[k].xc, wave_source[k].yc, phi, psi, xy_in); - printf("Adding wave at (%.2lg, %.2lg)\n", wave_source[k].xc, wave_source[k].yc); - wave_source[k].phase -= 1.0; - wave_source[k].sign *= -1; - } - } +// if (first_source) for (k=0; k<25; k++) +// { +// omega = DPI/25.0; +// wave_source[k].xc = 0.05*cos((double)k*omega);; +// wave_source[k].yc = 0.05*sin((double)k*omega); +// wave_source[k].phase = 0.99 - 1.4*sin(0.7*(1.0 + wave_source[k].xc/0.05)); +// wave_source[k].amp = 1.0; +// if (wave_source[k].phase) wave_source[k].sign = 1; +// else wave_source[k].sign = -1; +// first_source = 0; +// } +// +// for (k=0; k<25; k++) +// wave_source[k].phase += 1.4*sin(0.7*(1.0 + wave_source[k].xc*cos(angle)/0.05 + wave_source[k].yc*sin(angle)/0.05)); +// +// angle = DPI*(double)i/(double)NSTEPS; +// +// for (k=0; k<25; k++) +// { +// wave_source[k].phase += 0.07; +// wave_source[k].phase -= 1.4*sin(0.7*(1.0 + wave_source[k].xc*cos(angle)/0.05 + wave_source[k].yc*sin(angle)/0.05)); +// +// if (wave_source[k].phase > 1.0) +// { +// add_circular_wave((double)wave_source[k].sign*wave_source[k].amp, wave_source[k].xc, wave_source[k].yc, phi, psi, xy_in); +// printf("Adding wave at (%.2lg, %.2lg)\n", wave_source[k].xc, wave_source[k].yc); +// wave_source[k].phase -= 1.0; +// wave_source[k].sign *= -1; +// } +// } // p = 3; // y = -1.0; @@ -893,7 +893,7 @@ void animation() draw_wave_highres_palette(2, phi, psi, total_energy, total_flux, xy_in, scale, i, PLOT_B, COLOR_PALETTE_B, 0, 1.0); else draw_wave_epalette(phi, psi, total_energy, total_flux, color_scale, xy_in, scale, i, PLOT_B, COLOR_PALETTE_B, 0, 1.0); draw_billiard(0, 1.0); - if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0); + if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, CIRC_COLORBAR_B, 0, 1.0); if (PRINT_SPEED) print_speed(speed, 0, 1.0); if (PRINT_FREQUENCY) print_frequency(phase_shift, 0, 1.0); glutSwapBuffers(); @@ -922,7 +922,7 @@ void animation() if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, total_flux, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 0, 1.0); else draw_wave_epalette(phi, psi, total_energy, total_flux, color_scale, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 0, 1.0); draw_billiard(0, 1.0); - if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0); + if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, 0, 1.0); if (PRINT_SPEED) print_speed(speed, 0, 1.0); if (PRINT_FREQUENCY) print_frequency(phase_shift, 0, 1.0); glutSwapBuffers(); @@ -935,7 +935,7 @@ void animation() draw_wave_highres_palette(2, phi, psi, total_energy, total_flux, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 1, fade_value); else draw_wave_epalette(phi, psi, total_energy, total_flux, color_scale, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 1, fade_value); draw_billiard(1, fade_value); - if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, 1, fade_value); + if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, 1, fade_value); if (PRINT_SPEED) print_speed(speed, 1, fade_value); if (PRINT_FREQUENCY) print_frequency(phase_shift, 1, fade_value); if (!NO_EXTRA_BUFFER_SWAP) glutSwapBuffers(); @@ -948,7 +948,7 @@ void animation() draw_wave_highres_palette(2, phi, psi, total_energy, total_flux, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 0, 1.0); else draw_wave_epalette(phi, psi, total_energy, total_flux, color_scale, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 0, 1.0); draw_billiard(0, 1.0); - if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0); + if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, CIRC_COLORBAR_B, 0, 1.0); if (PRINT_SPEED) print_speed(speed, 0, 1.0); if (PRINT_FREQUENCY) print_frequency(phase_shift, 0, 1.0); glutSwapBuffers(); @@ -961,7 +961,7 @@ void animation() draw_wave_highres_palette(2, phi, psi, total_energy, total_flux, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 1, fade_value); else draw_wave_epalette(phi, psi, total_energy, total_flux, color_scale, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 1, fade_value); draw_billiard(1, fade_value); - if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 1, fade_value); + if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, CIRC_COLORBAR_B, 1, fade_value); if (PRINT_SPEED) print_speed(speed, 1, fade_value); if (PRINT_FREQUENCY) print_frequency(phase_shift, 0, 1.0); glutSwapBuffers(); diff --git a/wave_common.c b/wave_common.c index 8eaf8da..a33d601 100644 --- a/wave_common.c +++ b/wave_common.c @@ -386,8 +386,7 @@ void draw_wave(double *phi[NX], double *psi[NX], short int *xy_in[NX], double sc glEnd (); } -void draw_wave_e(double *phi[NX], double *psi[NX], double *total_energy[NX], double *color_scale[NX], short int *xy_in[NX], - double scale, int time, int plot) +void draw_wave_e(double *phi[NX], double *psi[NX], double *total_energy[NX], double *color_scale[NX], short int *xy_in[NX], double scale, int time, int plot) /* draw the field, new version with total energy option */ { int i, j, iplus, iminus, jplus, jminus; diff --git a/wave_comparison.c b/wave_comparison.c index a1143d8..5118942 100644 --- a/wave_comparison.c +++ b/wave_comparison.c @@ -41,40 +41,47 @@ #include #include /* Sam Leffler's libtiff library. */ #include +#include -#define MOVIE 0 /* set to 1 to generate movie */ -#define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */ -#define SAVE_MEMORY 0 /* set to 1 to save memory when writing tiff images */ +#define MOVIE 1 /* set to 1 to generate movie */ +#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */ +#define SAVE_MEMORY 1 /* set to 1 to save memory when writing tiff images */ +#define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */ #define TIME_LAPSE 0 /* set to 1 to add a time-lapse movie at the end */ #define TIME_LAPSE_FACTOR 4 /* factor of time-lapse movie */ -// #define WINWIDTH 1920 /* window width */ -// #define WINHEIGHT 1000 /* window height */ +#define WINWIDTH 1920 /* window width */ +#define WINHEIGHT 1150 /* window height */ // #define NX 1920 /* number of grid points on x axis */ // #define NY 1000 /* number of grid points on y axis */ // #define YMID 500 /* mid point of display */ -// #define XMIN -0.5 -// #define XMAX 3.5 /* x interval */ -// #define YMIN -1.041666667 -// #define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */ +#define NX 3840 /* number of grid points on x axis */ +#define NY 2300 /* number of grid points on y axis */ +#define YMID 1150 /* mid point of display */ +#define XMIN -2.0 +#define XMAX 2.0 /* x interval */ +#define YMIN -1.197916667 +#define YMAX 1.197916667 /* y interval for 9/16 aspect ratio */ -#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 -0.5 -#define XMAX 3.5 /* x interval */ -#define YMIN -1.125 -#define YMAX 1.125 /* y interval for 9/16 aspect ratio */ +// #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 -2.0 +// #define XMAX 2.0 /* x interval */ +// #define YMIN -1.125 +// #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ + +#define HIGHRES 1 /* set to 1 if resolution of grid is double that of displayed image */ #define JULIA_SCALE 1.0 /* scaling for Julia sets */ /* Choice of the billiard table */ -#define B_DOMAIN 43 /* choice of domain shape, see list in global_pdes.c */ -#define B_DOMAIN_B 47 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 999 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN_B 999 /* choice of domain shape, see list in global_pdes.c */ #define CIRCLE_PATTERN 13 /* pattern of circles, see list in global_pdes.c */ #define CIRCLE_PATTERN_B 13 /* pattern of circles, see list in global_pdes.c */ @@ -123,13 +130,13 @@ /* Physical parameters of wave equation */ #define TWOSPEEDS 1 /* set to 1 to replace hardcore boundary by medium with different speed */ -#define OSCILLATE_LEFT 1 /* set to 1 to add oscilating boundary condition on the left */ +#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.004 /* frequency of periodic excitation */ #define AMPLITUDE 1.0 /* amplitude of periodic excitation */ -#define COURANT 0.02 /* Courant number */ -#define COURANTB 0.01154 /* Courant number in medium B */ +#define COURANT 0.045 /* Courant number */ +#define COURANTB 0.045 /* Courant number in medium B */ // #define COURANTB 0.005 /* Courant number in medium B */ // #define COURANTB 0.008 /* Courant number in medium B */ #define GAMMA 0.0 /* damping factor in wave equation */ @@ -148,17 +155,22 @@ /* Increasing COURANT speeds up the simulation, but decreases accuracy */ /* For similar wave forms, COURANT^2*GAMMA should be kept constant */ +#define ADD_OSCILLATING_SOURCE 1 /* set to 1 to add an oscillating wave source */ +#define OSCILLATING_SOURCE_PERIOD 2 /* period of oscillating source */ +#define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */ +#define NSOURCES 48 /* number of sources */ + /* Boundary conditions, see list in global_pdes.c */ -#define B_COND 0 +#define B_COND 4 /* Parameters for length and speed of simulation */ -#define NSTEPS 2500 /* number of frames of movie */ -// #define NSTEPS 3300 /* number of frames of movie */ +// #define NSTEPS 500 /* number of frames of movie */ +#define NSTEPS 2600 /* 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 20 /* time after which to start saving frames */ +#define INITIAL_TIME 0 /* time after which to start saving frames */ #define COMPUTE_ENERGIES 0 /* set to 1 to compute and print energies */ #define BOUNDARY_WIDTH 2 /* width of billiard boundary */ @@ -168,24 +180,25 @@ #define SLEEP2 1 /* final sleeping time */ #define MID_FRAMES 20 /* number of still frames between movies */ #define END_FRAMES 100 /* number of still frames at end of movie */ +#define FADE 1 /* set to 1 to fade at end of movie */ /* Parameters of initial condition */ -#define INITIAL_AMP 0.75 /* amplitude of initial condition */ -// #define INITIAL_VARIANCE 0.0003 /* variance of initial condition */ -// #define INITIAL_WAVELENGTH 0.015 /* wavelength of initial condition */ -#define INITIAL_VARIANCE 0.0003 /* variance of initial condition */ -#define INITIAL_WAVELENGTH 0.02 /* wavelength of initial condition */ +#define INITIAL_AMP 0.002 /* amplitude of initial condition */ +// #define INITIAL_AMP 0.002 /* amplitude of initial condition */ +#define INITIAL_VARIANCE 0.0003 /* variance of initial condition */ +#define INITIAL_WAVELENGTH 0.015 /* wavelength of initial condition */ /* Plot type, see list in global_pdes.c */ #define PLOT 0 -#define PLOT_B 1 +#define PLOT_B 3 /* Color schemes */ #define COLOR_PALETTE 18 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE_B 13 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* background */ #define BLACK_TEXT 1 /* set to 1 to write text in black instead of white */ @@ -198,7 +211,7 @@ #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 200.0 /* scaling factor for energy representation */ +#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 FLUX_SCALE 1.0e4 /* scaling factor for enegy flux represtnation */ @@ -213,7 +226,7 @@ #define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */ #define COLORBAR_RANGE 1.5 /* scale of color scheme bar */ -#define COLORBAR_RANGE_B 2.5 /* scale of color scheme bar for 2nd part */ +#define COLORBAR_RANGE_B 12.5 /* scale of color scheme bar for 2nd part */ #define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ @@ -238,16 +251,16 @@ #define POT_MAZE 7 #define POTENTIAL 0 #define MAZE_WIDTH 0.02 /* half width of maze walls */ -#define VARIABLE_IOR 0 /* set to 1 for a variable index of refraction */ +#define VARIABLE_IOR 1 /* set to 1 for a variable index of refraction */ #define IOR 7 /* choice of index of refraction, see list in global_pdes.c */ #define IOR_TOTAL_TURNS 1.5 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */ #define MANDEL_IOR_SCALE -0.05 /* parameter controlling dependence of IoR on Mandelbrot escape speed */ -#define COURANT 0.04 /* Courant number */ -#define COURANTB 0.0 /* Courant number in medium B */ -#define INITIAL_AMP 0.5 /* amplitude of initial condition */ -#define INITIAL_VARIANCE 0.0003 /* variance of initial condition */ -#define INITIAL_WAVELENGTH 0.015 /* wavelength of initial condition */ -#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */ +// #define COURANT 0.04 /* Courant number */ +// #define COURANTB 0.0 /* Courant number in medium B */ +// #define INITIAL_AMP 0.5 /* amplitude of initial condition */ +// #define INITIAL_VARIANCE 0.0003 /* variance of initial condition */ +// #define INITIAL_WAVELENGTH 0.015 /* wavelength of initial condition */ +// #define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */ #define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */ #define N_WAVE_PACKETS 15 /* number of wave packets */ /* end of constants only used by sub_wave and sub_maze */ @@ -261,6 +274,7 @@ double courant2, courantb2; /* Courant parameters squared */ +FILE *monitor_sources; /*********************/ /* animation part */ @@ -677,20 +691,32 @@ void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *x } -void draw_color_bar(int plot, double range) +// void draw_color_bar(int plot, double range) +// { +// if (ROTATE_COLOR_SCHEME) draw_color_scheme(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range); +// else draw_color_scheme(XMAX - 0.3, YMIN + 0.1, XMAX - 0.1, YMAX - 0.1, plot, -range, range); +// } + +void draw_color_bar_palette(int plot, double range, int palette, int fade, double fade_value) { - if (ROTATE_COLOR_SCHEME) draw_color_scheme(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range); - else draw_color_scheme(XMAX - 0.3, YMIN + 0.1, XMAX - 0.1, YMAX - 0.1, plot, -range, range); + double width = 0.14; +// double width = 0.2; + + if (ROTATE_COLOR_SCHEME) + draw_color_scheme_palette_fade(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette, fade, fade_value); + else + draw_color_scheme_palette_fade(XMAX - 1.5*width, YMIN + 0.1, XMAX - 0.5*width, YMAX - 0.1, plot, -range, range, palette, fade, fade_value); } - - void animation() { - double time, scale, energies[6], top_energy, bottom_energy; - double *phi[NX], *psi[NX], *tmp[NX]; + double time, scale, energies[6], top_energy, bottom_energy, omega, angle, fade_value; + double *phi[NX], *psi[NX], *tmp[NX], *total_energy[NX]; short int *xy_in[NX]; - int i, j, s, counter = 0; + int i, j, s, counter = 0, k, first_source = 1, fade, resol = HIGHRES + 1; + t_wave_source wave_source[NSOURCES]; + + monitor_sources = fopen("monitor_sources.dat", "w"); /* Since NX and NY are big, it seemed wiser to use some memory allocation here */ for (i=0; i 0.0) */ + { + wave_source[k].phase += 0.06; +// if (k==1) printf("x = %.3lg, phase = %.3lg\n", wave_source[k].xc, wave_source[k].phase); +// fprintf(monitor_sources, "x = %.3lg, y = %.3lg, phase = %.3lg\n", wave_source[k].xc, wave_source[k].yc, wave_source[k].phase); + if (wave_source[k].phase > 1.0) + { + add_circular_wave_comp((double)wave_source[k].sign*wave_source[k].amp, wave_source[k].xc, wave_source[k].yc, phi, psi, xy_in, (wave_source[k].yc > 0)); + fprintf(monitor_sources, "x = %.3lg, y = %.3lg, phase = %.3lg\n", wave_source[k].xc, wave_source[k].yc, wave_source[k].phase); + printf("Adding wave at (%.2lg, %.2lg)\n", wave_source[k].xc, wave_source[k].yc); + fprintf(monitor_sources, "Adding wave at (%.2lg, %.2lg)\n", wave_source[k].xc, wave_source[k].yc); + wave_source[k].phase -= 1.0; + wave_source[k].sign *= -1; + } + } + } + draw_billiard_comp(); - - + if (COMPUTE_ENERGIES) { compute_energy_tblr(phi, psi, xy_in, energies); @@ -785,37 +874,39 @@ void animation() print_energies(energies, top_energy, bottom_energy); } - for (j=0; j= INITIAL_TIME) + if (i >= INITIAL_TIME) save_frame(); + + if ((i >= INITIAL_TIME)&&(DOUBLE_MOVIE)) { - save_frame(); +// save_frame(); if ((TIME_LAPSE)&&((i - INITIAL_TIME)%TIME_LAPSE_FACTOR == 0)) { save_frame_counter(NSTEPS + END_FRAMES + (i - INITIAL_TIME)/TIME_LAPSE_FACTOR); + if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0); counter++; } - else if (DOUBLE_MOVIE) + else /*if (DOUBLE_MOVIE)*/ { - draw_wave_comp(phi, psi, xy_in, scale, i, PLOT_B); +// draw_wave_comp(phi, psi, xy_in, scale, i, PLOT_B); + draw_wave_comp_highres_palette(resol, phi, psi, total_energy, xy_in, scale, i, PLOT_B, COLOR_PALETTE_B, 0, 1.0); draw_billiard_comp(); - if (DRAW_COLOR_SCHEME) draw_color_bar(PLOT_B, COLORBAR_RANGE_B); + if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0); glutSwapBuffers(); save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter); counter++; } } - else printf("Initial phase time %i of %i\n", i, INITIAL_TIME); + else if (NO_EXTRA_BUFFER_SWAP) glutSwapBuffers(); + +// else printf("Initial phase time %i of %i\n", i, INITIAL_TIME); @@ -835,20 +926,42 @@ void animation() { if (DOUBLE_MOVIE) { - draw_wave_comp(phi, psi, xy_in, scale, i, PLOT); +// draw_wave_comp(phi, psi, xy_in, scale, i, PLOT); + draw_wave_comp_highres_palette(resol, phi, psi, total_energy, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 0, 1.0); draw_billiard_comp(); - if (DRAW_COLOR_SCHEME) draw_color_bar(PLOT, COLORBAR_RANGE); + if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0); glutSwapBuffers(); } - for (i=0; i