21 #include "grass/N_pde.h"
24 static int make_les_entry_2d(
int i,
int j,
int offset_i,
int offset_j,
28 double entry,
int cell_type);
30 static int make_les_entry_3d(
int i,
int j,
int k,
int offset_i,
int offset_j,
31 int offset_k,
int count,
int pos,
N_les * les,
34 double entry,
int cell_type);
137 G_debug(5,
"N_create_5star: w %g e %g n %g s %g c %g v %g\n", star->
W,
138 star->
E, star->
N, star->
S, star->
C, star->
V);
160 double S,
double T,
double B,
double V)
175 G_debug(5,
"N_create_7star: w %g e %g n %g s %g t %g b %g c %g v %g\n",
176 star->
W, star->
E, star->
N, star->
S, star->
T, star->
B, star->
C,
201 double S,
double NW,
double SW,
double NE,
220 "N_create_9star: w %g e %g n %g s %g nw %g sw %g ne %g se %g c %g v %g\n",
221 star->
W, star->
E, star->
N, star->
S, star->
NW, star->
SW, star->
NE,
222 star->
SE, star->
C, star->
V);
264 double NW,
double SW,
double NE,
double SE,
265 double T,
double W_T,
double E_T,
double N_T,
266 double S_T,
double NW_T,
double SW_T,
267 double NE_T,
double SE_T,
double B,
double W_B,
268 double E_B,
double N_B,
double S_B,
double NW_B,
269 double SW_B,
double NE_B,
double SE_B,
double V)
309 "N_create_27star: w %g e %g n %g s %g nw %g sw %g ne %g se %g c %g v %g\n",
310 star->
W, star->
E, star->
N, star->
S, star->
NW, star->
SW, star->
NE,
311 star->
SE, star->
C, star->
V);
314 "N_create_27star: w_t %g e_t %g n_t %g s_t %g nw_t %g sw_t %g ne_t %g se_t %g t %g \n",
319 "N_create_27star: w_b %g e_b %g n_b %g s_b %g nw_b %g sw_b %g ne_b %g se_B %g b %g\n",
427 star->
E = 1 / geom->
dx;
428 star->
W = 1 / geom->
dx;
429 star->
N = 1 / geom->
dy;
430 star->
S = 1 / geom->
dy;
431 star->
T = 1 / geom->
dz;
432 star->
B = 1 / geom->
dz;
433 star->
C = -1 * (2 / geom->
dx + 2 / geom->
dy + 2 / geom->
dz);
437 "N_callback_template_3d: w %g e %g n %g s %g t %g b %g c %g v %g\n",
438 star->
W, star->
E, star->
N, star->
S, star->
T, star->
B, star->
C,
466 star->
E = 1 / geom->
dx;
467 star->
NE = 1 / sqrt(geom->
dx * geom->
dx + geom->
dy * geom->
dy);
468 star->
SE = 1 / sqrt(geom->
dx * geom->
dx + geom->
dy * geom->
dy);
469 star->
W = 1 / geom->
dx;
470 star->
NW = 1 / sqrt(geom->
dx * geom->
dx + geom->
dy * geom->
dy);
471 star->
SW = 1 / sqrt(geom->
dx * geom->
dx + geom->
dy * geom->
dy);
472 star->
N = 1 / geom->
dy;
473 star->
S = 1 / geom->
dy;
475 -1 * (star->
E + star->
NE + star->
SE + star->
W + star->
NW + star->
SW +
564 int cell_type_count = 0;
570 "N_assemble_les_2d: starting to assemble the linear equation system");
581 for (j = 0; j < geom->
rows; j++) {
582 for (i = 0; i < geom->
cols; i++) {
592 for (j = 0; j < geom->
rows; j++) {
593 for (i = 0; i < geom->
cols; i++) {
601 G_debug(2,
"N_assemble_les_2d: number of used cells %i\n",
604 if (cell_type_count == 0)
606 (
"Not enough cells [%i] to create the linear equation system. Check the cell status. Only active cells (value = 1) are used to create the equation system.",
611 index_ij = (
int **)G_calloc(cell_type_count,
sizeof(
int *));
612 for (i = 0; i < cell_type_count; i++)
613 index_ij[i] = (
int *)G_calloc(2,
sizeof(
int));
621 for (j = 0; j < geom->
rows; j++) {
622 for (i = 0; i < geom->
cols; i++) {
628 index_ij[count][0] = i;
629 index_ij[count][1] = j;
632 "N_assemble_les_2d: non-inactive cells count %i at pos x[%i] y[%i]\n",
639 index_ij[count][0] = i;
640 index_ij[count][1] = j;
643 "N_assemble_les_2d: active cells count %i at pos x[%i] y[%i]\n",
649 G_debug(2,
"N_assemble_les_2d: starting the parallel assemble loop");
652 #pragma omp parallel for private(i, j, pos, count) schedule(static)
653 for (count = 0; count < cell_type_count; count++) {
654 i = index_ij[count][0];
655 j = index_ij[count][1];
671 les->
b[count] = items->
V;
682 les->
A[count][count] = items->
C;
686 pos = make_les_entry_2d(i, j, -1, 0, count,
pos, les, spvect,
687 cell_count, status, start_val, items->
W,
691 if (i < geom->
cols - 1) {
692 pos = make_les_entry_2d(i, j, 1, 0, count,
pos, les, spvect,
693 cell_count, status, start_val, items->
E,
699 make_les_entry_2d(i, j, 0, -1, count,
pos, les, spvect,
700 cell_count, status, start_val, items->
N,
704 if (j < geom->rows - 1) {
705 pos = make_les_entry_2d(i, j, 0, 1, count,
pos, les, spvect,
706 cell_count, status, start_val, items->
S,
712 if (i > 0 && j > 0) {
713 pos = make_les_entry_2d(i, j, -1, -1, count,
pos, les, spvect,
714 cell_count, status, start_val,
715 items->
NW, cell_type);
718 if (i < geom->
cols - 1 && j > 0) {
719 pos = make_les_entry_2d(i, j, 1, -1, count,
pos, les, spvect,
720 cell_count, status, start_val,
721 items->
NE, cell_type);
724 if (i > 0 && j < geom->rows - 1) {
725 pos = make_les_entry_2d(i, j, -1, 1, count,
pos, les, spvect,
726 cell_count, status, start_val,
727 items->
SW, cell_type);
730 if (i < geom->
cols - 1 && j < geom->rows - 1) {
731 pos = make_les_entry_2d(i, j, 1, 1, count,
pos, les, spvect,
732 cell_count, status, start_val,
733 items->
SE, cell_type);
750 for (i = 0; i < cell_type_count; i++)
789 int i, j, x,
y, stat;
794 "N_les_integrate_dirichlet_2d: integrating the dirichlet boundary condition");
800 dvect1 = (
double *)G_calloc(les->
cols,
sizeof(
double));
801 dvect2 = (
double *)G_calloc(les->
cols,
sizeof(
double));
805 for (y = 0; y < rows; y++) {
806 for (x = 0; x <
cols; x++) {
819 #pragma omp parallel default(shared)
826 #pragma omp for schedule (static) private(i)
827 for (i = 0; i < les->
cols; i++)
828 les->
b[i] = les->
b[i] - dvect2[i];
834 for (y = 0; y < rows; y++) {
835 for (x = 0; x <
cols; x++) {
840 for (i = 0; i < les->
Asp[count]->
cols; i++)
843 for (i = 0; i < les->
rows; i++) {
844 for (j = 0; j < les->
Asp[i]->
cols; j++) {
856 for (i = 0; i < les->
cols; i++)
857 les->
A[count][i] = 0.0;
859 for (i = 0; i < les->
rows; i++)
860 les->
A[i][count] = 0.0;
863 les->
A[count][count] = 1.0;
878 int make_les_entry_2d(
int i,
int j,
int offset_i,
int offset_j,
int count,
881 N_array_2d * start_val,
double entry,
int cell_type)
899 if ((count + K) >= 0 && (count + K) < les->
cols) {
901 " make_les_entry_2d: (N_CELL_ACTIVE) create matrix entry at row[%i] col[%i] value %g\n",
902 count, count + K, entry);
909 les->
A[count][count + K] = entry;
919 if ((count + K) >= 0 && (count + K) < les->
cols) {
921 " make_les_entry_2d: (N_CELL_DIRICHLET) create matrix entry at row[%i] col[%i] value %g\n",
922 count, count + K, entry);
929 les->
A[count][count + K] = entry;
1015 int i, j, k, count = 0, pos = 0;
1016 int cell_type_count = 0;
1022 "N_assemble_les_3d: starting to assemble the linear equation system");
1033 for (k = 0; k < geom->
depths; k++) {
1034 for (j = 0; j < geom->
rows; j++) {
1035 for (i = 0; i < geom->
cols; i++) {
1048 for (k = 0; k < geom->
depths; k++) {
1049 for (j = 0; j < geom->
rows; j++) {
1050 for (i = 0; i < geom->
cols; i++) {
1062 "N_assemble_les_3d: number of used cells %i\n", cell_type_count);
1064 if (cell_type_count == 0.0)
1066 (
"Not enough active cells [%i] to create the linear equation system. Check the cell status. Only active cells (value = 1) are used to create the equation system.",
1073 index_ij = (
int **)G_calloc(cell_type_count,
sizeof(
int *));
1074 for (i = 0; i < cell_type_count; i++)
1075 index_ij[i] = (
int *)G_calloc(3,
sizeof(
int));
1080 for (k = 0; k < geom->
depths; k++) {
1081 for (j = 0; j < geom->
rows; j++) {
1082 for (i = 0; i < geom->
cols; i++) {
1089 index_ij[count][0] = i;
1090 index_ij[count][1] = j;
1091 index_ij[count][2] = k;
1094 "N_assemble_les_3d: non-inactive cells count %i at pos x[%i] y[%i] z[%i]\n",
1101 index_ij[count][0] = i;
1102 index_ij[count][1] = j;
1103 index_ij[count][2] = k;
1106 "N_assemble_les_3d: active cells count %i at pos x[%i] y[%i] z[%i]\n",
1113 G_debug(2,
"N_assemble_les_3d: starting the parallel assemble loop");
1115 #pragma omp parallel for private(i, j, k, pos, count) schedule(static)
1116 for (count = 0; count < cell_type_count; count++) {
1117 i = index_ij[count][0];
1118 j = index_ij[count][1];
1119 k = index_ij[count][2];
1134 les->
b[count] = items->
V;
1145 les->
A[count][count] = items->
C;
1150 make_les_entry_3d(i, j, k, -1, 0, 0, count, pos, les, spvect,
1151 cell_count, status, start_val, items->
W,
1155 if (i < geom->
cols - 1) {
1156 pos = make_les_entry_3d(i, j, k, 1, 0, 0, count, pos, les, spvect,
1157 cell_count, status, start_val, items->
E,
1163 make_les_entry_3d(i, j, k, 0, -1, 0, count, pos, les, spvect,
1164 cell_count, status, start_val, items->
N,
1168 if (j < geom->rows - 1) {
1169 pos = make_les_entry_3d(i, j, k, 0, 1, 0, count, pos, les, spvect,
1170 cell_count, status, start_val, items->
S,
1176 if (k < geom->depths - 1) {
1178 make_les_entry_3d(i, j, k, 0, 0, 1, count, pos, les,
1179 spvect, cell_count, status, start_val,
1180 items->
T, cell_type);
1185 make_les_entry_3d(i, j, k, 0, 0, -1, count, pos, les,
1186 spvect, cell_count, status, start_val,
1187 items->
B, cell_type);
1193 spvect->
cols = pos + 1;
1203 for (i = 0; i < cell_type_count; i++)
1240 int rows,
cols, depths;
1242 int i, j, x,
y, z, stat;
1247 "N_les_integrate_dirichlet_3d: integrating the dirichlet boundary condition");
1254 dvect1 = (
double *)G_calloc(les->
cols,
sizeof(
double));
1255 dvect2 = (
double *)G_calloc(les->
cols,
sizeof(
double));
1259 for (z = 0; z < depths; z++) {
1260 for (y = 0; y < rows; y++) {
1261 for (x = 0; x <
cols; x++) {
1269 dvect1[count] = 0.0;
1276 #pragma omp parallel default(shared)
1283 #pragma omp for schedule (static) private(i)
1284 for (i = 0; i < les->
cols; i++)
1285 les->
b[i] = les->
b[i] - dvect2[i];
1291 for (z = 0; z < depths; z++) {
1292 for (y = 0; y < rows; y++) {
1293 for (x = 0; x <
cols; x++) {
1298 for (i = 0; i < les->
Asp[count]->
cols; i++)
1301 for (i = 0; i < les->
rows; i++) {
1302 for (j = 0; j < les->
Asp[i]->
cols; j++) {
1303 if (les->
Asp[i]->
index[j] == count)
1314 for (i = 0; i < les->
cols; i++)
1315 les->
A[count][i] = 0.0;
1317 for (i = 0; i < les->
rows; i++)
1318 les->
A[i][count] = 0.0;
1321 les->
A[count][count] = 1.0;
1336 int make_les_entry_3d(
int i,
int j,
int k,
int offset_i,
int offset_j,
1337 int offset_k,
int count,
int pos,
N_les * les,
1340 double entry,
int cell_type)
1360 if ((count + K) >= 0 && (count + K) < les->
cols) {
1362 " make_les_entry_3d: (N_CELL_ACTIVE) create matrix entry at row[%i] col[%i] value %g\n",
1363 count, count + K, entry);
1370 les->
A[count][count + K] = entry;
1378 if ((count + K) >= 0 && (count + K) < les->
cols) {
1380 " make_les_entry_3d: (N_CELL_DIRICHLET) create matrix entry at row[%i] col[%i] value %g\n",
1381 count, count + K, entry);
1388 les->
A[count][count + K] = entry;