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
| Point start_point(scalar f,double xs, double ys){
double dist = sq(L0*pow(2.,0.5));
foreach(reduction(min:dist)){
if (f[]>0. && f[]<1.&&fabs(x-X0)>Delta&&fabs(x-(X0+L0))>Delta&&fabs(y-(Y0+L0))>Delta&&fabs(y-Y0)>Delta){
double dis = (sq(x-xs)+sq(y-ys));
if (dis<dist)
dist = dis;
}
}
double xp=(2*L0)+X0;// These coordinates cannot be located.
double yp=(2*L0)+Y0;
foreach(){
if (sq(x-xs)+sq(y-ys)==dist && f[]>0 && f[]<1.){ // Only true on a single PID.
xp=x;
yp=y;
}
}
return locate(xp,yp);
}
void find_facets(Point point,scalar f,double xy[4]){
coord n;
n = mycs (point, f);
double alpha = plane_alpha (f[], n);
coord segment[2];
if (facets (n, alpha, segment) == 2){
xy[0] = x + segment[0].x*Delta;
xy[1] = y + segment[0].y*Delta;
xy[2] = x + segment[1].x*Delta;
xy[3] = y + segment[1].y*Delta;
}else{
fprintf(stdout,"Warning:\nCould not find facets; expect unexpected behaviour.\n");
}
}
void next_point(Point point,scalar f,double xp,double yp,int dir[2],int * iter,double * xu,double * yu){
int nb=0;
foreach_dimension()
nb+=((f[-1]>0.)*(f[-1]<1.))+((f[1]>0.)*(f[1]<1.));
if (nb==2&&(fabs(dir[0])+fabs(dir[1])==1)){ // Only two neighbouring interfacial cells could be easy...
for(int ii = -1;ii<=1;ii++){
for(int jj = -1;jj<=1;jj++){
if ((abs(ii)+(abs(jj))==1) && (dir[0]!=-ii || dir[1]!=-jj) && f[ii,jj]>0. && f[ii,jj]<1.){// Simple neighbour
dir[0]=ii;
dir[1]=jj;
if (!is_refined(neighbor(ii,jj))){//Special care for refined next cells
*xu=(5*xp+x)/6.+((double)ii*Delta/3.);
*yu=(5*yp+y)/6.+((double)jj*Delta/3.);
return;
}else if(fine(f,(xp>x)+ii,(yp>y)+jj)>0 && fine(f,(xp>x)+ii,(yp>y)+jj)<1){//still OK
*xu=(5*xp+x)/6.+((double)ii*Delta/3.);
*yu=(5*yp+y)/6.+((double)jj*Delta/3.);
return;
}else if(ii==0){//y-dir
if(fine(f,(xp<x),((yp>y))+jj)>0 && fine(f,(xp<x),(yp>y)+jj)<1){ //poorly reconstructed, do not follow facet
*xu=x+(0.5*Delta*((xp<x)-0.5));
*yu=y+(0.75*Delta*jj) ;
return;
}
}else if(fine(f,(xp>x)+ii,(yp<y))>0 && fine(f,(xp>x)+ii,(yp<y))<1){ //x-dir, poorly reconstructed, do not follow facet
*yu=y+(0.5*Delta*((yp<y)-0.5));
*xu=x+(0.75*Delta*ii);
return;
}
}
}
}
}
if (fabs(xp-X0)<Delta||fabs(xp-(X0+L0))<=Delta||fabs(yp-Y0)<Delta||fabs(yp-(Y0+L0))<Delta){//domain boundary cell: stop tracing
*iter += (1<<25);// add large number
*xu=X0+L0/2.;
*yu=Y0+L0/2.;
return;
}
// If previous efforts have not worked, tryo to follow the interface
if (xp==x+Delta/2&&((f[1]<1.&&f[1]>0.)||(!is_leaf(neighbor(1,0))&&coarse(f,1,0)<1&&coarse(f,1,0)>0))){
dir[0]=1.;
dir[1]=0;
*xu=xp+Delta/10.;
*yu=yp;
return;
}
else if (xp==x-Delta/2.&&((f[-1]<1.&&f[-1]>0.)||(!is_leaf(neighbor(-1,0))&&coarse(f,-1,0)<1&&coarse(f,-1,0)>0))){
dir[0]=-1;
dir[1]=0;
*xu=xp-Delta/10.;
*yu=yp;
return;
}
else if (yp==y+Delta/2.&&((f[0,1]<1.&&f[0,1]>0.)||(!is_leaf(neighbor(0,1))&&coarse(f,0,1)<1&&coarse(f,0,1)>0))){
dir[0]=0;
dir[1]=1;
*xu=xp;
*yu=yp+Delta/10.;
return;
}
else if (yp==y-Delta/2.&&((f[0,-1]<1.&&f[0,-1]>0.)||(!is_leaf(neighbor(0,-1))&&coarse(f,0,-1)<1&&coarse(f,0,-1)>0))) {
dir[0]=0;
dir[1]=-1;
*xu=xp;
*yu=yp-Delta/10;
return;
}
//As a last resort, try a diagnonal neighbour, blindly. If this does not go well, we are in in trouble
double fac[4];
find_facets(point,f,fac);
double meanx = (fac[0]+fac[2])/2.;
double meany = (fac[1]+fac[3])/2.;
double dx,dy;
if (yp>meany){
dy= 0.3;
}else
dy = -0.3;
if (xp>meanx){
dx= 0.3;
}else
dx = -0.3;
dir[0]=0; dir[1]=0;
*xu=xp+Delta*dx;
*yu=yp+Delta*dy;
return;
}
void update_new_facet(Point point, double xyn[2],double xyf[4],int dir[2]){
if ((fabs(dir[0])+fabs(dir[1])==1)){ // Try to make use of the direction.
for (int g=-1;g<=1;g+=2){
if (dir[0]==g){ //x-dir
double dg = (double)g;
if(dg*xyf[0]==min(dg*xyf[0],dg*xyf[2])){
xyn[0]=xyf[2];xyn[1]=xyf[3];
return;
}else if(dg*xyf[2]==min(dg*xyf[0],dg*xyf[2])){
xyn[0]=xyf[0];xyn[1]=xyf[1];
return;
}
}else if(dir[1]==g){ //y-dir
double dg = (double)g;
if(dg*xyf[1]==min(dg*xyf[1],dg*xyf[3])){
xyn[0]=xyf[2];xyn[1]=xyf[3];
return;
}else if(dg*xyf[3]==min(dg*xyf[1],dg*xyf[3])){
xyn[0]=xyf[0];xyn[1]=xyf[1];
return;
}
}
}
}
// Else choose the one farest away from the previous. This is not very robust and is the main cause of 'bounce backs'.
if ((sq(xyf[0]-xyn[0])+sq(xyf[1]-xyn[1])) <= (sq(xyf[2]-xyn[0])+sq(xyf[3]-xyn[1]))){
xyn[0]=xyf[2]; xyn[1]=xyf[3];
}else{
xyn[0]=xyf[0]; xyn[1]=xyf[1];
}
}
void loop_interfacial_cells(FILE * fp, scalar f,scalar tag,double xynn[2]){
if (dimension!=2){
fprintf(stdout,"2D only!\n");
if (pid()==0)
fprintf(fp,"# 2D only!\n");
return;
}
boundary({f});//neighbors on trees etc.
double xyn[2];
xyn[0]=xynn[0]; xyn[1]=xynn[1];
int ifc= 0;
foreach(reduction(+:ifc)){
if ( f[]>0. && f[]<1.)
ifc++;
}
double xu,yu;
Point point;
point = start_point(f,xyn[0],xyn[1]);
// int i = 1;
//static FILE * fpp = popen ("gfsview2D tag.gfv -s","w"); //For trouble shooting
int igg=0;
double l =0,cur;
double xyf[4];
int dir[2]={0,0};
while (igg <= ifc){
xu=X0-L0;
yu=Y0-L0;
if (point.level>0){ // Local pid()
find_facets(point,f,xyf);
l+=pow(sq(xyf[0]-xyf[2])+sq(xyf[1]-xyf[3]),0.5);
cur = tag[];
//tag[]+=100.; // For trouble shooting
update_new_facet(point,xyn,xyf,dir); // For tracing the interface in the correct direction
next_point(point,f,xyn[0],xyn[1],dir,&igg,&xu,&yu); // Aims to Update xu and yu to be in next leaf cell along the interface
}else{ // Set small values for the tracing variables on the other threads
cur=-HUGE;
xyf[0]=xu; xyf[1]=yu; xyf[2]=xu; xyf[3]=yu;
xyn[0]=xu;xyn[1]=yu;
dir[0]=-2;dir[1]=-2;
}
@ if _MPI
// Set send buffers, so reduced variables can be recieved in the original ones.
double oldx=xu;double oldy=yu;double oldl = l; double xyfold[4];
xyfold[0]=xyf[0];xyfold[1]=xyf[1];xyfold[2]=xyf[2];xyfold[3]=xyf[3];
double oldcur=cur; double oldxyn[2];
oldxyn[0]=xyn[0];oldxyn[1]=xyn[1];
int olddir[2];
olddir[1]=dir[1];olddir[0]=dir[0];
int oldi=igg;
// Update all workers, just in case one has to take over, this should be improved for better performance.
MPI_Allreduce(&oldx, &xu, 1, MPI_DOUBLE, MPI_MAX,MPI_COMM_WORLD);
MPI_Allreduce(&oldy, &yu, 1, MPI_DOUBLE, MPI_MAX,MPI_COMM_WORLD);
MPI_Allreduce(&xyfold, &xyf, 4, MPI_DOUBLE, MPI_MAX,MPI_COMM_WORLD);
MPI_Allreduce(&oldxyn, &xyn, 2, MPI_DOUBLE, MPI_MAX,MPI_COMM_WORLD);
MPI_Allreduce(&oldl, &l, 1, MPI_DOUBLE, MPI_MAX,MPI_COMM_WORLD);
MPI_Allreduce(&oldcur, &cur, 1, MPI_DOUBLE, MPI_MAX,MPI_COMM_WORLD);
MPI_Allreduce(&olddir, &dir, 2, MPI_INT, MPI_MAX,MPI_COMM_WORLD);
MPI_Allreduce(&oldi, &igg, 1, MPI_INT, MPI_MAX,MPI_COMM_WORLD);
@ endif
if (pid()==0&&igg<ifc){ // Our favorite worker writes in the file
fprintf(fp,"%d\t%g\t%g\t%g\t%g\n",igg,(xyf[0]+xyf[2])/2,(xyf[1]+xyf[3])/2,l,cur);
fflush(fp);
}
// All threads try to find the new location.
point = locate(xu,yu);
igg++;
// i = igg;
//output_gfs(fpp); //For trouble shooting purposes:
}
}
|
