API tools faq
paste
Advertisement
phystota

DRR

phystota
May 30th, 2025
380
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 4.15 KB | None | 0 0
raw download clone embed print report
  1. #include <iostream>
  2. #include <random>
  3. #include <fstream>
  4. #include <assert.h>
  5.  
  6. #include <math.h>
  7. #include <iomanip>  // For std::fixed and std::setprecision
  8.  
  9. #include <vector>
  10. #include <cmath>
  11.  
  12.  
  13. int main() {
  14.  
  15. //    double D = 1.55986;
  16.     double D = 0.0;
  17.     double dx = 0.000523438;
  18.     double L = 0.067;
  19.     int N_G = 129;
  20.    
  21.     double dt = 1.0e-2;
  22.     double alpha = dt*D/(dx*dx);
  23.     std::cout << "alpha stability: " << alpha << "\n";
  24.     std::cout << "max dt for the diffusion: " << 0.5*dx*dx/D << "\n";
  25.  
  26.     double T_diff = L*L/(M_PI*M_PI*D); std::cout << T_diff << "\n";
  27.  
  28.     int max_iter = 5000000;    
  29.  
  30.     double GAS_DENSITY = 9.64e+20;
  31.     double R = sqrt(1.6e-4/M_PI);
  32.     double nu_loss = D/(R*R);
  33.  
  34.  
  35.     double accuracy = 1.0E-4;
  36.     double cycle = 0;                   // to sace the moment of iteraion
  37.  
  38.     double S0, S1, S2;
  39.     double cylindric_diff1 = 0.0, cylindric_diff2 = 0.0;  
  40.  
  41.     std::vector<double> e1_density_fw(N_G, 0);
  42.     std::vector<double> e1_density_bw(N_G, 0);
  43.     std::vector<double> e2_density_fw(N_G, 0);
  44.     std::vector<double> e2_density_bw(N_G, 0);    
  45.     std::vector<double> gas_dens_local(N_G, GAS_DENSITY);
  46.  
  47.     std::ofstream file0("initial.dat");
  48.     file0 << std::scientific << std::setprecision(6);    
  49.  
  50.     std::ofstream file1("final.dat");
  51.     file1 << std::scientific << std::setprecision(6);      
  52.  
  53.  
  54.     std::ifstream infile("dat/rates.dat");  
  55.     std::vector<double> avg_rate1f; std::vector<double> avg_rate1b;
  56.     std::vector<double> avg_rate2f; std::vector<double> avg_rate2b;
  57.  
  58.  
  59.     double a, b, c, d;
  60.  
  61.     while (infile >> a >> b >> c >> d){
  62.         avg_rate1f.push_back(a); avg_rate1b.push_back(b);
  63.         avg_rate2f.push_back(c); avg_rate2b.push_back(d);
  64.     }
  65.  
  66.     // initial condition:
  67.     for (int j = 0;  j < N_G; j++){
  68.         e1_density_bw[j] = 0.0;
  69.         e2_density_bw[j] = 0.0;
  70.     }
  71.  
  72.  
  73.     for (int k = 0; k < max_iter; k++){
  74.          
  75.         e1_density_bw[0] = e1_density_bw[1]; e1_density_bw[N_G-1] = e1_density_bw[N_G-2];    
  76.         e2_density_bw[0] = e2_density_bw[1]; e2_density_bw[N_G-1] = e2_density_bw[N_G-2];  
  77.         gas_dens_local[0] = gas_dens_local[1]; gas_dens_local[N_G-1] = gas_dens_local[N_G-2];
  78.  
  79.         for (int j = 1; j < N_G-1; j++){                      
  80.  
  81.             S0 = dt*(-(avg_rate1f[j]+avg_rate2f[j])*gas_dens_local[j] + avg_rate1b[j]*e1_density_bw[j] + avg_rate2b[j]*e2_density_bw[j]);
  82.             S1 = dt*(avg_rate1f[j]*gas_dens_local[j] - avg_rate1b[j]*e1_density_bw[j]);
  83.             S2 = dt*(avg_rate2f[j]*gas_dens_local[j] - avg_rate2b[j]*e2_density_bw[j]);    
  84.  
  85.             e1_density_fw[j] = e1_density_bw[j]*(1.0 - 2*alpha) + alpha*(e1_density_bw[j+1] + e1_density_bw[j-1]) + S1;      // step of calculation    
  86.             e2_density_fw[j] = e2_density_bw[j]*(1.0 - 2*alpha) + alpha*(e2_density_bw[j+1] + e2_density_bw[j-1]) + S2;
  87.             gas_dens_local[j] += S0;
  88.        
  89.         }
  90.  
  91.         // boundary conditions:
  92.                                    
  93.         e1_density_fw[0] = e1_density_fw[1]; e1_density_fw[N_G-1] = e1_density_fw[N_G-2];
  94.                                      
  95.         e2_density_fw[0] = e2_density_fw[1]; e2_density_fw[N_G-1] = e2_density_fw[N_G-2];              
  96.  
  97.  
  98.         double err1 = 0.0, err2 = 0.0;
  99.         for(int j = 1; j < N_G-1; ++j) {
  100.             err1 = std::max(err1, std::abs(e1_density_fw[j] - e1_density_bw[j])/dt);
  101.             err2 = std::max(err2, std::abs(e2_density_fw[j] - e2_density_bw[j])/dt);
  102.         }
  103.  
  104.         if(err1 < accuracy && err2 < accuracy) {
  105.             cycle = k;
  106.             std::cout << "Converged in " << k << " iterations" << "\n";
  107.             break;
  108.         }        
  109.  
  110.         for (int j = 0; j < N_G; j++){
  111.             e1_density_bw[j] = e1_density_fw[j];                                                      // swap of arrays each step
  112.             e2_density_bw[j] = e2_density_fw[j];
  113.         }
  114.     }
  115.  
  116.     for (int j = 0;  j < N_G; j++){
  117.         double x = j*dx;
  118.         file1 << x << " " << e1_density_bw[j] << " " << e2_density_bw[j] << " " << gas_dens_local[j]<< "\n";
  119.     }    
  120.  
  121.     file0.close();
  122.     file1.close();
  123.  
  124.     return 0;
  125.  
  126. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!