-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathSIMPLE_TS_MPI.cpp
12943 lines (10417 loc) · 396 KB
/
SIMPLE_TS_MPI.cpp
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
/*
============================================================================
Name : SIMPLE_TS_MPI
Author : Kiril S. Shterev
Version : v.1.1
Copyright : All rights reserved. The source code is freely distributed for non-commercial use.
Non-commercial use: Developers or distributors can compile, use all code or any part of the code, redistribute, sell results calculated using code or part of it and not-only, but except commercial use.
Commercial use : It is consider any use of the code or part of it in any way related to software, which is sold. In this case has to be contacted to Kiril Shterev to negotiate terms and conditions of use.
In any usage of the code, the derivatives has to include the following statement: "This software contains source code provided by Kiril Shterev."
In any case, that is used algorithm SIMPLE-TS has to be cited the main paper presented the algorithm [1] and any other related paper presented further development of the algorithm. The list of papers related to the algorithm are on web site contains source code of the algorithm or on the web page of the Kiril Shterev:
http://www.imbm.bas.bg/index.php/en_US/pressure-based-finite-volume-method
http://www.imbm.bas.bg/index.php/en_US/kiril-stoyanov-shterev
No Support : Kiril Shterev has no obligation to support or to continue providing or updating any of Materials.
No Warranties : This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
Description : Computational Fluid Dynamic using method unsteady SIMPLE-TS, MPI 1.2, C++
References:
1. K. Shterev and S. Stefanov, Pressure based finite volume method for calculation of compressible viscous gas flows, Journal of Computational Physics 229 (2010) pp. 461-480, doi:10.1016/j.jcp.2009.09.042
============================================================================
*/
//If is defined DEBUG_PROGRAM the program will write a lot of arrays to files
//to make debug easier.
//#define DEBUG_PROGRAM
//If is defined COMPILE_MPI the MPI will be compiled,
//otherwise no MPI will be compiled;).
//Just make or not #define COMPILE_MPI as comment.
#define COMPILE_MPI
#ifdef COMPILE_MPI
//Include 4 MPI_Barrier in loop2 and loop3 to stabilize the iterative process,
//when are used more processes (for 150 give error if are not included)
//#define USE_MPI_Barrier_in_iterative_process_to_be_more_stable
#define USE_MPI_Barrier_in_iterative_process_to_be_more_stable_loop2
//#define USE_MPI_Barrier_in_iterative_process_to_be_more_stable_loop3
#include <mpi.h>
#include "DomainDecomposition2D.h"
#endif //End of COMPILE_MPI
//Use fixed number of iteration in loop 3 (internal loop for calculation pressure and temperature)
//For unsteady pressure driven fluid flow in a 2D channel the results show that 2 or 4 iterations
//are good choice for good accuracy for a computational time
#define FIXED_NUMBER_OF_ITERATIONS_LOOP_3
// The solving of p and Temper have to be soved using Jacoby iterative method,
// when is used domain decomposition method (parralel - MPI).
#define Calculate_p_and_Temper_using_Jakoby_method
#ifdef Calculate_p_and_Temper_using_Jakoby_method
//Substitute P with p_old and TEMPER with Temper_old, when is used Jakoby iterative method
#define TEMPER Temper_old
#define P p_old
#else
//Substitute P with p and TEMPER with Temper, when is used Gauss-Seidel iterative method
#define TEMPER Temper
#define P p
#endif
#ifndef COMPILE_MPI
/*
Method of calculation p and Temper
Calculate p and Temper using Jakoby method
If variable Calculate_p_and_Temper_using_Jakoby_method is defined p and Temper
are calculated using method of Jacoby (p[ij] = a1p * p_old[i_1j] + ....).
If for loop for p and Temper is made one iteration, when is used Jacoby method, the
iteration process in some cases is not convergent. To make process converge can be made
some iterations, for example N_I_p_c = 5.
If variable Calculate_p_and_Temper_using_Jakoby_method is NOT defined p and Temper
are calculated using method of Gouse-Zaidel (p[ij] = a1p * p[i_1j] + ....)
*/
//#define Calculate_p_and_Temper_using_Jakoby_method
#endif //End of COMPILE_MPI
//From practical expiriance, when is used Jakoby method of calculation p and Temper
//have to make minimum 2 iteration, to garantee convergentbility and stablelity of this loop.
#ifdef Calculate_p_and_Temper_using_Jakoby_method
const int MinimalNumberOfIterationInLoopFor_p_and_Temper = 2;
#else
const int MinimalNumberOfIterationInLoopFor_p_and_Temper = 1;
#endif
#include "Objects_for_Solving.h"
#include "WriteReadData.h"
#include "stringConv.h"
#include "iostream"
#include <time.h>
#define _USE_MATH_DEFINES
#include <math.h>
#include <stdio.h>
#include <cstring>
//define all function which MUST be inline
//return maximum value
#define maximum(a, b) (((a) > (b)) ? (a) : (b))
//return maximum value
#define minimum(a, b) (((a) < (b)) ? (a) : (b))
//upwind shame - return rho * v
#define upwind(r0, r1, v1) (((v1) > 0) ? ((r0) * (v1)) : ((r1) * (v1)))
//upwind shame - return only rho
#define upwind_rho(r0, r1, v1) (((v1) == 0) ? (0.5 * ((r0) + (r1))) : (((v1) > 0) ? (r0) : (r1)))
//linear interpolation between two given values and steps in mesh
#define Li(fi_1, fi_2, h_1, h_2) (((fi_1) * (h_2) + (fi_2) * (h_1)) / ((h_1) + (h_2)))
//average harmonic between two given values and steps in mesh
#define Hi(fi_1, fi_2, h_1, h_2) (((h_1) + (h_2)) * (fi_1) * (fi_2) / ((h_1) * (fi_2) + (h_2) * (fi_1)))
//Solve index for points i+1, i-1, j+1, j-1
//This MACROS is NO need to be changed to be used in MPI.
#define ij_function(node1) \
ij = (node1);\
i_1 = i - 1 * (0 < i);\
i_2 = i - 1 * (0 < i) - 1 * (1 < i);\
i1 = i + 1 * (i < (Nx - 1));\
i_1j = ij - 1 * (0 < i);\
i1j = ij + 1 * (i < (Nx - 1));\
ij_1 = ij - Nx * (0 < j);\
ij1 = ij + Nx * (j < (Ny - 1));\
i_1j_1 = ij - 1 * (0 < i) - Nx * (0 < j);\
i_1j1 = ij - 1 * (0 < i) + Nx * (j < (Ny - 1));\
i1j_1 = ij + 1 * (i < (Nx - 1)) - Nx * (0 < j);\
i1j1 = ij + 1 * (i < (Nx - 1)) + Nx * (j < (Ny - 1));\
i_2j = i_2 + j * Nx;
//Solve index for points i+1, i-1, j+1, j-1,
//for loop which start from i = 0, j = 0 and are to the i = Nx - 1, j = Ny -1.
//This MACROS is NEEDED to be changed to be used in MPI.
#ifndef COMPILE_MPI
#define ij_function_i_and_j_to_the_boundaries(node1) \
ij = (node1);\
i_1 = i - 1 * (0 < i);\
i_2 = i - 1 * (0 < i) - 1 * (1 < i);\
i1 = i + 1 * (i < (Nx - 1));\
i_1j = ij - 1 * (0 < i);\
i1j = ij + 1 * (i < (Nx - 1));\
ij_1 = ij - Nx * (0 < j);\
ij1 = ij + Nx * (j < (Ny - 1));\
i_1j_1 = ij - 1 * (0 < i) - Nx * (0 < j);\
i_1j1 = ij - 1 * (0 < i) + Nx * (j < (Ny - 1));\
i1j_1 = ij + 1 * (i < (Nx - 1)) - Nx * (0 < j);\
i1j1 = ij + 1 * (i < (Nx - 1)) + Nx * (j < (Ny - 1));\
i_2j = i_2 + j * Nx;
#endif //End of COMPILE_MPI
#ifdef COMPILE_MPI
#define ij_function_i_and_j_to_the_boundaries(node1) \
ij = (node1);\
i_1 = i - 1 * (0 < i);\
i_2 = i - 1 * (0 < i) - 1 * (1 < i);\
i1 = i + 1 * (i < (Nx - 1));\
i_1j = ij - 1 * (0 < i);\
i1j = ij + 1 * (i < (Nx - 1));\
ij_1 = ij - Nx * (0 < j);\
ij1 = ij + Nx * (j < (Ny - 1));\
i_1j_1 = ij - 1 * (0 < i) - Nx * (0 < j);\
i_1j1 = ij - 1 * (0 < i) + Nx * (j < (Ny - 1));\
i1j_1 = ij + 1 * (i < (Nx - 1)) - Nx * (0 < j);\
i1j1 = ij + 1 * (i < (Nx - 1)) + Nx * (j < (Ny - 1));\
i_2j = i_2 + j * Nx;
#endif //End of COMPILE_MPI
//Solve index for points i+1, i-1, i-2, j+1, j-1 for Periodic Boundary Conditions
//This MACROS have to be checked to be used in MPI.
#define ij_function_PeriodicBC(ipBC, jpBC) \
if((ipBC) == 0) {i_1 = Nx - 1;} else {i_1 = (ipBC) - 1;}\
if((ipBC) == 1) {i_2 = Nx - 1;} else if((ipBC) == 0) {i_2 = Nx - 2;} else {i_2 = (ipBC) - 2;}\
if((ipBC) == (Nx - 1)) {i1 = 0;} else {i1 = (ipBC) + 1;}\
ij = (ipBC) + (jpBC) * Nx;\
i_1j = i_1 + (jpBC) * Nx;\
i1j = i1 + (jpBC) * Nx;\
ij_1 = (ipBC) + (jpBC - 1 * (0 < j)) * Nx;\
ij1 = (ipBC) + (jpBC + 1 * (j < (Ny - 1))) * Nx;\
i_1j_1 = i_1 + (jpBC - 1 * (0 < j)) * Nx;\
i_1j1 = i_1 + (jpBC + 1 * (j < (Ny - 1))) * Nx;\
i1j_1 = i1 + (jpBC - 1 * (0 < j)) * Nx;\
i1j1 = i1 + (jpBC + 1 * (j < (Ny - 1))) * Nx;\
i_2j = i_2 + (jpBC) * Nx;
//Solve index for points i+1, i-1, i-2, j+1, j-1 for Periodic Boundary Conditions
//for loop which start from i = 0, j = 0 and are to the i = Nx - 1, j = Ny -1.
//This MACROS have to be checked to be used in MPI.
#define ij_function_PeriodicBC_i_and_j_to_the_boundaries(ipBC, jpBC) \
if((ipBC) == 0) {i_1 = Nx - 1;} else {i_1 = (ipBC) - 1;}\
if((ipBC) == 1) {i_2 = Nx - 1;} else if((ipBC) == 0) {i_2 = Nx - 2;} else {i_2 = (ipBC) - 2;}\
if((ipBC) == (Nx - 1)) {i1 = 0;} else {i1 = (ipBC) + 1;}\
ij = (ipBC) + (jpBC) * Nx;\
i_1j = i_1 + (jpBC) * Nx;\
i1j = i1 + (jpBC) * Nx;\
ij_1 = (ipBC) + (jpBC - 1 * (0 < j)) * Nx;\
ij1 = (ipBC) + (jpBC + 1 * (j < (Ny - 1))) * Nx;\
i_1j_1 = i_1 + (jpBC - 1 * (0 < j)) * Nx;\
i_1j1 = i_1 + (jpBC + 1 * (j < (Ny - 1))) * Nx;\
i1j_1 = i1 + (jpBC - 1 * (0 < j)) * Nx;\
i1j1 = i1 + (jpBC + 1 * (j < (Ny - 1))) * Nx;\
i_2j = i_2 + (jpBC) * Nx;
//Macros to change the values of two pointers to double arrays pointer_double_tmp.
//It is used when have to change the pointers p and p_old instead copy of arrays.
//Have to be defined pointer to double
#define Change_two_pointers_double(pointer_0, pointer_1) \
pointer_double_tmp = pointer_1; \
pointer_1 = pointer_0; \
pointer_0 = pointer_double_tmp;
double * pointer_double_tmp;
using namespace std;
//const double pi = 3.14159265358979; //pi number to 15 singht
Objects_for_Solving * OforS; //Objects for Solving
//#ifndef SIMPLEST_ANL_compr_T_h_var_NoDim_Algorithm_H
//#define SIMPLEST_ANL_compr_T_h_var_NoDim_Algorithm_H
//#pragma once
//Nimber of Dimention of rule vectors for begind and end of area - begin end = 2
static const unsigned short ND = 2;
//template <typename T_a, typename T_b>
//inline double maximum(const T_a& a, const T_b& b)
//{
// return (a > b ? a : b);
//};
//it, it_begin, Nt, ht - this is absolutely time in dimensionless values.
//ht - step by time
double it, it_begin, Nt, ht;
unsigned int it_counter, Nx, Ny, Na, N_I, Iter, Iter_pr1, Iter_pr2;
/*Criterial to continue iteration process.
continue_iter == false -> stop iteration in this time step,
because the convergence criterions are not satisfyed
continue_iter == true -> continue iteration in this time step,
because the convergence criterions are satisfyed
continue_iter_p_c == false -> stop iteration in loop, where are calculated p and Temper,
because the convergence criterions are not satisfyed
continue_iter_p_c == true -> continue iteration in loop, where are calculated p and Temper,
because the convergence criterions are satisfyed
*/
bool continue_iter, continue_iter_p_c;
//countinue_time_step == true -> continue calculating next time step
//countinue_time_step == false -> stop calculating time step, stop program
bool countinue_time_step;
void Write_it_ToFile(void);
void Read_it_FromFile(void);
//Change u_gas_nd in time of solve.
//use_interpolation_for_u_BC_by_time == true; will use interpolation for u_BC
//use_interpolation_for_u_BC_by_time == false; will NOT use interpolation for u_BC
bool use_interpolation_for_u_BC_by_time;
//u_gas_nd will be solved using linear interpolation.
//u_BC_xb_tb - BC of u_gas_nd in time t_BC_b
//u_BC_xb_te - BC of u_gas_nd in time t_BC_e
double u_BC_xb_t_BC_b, u_BC_xb_t_BC_e;
//Time for what is given u_gas_nd. Time is equivalent to nondimensionall time.
//if(it < t_BC_b) then u_gas_nd = u_BC_xb_t_BC_b;
//if(t_BC_e < it) then u_gas_nd = u_BC_xb_t_BC_e;
double t_BC_b, t_BC_e;
//Apply Boundary Conditions which are function of time.
inline void BoundaryConditions_function_of_time(void);
unsigned int i, j, ij, i_1j, i1j, ij_1, ij1, i_1j_1, i_1j1, i1j_1, i1j1, i_2j,
i_1, i1, i_2,
node, counter;
//double rho_gas, mu_gas;
//variables for mesh
bool is_defined_array_for_mesh;
double * x_f, * y_f, * hx, * hy, * x_v, * y_v;
double h_tmp;
//kind of mesh
// _____ __ ______
//kind_of_mesh = Li_mesh - mesh is with linear interpolation step of kind: \/ \/
static const unsigned int Li_mesh = 1;
// _____ __ ______
//kind_of_mesh = Par_mesh - mesh is with parabolik interpolation step: \/ \/
static const unsigned int Par_mesh = 2;
// ______ _______
//kind_of_mesh = Li_flat_mesh - mesh is with linear interpolation step: \__/
static const unsigned int Li_mesh_flat = 3;
//kind_of_mesh = Polynom3_flat_mesh - mesh is with polynom from 3 degree.
//That mean that derivatiors in two given points are zero, ______ _______
//smooth on the two ends of interpoletion: \__/
static const unsigned int Par_mesh_flat = 4;
// ______ _______
//kind_of_mesh = Par_flat_mesh - mesh is with parabolik interpolation step: \__/
static const unsigned int Polynom3_flat_mesh = 5;
unsigned int kind_of_mesh;
//variable for mesh:
double hxmin, hxmax, hxmid, xfmin, xfmax, xbmin, xbmax, xmid,
axf3, axf2, axf1, axf0, axb3, axb2, axb1, axb0;
double hxfront(double& x);
double hxback(double& x);
double hxfront_Polynom3(double& x);
double hxback_Polynom3(double& x);
double hx_all(double& x);
double hx_all_par(double& x);
double hx_all_par_flat(double& x);
double hx_all_linear(double& x);
double hx_all_linear_flat(double& x);
double hx_all_Polynom3_flat_mesh(double& x);
double hymin, hymax, hymid, ytmin, ytmax, ybmin, ybmax, ymid,
ayt3, ayt2, ayt1, ayt0, ayb3, ayb2, ayb1, ayb0;
double hytop(double& y);
double hybottom(double& y);
double hytop_Polynom3(double& y);
double hybottom_Polynom3(double& y);
double hy_all(double& y);
double hy_all_par(double& y);
double hy_all_par_flat(double& y);
double hy_all_linear(double& y);
double hy_all_linear_flat(double& y);
double hy_all_Polynom3_flat_mesh(double& y);
double x_b, x_e, y_b, y_e; //this are area boundaries
//define rule vectors
//The target of rule vectors is to make loop without any check.
//define rule vectors.
void define_rule_vectors(void);
//is_rule_vectors_defined == true - rule vectors are already defined
//is_rule_vectors_defined == false - rule vectors are NOT defined
bool is_rule_vectors_defined;
//vol_inf_fb_bool == true -> if in volume is fluid
//vol_inf_fb_bool == false -> if in volume is body
//The target is to eliminate the check in loop.
//Example 1:
//if(vol_inf_fb[ij] > fluid) fi = 0;
//else fi = something;
//equivalence: fi = vol_inf_fb_bool[ij] * something;
//Example 2:
//if(vol_inf_fb[ij] > fluid || vol_inf_fb[i1j] > fluid) fi = 0;
//else fi = something;
//equivalence: fi = vol_inf_fb_bool[ij] * vol_inf_fb_bool[i1j] * something;
bool * vol_inf_fb_bool;
//define rule vectors and variables for u
//That is posible because the result from checks are the same EVERY LOOP!
//Define vectors and vatriables for area where is no check
//number of all elements of sloving for u of No Check - Na_nch_u
unsigned int Na_nch_u;
//vector contain ij of all elements of sloving for u of No Check - rule_vector_ij_nch_u
unsigned int * rule_vector_ij_nch_u;
//vector contain j of all elements of sloving for u of No Check - rule_vector_j_nch_u
//becouse ij = i + j * Nx => i = ij - j * Nx;
unsigned int * rule_vector_j_nch_u;
//Define vectors and vatriables for area where is Boundary Conditions - here we have some checks
//number of all elements of sloving for u of Boundary Conditions - Na_bc_u
unsigned int Na_bc_u;
//vector contain ij of all elements of sloving for u of Boundary Conditions - rule_vector_ij_bc_u
unsigned int * rule_vector_ij_bc_u;
//vector contain j of all elements of sloving for u of Boundary Conditions - rule_vector_j_bc_u
//becouse ij = i + j * Nx => i = ij - j * Nx;
unsigned int * rule_vector_j_bc_u;
//define rule vectors and variables for v
//That is posible because the result from checks are the same EVERY LOOP!
//Define vectors and vatriables for area where is no check
//number of all elements of sloving for v of No Check - Na_nch_v
unsigned int Na_nch_v;
//vector contain ij of all elements of sloving for v of No Check - rule_vector_ij_nch_v
unsigned int * rule_vector_ij_nch_v;
//vector contain j of all elements of sloving for v of No Check - rule_vector_j_nch_v
//becouse ij = i + j * Nx => i = ij - j * Nx;
unsigned int * rule_vector_j_nch_v;
//Define vectors and vatriables for area where is Boundary Conditions - here we have some checks
//number of all elements of sloving for v of Boundary Conditions - Na_bc_v
unsigned int Na_bc_v;
//vector contain ij of all elements of sloving for v of Boundary Conditions - rule_vector_ij_bc_v
unsigned int * rule_vector_ij_bc_v;
//vector contain j of all elements of sloving for v of Boundary Conditions - rule_vector_j_bc_v
//becouse ij = i + j * Nx => i = ij - j * Nx;
unsigned int * rule_vector_j_bc_v;
//define rule vectors and variables for others - p, Temper, rho in center of volume
//That is posible because the result from checks are the same EVERY LOOP!
//Define vectors and vatriables for area where is no check
//number of all elements of sloving for others of No Check - Na_nch_others
unsigned int Na_nch_others;
//vector contain ij of all elements of sloving for others of No Check - rule_vector_ij_nch_others
unsigned int * rule_vector_ij_nch_others;
//vector contain j of all elements of sloving for others of No Check - rule_vector_j_nch_others
//becouse ij = i + j * Nx => i = ij - j * Nx;
unsigned int * rule_vector_j_nch_others;
//Define vectors and vatriables for area where is Boundary Conditions - here we have some checks
//number of all elements of sloving for others of Boundary Conditions - Na_bc_others
unsigned int Na_bc_others;
//vector contain ij of all elements of sloving for others of Boundary Conditions - rule_vector_ij_bc_others
unsigned int * rule_vector_ij_bc_others;
//vector contain j of all elements of sloving for others of Boundary Conditions - rule_vector_j_bc_others
//becouse ij = i + j * Nx => i = ij - j * Nx;
unsigned int * rule_vector_j_bc_others;
//sqrt_T_ff_V - this is diffusion coefficient (dynamic viscosity) for velocity equations (u and v).
//sqrt_T_ff_V is interpolation of sqrt_T in points x_f[i], y_f[j].
//double * sqrt_T_ff_V;
//This is diffusion coefficient in control surface yf, from integration of diffusion terms of equation for u.
double * Gamma_yf;
//This is diffusion coefficient in control surface xf, from integration of diffusion terms of equation for v.
double * Gamma_xf;
//The fluxes in all area. They are the same for u and v.
double * Fx, * Fy;
//double au33, au323;
//double au44, au424;
//Data are stored in storage arrays.
//Pointer u and u_pr change the array.
//At this way is prevented copy of u_pr = u at the beginning of loop 1.
double * storage_u_0, * storage_u_1,
* u, //* u_old,
* u_pr,
* u_pseudo, u_gas,
* D_ux, * D_uy,
//D_u1, D_u2, D_u3, D_u4,
//Fxu0, Fxu1, Fxu2, Fxu3, Fxu4,
//Fyu3, Fyu4, Fyu23, Fyu24,
F_xu1, F_xu2,
au0, au1, au2, au3, au4, bu3,
Scu, Spu;
//This is temporary coefficient a:
double a_tmp, a_tmp_ij, a_tmp_i1j, a_tmp_ij1;
//Variables to apply velocity slip and temperature jump boundary conditions
//Number of body in control volume, when apply boundary conditions
int Number_of_body;
//Coefficients to apply temperature jump boundary condition
bool is_TemperatureJumpBC_body_i_1_tmp, is_TemperatureJumpBC_body_i1_tmp,
is_TemperatureJumpBC_body_j_1_tmp, is_TemperatureJumpBC_body_j1_tmp;
bool dTdn_on_wall_0_i_1_tmp, dTdn_on_wall_0_i1_tmp,
dTdn_on_wall_0_j_1_tmp, dTdn_on_wall_0_j1_tmp;
double T_body_i_1_tmp, T_body_i1_tmp, T_body_j_1_tmp, T_body_j1_tmp;
double aTBC_i_1, aTBC_i1, aTBC_j_1, aTBC_j1;
////contain information for u. If u > 0 is true, if u <= 0 is false.
//bool * positive_u;
double max_residual_in_u, max_residual_in_u_Iter_0;
//double au_c0, au_c1, au_c2, au_c3, au_c4, bu_c;
//double au_c33, au_c323, au_c44, au_c424;
//double aut0, aut1, aut2, aut3, aut4, but;
//double autce0, autce1, autce2, autce3, autce4, btuce;
//double aFu_c1, aFu_c2;
//double Fxu333, Fxu23323, Fxu444, Fxu24424;
//double av11, av141;
//double av22, av242;
//Data are stored in storage arrays.
//Pointer v and v_pr change the array.
//At this way is prevented copy of v_pr = v at the beginning of loop 1.
double * storage_v_0, * storage_v_1,
* v, //* v_old,
* v_pr,
* v_pseudo, v_gas,
* D_vx, * D_vy,
//D_v1, D_v2, D_v3, D_v4,
//* Fxv, * Fyv,
//Fxv1, Fxv2, Fxv41, Fxv42,
//Fyv0, Fyv1, Fyv2, Fyv3, Fyv4,
F_yv3, F_yv4,
av0, av1, av2, av3, av4, bv3,
Scv, Spv;
//contain information for v. If v > 0 is true, if v <= 0 is false.
//bool * positive_v;
double max_residual_in_v, max_residual_in_v_Iter_0;
//double av_c0, av_c1, av_c2, av_c3, av_c4, bv_c;
//double av_c11, av_c141, av_c22, av_c242;
//double avt0, avt1, avt2, avt3, avt4, bvt;
//double avtce0, avtce1, avtce2, avtce3, avtce4, btvce;
//double aFv_c3, aFv_c4;
//double Fyv111, Fyv41141, Fyv222, Fyv42242;
//Data are stored in storage arrays.
//Pointer p and p_pr change the array.
//At this way is prevented copy of p_pr = p at the beginning of loop 1 and p_old = p in loop 3.
double * storage_p_0, * storage_p_1,
//* p_c, * p_c_old,
* p, * p_pr,
* d_p_c_12, * d_p_c_34,
ap0,
//ap1, ap2, ap3, ap4,
* apx, * apy,
bp0,
MaxError_p_c;
//double fabs_max_p;
#ifdef Calculate_p_and_Temper_using_Jakoby_method
double * storage_p_2, * p_old;
#endif //End of Calculate_p_and_Temper_using_Jakoby_method
//double bp01, bp02, bp03, bp04;
double * bpx, * bpy;
double hyj_upwind_rho_1_0, hxi_upwind_rho_3_0;
//double rho_0_p, rho_1_p, rho_2_p, rho_3_p, rho_4_p;
//double max_residual_in_p_c;
//double max_residual_in_p_c_Itre0;
//double sum_ap_c_and_bpc;
unsigned int Iter_p_c, N_I_p_c;
//density - rho
//Data are stored in storage arrays.
//Pointer rho and rho_pr change the array.
//At this way is prevented copy of rho_pr = rho at the beginning of loop 1.
double * storage_rho_0, * storage_rho_1, * rho, * rho_pr;
////rho_1[ij] = 1 / rho[ij];
////Many coefficient must be zero in body and in this way we do not need to check.
//double * rho_1;
//rho_u[ij] = upwind_rho(rho[i_1j], rho[ij] ,u[ij])
//rho_u[ij] - contain rho in point u[ij] (x_f[i], y_v[j])
double * rho_u;
//rho_v[ij] = upwind_rho(rho[ij_1], rho[ij] ,v[ij])
//rho_v[ij] - contain rho in point v[ij] (x_v[i], y_f[j])
double * rho_v;
//define_OforS
void define_OforS(void);
//Number of matrixes of vol_inf
unsigned int N_vol_inf;
//vol_inf - matrixes witch give us information aboit point
void define_vol_inf(void);
//This variable contain information about fluid and body in point
static const unsigned int fluid = 0; //if vol_inf_fb[ij] is 0 that mean that there is fluid for solving
//Number of bodies circles and polyhedrons which can be counted in vil_inf is greater then fluid
//Ncp_beg_b - Number of circle or polyhedron begin body
static const unsigned int Ncp_beg_b = fluid + 1;
//If vol_inf_fb[ij] == fluid - in ij poin is fluid
//If vol_inf_fb[ij] == 3 - in ij poin is body number 3 - 1 = 2
unsigned int * vol_inf_fb;
//This variable contain information, where in domain is GIVEN p
static const unsigned int no_gp = 0; //if vol_inf_p[ij] is 0 that mean that there is no given pressure
//Number of bodies circles and polyhedrons which can be counted in vil_inf is greater then fluid
//Ncp_beg_p - Number of circle or polyhedron bedin pressure
static const unsigned int Ncp_beg_p = no_gp + 1;
//If vol_inf_p[ij] == flase - in ij poin is NOT GIVEN p
//If vol_inf_p[ij] == true - in ij poin is GIVEN p
unsigned int * vol_inf_p;
void ApplyPressureConditionsToPressure(void);
//This variable contain information, where in domain is GIVEN V(u,v)
static const unsigned int no_gV = 0; //if vol_inf_V[ij] is 0 that mean that there is no given Velocity
//Number of bodies circles and polyhedrons which can be counted in vil_inf is greater then fluid
//Ncp_beg_V - Number of circle or polyhedron bedin pressure
static const unsigned int Ncp_beg_V = no_gV + 1;
//If vol_inf_V[ij] == flase - in ij poin is NOT GIVEN V
//If vol_inf_V[ij] == true - in ij poin is GIVEN V
unsigned int * vol_inf_V;
void ApplyVelocityConditionsToVelocities(void);
//data for body
bool ToImport_data_for_circles_from_file_b,
ToImport_data_for_polyhedrons_from_file_b;
//data for given pressure
bool ToImport_data_for_circles_from_file_p,
ToImport_data_for_polyhedrons_from_file_p;
//data for given velocity
bool ToImport_data_for_circles_from_file_V,
ToImport_data_for_polyhedrons_from_file_V;
//variable which contain where in matrix OforS is information for
//body, pressure and Velocity
static const int gb = 0;
static const int gp = 1;
static const int gV = 2;
void Solve(void);
void define_Nx_Ny(void);
void WriteEnteredDataToFile(void);
void ReadEnteredDataFromFile(void);
bool ToReadSolvedDataFromFile;
//ToReadWriteSolvedDataToBinaryFile == true --> read solved data from binary file
//ToReadWriteSolvedDataToBinaryFile == false --> read solved data from text file
bool ToReadSolvedDataFromBinaryFile;
//ToWriteSolvedDataToBinaryFile == true --> write solved data to binary file
//ToWriteSolvedDataToBinaryFile == false --> write solved data to text file
bool ToWriteSolvedDataToBinaryFile;
//ToWriteSolvedDataToNewFiles == true --> write solved data to new file.
//Example for u: u.bin.100 - this is velocity u, solved in binary mode, it = 100
//To write date to binary ot text mode depend from ToWriteSolvedDataToBinaryFile varieble.
//ToWriteSolvedDataToNewFiles == false --> do not write solved data in new files.
bool ToWriteSolvedDataToNewFiles;
//This is true when we continue calculation readind data which are interpoleted
//from croase mesh.
bool ToContinueFromInterpolation;
void WriteSolvedDataToFile(void);
void ReadSolvedDataFromFile(void);
void WriteTempDataToFile(void);
//ToCheckForFileIsNow_Tmp == true - for first if file exist and the solve is begin will delete exist file
//ToCheckForFileIsNow_Tmp == false - if solve is begin file will not be deleted
bool ToCheckForFileIsNow_Tmp;
void WriteIterationsInLoopsToFile(void);
//ToCheckForFileIsNow_Iterations_in_loops == true - for first if file exist and the solve is begin will delete exist file
//ToCheckForFileIsNow_Iterations_in_loops == false - if solve is begin file will not be deleted
bool ToCheckForFileIsNow_Iterations_in_loops;
//Write data in short format for fast reading.
//The kind of format of Short Data:
//if WriteDataInShortFormat_kind == 0; ==> to not write data in short format
//if WriteDataInShortFormat_kind == 1; ==> write: Re Kn CD CD_p u_gas_nd q_xb q_xe p_BC_xb p_BC_xe
int WriteDataInShortFormat_kind;
//Convergents criterions
double MaxError_Velocities;
////Epsilon1 = 1e-4;
//double Epsilon1;
////Epsilon2 = 1e-6;
//double Epsilon2;
////Epsilon3 = 5 * 1e-5;
// double Epsilon3;
//// fabs_rho_u_h0_max_Iter0 = fabs(gas.rho_gas * u[ij] / h[0]) maximum after Iter = 0
//double fabs_rho_u_h0_max_Iter0;
////Full velocity after Iter = 0. V = sqrt(u^2 + v^2)
//double Vmax_Iter0;
//// fabs_u_u_old_max_Iter0 = fabs(u[ij] - u_old[ij]) maximum after Iter = 0
//double fabs_u_u_old_max_Iter0;
//// fabs_v_v_old_max_Iter0 = fabs(v[ij] - v_old[ij]) maximum after Iter = 0
//double fabs_v_v_old_max_Iter0;
double u_new, v_new;
double u_tmp, v_tmp, p_tmp, Temper_tmp;
//Temporary variables used for applying of Velocity Slip BC
double u_body_tmp, v_body_tmp, w_VelocitySlipBC_tmp, F_VelocitySlip_tmp;
bool is_VelocitySlipBC_body_tmp;
//Temporary variables used for applying of Temperature Jump BC
double T_body_tmp, dTdn_on_wall_0_tmp, F_TemperatureJump_tmp,
F_TemperatureJump_Kn_for_this_BC_tmp;
bool is_TemperatureJumpBC_body_tmp;
// fabs_rho_u_h0_max_Iter0 = fabs(gas.rho_gas * v[ij] / h[1]) maximum after Iter = 0
double fabs_rho_v_h1_max_Iter0;
// fabs_rho_u_h0_max = fabs(gas.rho_gas * u[ij] / h[0]) maximum after check
double fabs_rho_u_h0_max;
// fabs_rho_v_h1_max = fabs(gas.rho_gas * v[ij] / h[1]) maximum after check
double fabs_rho_v_h1_max;
//Full velocity after check. V = sqrt(u^2 + v^2)
double Vmax;
// fabs_u_u_old_max = fabs(u[ij] - u_old[ij]) maximum after check
double fabs_u_u_old_max;
// fabs_v_v_old_max = fabs(v[ij] - v_old[ij]) maximum after check
double fabs_v_v_old_max;
int Nt_save_solved_data, Nt_save_solved_data_counter;
int Nt_save_temp_data, Nt_save_temp_data_counter;
//Nimber of iterations, which are passed before check process for convergence
int N_I_check, N_I_check_counter;
void define_area_for_solving(void);
//Varieble which rules area for solving for velocities, pressure and etc.
//Nx_u[0] begin area for solving for u
//Nx_u[1] end area for solving for u
unsigned int * Nx_u;
//Ny_u[0] begin area for solving for u
//Ny_u[1] end area for solving for u
unsigned int * Ny_u;
//Nx_v[0] begin area for solving for v
//Nx_v[1] end area for solving for v
unsigned int * Nx_v;
//Ny_v[0] begin area for solving for v
//Ny_v[1] end area for solving for v
unsigned int * Ny_v;
//Nx_others[0] begin area for solving for others
//Nx_others[1] end area for solving for others
unsigned int * Nx_others;
//Ny_others[0] begin area for solving for others
//Ny_others[1] end area for solving for others
unsigned int * Ny_others;
//variables for energy equation
//Data are stored in storage arrays.
//Pointer Temper and Temper_pr change the array.
//At this way is prevented copy of Temper_pr = Temper at the beginning of loop 1 and Temper_old = Temper in loop 3.
double * storage_Temper_0, * storage_Temper_1,
* sqrt_T, cT, * Temper,
* Temper_pr,
cT1, cT2,
DT1, DT2, DT3, DT4, FT1, FT2, FT3, FT4,
* DTx, * DTy,
aT0 , aT1, aT2, aT3, aT4, bT0,
* ScT1,
FiT,
dudx_T, dvdy_T, dudy_T, dvdx_T,
//sqrt on boundary of area
sqrt_T_T_1b, sqrt_T_T_2b, sqrt_T_T_3b, sqrt_T_T_4b,
MaxError_T;
//This is Temperature of the fluid on the surface of the body.
//If There is no temperature jump Temperature of the fluid on the surface of the body
//is equal to thetemperature of the body.
//If there is temperature jump the temperature of fluid on the surface of the body will be
//calculated.
//Temper_of_fluid_on_the_wall_u[ij] - this is where is u velocity coordinates (x_f[i], y_v[j])
//Temper_of_fluid_on_the_wall_v[ij] - this is where is v velocity coordinates (x_v[i], y_f[j])
double * Temper_of_fluid_on_the_wall_u;
double * Temper_of_fluid_on_the_wall_v;
#ifdef Calculate_p_and_Temper_using_Jakoby_method
double * storage_Temper_2, * Temper_old;
//For iteration mathod of Jakoby have to be made minimum two iteration on p and Temper loop
//to be calculation process convergent. This is accepted from calculated problems.
const int minimum_iteration_of_p_and_Temper = 2;
#else
const int minimum_iteration_of_p_and_Temper = 1;
#endif //End of Calculate_p_and_Temper_using_Jakoby_method
//double max_residual_in_T;
//double max_residual_in_T_Itre0;
unsigned int Iter_T, N_I_T;
//MaxError_to_stop_program - if we have lower valies of
//MaxError_Velocities, MaxError_p_c, MaxError_T and current Iter == 0
//the program will stop, becouse that will be equivalent
//of stationary flow and final rezult.
//MaxError_to_stop_program == 1e-14 is excellent for 32-bit computers.
double MaxError_to_stop_program;
//max_time_for_solve - maximum time to solve hours. After that time program will be stoped.
//For GRID max time foe solve is 48 hours. Give 45 hours.
double max_time_for_solve;
//time_solve_stop - when program must stop
double time_solve_stop;
//When solve start. Get the number of seconds elapsed since 00:00 hours, Jan 1, 1970 UTC from the system clock.
int time_solve_start;
//solve_is_finished == true --> solve is finished
//solve_is_finished == false --> solve was stopped from some reson. Maybe time_solve_stop is reched.
bool solve_is_finished;
//double rho_0_T, rho_1_T, rho_2_T, rho_3_T, rho_4_T;
//in points
double p_T_0, p_T_1, p_T_2, p_T_3, p_T_4;
//on boundaries
double p_T_1b, p_T_2b, p_T_3b, p_T_4b;
//dimentionless parametars
double Ma, Pr, Fr, gamma1, c_mu, Kn, Re;
double M_local;
//use diferent type of non domentional volues
unsigned int TypeOfNonDimensiolization;
//NonDimentionalGiven_Re_pch_rho_V2_2 - that mean that system depend from
//Reynolds number:
//Re = Re_infinity = (rho_ch * V_ch * L / mu_ch) * (rho_non_dimentional_infinity * V_non_dimentional_infinity / mu_non_dimentional_infinity)
//therefore: cT = (1 / Re) * (rho_non_dimentional_infinity * V_non_dimentional_infinity / mu_non_dimentional_infinity);
//in program: cT = (1.0 / Re) * (u_given_xb * p_BC_xb / T_BC_xb);
//and nondimensiolization for pressure is: pch = rho_ch * V_ch * V_ch / 2;
//therefore: cPch = 0.5;
static const unsigned int NonDimentionalGiven_Re_pch_rho_V2_2 = 1;
//NonDimentionalGiven_Re_pch_rho_R_T - that mean that system depend from
//Reynolds number:
//Re = Re_infinity = (rho_ch * V_ch * L / mu_ch) * (rho_non_dimentional_infinity * V_non_dimentional_infinity / mu_non_dimentional_infinity)
//therefore: cT = (1 / Re) * (rho_non_dimentional_infinity * V_non_dimentional_infinity / mu_non_dimentional_infinity);
//in program: cT = (1.0 / Re) * (u_given_xb * p_BC_xb / T_BC_xb);
//and nondimensiolization for pressure is: pch = rho_ch * R * T_ch;
//therefore: cPch = 1 / (gamma1 * Ma * Ma);
//where:
//Ma = a / V_ch; Mach number, V_ch - character velocity;
static const unsigned int NonDimentionalGiven_Re_pch_rho_R_T = 2;
//NonDimentionalGiven_Kn_pch_rho_R_T - that mean that system depend from
//Knudsen number:
//Kn = l_0 / L;
//where: l_0 - mean-free molecolar path, L - character dimention.
//therefore: cT = c_mu * Kn / Ma;
//V_ch = V_th - thermal velocity; V_th = sqrt(2 * R * T);
//and nondimensiolization for pressure is: pch = rho_ch * R * T_ch;
//therefore: cPch = 1 / (gamma1 * Ma * Ma);
//where:
//Ma = a / V_ch; Mach number, V_ch - character velocity;
static const unsigned int NonDimentionalGiven_Kn_pch_rho_R_T = 3;
//gravity coefficient
//if(SolveGravityField) c_Fr_v = 2.0 / (gamma1 * Fr * Ma * Ma); - force from gravity field
//else c_Fr_v = 0; - force from gravity field
double c_Fr_v;
bool SolveGravityField;
//gas constant no dimention (nd)
double u_gas_nd, v_gas_nd, V_gas_nd;
double cPch, p_gas_nd, rho_gas_nd, T_gas_nd;
void ApplyBodiesConditionsToVariables(void);
inline void BoundaryConditions_Velocities(void);
//Because of the specific of the program velocity slip BC for u are in separate function.
//This BC can be applied right after calculating rule_vector_ij_bc_u, because
//in this calculation are included volumes to the wall.
//This procedure is used only when write data to output files.
inline void BoundaryConditions_Velocities_SlipBC_u(void);
//Vector containing the only the indexes if control volumes for field variables (others),
//which are to the wall and have to be calculated velocity slip boundary condition for
//horizontal component of velocity u.
int * rule_vector_ij_BoundaryConditions_Velocities_SlipBC_u,
* rule_vector_j_BoundaryConditions_Velocities_SlipBC_u,
Na_BoundaryConditions_Velocities_SlipBC_u;
//is_rule_vectors_BoundaryConditions_Velocities_SlipBC_u_defined == true - the rule vectors rule_vector_ij_BoundaryConditions_Velocities_SlipBC_u and rule_vector_j_BoundaryConditions_Velocities_SlipBC_u are defined
//is_rule_vectors_BoundaryConditions_Velocities_SlipBC_u_defined == false - the rule vectors rule_vector_ij_BoundaryConditions_Velocities_SlipBC_u and rule_vector_j_BoundaryConditions_Velocities_SlipBC_u are NOT defined
bool is_rule_vectors_BoundaryConditions_Velocities_SlipBC_u_defined;
//Because of the specific of the program velocity slip BC for v are in separate function.
//This BC can be applied right after calculating rule_vector_ij_bc_v, because
//in this calculation are included volumes to the wall.
//This procedure is used only when write data to output files.
inline void BoundaryConditions_Velocities_SlipBC_v(void);
//Vector containing the only the indexes if control volumes for field variables (others),
//which are to the wall and have to be calculated velocity slip boundary condition for
//vertical component of velocity v.
int * rule_vector_ij_BoundaryConditions_Velocities_SlipBC_v,
* rule_vector_j_BoundaryConditions_Velocities_SlipBC_v,
Na_BoundaryConditions_Velocities_SlipBC_v;
//is_rule_vectors_BoundaryConditions_Velocities_SlipBC_v_defined == true - the rule vectors rule_vector_ij_BoundaryConditions_Velocities_SlipBC_v and rule_vector_j_BoundaryConditions_Velocities_SlipBC_v are defined
//is_rule_vectors_BoundaryConditions_Velocities_SlipBC_v_defined == false - the rule vectors rule_vector_ij_BoundaryConditions_Velocities_SlipBC_v and rule_vector_j_BoundaryConditions_Velocities_SlipBC_v are NOT defined
bool is_rule_vectors_BoundaryConditions_Velocities_SlipBC_v_defined;
//boundary conditions for u dudx
//dudx_0_BC_xb give: du/dx = 0, on x = x_b
//dudx_0_BC_xe give: du/dx = 0, on x = x_e
bool dudx_0_BC_xb, dudx_0_BC_xe;
//durhodx_0_BC_xb - that give: d(u * rho) / dx = 0, on x = x_b;
//durhodx_0_BC_xe - that give: d(u * rho) / dx = 0, on x = x_b;
bool durhodx_0_BC_xb, durhodx_0_BC_xe;
//drhodt_durhodx_dvrhody_0_BC_xe - that give: drho/dt + d(u * rho)/dx + d(v * rho)/dy = 0; - continuity equation over the x_e boundary
//drhodt_durhodx_dvrhody_0_BC_xb - that give: drho/dt + d(u * rho)/dx + d(v * rho)/dy = 0; - continuity equation over the x_b boundary
bool drhodt_durhodx_dvrhody_0_BC_xb, drhodt_durhodx_dvrhody_0_BC_xe;
//u_uMin_Mout_0_BC_xe give: u[i1j] - u[ij] * Min / Mout = 0, on x = x_e
//--> u[i1j] - u[ij] * Min / Mout, on x = x_e
//where:
// Min - mass flux coming into domain
// Mout - mass flux leaveing domain
bool u_uMin_Mout_0_BC_xe;
double Min, Mout;
//This is calculation u[i_1j] from equation for u for node ij.
//The target for this kind of calculation is to keep the momentum on OX
//at the x = x_b.
bool ToCalculate_ui_1j_from_equation_for_uij_at_x_b;
//This is calculation u[i1j] from equation for u for node ij.
//The target for this kind of calculation is to keep the momentum on OX
//at the x = x_e.