-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathlinear_p_equation_historic.cpp
129 lines (77 loc) · 3.15 KB
/
linear_p_equation_historic.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
//#include"pParticles.h"
#include"linear.h"
#include"fields_enum.h"
#include"simu.h"
// Solve for pressure
// The famous pressure Poisson equation
void linear::p_equation(const FT dt ) {
cout << "Solving pressure equation " << endl;
// fill_Delta_DD(); // This may be important -- or not
FT ddt = dt;
if( dt < 1e-10 ) ddt = 1; // for debugging, mainly
// A
// Approximate Laplacian ~ Delta
// VectorXd divUstar = DD_scalar_vfield( vfield_list::Ustar );
// VectorXd p = Delta_solver.solve( divUstar );
// times (-0.5), because the Laplacian is approximated by -2 Delta / V
// vctr_to_field( -0.5 * p / ddt , sfield_list::p ) ;
// B
// Laplacian as div of grad :
// VectorXd divUstar = DD_scalar_vfield( vfield_list::Ustar );
// VectorXd p = LL_solver.solve( divUstar );
// vctr_to_field( p / ddt , sfield_list::p ) ;
// C
// As B, but Dvol source. This is an _iterative_ procedure,
// yielding a pressure change
// volumes( T );
VectorXd vol = field_to_vctr( sfield_list::vol ) ;
//FT target_vol_val = simu.meanV() ;
// FT target_vol_val = vol.array().sum() / FT( vol.size() );
//VectorXd Dvol = vol.array() - target_vol_val ;
VectorXd vol0 = field_to_vctr( sfield_list::vol0 ) ;
VectorXd Dvol = vol.array() - vol0.array() ;
int N = vol.size();
FT Dvol_sigma = Dvol.array().square().sum() / N ; // / FT( vol.size() );
FT Dvol_mean = vol.array().sum() / N ; // / FT( vol.size() );
cout << "Pressure "
<< " rel Dvol std dev: " << sqrt( Dvol_sigma ) / Dvol_mean
<< endl;
VectorXd Dp = LL_solver.solve( Dvol );
VectorXd p0 = field_to_vctr( sfield_list::p ) ;
//vctr_to_field( p0 + Dp / ( ddt * ddt) , sfield_list::p ) ;
vctr_to_field( Dp / ( ddt * ddt) , sfield_list::p ) ;
// vctr_to_field( vol , sfield_list::vol0 );
return;
}
void linear::u_add_press_grad( const FT dt ) {
VectorXd gradPx,gradPy;
DD_times_sfield( sfield_list::p , gradPx, gradPy);
VectorXd vol = field_to_vctr( sfield_list::vol );
// perhaps mean vol would be just fine
VectorXd Ustar_x, Ustar_y;
vfield_to_vctrs( vfield_list::Ustar , Ustar_x, Ustar_y );
VectorXd U_x, U_y;
FT ddt = dt;
// if( dt < 1e-10 ) ddt = 1; // for debugging, mainly
U_x = Ustar_x.array() - ddt * gradPx.array() / vol.array() ;
U_y = Ustar_y.array() - ddt * gradPy.array() / vol.array() ;
vctrs_to_vfield( U_x, U_y , vfield_list::U );
}
void linear::u_add_grads( const FT dt ) {
VectorXd gradPx,gradPy;
VectorXd gradsx,gradsy;
DD_times_sfield( sfield_list::p , gradPx, gradPy);
MM_times_sfield( sfield_list::s , gradsx, gradsy);
VectorXd vol = field_to_vctr( sfield_list::vol );
// perhaps mean vol would be just fine
VectorXd Ustar_x, Ustar_y;
vfield_to_vctrs( vfield_list::Ustar , Ustar_x, Ustar_y );
VectorXd U_x, U_y;
FT ddt = dt;
// if( dt < 1e-10 ) ddt = 1; // for debugging, mainly
VectorXd grad_x = gradPx + gradsx ;
VectorXd grad_y = gradPy + gradsy ;
U_x = Ustar_x.array() - ddt * grad_x.array() / vol.array() ;
U_y = Ustar_y.array() - ddt * grad_y.array() / vol.array() ;
vctrs_to_vfield( U_x, U_y , vfield_list::U );
}