Отчёт_неявная_схема (Весенний семестр), страница 5

They're calculated as field_param * processcoord + local cell coord.#pragma omp parallel forfor (int length = 0; length < field_length; ++length) {for (int width = 0; width < field_width; ++width) {array_data[field_width * length + width] =// Build test surface.exponential_coeff * sin((double)(x_offset + length) *x_step) ++ exponential_coeff * sin((double)(y_offset + width) *y_step);}}}// Perform a single iteration.void Trans_solver::Iterate(const int iter_num) {// Global offsets.int x_offset = field_length * coords[0];int y_offset = field_width * coords[1];//Set borders if process is responsible for the edge rectangle.exponential_coeff = pow(M_E, -(iter_num + 1) * DEFAULT_T_STEP);if (coords[1] == 0) { for (int i = 0; i < field_length; ++i) left_border[i] =exponential_coeff * sin(double(x_offset + i) * x_step); } // Set left edgeif (coords[1] == dim_sizes[1] - 1) { for (int i = 0; i < field_length; ++i)right_border[i] = exponential_coeff * sin(double(x_offset + i) * x_step); } // Set rightedgeif (coords[0] == 0) { for (int i = 0; i < field_width; ++i) upper_border[i] =exponential_coeff * sin(double(y_offset + i) * y_step); } // Set upper edgeif (coords[0] == dim_sizes[0] - 1) { for (int i = 0; i < field_width; ++i)lower_border[i] = exponential_coeff * sin(double(y_offset + i) * y_step); } // Set loweredge// This coeddicient is used in string subtraction process.double subtract_coeff_y_axis = 0;// First step is "width" solution// --------------------------------------------------------------// --------------------------------------------------------------for (int i = 0; i < field_length; ++i) { // Solution is done separately for everyhorizontal line.// Don't forget to reset appendix from previous step.appendix = 0;auxillary = 0;// Initial action is calculating all coefficients of the tridiagonalequation.

Note that A,B,C are same for all equations and only F differs.#pragma omp parallel num_threads(4){29#pragma omp forfor (int j = 0; j < field_width; ++j) {// Pass A, B, C coeffs to the tridiagonal equation.alpha_coeff_y_axis[j] = A_coeff_y_axis;charlie_coeff_y_axis[j] = C_coeff_y_axis;beta_coeff_y_axis[j] = B_coeff_y_axis;// Calculate F_coeff.F_coeff_y_axis[j] =- array_data[i * field_width + j] * (DEFAULT_RO* DEFAULT_C / (DEFAULT_T_STEP)) generate_heat(iter_num, i, j);// Left and right data buffers contain only zeros bydefault.left_data_buffer[j] = 0;right_data_buffer[j] = 0;}}// Prepare left data buffer for storing "defective" lines.left_data_buffer[0] = alpha_coeff_y_axis[0];// In case our processor is responsible for the edge rectangle,// we can adjust F_coefficient right now according to REAL border values.// This also means that left data buffer will not be actually used.

Westill keep it, since other processors will be busy with them anyway.if (coords[1] == 0) {F_coeff_y_axis[0] -= A_coeff_y_axis * left_border[i];alpha_coeff_y_axis[0] = 0;left_data_buffer[0] = 0;}if (coords[1] == dim_sizes[1] - 1) {F_coeff_y_axis[field_width - 1] -= B_coeff_y_axis *right_border[i];beta_coeff_y_axis[field_width - 1] = 0;}// Now we have to solve the equation system, which is done in severalbasic steps.//------------------------------------------------------------------------------// At first step, we subtract equation lines top-down; this process makesall alpha coeffs = 0 and adjusts f_coeffs.

However,// it also causes left data buffers to become non-zero and thus makingthose "defective" lines.for (int j = 1; j < field_width; ++j) {// Subtraction coefficient calculation.subtract_coeff_y_axis = alpha_coeff_y_axis[j] /charlie_coeff_y_axis[j - 1];// Solution vector.F_coeff_y_axis[j] -= subtract_coeff_y_axis * F_coeff_y_axis[j -301];// Alpha and charlie vectors. Beta vector doesn't change duringthis step.alpha_coeff_y_axis[j] = 0;charlie_coeff_y_axis[j] -= beta_coeff_y_axis[j - 1] *subtract_coeff_y_axis;// Left data buffer gets filled with defective data :(left_data_buffer[j] -= left_data_buffer[j - 1] *subtract_coeff_y_axis;}// At second step same sequence of operations is done in reversal.// This nullifies beta coefficient and fills right data buffers withunnecessary data.// Right data buffer values are initialized here;if (field_width >= 2) {right_data_buffer[field_width - 1] =charlie_coeff_y_axis[field_width - 1];right_data_buffer[field_width - 2] =beta_coeff_y_axis[field_width - 2];}// This also changes data in left data buffer !!!for (int j = field_width - 3; j >= 0; --j) {// Subtraction coefficient calculation.subtract_coeff_y_axis = beta_coeff_y_axis[j] /charlie_coeff_y_axis[j + 1];// Solution vector.F_coeff_y_axis[j] -= subtract_coeff_y_axis * F_coeff_y_axis[j +1];// Beta vector becomes zero, charlie vector is not affected( alpha is zero ).beta_coeff_y_axis[j] -= subtract_coeff_y_axis *charlie_coeff_y_axis[j + 1];charlie_coeff_y_axis[j] -= subtract_coeff_y_axis *alpha_coeff_y_axis[j + 1];// Right data buffer is filled with defective data.

Left databuffer is adjusted.left_data_buffer[j] -= left_data_buffer[j + 1] *subtract_coeff_y_axis;right_data_buffer[j] -= right_data_buffer[j + 1] *subtract_coeff_y_axis;}// Now the tricky part: we take the upper line, take CHARLIE, LEFT, RIGHTand SOLUTION// values from there and send them to the previous process, if there'sone.31send_buf[0]send_buf[1]send_buf[2]send_buf[3]====left_data_buffer[0];charlie_coeff_y_axis[0];right_data_buffer[0];F_coeff_y_axis[0];if (dim_sizes[1] > 1) { // Send only there is somebody to receive.MPI_Cart_shift(cart_comm, 1, -1, &shift_source, &shift_dest);MPI_Sendrecv(send_buf, 4, MPI_DOUBLE, shift_dest, 0, recv_buf,4, MPI_DOUBLE, shift_source, 0, cart_comm, &status);}else {for (int k = 0; k <= 3; ++k) { recv_buf[k] = send_buf[k]; }}// Then we adjust actual values according to what we received and fillONE line of a global matrix.if (coords[1] != dim_sizes[1] - 1) { // It's not the "rightest" process.// Calculate subtraction coefficient.subtract_coeff_y_axis = beta_coeff_y_axis[field_width - 1] /recv_buf[1];// Calculate actual new values.charlie_coeff_y_axis[field_width - 1] -= subtract_coeff_y_axis *recv_buf[0];beta_coeff_y_axis[field_width - 1] = 0;appendix = -subtract_coeff_y_axis * recv_buf[2];F_coeff_y_axis[field_width - 1] -= subtract_coeff_y_axis *recv_buf[3];}// Form a line for the global matrix.send_buf[0] = left_data_buffer[field_width - 1];send_buf[1] = charlie_coeff_y_axis[field_width - 1];send_buf[2] = appendix;send_buf[3] = F_coeff_y_axis[field_width - 1];if (dim_sizes[1] > 1) { // Send only there is somebody to receive.if (coords[1] == 0) { // Now processes with coords[1] == 0receive data to global matrix.dub_coords[0] = coords[0];for (int j = 0; j <= 3; ++j) { global_y_mtrx[j] =send_buf[j]; }for (int j = 1; j < dim_sizes[1]; ++j) {dub_coords[1] = j;MPI_Cart_rank(cart_comm, dub_coords,&shift_source);MPI_Recv(&global_y_mtrx[4 * j], 4, MPI_DOUBLE,shift_source, 0, cart_comm, &status);}}else { // And other processes send data to them.dub_coords[0] = coords[0];dub_coords[1] = 0;32MPI_Cart_rank(cart_comm, dub_coords, &shift_dest);MPI_Send(send_buf, 4, MPI_DOUBLE, shift_dest, 0,cart_comm);}}else {for (int k = 0; k <= 3; ++k) { global_y_mtrx[k] = send_buf[k]; }}// Simple version of the method is used for solving the equation withglobal matrix.simple_solve(global_y_mtrx, global_y_array, dim_sizes[1]);// Results are sent to processes, directly into appopriate locations oftheir dub_arrays.if (dim_sizes[1] > 1) { // There is somebody to receive.if (coords[1] == 0) {dub_array[i * field_width + (field_width - 1)] =global_y_array[0];dub_coords[0] = coords[0];for (int j = 1; j < dim_sizes[1]; ++j) {dub_coords[1] = j;MPI_Cart_rank(cart_comm, dub_coords,&shift_source);MPI_Send(&global_y_array[j], 1, MPI_DOUBLE,shift_source, 0, cart_comm);MPI_Send(&global_y_array[j - 1], 1, MPI_DOUBLE,shift_source, 0, cart_comm);}}else {MPI_Recv(&dub_array[i * field_width + (field_width 1)], 1, MPI_DOUBLE, shift_dest,0, cart_comm, &status);MPI_Recv(&auxillary, 1, MPI_DOUBLE, shift_dest,0, cart_comm, &status);}}else {dub_array[i * field_width + (field_width - 1)] =global_y_array[0];}// Finally, every process locally fills the remaining dub_array cells.for (int j = field_width - 2; j >= 0; --j) {dub_array[i * field_width + j] =(F_coeff_y_axis[j] - left_data_buffer[j] * auxillary dub_array[i * field_width + (field_width - 1)] * right_data_buffer[j])/ charlie_coeff_y_axis[j];}}// --------------------------------------------------------------// ---------------------------------------------------------------33// This coeddicient is used in string subtraction process.double subtract_coeff_x_axis = 0;// Second step is "length" solution// --------------------------------------------------------------// --------------------------------------------------------------for (int i = 0; i < field_width; ++i) { // Solution is done separately for everyvertical line.// Don't forget to reset appendix from previous step.appendix = 0;auxillary = 0;// Initial action is calculating all A,B,C,F coefficients of thetridiagonal equation.

Note that A,B,C are same for all equations and only F differs.// So we only need to calculate F coefficient.#pragma omp parallel{#pragma omp forfor (int j = 0; j < field_length; ++j) {// Pass A, B, C coeffs to the tridiagonal equation.alpha_coeff_x_axis[j] = A_coeff_x_axis;charlie_coeff_x_axis[j] = C_coeff_x_axis;beta_coeff_x_axis[j] = B_coeff_x_axis;// Calculate F_coeff.F_coeff_x_axis[j] =- dub_array[j * field_length + i] * (DEFAULT_RO* DEFAULT_C / (DEFAULT_T_STEP))- generate_heat(iter_num, j, i);// Left and right data buffers contain only zeros bydefault.up_data_buffer[j] = 0;down_data_buffer[j] = 0;}}// Prepare up data buffer for storing "defective" lines.up_data_buffer[0] = alpha_coeff_x_axis[0];// In case our processor is responsible for the edge rectangle,// we can adjust F_coefficient right now according to REAL border values.// This also means that up data buffer will not be actually used.