Peano 4
Loading...
Searching...
No Matches
Poisson.cpp
Go to the documentation of this file.
2
5
6
8
10
12 _currentMaxGradientDelta(1.0),
13 _previousMaxGradientDelta(1.0) {}
14
15
17 double globalMinTimeStamp, double globalMaxTimeStamp, double globalMinTimeStepSize, double globalMaxTimeStepSize
18) {
19 AbstractPoisson::startTimeStep(globalMinTimeStamp, globalMaxTimeStamp, globalMinTimeStepSize, globalMaxTimeStepSize);
20 _previousMaxGradientDelta = _currentMaxGradientDelta;
21 _currentMaxGradientDelta = 0.0;
22 if (_previousMaxGradientDelta > 1e-8) {
23 logInfo(
24 "startTimeStep(...)",
25 "max gradient difference="
26 << _previousMaxGradientDelta
27 << ". Solver has not converged, i.e. do not use |dt|_max=" << AbstractPoisson::getAdmissibleTimeStepSize()
28 );
29 } else {
30 logInfo(
31 "startTimeStep(...)",
32 "max gradient difference="
33 << _previousMaxGradientDelta
34 << ". Solver has converged. Use |dt|_max=" << AbstractPoisson::getAdmissibleTimeStepSize()
35 );
36 }
37}
38
39
41 double originalAdmissibleTimeStepSize = AbstractPoisson::getAdmissibleTimeStepSize();
42
43 if (_previousMaxGradientDelta > 1.0e-8) {
44 return tarch::la::NUMERICAL_ZERO_DIFFERENCE * 2.0; // not converged yet
45 } else {
46 return originalAdmissibleTimeStepSize;
47 }
48}
49
50
52 tarch::multicore::Lock lock(_maxGradientDeltaSemaphore);
53 _currentMaxGradientDelta = std::max(_currentMaxGradientDelta, delta);
54}
55
56
58 double* __restrict__ Q,
59 const tarch::la::Vector<Dimensions, double>& volumeCentre,
61 bool gridIsConstructed
62) {
63 logTraceInWith3Arguments("initialCondition(...)", volumeCentre, volumeH, gridIsConstructed);
64
65 // Trivial initial conditions
66 Q[0] = 0.0;
67 Q[1] = 0.0;
68 Q[2] = 0.0;
69 Q[3] = 0.0;
70 Q[4] = 1.0; // we need a proper initial density for the very first time step
71 // just before the actual Poisson solve has already started.
72
73 logTraceOut("initialCondition(...)");
74}
75
76
78 const double* __restrict__ Qinside, // Qinside[4+0]
79 double* __restrict__ Qoutside, // Qoutside[4+0]
82 double t,
83 int normal
84) {
85 logTraceInWith4Arguments("boundaryConditions(...)", faceCentre, volumeH, t, normal);
86
87 // Neumann conditions
88 Qoutside[0] = Qinside[0];
89 Qoutside[1] = Qinside[1];
90 Qoutside[2] = Qinside[2];
91 Qoutside[3] = Qinside[3];
92
93 logTraceOut("boundaryConditions(...)");
94}
95
96
98 const double* __restrict__ Q, // Q[4+0],
101 double t,
102 double dt,
103 int normal
104) {
105 logTraceInWith4Arguments("maxEigenvalue(...)", faceCentre, volumeH, t, normal);
106
107 // @todo needs cross-checking
108 const double T_tau = 1.0 / 4.0 / tarch::la::PI / tarch::la::PI;
109
110 double result = std::sqrt(1 / T_tau);
111
112 logTraceOutWith1Argument("maxEigenvalue(...)", result);
113 return result;
114}
115
116
118 const double* __restrict__ Q, // Q[4+0],
121 double t,
122 double dt,
123 int normal,
124 double* __restrict__ F // F[4]
125) {
126 logTraceInWith4Arguments("flux(...)", faceCentre, volumeH, t, normal);
127
128 const double T_tau = 1.0 / 4.0 / tarch::la::PI / tarch::la::PI;
129 F[0] = Q[normal + 1];
130 F[1] = normal == 0 ? -Q[0] / T_tau : 0.0;
131 F[2] = normal == 1 ? -Q[0] / T_tau : 0.0;
132 F[3] = normal == 2 ? -Q[0] / T_tau : 0.0;
133
134 logTraceOut("flux(...)");
135}
136
137
139 const double* __restrict__ Q, // Q[4+1]
142 double t,
143 double dt,
144 double* __restrict__ S // S[4
145) {
146 logTraceInWith4Arguments("sourceTerm(...)", volumeX, volumeH, t, dt);
147
148 const double T_tau = 1.0 / 4.0 / tarch::la::PI / tarch::la::PI;
149 // const double G = 6.67e-11;
150 const double G = 6.67e-3; // @todo No clue what normalised constant belong in here. 1e-11 is definitely too small
151 const double rho = Q[4]; // This is the input term
152
153 S[0] = -tarch::la::PI * G / rho;
154 S[1] = -Q[1] / T_tau;
155 S[2] = -Q[2] / T_tau;
156 S[3] = -Q[3] / T_tau;
157
158 logTraceOut("sourceTerm(...)");
159}
#define logTraceOutWith1Argument(methodName, argument0)
Definition Log.h:380
#define logTraceOut(methodName)
Definition Log.h:379
#define logTraceInWith4Arguments(methodName, argument0, argument1, argument2, argument3)
Definition Log.h:373
#define logTraceInWith3Arguments(methodName, argument0, argument1, argument2)
Definition Log.h:372
#define logInfo(methodName, logMacroMessageStream)
Wrapper macro around tarch::tarch::logging::Log to improve logging.
Definition Log.h:411
void startTimeStep(double globalMinTimeStamp, double globalMaxTimeStamp, double globalMinTimeStepSize, double globalMaxTimeStepSize) override
Definition Poisson.cpp:16
void sourceTerm(const double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &volumeCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t, double dt, double *__restrict__ S) override
Definition Poisson.cpp:138
virtual double maxEigenvalue(const double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &faceCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t, double dt, int normal) override
Determine max eigenvalue over Jacobian in a given point with solution values (states) Q.
Definition Poisson.cpp:97
double getAdmissibleTimeStepSize() const override
overwritten
Definition Poisson.cpp:40
static tarch::multicore::BooleanSemaphore _maxGradientDeltaSemaphore
Definition Poisson.h:15
void initialCondition(double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &volumeCentre, const tarch::la::Vector< Dimensions, double > &volumeH, bool gridIsConstructed) override
Definition Poisson.cpp:57
virtual void flux(const double *__restrict__ Q, const tarch::la::Vector< Dimensions, double > &faceCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t, double dt, int normal, double *__restrict__ F) override
Definition Poisson.cpp:117
virtual void boundaryConditions(const double *__restrict__ Qinside, double *__restrict__ Qoutside, const tarch::la::Vector< Dimensions, double > &faceCentre, const tarch::la::Vector< Dimensions, double > &volumeH, double t, int normal) override
Definition Poisson.cpp:77
Log Device.
Definition Log.h:516
Create a lock around a boolean semaphore region.
Definition Lock.h:19
constexpr double NUMERICAL_ZERO_DIFFERENCE
Definition Scalar.h:17
constexpr double PI
Definition Scalar.h:12
Simple vector class.
Definition Vector.h:134