20 #include <grass/N_solute_transport.h> 32 double Df_e = 0, Df_w = 0, Df_n = 0, Df_s = 0, Df_t = 0, Df_b = 0;
33 double dx, dy, dz, Az;
34 double diff_x, diff_y, diff_z;
35 double diff_xw, diff_yn;
36 double diff_xe, diff_ys;
37 double diff_zt, diff_zb;
38 double cin = 0, cg, cg_start;
40 double C,
W, E,
N, S, T, B, V;
41 double vw = 0, ve = 0, vn = 0, vs = 0, vt = 0, vb = 0;
42 double Ds_w = 0, Ds_e = 0, Ds_n = 0, Ds_s = 0, Ds_t = 0, Ds_b = 0;
43 double Dw = 0, De = 0, Dn = 0, Ds = 0, Dt = 0, Db = 0;
44 double rw = 0.5, re = 0.5, rn = 0.5, rs = 0.5, rt = 0.5, rb = 0.5;
95 Dw = ((Df_w + Ds_w)) / dx;
96 De = ((Df_e + Ds_e)) / dx;
97 Dn = ((Df_n + Ds_n)) / dy;
98 Ds = ((Df_s + Ds_s)) / dy;
99 Dt = ((Df_t + Ds_t)) / dz;
100 Db = ((Df_b + Ds_b)) / dz;
110 W = -1 * (Dw) * dy * dz - vw * (1 - rw) * dy * dz;
112 E = -1 * (De) * dy * dz + ve * (1 - re) * dy * dz;
114 S = -1 * (Ds) * dx * dz - vs * (1 - rs) * dx * dz;
116 N = -1 * (Dn) * dx * dz + vn * (1 - rn) * dx * dz;
118 B = -1 * (Db) * Az - vb * (1 - rb) * Az;
120 T = -1 * (Dt) * Az + vt * (1 - rt) * Az;
134 C = ((Dw - vw) * dy * dz +
135 (De + ve) * dy * dz +
136 (Ds - vs) * dx * dz +
137 (Dn + vn) * dx * dz +
138 (Db - vb) * Az + (Dt + vt) * Az + Az * dz * R / data->
dt - q / nf);
141 V = (cs + cg_start * Az * dz * R / data->
dt - q / nf * cin);
155 G_debug(6,
"N_callback_solute_transport_3d: called [%i][%i][%i]", row,
188 double Df_e = 0, Df_w = 0, Df_n = 0, Df_s = 0;
189 double z_e = 0, z_w = 0, z_n = 0, z_s = 0;
191 double diff_x, diff_y;
192 double disp_x, disp_y;
194 double diff_xw, diff_yn;
195 double disp_xw, disp_yn;
197 double diff_xe, diff_ys;
198 double disp_xe, disp_ys;
200 double cin = 0, cg, cg_start;
203 double vw = 0, ve = 0, vn = 0, vs = 0;
204 double Ds_w = 0, Ds_e = 0, Ds_n = 0, Ds_s = 0;
205 double Dw = 0, De = 0, Dn = 0, Ds = 0;
206 double rw = 0.5, re = 0.5, rn = 0.5, rs = 0.5;
306 Dw = ((Df_w + Ds_w)) / dx;
307 De = ((Df_e + Ds_e)) / dx;
308 Ds = ((Df_s + Ds_s)) / dy;
309 Dn = ((Df_n + Ds_n)) / dy;
330 W = -1 * (Dw) * dy * z_w + vw * (1 - rw) * dy * z_w;
332 E = -1 * (De) * dy * z_e + ve * (1 - re) * dy * z_e;
334 S = -1 * (Ds) * dx * z_s + vs * (1 - rs) * dx * z_s;
336 N = -1 * (Dn) * dx * z_n + vn * (1 - rn) * dx * z_n;
355 C = (Dw + vw * rw) * dy * z_w +
356 (De + ve * re) * dy * z_e +
357 (Ds + vs * rs) * dx * z_s +
358 (Dn + vn * rn) * dx * z_n + Az * z * R / data->
dt - q / nf;
361 V = (cs + cg_start * Az * z * R / data->
dt + q / nf * cin);
376 G_debug(6,
"N_callback_solute_transport_2d: called [%i][%i]", row, col);
589 "N_calc_solute_transport_transmission_2d: calculating transmission boundary");
591 for (j = 0; j < rows; j++) {
592 for (i = 0; i < cols; i++) {
626 c = c / (double)count;
628 if (c > 0 || c == 0 || c < 0)
656 double disp_xx, disp_yy, disp_xy;
663 "N_calc_solute_transport_disptensor_2d: calculating the dispersivity tensor");
665 for (j = 0; j < rows; j++) {
666 for (i = 0; i < cols; i++) {
674 vx = (grad.
WC + grad.
EC) / 2;
675 vy = (grad.
NC + grad.
SC) / 2;
676 vv = sqrt(vx * vx + vy * vy);
679 disp_xx = data->
al * vx * vx / vv + data->
at * vy * vy / vv;
680 disp_yy = data->
at * vx * vx / vv + data->
al * vy * vy / vv;
681 disp_xy = (data->
al - data->
at) * vx * vy / vv;
685 "N_calc_solute_transport_disptensor_2d: [%i][%i] disp_xx %g disp_yy %g disp_xy %g",
686 i, j, disp_xx, disp_yy, disp_xy);
713 int cols, rows, depths;
714 double vx, vy, vz, vv;
715 double disp_xx, disp_yy, disp_zz, disp_xy, disp_xz, disp_yz;
723 "N_calc_solute_transport_disptensor_3d: calculating the dispersivity tensor");
725 for (k = 0; k < depths; k++) {
726 for (j = 0; j < rows; j++) {
727 for (i = 0; i < cols; i++) {
737 vx = (grad.
WC + grad.
EC) / 2;
738 vy = (grad.
NC + grad.
SC) / 2;
739 vz = (grad.
BC + grad.
TC) / 2;
740 vv = sqrt(vx * vx + vy * vy + vz * vz);
744 data->
al * vx * vx / vv + data->
at * vy * vy / vv +
745 data->
at * vz * vz / vv;
747 data->
at * vx * vx / vv + data->
al * vy * vy / vv +
748 data->
at * vz * vz / vv;
750 data->
at * vx * vx / vv + data->
at * vy * vy / vv +
751 data->
al * vz * vz / vv;
752 disp_xy = (data->
al - data->
at) * vx * vy / vv;
753 disp_xz = (data->
al - data->
at) * vx * vz / vv;
754 disp_yz = (data->
al - data->
at) * vy * vz / vv;
758 "N_calc_solute_transport_disptensor_3d: [%i][%i][%i] disp_xx %g disp_yy %g disp_zz %g disp_xy %g disp_xz %g disp_yz %g ",
759 i, j, k, disp_xx, disp_yy, disp_zz, disp_xy, disp_xz,
N_gradient_2d * N_get_gradient_2d(N_gradient_field_2d *field, N_gradient_2d *gradient, int col, int row)
Return a N_gradient_2d structure calculated from the input gradient field at position [row][col]...
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.
double N_calc_geom_mean(double a, double b)
Calculate the geometrical mean of values a and b.
#define N_CELL_TRANSMISSION
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.
Gradient between the cells in X and Y direction.
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_gradient_field_3d * N_alloc_gradient_field_3d(int cols, int rows, int depths)
Allocate a N_gradient_field_3d.
void G_free(void *)
Free allocated memory.
void N_calc_solute_transport_transmission_2d(N_solute_transport_data2d *data)
Compute the transmission boundary condition in 2d.
N_gradient_field_3d * grad
void N_calc_solute_transport_disptensor_3d(N_solute_transport_data3d *data)
Compute the dispersivity tensor based on the solute transport data in 3d.
Geometric information about the structured grid.
void N_free_array_2d(N_array_2d *data)
Release the memory of a N_array_2d 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.
double N_exp_upwinding(double sprod, double distance, double D)
exponential upwinding stabilization algorithm
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
int N_is_array_2d_value_null(N_array_2d *data, int col, int row)
Returns 1 if the value of N_array_2d struct at position col, row is of type null, otherwise 0...
Gradient between the cells in X, Y and Z direction.
double N_get_geom_data_area_of_cell(N_geom_data *geom, int row)
Get the areay size in square meter of one cell (x*y) at row.
void N_calc_solute_transport_disptensor_2d(N_solute_transport_data2d *data)
Compute the dispersivity tensor based on the solute transport data in 2d.
void N_free_solute_transport_data2d(N_solute_transport_data2d *data)
Release the memory of the solute transport data structure in two dimensions.
void N_free_solute_transport_data3d(N_solute_transport_data3d *data)
Release the memory of the solute transport data structure in three dimensions.
N_gradient_field_2d * grad
N_gradient_field_2d * N_alloc_gradient_field_2d(int cols, int rows)
Allocate a N_gradient_field_2d.
N_data_star * N_callback_solute_transport_2d(void *solutedata, N_geom_data *geom, int col, int row)
This callback function creates the mass balance of a 5 point star.
void N_put_array_2d_d_value(N_array_2d *data, int col, int row, DCELL value)
Writes a DCELL value to the N_array_2d struct at position col, row.
double N_calc_harmonic_mean(double a, double b)
Calculate the harmonical mean of values a and b.
N_data_star * N_callback_solute_transport_3d(void *solutedata, N_geom_data *geom, int col, int row, int depth)
This is just a placeholder.
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 G_debug(int, const char *,...) __attribute__((format(printf
double N_full_upwinding(double sprod, double distance, double D)
full upwinding stabilization algorithm
N_gradient_3d * N_get_gradient_3d(N_gradient_field_3d *field, N_gradient_3d *gradient, int col, int row, int depth)
Return a N_gradient_3d structure calculated from the input gradient field at position [depth][row][co...
N_solute_transport_data3d * N_alloc_solute_transport_data3d(int cols, int rows, int depths)
Alllocate memory for the solute transport data structure in three dimensions.
N_solute_transport_data2d * N_alloc_solute_transport_data2d(int cols, int rows)
Alllocate memory for the solute transport data structure in two dimensions.