Peano 4
Loading...
Searching...
No Matches
ContextDiffuse.h
Go to the documentation of this file.
1#ifndef EXASEIS_CONTEXT_DIFFUSE_HEADER
2#define EXASEIS_CONTEXT_DIFFUSE_HEADER
3
4#include "Context.h"
8#include "../Topography/TopographyFactory.h"
9#include "../../../ExaHyPE/kernels/GaussLegendreBasis.h"
10#include "../../../../ExaHyPE/kernels/KernelUtils.h"
11
12template <class Shortcuts, int basisSize>
13class ContextDiffuse: public Context<Shortcuts, basisSize> {
14
15public:
17 std::string& scenario_string,
18 std::string& topography_string,
19 DomainInformation* a_domainInfo,
20 SolverInformation* a_solverInfo,
21 double a_interface_factor,
22 double a_max_distance_to_direct_surface
23 ):
24 Context<Shortcuts, basisSize>(scenario_string, a_domainInfo) {
25
26 domainInfo = a_domainInfo;
27 solverInfo = a_solverInfo;
28 Topography* topography = TopographyFactory::createTopography(topography_string, a_domainInfo);
29 interface = new DiffuseInterface<Shortcuts, basisSize>(
30 topography, a_domainInfo, a_solverInfo, a_interface_factor, a_max_distance_to_direct_surface
31 );
32 };
34 delete interface;
35 delete domainInfo;
36 delete solverInfo;
37 };
38
39 virtual void initUnknownsPatch(
40 double* luh,
43 double t,
44 double dt
45 ) override {
46 Shortcuts s;
47 int numberOfData = s.SizeVariables + s.SizeParameters;
48
49 kernels::idx3 id_xyz(basisSize, basisSize, basisSize);
50 kernels::idx4 id_xyzn(basisSize, basisSize, basisSize, Dimensions);
51 kernels::idx4 id_xyzf(basisSize, basisSize, basisSize, numberOfData);
52
53 double offset_x = center[0] - 0.5 * dx[0];
54 double offset_y = center[1] - 0.5 * dx[1];
55 double offset_z = center[2] - 0.5 * dx[2];
56
57 double width_x = dx[0];
58 double width_y = dx[1];
59 double width_z = dx[2];
60
61 double nodes[basisSize * basisSize * basisSize * Dimensions];
62 double alpha[basisSize * basisSize * basisSize];
63 double quad_nodes[basisSize];
64
65 if (solverInfo->isDG()) {
66 std::copy_n(kernels::legendre::nodes[basisSize - 1], basisSize, quad_nodes);
67 } else {
68 for (int i = 1; i < basisSize + 1; i++) {
69 quad_nodes[i - 1] = i / static_cast<double>(basisSize + 1);
70 }
71 }
72
73 for (int k = 0; k < basisSize; k++) {
74 double z = (offset_z + width_z * quad_nodes[k]);
75 for (int j = 0; j < basisSize; j++) {
76 double y = (offset_y + width_y * quad_nodes[j]);
77 for (int i = 0; i < basisSize; i++) {
78 double x = (offset_x + width_x * quad_nodes[i]);
79 nodes[id_xyzn(k, j, i, 0)] = x;
80 nodes[id_xyzn(k, j, i, 1)] = y;
81 nodes[id_xyzn(k, j, i, 2)] = z;
82 }
83 }
84 }
85
86 interface->getAlphaPatch(nodes, alpha);
87
88 for (int k = 0; k < basisSize; k++) {
89 for (int j = 0; j < basisSize; j++) {
90 for (int i = 0; i < basisSize; i++) {
91 luh[id_xyzf(k, j, i, s.alpha)] = alpha[id_xyz(k, j, i)];
92 this->scenario->initUnknownsPointwise(nodes + id_xyzn(k, j, i, 0), center, t, dt, luh + id_xyzf(k, j, i, 0));
93 assertion(luh[id_xyzf(k, j, i, s.rho)] > 0.0);
94 }
95 }
96 }
97#ifdef Asserts
98 for (int k = 0; k < basisSize; k++) {
99 for (int j = 0; j < basisSize; j++) {
100 for (int i = 0; i < basisSize; i++) {
101 for (int m = 0; m < numberOfData; m++) {
102 assertion5(std::isfinite(luh[id_xyzf(k, j, i, m)]), k, j, i, m, domainInfo->meshLevel);
103 }
104 assertion(luh[id_xyzf(k, j, i, s.rho)] > 0.0);
105 }
106 }
107 }
108#endif
109 };
110
111 void initUnknownsPointwise(const double* const x, const double t, const double dt, double* Q) {
112 Shortcuts s;
113 const int numberOfData = s.SizeVariables + s.SizeParameters;
114 Q[s.alpha] = interface->getAlpha(x);
116 center[0] = x[0];
117 center[1] = x[1];
118 center[2] = x[2];
119
120 this->scenario->initUnknownsPointwise(x, center, t, dt, Q);
121#ifdef Asserts
122 for (int m = 0; m < numberOfData; m++) {
123 assertion2(std::isfinite(Q[m]), m, domainInfo->meshLevel);
124 }
125#endif
126 };
127
128private:
132};
133
134#endif
#define assertion2(expr, param0, param1)
#define assertion(expr)
#define assertion5(expr, param0, param1, param2, param3, param4)
And from this we can write down f$ nabla phi_i nabla phi_i dx but since we are constructing matrix let s investigate the f$ j
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ k
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ so we are integrating with f$ phi_k phi_k dx
virtual void initUnknownsPatch(double *luh, const tarch::la::Vector< Dimensions, double > &center, const tarch::la::Vector< Dimensions, double > &dx, double t, double dt) override
DomainInformation * domainInfo
SolverInformation * solverInfo
DiffuseInterface< Shortcuts, basisSize > * interface
ContextDiffuse(std::string &scenario_string, std::string &topography_string, DomainInformation *a_domainInfo, SolverInformation *a_solverInfo, double a_interface_factor, double a_max_distance_to_direct_surface)
void initUnknownsPointwise(const double *const x, const double t, const double dt, double *Q)
virtual bool isDG()=0
examples::exahype2::elastic::VariableShortcuts s
Definition loh.cpp:10
Simple vector class.
Definition Vector.h:134