BookGPU/Chapters/chapter5/src/chapter5/heat_equation/main.cu

   1 ////////////////////////////////////////////////////////////////////////////////\r
   2 // HEAT EQUATION TEST\r
   3 ////////////////////////////////////////////////////////////////////////////////\r
   4 \r
   5 #include <gpulab/gpulab.h>\r
   6 #include <gpulab/grid.h>\r
   7 #include <gpulab/FD/stencil.h>\r
   8 #include <gpulab/integration/ode_solver.h>\r
   9 #include <gpulab/integration/forward_euler.h>\r
  10 #include <gpulab/io/print.h>\r
  11 \r
  12 #include "../../../code/laplacian.cu"\r
  13 \r
  14 #define PI 3.141592653589793\r
  15 \r
  16 // Define and assmble heat solver components here\r
  17 typedef float                                                                                           value_type;\r
  18 typedef laplacian<value_type>                                                           rhs_type;\r
  19 typedef gpulab::grid<value_type,gpulab::device_memory>          vector_type;\r
  20 typedef gpulab::grid<value_type,gpulab::host_memory>            host_vector_type;\r
  21 typedef vector_type::property_type                                                      property_type;\r
  22 typedef gpulab::integration::forward_euler                                      time_integrator_type;\r
  23 \r
  24 int main(int argc, char *argv[])\r
  25 {\r
  26         // Initialize gpulab\r
  27         gpulab::init(argc,argv);\r
  28         GPULAB_LOG_INF("=== BEGIN HEAT EQUATION SOLVER ===\n");\r
  29 \r
  30         // Get grid dimension from user unput\r
  31         // TODO: Should we use cofig_reader instead?\r
  32         int alpha               = (argc > 1) ? atoi(argv[1]) : 1;\r
  33         int N                   = (argc > 2) ? atoi(argv[2]) : 16;\r
  34         value_type tend = (argc > 3) ? atof(argv[3]) : 1.;\r
  35         \r
  36         // Setup grids\r
  37         gpulab::grid_dim<int>        dim(N,N,1);       // Discrete dimension\r
  38         gpulab::grid_dim<value_type> p0(0,0);          // Physical dimension start values\r
  39         gpulab::grid_dim<value_type> p1(1,1);          // Physical dimension end values\r
  40         property_type                props(dim,p0,p1); // Gather in property class\r
  41         host_vector_type             u0(props);        // Initialize host vector u0 for initial values\r
  42         host_vector_type             ut(props);        // Initialize host vector ut for exact solution\r
  43         \r
  44         value_type dx = min(u0.local_props().delta.x,u0.local_props().delta.y);\r
  45         value_type dt = 1./(alpha*alpha*10.0)*dx*dx;            // Time step. For explicit Euler and alpha=1, dt <= 0.5*dx^2\r
  46         GPULAB_LOG_INF("Using dt=%e\n",dt);\r
  47 \r
  48         // Initialize u_h\r
  49         for(int j=0; j<u0.Ny(); ++j)\r
  50         {\r
  51                 for(int i=0; i<u0.Nx(); ++i)\r
  52                 {\r
  53                         value_type x = i/(value_type)(u0.Nx()-1);\r
  54                         value_type y = j/(value_type)(u0.Ny()-1);\r
  55                         u0[j*u0.Nx()+i] = sin(PI*x)*sin(PI*y);\r
  56                         ut[j*ut.Nx()+i] = exp(-2.0*PI*PI*tend)*sin(PI*x)*sin(PI*y);\r
  57                 }\r
  58         }\r
  59         vector_type u(u0);                                                                                      // Initialize device vector u from host vector u0\r
  60 \r
  61         // Create time integrator and take time steps\r
  62         time_integrator_type solver;                                                            // Create time integrator\r
  63         rhs_type rhs(alpha);                                                                            // Create right hand side operator\r
  64         solver(&rhs,u,0.0f,tend,dt);                                                            // Integrate t=0 to tend using dt\r
  65 \r
  66         gpulab::io::print(ut,gpulab::io::TO_BINARY_FILE,1,"ut.bin");// Print initial state\r
  67         \r
  68         ut.axpy(-1.0,u);\r
  69         GPULAB_LOG_INF("||e||_inf = %e\n", ut.nrmi());\r
  70 \r
  71         value_type enrm = ut.nrmi();\r
  72         gpulab::io::print(ut,gpulab::io::TO_BINARY_FILE,1,"e.bin");     // Print error\r
  73         gpulab::io::print(u0,gpulab::io::TO_BINARY_FILE,1,"u0.bin");// Print initial state\r
  74         gpulab::io::print(u, gpulab::io::TO_BINARY_FILE,1,"u.bin");     // Print final state\r
  75 \r
  76         gpulab::finalize();\r
  77     return 0;\r
  78 }