-
Notifications
You must be signed in to change notification settings - Fork 44
/
Copy pathrde.c
2581 lines (2257 loc) · 107 KB
/
rde.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/*********************************************************************************/
/* */
/* Animation of reaction-diffusion equation in a planar domain */
/* */
/* N. Berglund, January 2022 */
/* */
/* Feel free to reuse, but if doing so it would be nice to drop a */
/* line to [email protected] - Thanks! */
/* */
/* compile with */
/* gcc -o rde rde.c */
/* -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -O3 -fopenmp */
/* */
/* OMP acceleration may be more effective after executing */
/* export OMP_NUM_THREADS=2 in the shell before running the program */
/* */
/* To make a video, set MOVIE to 1 and create subfolder tif_bz */
/* It may be possible to increase parameter PAUSE */
/* */
/* create movie using */
/* ffmpeg -i wave.%05d.tif -vcodec libx264 wave.mp4 */
/* */
/*********************************************************************************/
/*********************************************************************************/
/* */
/* NB: The algorithm used to simulate the wave equation is highly paralellizable */
/* One could make it much faster by using a GPU */
/* */
/*********************************************************************************/
#include <math.h>
#include <string.h>
#include <GL/glut.h>
#include <GL/glu.h>
#include <unistd.h>
#include <sys/types.h>
#include <tiffio.h> /* Sam Leffler's libtiff library. */
#include <omp.h>
#include <time.h>
#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 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 */
/* General geometrical parameters */
#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 HRES 6 /* factor for high resolution plots */
#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 */
/* Choice of simulated equation */
#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 */
#define SPHERE 1 /* set to 1 to simulate equation on sphere */
#define DPOLE 0 /* safety distance to poles */
#define DSMOOTH 1 /* size of neighbourhood of poles that are smoothed */
#define SMOOTHPOLE 0.05 /* smoothing coefficient at poles */
#define SMOOTHCOTPOLE 0.05 /* smoothing coefficient of cotangent at poles */
#define PHISHIFT 0.0 /* shift of phi in 2D plot (in degrees) */
#define SMOOTHBLOCKS 1 /* set to 1 to use blocks of points near the poles */
#define BLOCKDIST 64 /* distance to poles where points are blocked */
#define ZERO_MERIDIAN 190.0 /* choice of zero meridian (will be at left/right boundary of 2d plot) */
#define POLE_NODRAW 2 /* distance around poles where wave is not drawn */
#define ADD_POTENTIAL 0 /* set to 1 to add a potential (for Schrodinger equation) */
#define ADD_MAGNETIC_FIELD 0 /* set to 1 to add a magnetic field (for Schrodinger equation) - then set POTENTIAL 1 */
#define ADD_FORCE_FIELD 1 /* 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 6 /* 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 VARIABLE_DEPTH 1 /* set to 1 for variable depth in shallow water equation */
#define SWATER_DEPTH 10 /* variable depth in shallow water equation */
#define ANTISYMMETRIZE_WAVE_FCT 0 /* set tot 1 to make wave function antisymmetric */
#define ADAPT_STATE_TO_BC 1 /* to smoothly adapt initial state to obstacles */
#define OBSTACLE_GEOMETRY 84 /* geometry of obstacles, as in B_DOMAIN */
#define BC_STIFFNESS 0.25 /* 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 */
#define JULIA_ROT -20.0 /* rotation of Julia set, in degrees */
#define JULIA_RE 0.5
#define JULIA_IM 0.462 /* parameters for Julia sets */
/* Choice of the billiard table */
#define B_DOMAIN 999 /* choice of domain shape, 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 PDISC_FACTOR 3.25 /* controls density of Poisson disc process (default: 3.25) */
#define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
#define LAMBDA 1.0 /* parameter controlling the dimensions of domain */
#define MU 1.0 /* 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 */
#define MRATIO 5 /* ratio defining Menger gasket */
#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */
#define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */
#define FOCI 1 /* set to 1 to draw focal points of ellipse */
#define NGRIDX 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 WALL_WIDTH 0.05 /* width of wall separating lenses */
#define RADIUS_FACTOR 0.3 /* controls inner radius for C_RING arrangements */
#define X_SHOOTER -0.2
#define Y_SHOOTER -0.6
#define X_TARGET 0.4
#define Y_TARGET 0.7 /* shooter and target positions in laser fight */
#define ISO_XSHIFT_LEFT -1.65
#define ISO_XSHIFT_RIGHT 0.4
#define ISO_YSHIFT_LEFT -0.05
#define ISO_YSHIFT_RIGHT -0.05
#define ISO_SCALE 0.85 /* coordinates for isospectral billiards */
/* You can add more billiard tables by adapting the functions */
/* xy_in_billiard and draw_billiard in sub_wave.c */
/* Physical parameters of wave equation */
#define DT 0.00000025
#define VISCOSITY 0.02
#define POISSON_STIFFNESS 1.0 /* stiffness of Poisson equation solver for incompressible Euler */
#define DISSIPATION 0.0
#define DISSIPATION_EXT 1.0e-3 /* dissipation of shallow water eq. outside domain */
#define RPSA 0.75 /* parameter in Rock-Paper-Scissors-type interaction */
#define RPSLZB 0.0 /* second parameter in Rock-Paper-Scissors-Lizard-Spock type interaction */
#define K_AC 0.1 /* force constant in Allen-Cahn equation */
#define EPSILON 0.8 /* time scale separation */
#define DELTA 0.1 /* time scale separation */
#define FHNA 1.0 /* parameter in FHN equation */
#define FHNC -0.01 /* parameter in FHN equation */
#define K_HARMONIC 1.0 /* spring constant of harmonic potential */
#define K_COULOMB 0.5 /* constant in Coulomb potential */
#define V_MAZE 0.4 /* potential in walls of maze */
#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.025 /* gravity/Coriolis force */
#define BC_FIELD 1.0e-6 /* 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 */
#define C_EULER_COMP 0.1 /* constant in compressible Euler equation */
#define SWATER_VARDEPTH_FACTOR 0.1 /* force constant in front of variable depth term */
#define SWATER_CORIOLIS_FORCE 0.002 /* Coriolis force for shallow water equation */
#define SMOOTHEN_VORTICITY 0 /* set to 1 to smoothen vorticity field in Euler equation */
#define SMOOTHEN_VELOCITY 1 /* set to 1 to smoothen velocity field in Euler equation */
#define SMOOTHEN_PERIOD 7 /* period between smoothenings */
#define SMOOTH_FACTOR 0.2 /* factor by which to smoothen */
#define ADD_OSCILLATING_SOURCE 0 /* set to 1 to add an oscillating wave source */
#define OSCILLATING_SOURCE_PERIOD 1 /* period of oscillating source */
#define OSCILLATING_SOURCE_OMEGA 0.2 /* frequency of oscillating source */
#define ADD_TRACERS 1 /* set to 1 tof add tracer particles (for Euler equations) */
#define N_TRACERS 2000 /* number of tracer particles */
#define TRACERS_STEP 0.1 /* step size in tracer evolution */
#define RESPAWN_TRACERS 1 /* set to 1 to randomly move tracer position */
#define RESPAWN_PROBABILTY 0.005 /* probability of moving tracer */
#define T_OUT 2.0 /* outside temperature */
#define T_IN 0.0 /* inside temperature */
#define SPEED 0.0 /* speed of drift to the right */
#define ADD_NOISE 0 /* set to 1 to add noise, set to 2 to add noise in right half */
#define NOISE_INTENSITY 0.01 /* noise intensity */
#define CHANGE_NOISE 0 /* set to 1 to increase noise intensity */
#define NOISE_FACTOR 40.0 /* factor by which to increase noise intensity */
#define NOISE_INITIAL_TIME 100 /* initial time during which noise remains constant */
#define CHANGE_VISCOSITY 0 /* set to 1 to change the viscosity in the course of the simulation */
#define ADJUST_INTSTEP 0 /* set to 1 to decrease integration step when viscosity increases */
#define VISCOSITY_INITIAL_TIME 10 /* initial time during which viscosity remains constant */
#define VISCOSITY_FACTOR 100.0 /* factor by which to change viscosity */
#define VISCOSITY_MAX 2.0 /* max value of viscosity beyond which NVID is increased and integration step is decrase, for numerical stability */
#define CHANGE_RPSLZB 0 /* set to 1 to change second parameter in Rock-Paper-Scissors-Lizard-Spock equation */
#define RPSLZB_CHANGE 0.75 /* factor by which to rpslzb parameter */
#define RPSLZB_INITIAL_TIME 0 /* initial time during which rpslzb remains constant */
#define RPSLZB_FINAL_TIME 500 /* final time during which rpslzb remains constant */
#define CHANGE_FLOW_SPEED 0 /* set to 1 to change speed of laminar flow */
#define IN_OUT_FLOW_BC 0 /* type of in-flow/out-flow boundary conditions for Euler equation, 0 for no b.c. */
#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.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.2 /* amplitude of pressure gradient for Euler equation */
#define SWATER_MIN_HEIGHT 0.05 /* min height of initial condition for shallow water equation */
#define DEPTH_FACTOR 0.00002 /* proportion of min height in variable depth */
#define TANH_FACTOR 5.0 /* steepness of variable depth */
#define EULER_GRADIENT_YSHIFT 0.0 /* y-shift in computation of gradient in Euler equation */
#define ADD_MOON_FORCING 1 /* set to 1 to simulate tidal forces from Moon */
#define FORCING_AMP 5.0e-6 /* amplitude of periodic forcing */
#define FORCING_CONST_AMP 0.0 /* amplitude of periodic forcing */
#define FORCING_PERIOD 1600 /* period of forcing */
#define FORCING_SHIFT 1.5 /* phase shift of forcing in units of Pi */
/* Boundary conditions, see list in global_pdes.c */
#define B_COND 1
#define B_COND_LEFT 0
#define B_COND_RIGHT 0
#define B_COND_TOP 0
#define B_COND_BOTTOM 0
/* Parameters for length and speed of simulation */
#define NSTEPS 1400 /* number of frames of movie */
#define NVID 80 /* 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 999 /* number of segments of boundary */
#define BOUNDARY_WIDTH 2 /* width of billiard boundary */
#define PAUSE 100 /* number of frames after which to pause */
#define PSLEEP 2 /* sleep time during pause */
#define SLEEP1 2 /* initial sleeping time */
#define SLEEP2 1 /* final sleeping time */
#define INITIAL_TIME 0 /* initial still time */
#define MID_FRAMES 100 /* number of still frames between parts of two-part movie */
#define END_FRAMES 250 /* number of still frames at end of movie */
#define FADE 1 /* set to 1 to fade at end of movie */
/* Visualisation */
#define PLOT_3D 0 /* controls whether plot is 2D or 3D */
#define PLOT_SPHERE 1 /* draws fields on a sphere */
#define ROTATE_VIEW 1 /* set to 1 to rotate position of observer */
#define ROTATE_ANGLE -45.0 /* total angle of rotation during simulation */
#define SHADE_3D 0 /* set to 1 to change luminosity according to normal vector */
#define SHADE_2D 1 /* set to 1 to change luminosity according to normal vector */
#define VIEWPOINT_TRAJ 1 /* type of viewpoint trajectory */
#define MAX_LATITUDE 45.0 /* maximal latitude for viewpoint trajectory VP_ORBIT2 */
#define DRAW_PERIODICISED 0 /* set to 1 to repeat wave periodically in x and y directions */
#define DRAW_MOON_POSITION 0 /* set to 1 to draw vertical lines for tidal sim */
/* Plot type - color scheme */
#define CPLOT 70
#define CPLOT_B 74
/* Plot type - height of 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 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 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 DRAW_DEPTH 0 /* 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 FLOODING 1 /* set to 1 for drawing water when higher than continents */
// #define FLOODING_VSHIFT 0.51 /* controls when wave is considered higher than land */
#define FLOODING_VSHIFT 0.56 /* controls when wave is considered higher than land */
#define FLOODING_VSHIFT_2D 0.56 /* controls when wave is considered higher than land */
// #define FLOODING_VSHIFT_2D 0.61 /* controls when wave is considered higher than land */
#define PLOT_SCALE_ENERGY 0.05 /* vertical scaling in energy plot */
#define PRINT_TIME 0 /* set to 1 to print running time */
#define PRINT_VISCOSITY 0 /* set to 1 to print viscosity */
#define PRINT_RPSLZB 0 /* set to 1 to print rpslzb parameter */
#define PRINT_PROBABILITIES 0 /* set to 1 to print probabilities (for Ehrenfest urn configuration) */
#define PRINT_NOISE 0 /* set to 1 to print noise intensity */
#define PRINT_FLOW_SPEED 0 /* set to 1 to print speed of flow */
#define PRINT_AVERAGE_SPEED 0 /* set to 1 to print average speed of flow */
#define PRINT_LEFT 1 /* set to 1 to print parameters at left side */
#define DRAW_FIELD_LINES 0 /* set to 1 to draw field lines */
#define FIELD_LINE_WIDTH 1 /* width of field lines */
#define N_FIELD_LINES 120 /* number of field lines */
#define FIELD_LINE_FACTOR 120 /* factor controlling precision when computing origin of field lines */
#define DRAW_BILLIARD 1 /* set to 1 to draw boundary */
#define DRAW_BILLIARD_FRONT 1 /* set to 1 to draw boundary */
#define FILL_BILLIARD_COMPLEMENT 1 /* set to 1 to fill complement of billiard (for certain shapes only) */
/* 3D representation */
#define REPRESENTATION_3D 1 /* choice of 3D representation */
#define REP_AXO_3D 0 /* linear projection (axonometry) */
#define REP_PROJ_3D 1 /* projection on plane orthogonal to observer line of sight */
/* Color schemes, see list in global_pdes.c */
#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */
#define COLOR_PALETTE_B 16 /* Color palette, see list in global_pdes.c */
#define BLACK 1 /* black background */
#define COLOR_OUT_R 1.0 /* color outside domain */
#define COLOR_OUT_G 1.0
#define COLOR_OUT_B 1.0
#define COLOR_SCHEME 3 /* choice of color scheme */
#define PHASE_SHIFT -0.25 /* phase shift of color scheme, in units of Pi (formerly COLOR_PHASE_SHIFT) */
#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 VSHIFT_AMPLITUDE 0.0 /* additional shift for wave amplitude */
#define VSCALE_AMPLITUDE 15.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */
#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */
#define CURL_SCALE 1.0 /* 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 10.0 /* sensitivity of luminosity on module, for color scheme Z_ARGUMENT */
#define MIN_SCHROD_LUM 0.1 /* 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 15.0 /* vertical scaling of water height */
#define ADD_HEIGHT_CONSTANT -0.016 /* constant added to wave height */
#define SHADE_SCALE_2D 0.25 /* controls "depth" of 2D shading */
#define COLORHUE 260 /* initial hue of water color for scheme C_LUM */
#define COLORDRIFT 0.0 /* how much the color hue drifts during the whole simulation */
#define LUMMEAN 0.5 /* amplitude of luminosity variation for scheme C_LUM */
#define LUMAMP 0.3 /* amplitude of luminosity variation for scheme C_LUM */
#define HUEMEAN 359.0 /* mean value of hue for color scheme C_HUE */
#define HUEAMP -359.0 /* amplitude of variation of hue for color scheme C_HUE */
#define E_SCALE 100.0 /* scaling factor for energy representation */
#define FLUX_SCALE 100.0 /* scaling factor for energy representation */
#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 50.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 1.0 /* shift for color scheme Z_EULER_DENSITY */
#define VSCALE_DENSITY 30.0 /* additional scaling factor for color scheme Z_EULER_DENSITY */
#define VSCALE_VORTICITY 15.0 /* additional scaling factor for color scheme Z_EULERC_VORTICITY */
#define VORTICITY_SHIFT 0.0 /* vertical shift of vorticity */
#define ZSCALE_SPEED 300.0 /* additional scaling factor for z-coord Z_EULER_SPEED and Z_SWATER_SPEED */
#define ZSHIFT_SPEED 0.0 /* additional shift of z-coord Z_EULER_SPEED and Z_SWATER_SPEED */
#define ZSCALE_NORMGRADIENT -0.0001 /* vertical scaling for Z_NORM_GRADIENT */
#define VSCALE_SWATER 100.0 /* additional scaling factor for color scheme Z_EULER_DENSITY */
#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.04 /* half width of maze walls */
#define DRAW_COLOR_SCHEME 0 /* set to 1 to plot the color scheme */
#define COLORBAR_RANGE 3.0 /* 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 */
#define VARIABLE_IOR 0 /* set to 1 for a variable index of refraction */
#define IOR 4 /* 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 OMEGA 0.005 /* frequency of periodic excitation */
#define OSCIL_YMAX 0.2 /* defines oscillation range */
#define COURANT 0.08 /* Courant number */
#define COURANTB 0.03 /* Courant number in medium B */
#define INITIAL_AMP 0.5 /* amplitude of initial condition */
#define INITIAL_VARIANCE 0.0002 /* variance of initial condition */
#define INITIAL_WAVELENGTH 0.1 /* wavelength of initial condition */
#define VSCALE_ENERGY 200.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 AVRG_E_FACTOR 0.99 /* controls time window size in P_AVERAGE_ENERGY scheme */
#define OSCILLATION_SCHEDULE 0 /* oscillation schedule, see list in global_pdes.c */
#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 COMPARISON 0 /* set to 1 to compare two different patterns (beta) */
#define B_DOMAIN_B 20 /* second domain shape, for comparisons */
#define CIRCLE_PATTERN_B 0 /* second pattern of circles or polygons */
#define FLUX_WINDOW 20 /* averaging window for energy flux */
#define ADD_WAVE_PACKET_SOURCES 0 /* set to 1 to add several sources emitting wave packets */
#define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */
#define N_SOURCES 1 /* number of wave sources */
#define N_WAVE_PACKETS 15 /* number of wave packets */
#define WAVE_PACKET_RADIUS 20 /* radius of wave packets */
#define OSCIL_LEFT_YSHIFT 25.0 /* y-dependence of left oscillation (for non-horizontal waves) */
#define INITIAL_SHIFT 20.0 /* time shift of initial wave packet (in oscillation periods) */
#define WAVE_PACKET_SHIFT 200.0 /* time shift between wave packets (in oscillation periods) */
#define DRAW_WAVE_PROFILE 0 /* set to 1 to draw a profile of the wave */
#define HORIZONTAL_WAVE_PROFILE 0 /* set to 1 to draw wave profile vertically */
#define VERTICAL_WAVE_PROFILE 0 /* set to 1 to draw wave profile vertically */
#define WAVE_PROFILE_X 1.9 /* value of x to sample wave profile */
#define WAVE_PROFILE_Y -1.0 /* value of y to sample wave profile */
#define PROFILE_AT_BOTTOM 1 /* draw wave profile at bottom instead of top */
#define AVERAGE_WAVE_PROFILE 1 /* set to 1 to draw time-average of wave profile squared*/
#define MU_B 1.0 /* parameter controlling the dimensions of domain */
#define GAMMA 0.0 /* damping factor in wave equation */
#define GAMMAB 0.0 /* damping factor in wave equation */
#define VERTICAL_WAVE_PROFILE 0 /* set to 1 to draw wave profile vertically */
#define DRAW_WAVE_TIMESERIES 0 /* set to 1 to draw a time series of the wave */
#define TIMESERIES_NVALUES 400 /* number of values plotted in time series */
#define DRAW_WAVE_SOURCE 0 /* set to 1 to draw source of wave at (wave_source_x, wave_source_y) */
#define MESSAGE_LDASH 1 /* length of dash for Morse code message */
#define MESSAGE_LDOT 1 /* length of dot for Morse code message */
#define MESSAGE_LINTERVAL 1 /* length of interval between dashes/dots for Morse code message */
#define MESSAGE_LINTERLETTER 1 /* length of interval between letters for Morse code message */
#define MESSAGE_LSPACE 1 /* length of space for Morse code message */
#define MESSAGE_INITIAL_TIME 1 /* initial time before starting message for Morse code message */
#define IMAGE_FILE 0 /* for option D_IMAGE */
/* end of constants added only for compatibility with wave_common.c */
double u_3d[2] = {0.75, -0.45}; /* projections of basis vectors for REP_AXO_3D representation */
double v_3d[2] = {-0.75, -0.45};
double w_3d[2] = {0.0, 0.015};
double light[3] = {0.40824829, -0.816496581, 0.40824829}; /* vector of "light" direction for P_3D_ANGLE color scheme */
double observer[3] = {0.0, -6.0, 2.5}; /* location of observer for REP_PROJ_3D representation */
int reset_view = 0; /* switch to reset 3D view parameters (for option ROTATE_VIEW) */
/* constants for simulations on planets */
#define ADD_DEM 1 /* add DEM (digital elevation model) */
#define ADD_NEGATIVE_DEM 1 /* add DEM with bathymetric data */
#define RSCALE_DEM 0.075 /* scaling factor of radial component for DEM */
#define SMOOTH_DEM 1 /* set to 1 to smoothen DEM (to make altitude less constant) */
#define DEM_SMOOTH_STEPS 10 /* number of smoothening steps */
#define DEM_SMOOTH_HEIGHT 2.0 /* relative height below which to smoothen */
#define DEM_MAXHEIGHT 9000.0 /* max height of DEM (estimated from Everest/Olympus Mons) */
#define DEM_MAXDEPTH -10000 /* max depth of DEM */
#define PLANET_SEALEVEL 0.0 /* sea level for flooded planet */
#define VENUS_NODATA_FACTOR 0.5 /* altitude to assign to DEM points without data (fraction of mean altitude) */
#define Z_SCALING_FACTOR 0.8 /* overall scaling factor of z axis for REP_PROJ_3D representation */
#define XY_SCALING_FACTOR 2.1 /* 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.0 /* 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 */
#define DRAW_ARROW 0 /* set to 1 to draw arrow above sphere */
#define RSCALE 0.15 /* scaling factor of radial component */
#define RSHIFT -0.075 /* shift in radial component */
#define RMAX 2.0 /* max value of radial component */
#define RMIN 0.5 /* min value of radial component */
#define COS_VISIBLE -0.35 /* limit on cosine of normal to shown facets */
/* For debugging purposes only */
#define FLOOR 1 /* set to 1 to limit wave amplitude to VMAX */
#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)||(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))
#define ASYM_SPEED_COLOR (VMEAN_SPEED == 0.0)
#define XYIN_INITIALISED (B_DOMAIN == D_IMAGE)
int block_sizes[NY]; /* table of block sizes for blocking around poles */
int block_numbers[NY]; /* table of block numbers for blocking around poles */
#include "global_pdes.c"
#include "global_3d.c" /* constants and global variables */
#include "sub_maze.c"
#include "sub_wave.c"
#include "wave_common.c" /* common functions for wave_billiard, wave_comparison, etc */
#include "sub_wave_3d_rde.c" /* should be later replaced by sub_wave_rde.c */
#include "sub_rde.c"
double f_aharonov_bohm(double r2)
/* radial part of Aharonov-Bohm vector potential */
{
double r02 = AB_RADIUS*AB_RADIUS;
if (r2 > r02) return(-0.25*r02/r2);
else return(0.25*(r2 - 2.0*r02)/r02);
// if (r2 > r02) return(1.0/r2);
// else return((2.0*r02 - r2)/(r02*r02));
}
double potential(int i, int j)
/* compute potential (e.g. for Schrödinger equation), or potential part if there is a magnetic field */
{
double x, y, xy[2], r, small = 1.0e-1, kx, ky, lx = XMAX - XMIN, r1, r2, r3, f;
int rect;
ij_to_xy(i, j, xy);
x = xy[0];
y = xy[1];
switch (POTENTIAL) {
case (POT_HARMONIC):
{
return (K_HARMONIC*(x*x + y*y));
}
case (POT_COULOMB):
{
// r = module2(x, y);
r = sqrt(x*x + y*y + small*small);
// if (r < small) r = small;
return (-K_COULOMB/r);
}
case (POT_PERIODIC):
{
kx = 4.0*DPI/(XMAX - XMIN);
ky = 2.0*DPI/(YMAX - YMIN);
return(-K_HARMONIC*cos(kx*x)*cos(ky*y));
}
case (POT_FERMIONS):
{
r = sqrt((x-y)*(x-y) + small*small);
return (-K_COULOMB/r);
}
case (POT_FERMIONS_PERIODIC):
{
r1 = sqrt((x-y)*(x-y) + small*small);
r2 = sqrt((x-lx-y)*(x-lx-y) + small*small);
r3 = sqrt((x+lx-y)*(x+lx-y) + small*small);
// r = r/3.0;
return (-0.5*K_COULOMB*(1.0/r1 + 1.0/r2 + 1.0/r3));
}
case (VPOT_CONSTANT_FIELD):
{
return (K_HARMONIC*(x*x + y*y)); /* magnetic field strength b is chosen such that b^2 = K_HARMONIC */
}
case (VPOT_AHARONOV_BOHM):
{
r2 = x*x + y*y;
f = f_aharonov_bohm(r2);
return (B_FIELD*B_FIELD*f*f*r2); /* magnetic field strength b is chosen such that b^2 = K_HARMONIC */
// return (K_HARMONIC*f); /* magnetic field strength b is chosen such that b^2 = K_HARMONIC */
}
case (POT_MAZE):
{
for (rect=0; rect<npolyrect; rect++)
if (ij_in_polyrect(i, j, polyrect[rect])) return(V_MAZE);
return(0.0);
}
default:
{
return(0.0);
}
}
}
void compute_vector_potential(int i, int j, double *ax, double *ay)
/* initialize the vector potential, for Schrodinger equation in a magnetic field */
{
double x, y, xy[2], r2, f;
ij_to_xy(i, j, xy);
x = xy[0];
y = xy[1];
switch (POTENTIAL) {
case (VPOT_CONSTANT_FIELD):
{
*ax = B_FIELD*y;
*ay = -B_FIELD*x;
break;
}
case (VPOT_AHARONOV_BOHM):
{
r2 = x*x + y*y;
f = f_aharonov_bohm(r2);
*ax = B_FIELD*y*f;
*ay = -B_FIELD*x*f;
break;
}
default:
{
*ax = 0.0;
*ay = 0.0;
}
}
}
void compute_gfield(int i, int j, double *gx, double *gy)
/* initialize the exterior field, for the compressible Euler equation */
{
double x, y, xy[2], r, f, a = 0.4, x1, y1, hx, hy, h;
ij_to_xy(i, j, xy);
x = xy[0];
y = xy[1];
switch (FORCE_FIELD) {
case (GF_VERTICAL):
{
*gx = 0.0;
*gy = -G_FIELD;
break;
}
case (GF_CIRCLE):
{
r = module2(x,y) + 1.0e-2;
f = 0.5*(1.0 - tanh(BC_STIFFNESS*(r - LAMBDA)));
*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 = BC_FIELD*f*x/(LAMBDA*LAMBDA);
*gy = BC_FIELD*f*y/(MU*MU);
break;
}
case (GF_AIRFOIL):
{
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 = BC_FIELD*f*(x/(LAMBDA*LAMBDA) + a*y1/(MU*MU));
*gy = BC_FIELD*f*y1/(MU*MU);
break;
}
case (GF_WING):
{
if (x >= LAMBDA)
{
*gx = 0.0;
*gy = 0.0;
}
else
{
x1 = 1.0 - x/LAMBDA;
if (x1 < 0.1) x1 = 0.1;
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 = 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);
// hy = y1/(MU*MU*x1*x1);
// h = module2(hx, hy) + 1.0e-2;
// *gx = G_FIELD*f*hx/h;
// *gy = G_FIELD*f*hy/h;
}
break;
}
default:
{
*gx = 0.0;
*gy = 0.0;
}
}
}
void initialize_potential(double potential_field[NX*NY])
/* initialize the potential field, e.g. for the Schrödinger equation */
{
int i, j;
#pragma omp parallel for private(i,j)
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
potential_field[i*NY+j] = potential(i,j);
}
}
}
void initialize_vector_potential(double vpotential_field[2*NX*NY])
/* initialize the potential field, e.g. for the Schrödinger equation */
{
int i, j;
#pragma omp parallel for private(i,j)
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
compute_vector_potential(i, j, &vpotential_field[i*NY+j], &vpotential_field[NX*NY+i*NY+j]);
}
}
}
void initialize_gfield(double gfield[2*NX*NY], double bc_field[NX*NY], double bc_field2[NX*NY], t_wave_sphere *wsphere, t_wave_sphere *wsphere_hr)
/* initialize the exterior field, e.g. for the compressible Euler equation */
{
int i, j;
double dx, dy;
switch (FORCE_FIELD)
{
case (GF_COMPUTE_FROM_BC):
{
dx = (XMAX - XMIN)/(double)NX;
dy = (YMAX - YMIN)/(double)NY;
#pragma omp parallel for private(i,j)
for (i=1; i<NX-1; i++){
for (j=1; j<NY-1; j++){
gfield[i*NY+j] = BC_FIELD*(bc_field2[(i+1)*NY+j] - bc_field2[(i-1)*NY+j])/dx;
gfield[NX*NY+i*NY+j] = BC_FIELD*(bc_field2[i*NY+j+1] - bc_field2[i*NY+j-1])/dy;
// printf("gfield at (%i,%i): (%.3lg, %.3lg)\n", i, j, gfield[i*NY+j], gfield[NX*NY+i*NY+j]);
}
}
/* boundaries */
for (i=0; i<NX; i++)
{
gfield[i*NY] = 0.0;
gfield[NX*NY+i*NY] = 0.0;
gfield[i*NY+NY-1] = 0.0;
gfield[NX*NY+i*NY+NY-1] = 0.0;
}
for (j=0; j<NY; j++)
{
gfield[j] = 0.0;
gfield[NX*NY+j] = 0.0;
gfield[(NX-1)*NY+j] = 0.0;
gfield[NX*NY+(NX-1)*NY+j] = 0.0;
}
break;
}
case (GF_EARTH):
{
dx = (XMAX - XMIN)/(double)NX;
dy = (YMAX - YMIN)/(double)NY;
init_earth_map_rde(wsphere, 1);
init_earth_map_rde(wsphere_hr, HRES);
#pragma omp parallel for private(i,j)
for (i=1; i<NX-1; i++){
for (j=1; j<NY-1; j++)
{
if (wsphere[i*NY+j].altitude > 0.0)
{
gfield[i*NY+j] = -BC_FIELD*(wsphere[(i+1)*NY+j].altitude - wsphere[(i-1)*NY+j].altitude)/dx;
gfield[NX*NY+i*NY+j] = -BC_FIELD*(wsphere[i*NY+j+1].altitude - wsphere[i*NY+j-1].altitude)/dy;
}
else
{
gfield[i*NY+j] = 0.0;
gfield[NX*NY+i*NY+j] = 0.0;
}
}
}
/* boundaries TODO */
for (i=0; i<NX; i++)
{
gfield[i*NY] = 0.0;
gfield[NX*NY+i*NY] = 0.0;
gfield[i*NY+NY-1] = 0.0;
gfield[NX*NY+i*NY+NY-1] = 0.0;
}
for (j=0; j<NY; j++)
{
gfield[j] = 0.0;
gfield[NX*NY+j] = 0.0;
gfield[(NX-1)*NY+j] = 0.0;
gfield[NX*NY+(NX-1)*NY+j] = 0.0;
}
break;
}
case (GF_MARS):
{
dx = (XMAX - XMIN)/(double)NX;
dy = (YMAX - YMIN)/(double)NY;
init_planet_map_rde(wsphere, D_SPHERE_MARS, 1);
init_planet_map_rde(wsphere_hr, D_SPHERE_MARS, HRES);
#pragma omp parallel for private(i,j)
for (i=1; i<NX-1; i++){
for (j=1; j<NY-1; j++){
gfield[i*NY+j] = BC_FIELD*(wsphere[(i+1)*NY+j].altitude - wsphere[(i-1)*NY+j].altitude)/dx;
gfield[NX*NY+i*NY+j] = BC_FIELD*(wsphere[i*NY+j+1].altitude - wsphere[i*NY+j-1].altitude)/dy;
}
}
/* boundaries TODO */
for (i=0; i<NX; i++)
{
gfield[i*NY] = 0.0;
gfield[NX*NY+i*NY] = 0.0;
gfield[i*NY+NY-1] = 0.0;
gfield[NX*NY+i*NY+NY-1] = 0.0;
}
for (j=0; j<NY; j++)
{
gfield[j] = 0.0;
gfield[NX*NY+j] = 0.0;
gfield[(NX-1)*NY+j] = 0.0;
gfield[NX*NY+(NX-1)*NY+j] = 0.0;
}
break;
}
case (GF_VENUS):
{
dx = (XMAX - XMIN)/(double)NX;
dy = (YMAX - YMIN)/(double)NY;
init_planet_map_rde(wsphere, D_SPHERE_VENUS, 1);
init_planet_map_rde(wsphere_hr, D_SPHERE_VENUS, HRES);
#pragma omp parallel for private(i,j)
for (i=1; i<NX-1; i++){
for (j=1; j<NY-1; j++){
gfield[i*NY+j] = BC_FIELD*(wsphere[(i+1)*NY+j].altitude - wsphere[(i-1)*NY+j].altitude)/dx;
gfield[NX*NY+i*NY+j] = BC_FIELD*(wsphere[i*NY+j+1].altitude - wsphere[i*NY+j-1].altitude)/dy;
}
}
/* boundaries TODO */
for (i=0; i<NX; i++)
{
gfield[i*NY] = 0.0;
gfield[NX*NY+i*NY] = 0.0;
gfield[i*NY+NY-1] = 0.0;
gfield[NX*NY+i*NY+NY-1] = 0.0;
}
for (j=0; j<NY; j++)
{
gfield[j] = 0.0;
gfield[NX*NY+j] = 0.0;
gfield[(NX-1)*NY+j] = 0.0;
gfield[NX*NY+(NX-1)*NY+j] = 0.0;
}
break;
}
default:
{
#pragma omp parallel for private(i,j)
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
compute_gfield(i, j, &gfield[i*NY+j], &gfield[NX*NY+i*NY+j]);
}
}
}
}
}
void compute_water_depth(int i, int j, t_rde *rde, t_wave_sphere wsphere[NX*NY], int ncircles)
/* initialize the vector potential, for Schrodinger equation in a magnetic field */
{
double x, y, xy[2], r0, r, z, h, h0, hh, x1, y1, z1, d, rmax, rmin, phi, vshift;
int n, nx, ny;
static double phi1, phi2, phi3, sq3, rico;
static int first = 1;
if (first)
{
phi = 0.5*(sqrt(5.0)+1.0);
phi1 = phi - 1.0;
phi2 = phi - 2.0;
phi3 = 2.0*phi - 3.0;
sq3 = 0.5/sqrt(3.0);
rico = 0.5/(2.0 + phi);
first = 0;
}
ij_to_xy(i, j, xy);
x = xy[0];
y = xy[1];
switch (SWATER_DEPTH) {
case (SH_CIRCLE):
{
r0 = module2(x,y);
r = r0/LAMBDA;
h = tanh(TANH_FACTOR*(r-1.0));
hh = 1.0 - h*h;
z = SWATER_MIN_HEIGHT*DEPTH_FACTOR*TANH_FACTOR*hh/(r0*LAMBDA);
rde->depth = 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;
}
case (SH_SPHERE_CUBE):
{
rmax = 0.0;
/* compute distance from cube to origin, given by 0.5/rmax */
if (vabs(wsphere[i*NY+j].x) > rmax) rmax = vabs(wsphere[i*NY+j].x);
if (vabs(wsphere[i*NY+j].y) > rmax) rmax = vabs(wsphere[i*NY+j].y);
if (vabs(wsphere[i*NY+j].z) > rmax) rmax = vabs(wsphere[i*NY+j].z);
h = 1.0 - 0.5/rmax;
// printf("h = %.3lg\n", h);
if (h < 0.0) h = 0.0;
rde->depth = SWATER_MIN_HEIGHT*DEPTH_FACTOR*h;
break;
}
case (SH_SPHERE_OCTAHEDRON):
{
rmax = 0.0;
/* compute distance from octahedron to origin, given by 0.5/rmax */
rmax = 1.0/(vabs(wsphere[i*NY+j].x) + vabs(wsphere[i*NY+j].y) + vabs(wsphere[i*NY+j].z));
h = 1.0 - 0.5/rmax;
printf("h = %.3lg\n", h);
if (h < 0.0) h = 0.0;
rde->depth = SWATER_MIN_HEIGHT*DEPTH_FACTOR*h;
break;
}
case (SH_SPHERE_DODECAHEDRON):
{
rmax = 0.0;
x1 = wsphere[i*NY+j].x;
y1 = wsphere[i*NY+j].y;
z1 = wsphere[i*NY+j].z;
/* compute distance from dodecahedron */
d = vabs(phi1*y1 - phi2*z1) - sq3;
if (d > rmax) rmax = d;
d = vabs(-phi2*x1 + phi1*z1) - sq3;
if (d > rmax) rmax = d;
d = vabs(phi1*x1 - phi2*y1) - sq3;
if (d > rmax) rmax = d;
d = vabs(phi1*y1 + phi2*z1) - sq3;
if (d > rmax) rmax = d;
d = vabs(phi2*x1 + phi1*z1) - sq3;
if (d > rmax) rmax = d;
d = vabs(phi1*x1 + phi2*y1) - sq3;
if (d > rmax) rmax = d;
h = rmax;
printf("h = %.3lg\n", h);
if (h < 0.0) h = 0.0;
rde->depth = SWATER_MIN_HEIGHT*DEPTH_FACTOR*h;
break;
}
case (SH_SPHERE_ICOSAHEDRON):
{
rmax = 0.0;
x1 = wsphere[i*NY+j].x;
y1 = wsphere[i*NY+j].y;