SciHPC/benchmarks_lsf_initialization.cpp at master · luisleo526/SciHPC · 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
84
85
86
87
88
89
90
91
#include <iostream>
#include <cmath>
#include <fstream>
#include <chrono>
#include <random>
#include <ctime>

#include "scihpc/wrapper.h"
#include "scihpc/vtkWriter.h"
#include "scihpc/structured_grid.h"
#include "scihpc/source.h"
#include "scihpc/runge_kutta.h"
#include "scihpc/flux.h"
#include "scihpc/wrapper_func.h"

int main() {

    auto geo = structured_grid(axis{-1.1, 1.1, 256},
                               axis{-1.1, 1.1, 256});

    auto phi = wrapper(true, &geo,
                       bc_info{NEUMANN}, bc_info{NEUMANN},
                       bc_info{NEUMANN}, bc_info{NEUMANN});
    auto vel = wrapper(false, &geo,
                       bc_info{NEUMANN}, bc_info{NEUMANN},
                       bc_info{NEUMANN}, bc_info{NEUMANN});
    auto solver = runge_kutta(phi.scalar->Nx, phi.scalar->Ny, phi.scalar->Nz);
    auto vtk = vtkWriter(&geo, "lsf_iris");
    auto param = new problem_parameters{};
    auto deri_solvers = SharedSolvers_alloc(phi.scalar, &geo);
    auto dummy = dummy_data_alloc(phi.scalar);
    phi.link_params(param);
    phi.link_solvers(deri_solvers);
    phi.link_dummy(dummy);

    param->dt = 0.01 * geo.h;
    param->rdt = 0.25 * geo.h;
    param->ls_width = 3.0 * geo.h;

    std::default_random_engine generator(time(NULL));
    // noise
    std::uniform_real_distribution<DataType> unif(-0.00, 0.00);
    std::uniform_real_distribution<DataType> unif2(0, 0.0);

    // Initialize phi
    for (int i = 0; i < phi.scalar->nx; ++i) {
        for (int j = 0; j < phi.scalar->ny; ++j) {
            auto index = phi.scalar->index_mapping(i + 1, j + 1, 1);
            if (geo.xc[i] > -1.0 && geo.xc[i] < 1.0) {
                if (geo.yc[j] > -0.85 * cos(0.5 * pi * geo.xc[i]) - unif(generator) and
                    geo.yc[j] < 0.85 * cos(0.5 * pi * geo.xc[i]) + unif(generator)) {
                    if (pow(geo.xc[i], 2) + pow(geo.yc[j], 2) < 0.25 * (1.0 + unif(generator))) {
                        phi.scalar->data[index.i][index.j][index.k] = 255.0 - unif2(generator);
                    } else {
                        phi.scalar->data[index.i][index.j][index.k] = -255.0 + unif2(generator);
                    }
                } else {
                    phi.scalar->data[index.i][index.j][index.k] = 255.0 - unif2(generator);
                }
            } else {
                phi.scalar->data[index.i][index.j][index.k] = 255.0 + unif2(generator);
            }

            phi.scalar->data[index.i][index.j][index.k] /= 255.0;
        }
    }
    phi.apply_scalar_bc();

    vtk.create(0);
    vtk.add_scalar(phi.scalar, "phi");
    vtk.close();

    find_sign(&phi);
    stabilized_upon_gradient(&phi);
    int step = 0;
    do {
        store_tmp(&phi);
        if (step % 100 == 0 and step < 500) {
            find_sign(&phi);
            stabilized_upon_gradient(&phi);
        }
        solver.tvd_rk3(&phi, &vel, &identity_flux, &lsf_redistance_lambda);
        std::cout << l2norm(&phi) << std::endl;
    } while (++step < 1500 and l2norm(&phi) > 1e-6);

    vtk.create(1);
    vtk.add_scalar(phi.scalar, "phi");
    vtk.close();

    return 0;
}