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

Onlydeclaration is here.// ------------------------------------------------------------------------------class Trans_solver{private:// Step parameters.double x_step;double y_step;double exponential_coeff;// Data field parameters.int field_length; // Data field length.int field_width; // Data field width.double* array_data; // A field of dots calculated.13double* dub_array; // A field of dots calculated duplicate.// Border parameters.double* upper_border; // Border data :)double* lower_border; // Border data :)double* left_border; // Border data :)double* right_border; // Border data :)double border_value; // Default border value.// Tridiagonal equation system coefficients for horizontal solution step.double* alpha_coeff_y_axis;double* charlie_coeff_y_axis;double* beta_coeff_y_axis;double* alpha_coeff_x_axis;double* charlie_coeff_x_axis;double* beta_coeff_x_axis;// Equation coefficients Ay(i-1) - Cy(i) + By(i+1) = - Fdouble A_coeff_x_axis;double B_coeff_x_axis;double C_coeff_x_axis;double A_coeff_y_axis;double B_coeff_y_axis;double C_coeff_y_axis;double* F_coeff_x_axis;double* F_coeff_y_axis;// This internal function is used to simulate heat generators.double generate_heat(const int iteration_number, const int x_pos, const inty_pos);public:// Constructor.Trans_solver();void Init(int data_length, int data_width, int p_rank, int p_size); // Initializesolver with processor cart.void Fill_data(const int iter_num); // Fill solver with data.void Iterate(const int iter_num); // Perform a local iteration.void Output(const int mode); // Output data.// Destructor.virtual ~Trans_solver();};// Project inclusions.#include "Trans_solver.h"// ------------------------------------------------------------------------------// This file contains functions used in trans_solver class// -------------------------------------------------------------------------------14// ------------------------------------------------------------------------------// ------------------------------------------------------------------------------// This internal function is used to simulate heat generators.double Trans_solver::generate_heat(const int iter_num, const int x_pos, const int y_pos) {// This function can give an additional heat source to any cell in the field.// As you can see, it depends on coordinates and iteration number.// So it's basically f(r,t).// For test purposes it's considered a simple given function, nondependant ofposition// This doesn't change any of calculation principles, since these additional heatsources always provide a known value.// Actual function.//double result = (double)iter_num * DEFAULT_T_STEP + DEFAULT_T_STEP / 2.0;//return 0.2 * cos(((double)result / DEFAULT_T_STEP) / 20.0);// Test variant for checking.return 0.0;}// ------------------------------------------------------------------------------// ------------------------------------------------------------------------------// Default constructor.Trans_solver::Trans_solver() {}// Initialize solver ( make a cart of processes ).void Trans_solver::Init(int data_length, int data_width, int p_rank, int p_size) {// Calculate field parameters and initialize field data.field_length = data_length;field_width = data_width;array_data = new double[field_length * field_width];dub_array = new double[field_length * field_width];// Allocate memory for border fields.// Border values will be set later during data filling.upper_border = new double[field_length];lower_border = new double[field_length];left_border = new double[field_width];right_border = new double[field_width];// Allocate memory for Tridiagonal coefficients.alpha_coeff_y_axis = new double[field_width];charlie_coeff_y_axis = new double[field_width];beta_coeff_y_axis = new double[field_width];alpha_coeff_x_axis = new double[field_length];charlie_coeff_x_axis = new double[field_length];beta_coeff_x_axis = new double[field_length];// Calculate step sizes.x_step = (2.0 * DEFAULT_X_STEP) / (double)(data_length - 1);15y_step = (2.0 * DEFAULT_Y_STEP) / (double)(data_width - 1);// Calculate static equation coefficients.A_coeff_x_axis = DEFAULT_LAMBDA / (x_step * x_step);B_coeff_x_axis = DEFAULT_LAMBDA / (x_step * x_step);C_coeff_x_axis = -2.0 * (DEFAULT_LAMBDA / (x_step * x_step)) - (DEFAULT_RO *DEFAULT_C) / DEFAULT_T_STEP;A_coeff_y_axis = DEFAULT_LAMBDA / (y_step * y_step);B_coeff_y_axis = DEFAULT_LAMBDA / (y_step * y_step);C_coeff_y_axis = -2.0 * (DEFAULT_LAMBDA / (y_step * y_step)) - (DEFAULT_RO *DEFAULT_C) / DEFAULT_T_STEP;// Allocate memory for equation coefficients.F_coeff_x_axis = new double[field_length];F_coeff_y_axis = new double[field_width];}// Fill solver with test source data.void Trans_solver::Fill_data(const int iter_num) {// Formula used is sqrt((x - global_length/2)(y - global_width/2)) +sqrt((global_length/2)(global_width/2)) + 1// X and Y are global coordinates.

They're calculated as field_param * processcoord + local cell coord.exponential_coeff = pow(M_E, -(iter_num) * DEFAULT_T_STEP);for (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)(length) * x_step) ++exponential_coeff * sin((double)(width) * y_step);}}}// Perform a single iteration.void Trans_solver::Iterate(const int iter_num) {// Set new border values ( we know a precise solution ).exponential_coeff = pow(M_E, -(iter_num + 1) * DEFAULT_T_STEP);for (int i = 0; i < field_length; ++i) left_border[i] = exponential_coeff *sin(double(i) * x_step); // Left edgefor (int i = 0; i < field_length; ++i) right_border[i] = exponential_coeff *sin(double(i) * x_step); // Right edgefor (int i = 0; i < field_width; ++i) upper_border[i] = exponential_coeff *sin(double(i) * y_step); // Upper edgefor (int i = 0; i < field_width; ++i) lower_border[i] = exponential_coeff *sin(double(i) * y_step); // Lower edge.// This coeddicient is used in string subtraction process.double subtract_coeff_y_axis = 0;// First step is "width solution".// --------------------------------------------------------------// ---------------------------------------------------------------16for (int i = 0; i < field_length; ++i) { // Solution is done separately for everyhorizontal line.// Initial action is calculating all coefficients of the tridiagonalequation.

Note that A,B,C are same for all equations and only F differs.for (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);}// 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.F_coeff_y_axis[0] -= A_coeff_y_axis * left_border[i];alpha_coeff_y_axis[0] = 0;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.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 1];// Alpha and charlie vectors.

Beta vector doesn't change duringthis step.alpha_coeff_y_axis[j] = 0; // -= subtract_coeff_y_axis *charlie_coeff_y_axis[j - 1];charlie_coeff_y_axis[j] -= beta_coeff_y_axis[j - 1] *subtract_coeff_y_axis;}// At second step same sequence of operations is done in reversal.for (int j = field_width - 2; j >= 0; --j) {17// 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 should bezero ).// charlie_coeff_y_axis[j] -= subtract_coeff_y_axis *alpha_coeff_y_axis[j + 1];beta_coeff_y_axis[j] = 0;}// Calculate actual values now.for (int j = field_width - 1; j >= 0; --j) {dub_array[i * field_width + j] = F_coeff_y_axis[j] /charlie_coeff_y_axis[j];}}// --------------------------------------------------------------// --------------------------------------------------------------// 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.// 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.for (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);}// In case our processor is responsible for the edge rectangle,// we can adjust F_coefficient right now according to REAL border values.F_coeff_x_axis[0] -= A_coeff_x_axis * upper_border[i];alpha_coeff_x_axis[0] = 0;18F_coeff_x_axis[field_length - 1] -= B_coeff_x_axis * lower_border[i];beta_coeff_x_axis[field_length - 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,for (int j = 1; j < field_length; ++j) {// Subtraction coefficient calculation.subtract_coeff_x_axis = alpha_coeff_x_axis[j] /charlie_coeff_x_axis[j - 1];// Solution vector.F_coeff_x_axis[j] -= subtract_coeff_x_axis * F_coeff_x_axis[j 1];// Alpha and charlie vectors. Beta vector doesn't change duringthis step.alpha_coeff_x_axis[j] = 0; // -= subtract_coeff_x_axis *charlie_coeff_x_axis[j - 1];charlie_coeff_x_axis[j] -= beta_coeff_x_axis[j - 1] *subtract_coeff_x_axis;}// At second step same sequence of operations is done in reversal.for (int j = field_length - 2; j >= 0; --j) {// Subtraction coefficient calculation.subtract_coeff_x_axis = beta_coeff_x_axis[j] /charlie_coeff_x_axis[j + 1];// Solution vector.F_coeff_x_axis[j] -= subtract_coeff_x_axis * F_coeff_x_axis[j +1];// Beta vector becomes zero, charlie vector is not affected( alpha is zero ).// charlie_coeff_x_axis[j] -= subtract_coeff_x_axis *alpha_coeff_x_axis[j + 1];beta_coeff_x_axis[j] = 0;}// Finally, every process locally fills the remaining dub_array cells.for (int j = field_length - 1; j >= 0; --j) {array_data[j * field_length + i] = F_coeff_x_axis[j] /charlie_coeff_x_axis[j];}}// As horizontal solution is over, we should have fully updated data ready forthe next step.