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
|
#include "Blues.h"
event movie(t += 0.1){
for (scalar s in {f}){
view(camera = "top");
box();
draw_vof("f");
save("top.mp4");
}
}
int dissipation_rate (double* rates){
double rateWater = 0.0;
double rateAir = 0.0;
double rateWater2 = 0.0;
foreach (reduction (+:rateWater), reduction (+:rateAir), reduction (+:rateWater2)) {
if (cs[] == 1) {
double dcdx = (c[1] - c[-1] )/(2.*Delta);
double dcdy = (c[0,1] - c[0,-1] )/(2.*Delta);
double dcdz = (c[0,0,1] - c[0,0,-1])/(2.*Delta);
double dudx = (u.x[1] - u.x[-1] )/(2.*Delta);
double dudy = (u.x[0,1] - u.x[0,-1] )/(2.*Delta);
double dudz = (u.x[0,0,1] - u.x[0,0,-1])/(2.*Delta);
double dvdx = (u.y[1] - u.y[-1] )/(2.*Delta);
double dvdy = (u.y[0,1] - u.y[0,-1] )/(2.*Delta);
double dvdz = (u.y[0,0,1] - u.y[0,0,-1])/(2.*Delta);
double dwdx = (u.z[1] - u.z[-1] )/(2.*Delta);
double dwdy = (u.z[0,1] - u.z[0,-1] )/(2.*Delta);
double dwdz = (u.z[0,0,1] - u.z[0,0,-1])/(2.*Delta);
double SDeformxx = dudx;
double SDeformxy = 0.5*(dudy + dvdx);
double SDeformxz = 0.5*(dudz + dwdx);
double SDeformyx = SDeformxy;
double SDeformyy = dvdy;
double SDeformyz = 0.5*(dvdz + dwdy);
double SDeformzx = SDeformxz;
double SDeformzy = SDeformyz;
double SDeformzz = dwdz;
double sqterm = 2.*dv()*(sq(SDeformxx) + sq(SDeformxy) + sq(SDeformxz) +
sq(SDeformyx) + sq(SDeformyy) + sq(SDeformyz) +
sq(SDeformzx) + sq(SDeformzy) + sq(SDeformzz)) ;
rateWater += mu1/rho[]*f[]*sqterm; //water
rateAir += mu2/rho[]*(1. - f[])*sqterm; //air
rateWater2 += c.D1*f[]*dv()*(sq(dcdx) + sq(dcdy) + sq(dcdz)); //water
}
}
rates[0] = rateWater;
rates[1] = rateAir;
rates[2] = rateWater2;
return 0;
}
#include "../output_fields/available_potential.h"
event time_series(t += 0.1){
scalar zref[];
double zmix = reference_height(zref, f, c, 0.0, 1.0, true, _H0/2.);
double ekin = 0., epot[3] = {0.}, bpot = 0.;
foreach (reduction(+ : ekin), reduction(+ : epot[:3]), reduction(+: bpot)){
epot[0] += dv() * rhov[] * x;
epot[1] += dv() * rhov[] * y;
epot[2] += dv() * rhov[] * z;
bpot += dv() * f[] * rhov[] * zref[];
foreach_dimension(){
ekin += dv() * 0.5 * rhov[] * sq(u.x[]);
}
}
double rates[3];
dissipation_rate(rates);
double disvWater = rates[0];
double disvAir = rates[1];
double discWater = rates[2];
if (pid() == 0){
FILE *fp = fopen("globals.asc", "a");
fprintf(fp, "%f %.9g %.9g %.9g %.9g %.9g %.9g %.9g %.9g %.9g %.9g %.9g \n", t, force.Gn, ekin, epot[0], epot[1], epot[2], bpot, zmix, disvWater, disvAir, discWater, force.ramp);
fclose(fp);
}
}
event save_profiles(t += 0.2){
scalar mfrac[];
foreach (){
mfrac[] = rho1*f[]/(f[]*(rho1 - rho2) + rho2);
}
scalar *list = {u, p, mfrac, c, zref};
int n = (1 << MAXLEVEL);
int m1 = (1 << MAXLEVEL) / (L0 / _H0);
int m2 = (1 << MAXLEVEL) / (L0 / _D0);
profile_foreach_region2(list, rhov, filename="profiles_favre.asc", xmin = -L0/2., xmax = L0/2., ymin = -_D0 / 2., ymax = _D0 / 2., hmin = -_H0 / 2. + _mindel, hmax = _H0 / 2. - _mindel, n = n, m1 = m1, m2 = m2);
}
#include "curvature.h"
event save_facets(t += 0.2){
scalar kappa[];
curvature(f, kappa);
char fname[99];
sprintf(fname, "interface_t%.5f", t);
output_facets_vtu(f, kappa, fname);
}
event save_slices(t += 0.7854){
reference_height(yref, f, c, 0.0, 1.0, true, _H0/2.);
foreach()
yref[] *= cs[]*f[];
char filename[99];
sprintf(filename, "slice_x_t%09.5f.vtkhdf", t);
output_vtkhdf_slice({f, p, c, rhov, yref}, {u}, filename, (coord){1, 0, 0}, 0);
sprintf(filename, "slice_y_t%09.5f.vtkhdf", t);
output_vtkhdf_slice({f, p, c, rhov, yref}, {u}, filename, (coord){0, 1, 0}, 0);
sprintf(filename, "slice_z_t%09.5f.vtkhdf", t);
output_vtkhdf_slice({f, p, c, rhov, yref}, {u}, filename, (coord){0, 0, 1}, 0);
}
event save_snapshots(t += 6.2832){
char filename[99];
sprintf(filename, "snapshot_t%09.5f.vtkhdf", t);
output_vtkhdf({f, p, c, rhov, yref}, {u}, filename);
}
|
