21 #include <grass/N_pde.h> 24 static int make_les_entry_2d(
int i,
int j,
int offset_i,
int offset_j,
31 static int make_les_entry_3d(
int i,
int j,
int k,
int offset_i,
int offset_j,
139 G_debug(5,
"N_create_5star: w %g e %g n %g s %g c %g v %g\n", star->
W,
140 star->
E, star->
N, star->
S, star->
C, star->
V);
162 double S,
double T,
double B,
double V)
177 G_debug(5,
"N_create_7star: w %g e %g n %g s %g t %g b %g c %g v %g\n",
178 star->
W, star->
E, star->
N, star->
S, star->
T, star->
B, star->
C,
203 double S,
double NW,
double SW,
double NE,
222 "N_create_9star: w %g e %g n %g s %g nw %g sw %g ne %g se %g c %g v %g\n",
223 star->
W, star->
E, star->
N, star->
S, star->
NW, star->
SW, star->
NE,
224 star->
SE, star->
C, star->
V);
266 double NW,
double SW,
double NE,
double SE,
267 double T,
double W_T,
double E_T,
double N_T,
268 double S_T,
double NW_T,
double SW_T,
269 double NE_T,
double SE_T,
double B,
double W_B,
270 double E_B,
double N_B,
double S_B,
double NW_B,
271 double SW_B,
double NE_B,
double SE_B,
double V)
311 "N_create_27star: w %g e %g n %g s %g nw %g sw %g ne %g se %g c %g v %g\n",
312 star->
W, star->
E, star->
N, star->
S, star->
NW, star->
SW, star->
NE,
313 star->
SE, star->
C, star->
V);
316 "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",
321 "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",
429 star->
E = 1 / geom->
dx;
430 star->
W = 1 / geom->
dx;
431 star->
N = 1 / geom->
dy;
432 star->
S = 1 / geom->
dy;
433 star->
T = 1 / geom->
dz;
434 star->
B = 1 / geom->
dz;
435 star->
C = -1 * (2 / geom->
dx + 2 / geom->
dy + 2 / geom->
dz);
439 "N_callback_template_3d: w %g e %g n %g s %g t %g b %g c %g v %g\n",
440 star->
W, star->
E, star->
N, star->
S, star->
T, star->
B, star->
C,
468 star->
E = 1 / geom->
dx;
469 star->
NE = 1 / sqrt(geom->
dx * geom->
dx + geom->
dy * geom->
dy);
470 star->
SE = 1 / sqrt(geom->
dx * geom->
dx + geom->
dy * geom->
dy);
471 star->
W = 1 / geom->
dx;
472 star->
NW = 1 / sqrt(geom->
dx * geom->
dx + geom->
dy * geom->
dy);
473 star->
SW = 1 / sqrt(geom->
dx * geom->
dx + geom->
dy * geom->
dy);
474 star->
N = 1 / geom->
dy;
475 star->
S = 1 / geom->
dy;
477 -1 * (star->
E + star->
NE + star->
SE + star->
W + star->
NW + star->
SW +
565 int i, j,
count = 0, pos = 0;
566 int cell_type_count = 0;
572 "N_assemble_les_2d: starting to assemble the linear equation system");
583 for (j = 0; j < geom->
rows; j++) {
584 for (i = 0; i < geom->
cols; i++) {
594 for (j = 0; j < geom->
rows; j++) {
595 for (i = 0; i < geom->
cols; i++) {
603 G_debug(2,
"N_assemble_les_2d: number of used cells %i\n",
606 if (cell_type_count == 0)
608 (
"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.",
613 index_ij = (
int **)
G_calloc(cell_type_count,
sizeof(
int *));
614 for (i = 0; i < cell_type_count; i++)
615 index_ij[i] = (
int *)
G_calloc(2,
sizeof(
int));
623 for (j = 0; j < geom->
rows; j++) {
624 for (i = 0; i < geom->
cols; i++) {
630 index_ij[
count][0] = i;
631 index_ij[
count][1] = j;
634 "N_assemble_les_2d: non-inactive cells count %i at pos x[%i] y[%i]\n",
641 index_ij[
count][0] = i;
642 index_ij[
count][1] = j;
645 "N_assemble_les_2d: active cells count %i at pos x[%i] y[%i]\n",
651 G_debug(2,
"N_assemble_les_2d: starting the parallel assemble loop");
654 #pragma omp parallel for private(i, j, pos, count) schedule(static) 655 for (count = 0; count < cell_type_count; count++) {
656 i = index_ij[
count][0];
657 j = index_ij[
count][1];
681 spvect->
values[pos] = items->
C;
688 pos = make_les_entry_2d(i, j, -1, 0, count, pos, les, spvect,
689 cell_count, status, start_val, items->
W,
693 if (i < geom->
cols - 1) {
694 pos = make_les_entry_2d(i, j, 1, 0, count, pos, les, spvect,
695 cell_count, status, start_val, items->
E,
701 make_les_entry_2d(i, j, 0, -1, count, pos, les, spvect,
702 cell_count, status, start_val, items->
N,
706 if (j < geom->
rows - 1) {
707 pos = make_les_entry_2d(i, j, 0, 1, count, pos, les, spvect,
708 cell_count, status, start_val, items->
S,
714 if (i > 0 && j > 0) {
715 pos = make_les_entry_2d(i, j, -1, -1, count, pos, les, spvect,
716 cell_count, status, start_val,
717 items->
NW, cell_type);
720 if (i < geom->
cols - 1 && j > 0) {
721 pos = make_les_entry_2d(i, j, 1, -1, count, pos, les, spvect,
722 cell_count, status, start_val,
723 items->
NE, cell_type);
726 if (i > 0 && j < geom->
rows - 1) {
727 pos = make_les_entry_2d(i, j, -1, 1, count, pos, les, spvect,
728 cell_count, status, start_val,
729 items->
SW, cell_type);
732 if (i < geom->
cols - 1 && j < geom->
rows - 1) {
733 pos = make_les_entry_2d(i, j, 1, 1, count, pos, les, spvect,
734 cell_count, status, start_val,
735 items->
SE, cell_type);
741 spvect->
cols = pos + 1;
752 for (i = 0; i < cell_type_count; i++)
791 int i, j,
x, y, stat;
796 "N_les_integrate_dirichlet_2d: integrating the dirichlet boundary condition");
807 for (y = 0; y <
rows; y++) {
808 for (x = 0; x <
cols; x++) {
821 #pragma omp parallel default(shared) 828 #pragma omp for schedule (static) private(i) 829 for (i = 0; i < les->
cols; i++)
830 les->
b[i] = les->
b[i] - dvect2[i];
836 for (y = 0; y <
rows; y++) {
837 for (x = 0; x <
cols; x++) {
845 for (i = 0; i < les->
rows; i++) {
846 for (j = 0; j < les->
Asp[i]->
cols; j++) {
858 for (i = 0; i < les->
cols; i++)
859 les->
A[count][i] = 0.0;
861 for (i = 0; i < les->
rows; i++)
862 les->
A[i][count] = 0.0;
865 les->
A[count][count] = 1.0;
880 int make_les_entry_2d(
int i,
int j,
int offset_i,
int offset_j,
int count,
883 N_array_2d * start_val,
double entry,
int cell_type)
901 if ((count + K) >= 0 && (count + K) < les->
cols) {
903 " make_les_entry_2d: (N_CELL_ACTIVE) create matrix entry at row[%i] col[%i] value %g\n",
904 count, count + K, entry);
907 spvect->
index[pos] = count + K;
908 spvect->
values[pos] = entry;
911 les->
A[
count][count + K] = entry;
921 if ((count + K) >= 0 && (count + K) < les->
cols) {
923 " make_les_entry_2d: (N_CELL_DIRICHLET) create matrix entry at row[%i] col[%i] value %g\n",
924 count, count + K, entry);
927 spvect->
index[pos] = count + K;
928 spvect->
values[pos] = entry;
931 les->
A[
count][count + K] = entry;
1017 int i, j, k, count = 0, pos = 0;
1018 int cell_type_count = 0;
1024 "N_assemble_les_3d: starting to assemble the linear equation system");
1035 for (k = 0; k < geom->
depths; k++) {
1036 for (j = 0; j < geom->
rows; j++) {
1037 for (i = 0; i < geom->
cols; i++) {
1050 for (k = 0; k < geom->
depths; k++) {
1051 for (j = 0; j < geom->
rows; j++) {
1052 for (i = 0; i < geom->
cols; i++) {
1064 "N_assemble_les_3d: number of used cells %i\n", cell_type_count);
1066 if (cell_type_count == 0.0)
1068 (
"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.",
1075 index_ij = (
int **)
G_calloc(cell_type_count,
sizeof(
int *));
1076 for (i = 0; i < cell_type_count; i++)
1077 index_ij[i] = (
int *)
G_calloc(3,
sizeof(
int));
1082 for (k = 0; k < geom->
depths; k++) {
1083 for (j = 0; j < geom->
rows; j++) {
1084 for (i = 0; i < geom->
cols; i++) {
1091 index_ij[
count][0] = i;
1092 index_ij[
count][1] = j;
1093 index_ij[
count][2] = k;
1096 "N_assemble_les_3d: non-inactive cells count %i at pos x[%i] y[%i] z[%i]\n",
1103 index_ij[
count][0] = i;
1104 index_ij[
count][1] = j;
1105 index_ij[
count][2] = k;
1108 "N_assemble_les_3d: active cells count %i at pos x[%i] y[%i] z[%i]\n",
1115 G_debug(2,
"N_assemble_les_3d: starting the parallel assemble loop");
1117 #pragma omp parallel for private(i, j, k, pos, count) schedule(static) 1118 for (count = 0; count < cell_type_count; count++) {
1119 i = index_ij[
count][0];
1120 j = index_ij[
count][1];
1121 k = index_ij[
count][2];
1144 spvect->
values[pos] = items->
C;
1152 make_les_entry_3d(i, j, k, -1, 0, 0, count, pos, les, spvect,
1153 cell_count, status, start_val, items->
W,
1157 if (i < geom->
cols - 1) {
1158 pos = make_les_entry_3d(i, j, k, 1, 0, 0, count, pos, les, spvect,
1159 cell_count, status, start_val, items->
E,
1165 make_les_entry_3d(i, j, k, 0, -1, 0, count, pos, les, spvect,
1166 cell_count, status, start_val, items->
N,
1170 if (j < geom->
rows - 1) {
1171 pos = make_les_entry_3d(i, j, k, 0, 1, 0, count, pos, les, spvect,
1172 cell_count, status, start_val, items->
S,
1178 if (k < geom->
depths - 1) {
1180 make_les_entry_3d(i, j, k, 0, 0, 1, count, pos, les,
1181 spvect, cell_count, status, start_val,
1182 items->
T, cell_type);
1187 make_les_entry_3d(i, j, k, 0, 0, -1, count, pos, les,
1188 spvect, cell_count, status, start_val,
1189 items->
B, cell_type);
1195 spvect->
cols = pos + 1;
1205 for (i = 0; i < cell_type_count; i++)
1244 int i, j,
x, y, z, stat;
1249 "N_les_integrate_dirichlet_3d: integrating the dirichlet boundary condition");
1256 dvect1 = (
double *)
G_calloc(les->
cols,
sizeof(
double));
1257 dvect2 = (
double *)
G_calloc(les->
cols,
sizeof(
double));
1261 for (z = 0; z <
depths; z++) {
1262 for (y = 0; y <
rows; y++) {
1263 for (x = 0; x <
cols; x++) {
1271 dvect1[
count] = 0.0;
1278 #pragma omp parallel default(shared) 1285 #pragma omp for schedule (static) private(i) 1286 for (i = 0; i < les->
cols; i++)
1287 les->
b[i] = les->
b[i] - dvect2[i];
1293 for (z = 0; z <
depths; z++) {
1294 for (y = 0; y <
rows; y++) {
1295 for (x = 0; x <
cols; x++) {
1303 for (i = 0; i < les->
rows; i++) {
1304 for (j = 0; j < les->
Asp[i]->
cols; j++) {
1305 if (les->
Asp[i]->
index[j] == count)
1316 for (i = 0; i < les->
cols; i++)
1317 les->
A[count][i] = 0.0;
1319 for (i = 0; i < les->
rows; i++)
1320 les->
A[i][count] = 0.0;
1323 les->
A[count][count] = 1.0;
1338 int make_les_entry_3d(
int i,
int j,
int k,
int offset_i,
int offset_j,
1339 int offset_k,
int count,
int pos,
N_les * les,
1342 double entry,
int cell_type)
1362 if ((count + K) >= 0 && (count + K) < les->
cols) {
1364 " make_les_entry_3d: (N_CELL_ACTIVE) create matrix entry at row[%i] col[%i] value %g\n",
1365 count, count + K, entry);
1368 spvect->
index[pos] = count + K;
1369 spvect->
values[pos] = entry;
1372 les->
A[
count][count + K] = entry;
1380 if ((count + K) >= 0 && (count + K) < les->
cols) {
1382 " make_les_entry_3d: (N_CELL_DIRICHLET) create matrix entry at row[%i] col[%i] value %g\n",
1383 count, count + K, entry);
1386 spvect->
index[pos] = count + K;
1387 spvect->
values[pos] = entry;
1390 les->
A[
count][count + K] = entry;
N_data_star *(* callback)()
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
N_data_star * N_create_9star(double C, double W, double E, double N, double S, double NW, double SW, double NE, double SE, double V)
allocate and initialize a 9 point star data structure
Matrix entries for a mass balance 5/7/9 star system.
N_les * N_alloc_les_Ax_b(int rows, int type)
Allocate memory for a quadratic linear equation system which includes the Matrix A, vector x and vector b.
N_les * N_assemble_les_2d_dirichlet(int les_type, N_geom_data *geom, N_array_2d *status, N_array_2d *start_val, void *data, N_les_callback_2d *call)
Assemble a linear equation system (les) based on 2d location data (raster) and active and dirichlet c...
G_math_spvector * G_math_alloc_spvector(int)
Allocate memory for a sparse vector.
The row vector of the sparse matrix.
int N_les_integrate_dirichlet_2d(N_les *les, N_geom_data *geom, N_array_2d *status, N_array_2d *start_val)
Integrate Dirichlet or Transmission boundary conditions into the les (2s)
double N_get_array_3d_d_value(N_array_3d *data, int col, int row, int depth)
This function returns the value of type float at position col, row, depth.
N_array_3d * N_alloc_array_3d(int cols, int rows, int depths, int offset, int type)
Allocate memory for a N_array_3d data structure.
N_array_2d * N_alloc_array_2d(int cols, int rows, int offset, int type)
Allocate memory for a N_array_2d data structure.
N_les_callback_3d * N_alloc_les_callback_3d(void)
Allocate the structure holding the callback function.
N_data_star * N_alloc_5star(void)
allocate a 5 point star data structure
void G_free(void *)
Free allocated memory.
#define N_MAX_CELL_STATE
the maximum number of available cell states (eg: boundary condition, inactiven active) ...
int N_les_integrate_dirichlet_3d(N_les *les, N_geom_data *geom, N_array_3d *status, N_array_3d *start_val)
Integrate Dirichlet or Transmission boundary conditions into the les (3d)
N_les * N_assemble_les_2d(int les_type, N_geom_data *geom, N_array_2d *status, N_array_2d *start_val, void *data, N_les_callback_2d *call)
Assemble a linear equation system (les) based on 2d location data (raster) and active cells...
void N_put_array_2d_c_value(N_array_2d *data, int col, int row, CELL value)
Writes a CELL value to the N_array_2d struct at position col, row.
callback structure for 2d matrix assembling
N_data_star * N_callback_template_2d(void *data, N_geom_data *geom, int col, int row)
A callback template creates a 9 point star structure.
void G_math_Ax_sparse(G_math_spvector **, double *, double *, int)
Compute the matrix - vector product of sparse matrix **Asp and vector x.
N_les * N_assemble_les_2d_param(int les_type, N_geom_data *geom, N_array_2d *status, N_array_2d *start_val, void *data, N_les_callback_2d *call, int cell_type)
Assemble a linear equation system (les) based on 2d location data (raster)
N_les * N_assemble_les_3d_dirichlet(int les_type, N_geom_data *geom, N_array_3d *status, N_array_3d *start_val, void *data, N_les_callback_3d *call)
Assemble a linear equation system (les) based on 3d location data (g3d) active and dirichlet cells...
Geometric information about the structured grid.
N_les * N_assemble_les_3d(int les_type, N_geom_data *geom, N_array_3d *status, N_array_3d *start_val, void *data, N_les_callback_3d *call)
Assemble a linear equation system (les) based on 3d location data (g3d) active cells.
void N_free_array_2d(N_array_2d *data)
Release the memory of a N_array_2d structure.
N_data_star * N_alloc_9star(void)
allocate a 9 point star data structure
N_les_callback_2d * N_alloc_les_callback_2d(void)
Allocate the structure holding the callback function.
N_data_star * N_create_27star(double C, double W, double E, double N, double S, double NW, double SW, double NE, double SE, double T, double W_T, double E_T, double N_T, double S_T, double NW_T, double SW_T, double NE_T, double SE_T, double B, double W_B, double E_B, double N_B, double S_B, double NW_B, double SW_B, double NE_B, double SE_B, double V)
allocate and initialize a 27 point star data structure
void N_put_array_3d_d_value(N_array_3d *data, int col, int row, int depth, double value)
Writes a double value to the N_array_3d struct at position col, row, depth.
N_data_star * N_alloc_7star(void)
allocate a 7 point star data structure
void N_set_les_callback_2d_func(N_les_callback_2d *data, N_data_star *(*callback_func_2d)())
Set the callback function which is called while assembling the les in 2d.
N_data_star * N_create_7star(double C, double W, double E, double N, double S, double T, double B, double V)
allocate and initialize a 7 point star data structure
N_data_star * N_alloc_27star(void)
allocate a 27 point star data structure
int depths
number of depths for 3D data
N_les * N_assemble_les_3d_active(int les_type, N_geom_data *geom, N_array_3d *status, N_array_3d *start_val, void *data, N_les_callback_3d *call)
Assemble a linear equation system (les) based on 3d location data (g3d) active cells.
int cols
Number of columns for 2D data.
CELL N_get_array_2d_c_value(N_array_2d *data, int col, int row)
Returns the value of type CELL at position col, row.
void N_set_les_callback_3d_func(N_les_callback_3d *data, N_data_star *(*callback_func_3d)())
Set the callback function which is called while assembling the les in 3d.
N_data_star *(* callback)()
N_les * N_assemble_les_2d_active(int les_type, N_geom_data *geom, N_array_2d *status, N_array_2d *start_val, void *data, N_les_callback_2d *call)
Assemble a linear equation system (les) based on 2d location data (raster) and active cells...
N_les * N_assemble_les_3d_param(int les_type, N_geom_data *geom, N_array_3d *status, N_array_3d *start_val, void *data, N_les_callback_3d *call, int cell_type)
Assemble a linear equation system (les) based on 3d location data (g3d)
void G_math_d_Ax(double **, double *, double *, int, int)
Compute the matrix - vector product of matrix A and vector x.
callback structure for 3d matrix assembling
DCELL N_get_array_2d_d_value(N_array_2d *data, int col, int row)
Returns the value of type DCELL at position col, row.
void N_free_array_3d(N_array_3d *data)
Release the memory of a N_array_3d.
int rows
Number of rows for 2D data.
int G_math_add_spvector(G_math_spvector **, G_math_spvector *, int)
Adds a sparse vector to a sparse matrix at position row.
int G_debug(int, const char *,...) __attribute__((format(printf
N_data_star * N_create_5star(double C, double W, double E, double N, double S, double V)
allocate and initialize a 5 point star data structure
N_data_star * N_callback_template_3d(void *data, N_geom_data *geom, int col, int row, int depth)
A callback template creates a 7 point star structure.
The linear equation system (les) structure.