TurbulentFlowSimulation-HPC/TurbulentSimulation.cpp at master · ludcila/TurbulentFlowSimulation-HPC · 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
92
93
94
95
96
97
98
99
100
101
#include "TurbulentSimulation.h"
#include "stencils/WallDistanceStencil.h"

TurbulentSimulation::TurbulentSimulation(Parameters &parameters, TurbulentFlowField &flowField):
	Simulation(parameters, flowField),
	_turbulentFlowField(flowField),
	_turbulentFghStencil(parameters),
	_turbulentFghIterator(_turbulentFlowField, parameters, _turbulentFghStencil),
	_turbulentViscosityStencil(parameters),
	_turbulentViscosityIterator(_turbulentFlowField, parameters, _turbulentViscosityStencil, 1, 0),
	_turbulentVtkStencil(parameters),
	_turbulentVtkIterator(_turbulentFlowField, parameters, _turbulentVtkStencil, 1, 0),
	_maxNuStencil(parameters),
    _maxNuFieldIterator(_turbulentFlowField,parameters,_maxNuStencil),
    _maxNuBoundaryIterator(_turbulentFlowField,parameters,_maxNuStencil),
	_turbulentViscosityBoundaryStencil(parameters),
	_turbulentViscosityBoundaryIterator(_turbulentFlowField, parameters, _turbulentViscosityBoundaryStencil, 1, 0),
	_parallelManagerTurbulent(_turbulentFlowField, parameters)
{
}

// TODO: Changes for turbulent simulation not implemented yet!
void TurbulentSimulation::solveTimestep(){

	// determine and set max. timestep which is allowed in this simulation
	setTimeStep();
	// compute fgh
	_turbulentFghIterator.iterate();
	// set global boundary values
	_wallFGHIterator.iterate();
	// compute the right hand side
	_rhsIterator.iterate();
	// solve for pressure
	_solver.solve();
	// TODO WS2: communicate pressure values
	_parallelManagerTurbulent.communicatePressure();
	// compute velocity
	_velocityIterator.iterate();
	// set obstacle boundaries
	_obstacleIterator.iterate();
	// TODO WS2: communicate velocity values
	_parallelManagerTurbulent.communicateVelocity();
	// Iterate for velocities on the boundary
	_wallVelocityIterator.iterate();

	_parallelManagerTurbulent.communicateCenterLineVelocity();
	_turbulentViscosityIterator.iterate();
	_parallelManagerTurbulent.communicateViscosity();
	_turbulentViscosityBoundaryIterator.iterate();
}

// TODO: Changes for turbulent simulation not implemented yet!
void TurbulentSimulation::setTimeStep(){

	const FLOAT cinematicviscosity=1.0/_parameters.flow.Re;
	FLOAT localMin, globalMin;
	assertion(_parameters.geometry.dim == 2 || _parameters.geometry.dim == 3);

	// determine maximum velocity
	_maxUStencil.reset();
	_maxUFieldIterator.iterate();
	_maxUBoundaryIterator.iterate();

	_maxNuStencil.reset();
	_maxNuFieldIterator.iterate();
	_maxNuBoundaryIterator.iterate();

	if (_parameters.geometry.dim == 3) {
		_parameters.timestep.dt = 1.0 / _maxUStencil.getMaxValues()[2];
	} else {
		_parameters.timestep.dt = 1.0 / _maxUStencil.getMaxValues()[0];
	}

	localMin = std::min(_parameters.timestep.dt,
		                            std::min(std::min(1.0/_maxNuStencil.getMaxValue(),
		                            1.0 / _maxUStencil.getMaxValues()[0]),
		                            1.0 / _maxUStencil.getMaxValues()[1]));

	// Here, we select the type of operation before compiling. This allows to use the correct
	// data type for MPI. Not a concern for small simulations, but useful if using heterogeneous
	// machines.

	globalMin = MY_FLOAT_MAX;
	MPI_Allreduce(&localMin, &globalMin, 1, MY_MPI_FLOAT, MPI_MIN, PETSC_COMM_WORLD);

	_parameters.timestep.dt = globalMin;
	_parameters.timestep.dt *= _parameters.timestep.tau;

}

void TurbulentSimulation::plotVTK(int timeStep, std::string foldername) {
	_turbulentVtkIterator.iterate();
	_turbulentVtkStencil.write(this->_turbulentFlowField, timeStep, foldername);
}

void TurbulentSimulation::initializeFlowField() {
	Simulation::initializeFlowField();
	WallDistanceStencil wds(_parameters);
	FieldIterator<TurbulentFlowField> it(_turbulentFlowField, _parameters, wds);
	it.iterate();
};