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
|
/**
* main.c
* ------
*
* executable file
*
* @author ilhan oezgen, wahyd, ilhan.oezgen@wahyd.tu-berlin.de
* @date 6 apr 17
* 10 apr 17 : switch to 1d vectors
* 2 may 17 : switch to structs
**/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "global.h"
#include "sysmat.h"
#include "linalg.h"
#include "solver.h"
int main (void)
{
// build the domain
int nx = 50;
int ny = 50;
int dim = nx * ny;
double xL = -1.0; // left boundary
double xR = 1.0; // right boundary
double yD = -1.0; // lower boundary
double yU = 1.0; // upper boundary
// boundary type
int north = 0;
int south = 1;
int east = 2;
int west = 2;
double val_north = 40.0;
double val_south = 1.8e-5;
double val_east = 0.0;
double val_west = 0.0;
double dx = (xR - xL) / nx;
double dy = (yU - yD) / ny;
double kf1 = 1.0e-4;
double kf2 = 1.0e-3;
int ids[5] = { - nx, -1, 0, 1, nx }; // system matrix
int pre_ids[1] = { 0 }; // preconditioner
// allocate memory
double* kf = (double*) calloc(dim, sizeof(double));
double* a = (double*) calloc(5 * dim, sizeof(double));
double* c = (double*) calloc(dim, sizeof(double));
double* r = (double*) calloc(dim, sizeof(double));
double* x = (double*) calloc(dim, sizeof(double));
double* vx = (double*) calloc(dim, sizeof(double));
double* vy = (double*) calloc(dim, sizeof(double));
// initialize permeability matrix
for (int i = 0; i < ny; i ++) {
for (int j = 0; j < nx; j ++) {
double xij = xL + j * dx + 0.5 * dx;
double yij = yD + i * dy + 0.5 * dy;
if (xij < 0.2 && xij > -0.2 && yij < 0.2 && yij > -0.2) {
kf[i * nx + j] = kf2;
} else {
kf[i * nx + j] = kf1;
}
}
}
// initialize system matrix
//
struct diag m;
m.values = a;
m.ids = ids;
m.len_c = dim;
m.len_r = dim;
m.len_d = dim;
m.len_id = 5;
system_matrix(&m, nx, ny, kf);
// initialize boundary conditions
//
boundary(dim, nx, ny, north, south, east, west, val_north, val_south,
val_east, val_west, &m, r, x);
// solve the system A*x = b
//
for (int i = 0; i < dim; i ++) {
x[i] = 1.0;
}
// define preconditioner
for (int i = 0; i < dim; i ++) {
c[i] = 1.0; // / get(&m, i, i);
}
struct diag pre = { c, pre_ids, dim, dim, dim, 1 };
pcg(&m, r, x, &pre, 1.0e-12);
darcy(x, kf, nx, ny, dx, dy, vx, vy);
// i/o
//
FILE *fh = fopen("result/head.dat", "w");
if (fh == NULL) {
printf("Error opening file.\n");
}
fprintf(fh, "\n>>> (app) :: solution:\n");
for (int i = 0; i < ny; i ++) {
//printf(">>>");
for (int j = 0; j < nx; j ++) {
fprintf(fh, " %6g ", x[i * nx + j]);
}
fprintf(fh, "\n");
}
fclose(fh);
FILE *fvx = fopen("result/velx.dat", "w");
if (fvx == NULL) {
printf("Error opening file.\n");
}
fprintf(fvx, "\n>>> (app) :: solution:\n");
for (int i = 0; i < ny; i ++) {
for (int j = 0; j < nx; j ++) {
fprintf(fvx, " %6g ", vx[i * nx + j]);
}
fprintf(fvx, "\n");
}
fclose(fvx);
FILE *fvy = fopen("result/vely.dat", "w");
if (fvy == NULL) {
printf("Error opening file.\n");
}
fprintf(fvy, "\n>>> (app) :: solution:\n");
for (int i = 0; i < ny; i ++) {
for (int j = 0; j < nx; j ++) {
fprintf(fvy, " %6g ", vy[i * nx + j]);
}
fprintf(fvy, "\n");
}
fclose(fvy);
// clean up
free(kf);
free(a);
free(c);
free(r);
free(x);
free(vx);
free(vy);
return 0;
}
|
