GasFunctionWithDiffusionCond/main.cpp at master · Nevsky12/GasFunctionWithDiffusionCond · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
#include <cmath>
#include <cstdio>
#include <fstream>
#include <iostream>
const int N = 50, V = 11, T = 2000;
const double dv = 1.0, dx = 10.0, dt = 1.0;

double get_V(int id) { return dv * (id - V / 2) - dv / 2; }

void set_initial_data(double other[N][V]) {

    for (int i = 0; i < N; i++) {
        for (int j = 0; j < V; j++) {
            auto v = get_V(j);
            auto x = (i - N / 2);
            other[i][j] = std::exp(-v * v / 2) * std::exp( -x * x / 2);
        }
    }
}

double find_const(double other[N][V]) {

    auto sum_positive=0;
    auto sum_negative=0;
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < V; j++) {
            auto v = get_V(j);
            auto x = (i - N / 2);
            if (v > 0) { sum_positive += v * std::exp(-v * v / 2); }
            else { sum_negative += other[i][j] * v; }
        }
    }
    if (sum_positive != 0) { return (double)-sum_negative / sum_positive; }
    else return 0;
}

double find_diffusion_condition(double v, double Const) { return Const * std::exp(- v * v / 2); }

int main() {

    char filename[100];
    double ar_current[N][V]{0};
    double ar_next[N][V]{0};
    double ar_concentration[N]{0};
    set_initial_data(ar_current);
    std::fstream file_out;
    for (int k = 0; k < T; k++) {
        for (int i = 0; i < N; i++) {

            auto normical_const = find_const(ar_current);
            for (int j = 0; j < V; j++) {
                auto v = get_V(j);
                double gamma = v * dt / dx;

                if (v > 0) {
                    if (i == 0) {
                        ar_next[0][j] = find_diffusion_condition(v,normical_const);
                    } else {
                        ar_next[i][j] = ar_current[i][j] - gamma * (ar_current[i][j] - ar_current[i - 1][j]);
                    }
                } else {
                    if (i == N - 1) {
                        ar_next[N - 1][j] = ar_current[0][j];
                    } else {
                        ar_next[i][j] = ar_current[i][j] - gamma * (ar_current[i + 1][j] - ar_current[i][j]);
                    }
                }
            }
        }
        std::swap(ar_current, ar_next);

        sprintf(filename, "data/out_%03d.txt", k);
        file_out.open(filename, std::fstream::out);
        for (int i = 0; i < N; i++) {
            ar_concentration[i] = 0;
            for (int j = 0; j < V; j++) {
                ar_concentration[i] += ar_current[i][j];
            }
            file_out << i << '\t' << ar_concentration[i] << std::endl;
        }
        file_out.close();
    }
}