Peano
Loading...
Searching...
No Matches
CCZ4Helper.py
Go to the documentation of this file.
1def get_body_of_enforceCCZ4constraint():
2 return """
3 {
4 #if Dimensions==2
5 constexpr int itmax = {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}} * {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}};
6 #endif
7
8 #if Dimensions==3
9 constexpr int itmax = {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}} * {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}} * {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}};
10 #endif
11
12 int index = 0;
13 for (int i=0;i<itmax;i++)
14 {
15 applications::exahype2::ccz4::enforceCCZ4constraints( newQ+index );
16 index += {{NUMBER_OF_UNKNOWNS}} + {{NUMBER_OF_AUXILIARY_VARIABLES}};
17 }
18 }
19"""
20
21def get_body_of_adm_constraints(patch_size, number_of_output_variable, so_flag):
22 if so_flag:
23 extra_derivative_assign="""
24 ::exahype2::enumerator::AoSLexicographicEnumerator enumeratorWithAuxiliaryVariablesOnReconstructedPatch( 1, {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}}, {{HALO_SIZE}}, {{NUMBER_OF_UNKNOWNS}}, {{NUMBER_OF_AUXILIARY_VARIABLES}});
25 ::exahype2::enumerator::AoSLexicographicEnumerator enumeratorWithoutAuxiliaryVariables( {{RK_STEPS}}, {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}}, 0, {{NUMBER_OF_UNKNOWNS}}, 0 );
26 int primaryVarsIndices[25]={ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,
27 10,11,12,13,14,15,16,17,18,19,
28 20,21,22,53,54};
29 int auxiliaryVarsIndices[34]={23,24,25,26,27,28,29,30,31,32,
30 33,34,35,36,37,38,39,40,41,42,
31 43,44,45,46,47,48,49,50,51,52,
32 55,56,57,58};
33
34 dfor( dof, 3 ) {
35 for (int unknown : auxiliaryVarsIndices) {
36 oldQWithHalo[ enumeratorWithAuxiliaryVariablesOnReconstructedPatch(0,dof,unknown) ] = //abuse of RHS, auxiliary variables are stored there.
37 fineGridCellCCZ4QRhsEstimates.value[ enumeratorWithoutAuxiliaryVariables(0,dof,unknown) ];
38 }
39 }
40 """
41 else:
42 extra_derivative_assign=""
43 return """
44 const int patchSize = """ + str( patch_size ) + """;
45 double volumeH = ::exahype2::fv::getVolumeLength(marker.h(),patchSize);
46
47 const int n_a_v="""+str(number_of_output_variable)+""";
48 const int overlap=3; //make sure you are using fd4 solver!
49
50 if (not marker.willBeRefined() and repositories::instanceOfCCZ4.getSolverState()!=CCZ4::SolverState::GridConstruction and repositories::instanceOfCCZ4.getSolverState()==CCZ4::SolverState::RungeKuttaSubStep0) {"""+extra_derivative_assign+"""
51
52 dfor(cell,patchSize) {
53 tarch::la::Vector<Dimensions,int> currentCell = cell + tarch::la::Vector<Dimensions,int>(overlap);
54
55 double gradQ[3*59]={ 0 };
56
57 /*for (int d=0; d<3; d++) {
58 tarch::la::Vector<Dimensions,int> leftCell = currentCell;
59 tarch::la::Vector<Dimensions,int> rightCell = currentCell;
60 leftCell(d) -= 1;
61 rightCell(d) += 1;
62 const int leftCellSerialised = peano4::utils::dLinearised(leftCell, patchSize + 2*overlap);
63 const int rightCellSerialised = peano4::utils::dLinearised(rightCell,patchSize + 2*overlap);
64 for(int i=0; i<59; i++) {
65 gradQ[d*59+i] = ( oldQWithHalo[rightCellSerialised*(59+n_a_v)+i] - oldQWithHalo[leftCellSerialised*(59+n_a_v)+i] ) / 2.0 / volumeH;
66 }
67 }*/
68
69 for (int d=0; d<3; d++) {
70 tarch::la::Vector<Dimensions,int> leftCell = currentCell;
71 tarch::la::Vector<Dimensions,int> rightCell = currentCell;
72 tarch::la::Vector<Dimensions,int> DouleftCell = currentCell;
73 tarch::la::Vector<Dimensions,int> DourightCell = currentCell;
74 leftCell(d) -= 1; DouleftCell(d) -= 2;
75 rightCell(d) += 1; DourightCell(d) += 2;
76 const int leftCellSerialised = peano4::utils::dLinearised(leftCell, patchSize + 2*overlap);
77 const int rightCellSerialised = peano4::utils::dLinearised(rightCell,patchSize + 2*overlap);
78 const int DouleftCellSerialised = peano4::utils::dLinearised(DouleftCell, patchSize + 2*overlap);
79 const int DourightCellSerialised = peano4::utils::dLinearised(DourightCell,patchSize + 2*overlap);
80 for(int i=0; i<59; i++) {
81 gradQ[d*59+i] = ( -1*oldQWithHalo[DourightCellSerialised*(59+n_a_v)+i] + 8*oldQWithHalo[rightCellSerialised*(59+n_a_v)+i] - 8*oldQWithHalo[leftCellSerialised*(59+n_a_v)+i] + 1*oldQWithHalo[DouleftCellSerialised*(59+n_a_v)+i]) / 12.0 / volumeH;
82 }
83 }
84
85 const int cellSerialised = peano4::utils::dLinearised(currentCell, patchSize + 2*overlap);
86
87 double constraints[n_a_v]={0};
88 admconstraints(constraints, oldQWithHalo+cellSerialised*(59+n_a_v), gradQ);
89 ::exahype2::enumerator::AoSLexicographicEnumerator enumeratorWithAuxiliaryVariables ( 1, patchSize, 0, 59, n_a_v);
90 for (int i=0;i<n_a_v;i++){
91 fineGridCellCCZ4Q.value[ enumeratorWithAuxiliaryVariables(0,cell,59+i) ] = constraints[i];
92 }
93 }
94 }
95"""
96
97def get_body_of_Psi_Calc(patch_size, number_of_output_variable, so_flag):
98 if so_flag:
99 extra_derivative_assign="""
100 ::exahype2::enumerator::AoSLexicographicEnumerator enumeratorWithAuxiliaryVariablesOnReconstructedPatch( 1, {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}}, {{HALO_SIZE}}, {{NUMBER_OF_UNKNOWNS}}, {{NUMBER_OF_AUXILIARY_VARIABLES}});
101 ::exahype2::enumerator::AoSLexicographicEnumerator enumeratorWithoutAuxiliaryVariables( {{RK_STEPS}}, {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}}, 0, {{NUMBER_OF_UNKNOWNS}}, 0 );
102 int primaryVarsIndices[25]={ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,
103 10,11,12,13,14,15,16,17,18,19,
104 20,21,22,53,54};
105 int auxiliaryVarsIndices[34]={23,24,25,26,27,28,29,30,31,32,
106 33,34,35,36,37,38,39,40,41,42,
107 43,44,45,46,47,48,49,50,51,52,
108 55,56,57,58};
109
110 dfor( dof, 3 ) {
111 for (int unknown : auxiliaryVarsIndices) {
112 oldQWithHalo[ enumeratorWithAuxiliaryVariablesOnReconstructedPatch(0,dof,unknown) ] = //abuse of RHS, auxiliary variables are stored there.
113 fineGridCellCCZ4QRhsEstimates.value[ enumeratorWithoutAuxiliaryVariables(0,dof,unknown) ];
114 }
115 }
116 """
117 else:
118 extra_derivative_assign=""
119 return """
120 const int patchSize = """ + str( patch_size ) + """;
121 double volumeH = ::exahype2::fv::getVolumeLength(marker.h(),patchSize);
122
123 const int n_a_v="""+str(number_of_output_variable)+""";
124 const int overlap=3; //make sure you are using fd4 solver!
125
126 if (not marker.willBeRefined() and repositories::instanceOfCCZ4.getSolverState()!=CCZ4::SolverState::GridConstruction and repositories::instanceOfCCZ4.getSolverState()==CCZ4::SolverState::RungeKuttaSubStep0) {"""+extra_derivative_assign+"""
127
128 dfor(cell,patchSize) {
129 tarch::la::Vector<Dimensions,int> currentCell = cell + tarch::la::Vector<Dimensions,int>(overlap);
130
131 double gradQ[3*59]={ 0 };
132
133 /*for (int d=0; d<3; d++) {
134 tarch::la::Vector<Dimensions,int> leftCell = currentCell;
135 tarch::la::Vector<Dimensions,int> rightCell = currentCell;
136 leftCell(d) -= 1;
137 rightCell(d) += 1;
138 const int leftCellSerialised = peano4::utils::dLinearised(leftCell, patchSize + 2*overlap);
139 const int rightCellSerialised = peano4::utils::dLinearised(rightCell,patchSize + 2*overlap);
140 for(int i=0; i<59; i++) {
141 gradQ[d*59+i] = ( oldQWithHalo[rightCellSerialised*(59+n_a_v)+i] - oldQWithHalo[leftCellSerialised*(59+n_a_v)+i] ) / 2.0 / volumeH;
142 }
143 }*/
144
145 for (int d=0; d<3; d++) {
146 tarch::la::Vector<Dimensions,int> leftCell = currentCell;
147 tarch::la::Vector<Dimensions,int> rightCell = currentCell;
148 tarch::la::Vector<Dimensions,int> DouleftCell = currentCell;
149 tarch::la::Vector<Dimensions,int> DourightCell = currentCell;
150 leftCell(d) -= 1; DouleftCell(d) -= 2;
151 rightCell(d) += 1; DourightCell(d) += 2;
152 const int leftCellSerialised = peano4::utils::dLinearised(leftCell, patchSize + 2*overlap);
153 const int rightCellSerialised = peano4::utils::dLinearised(rightCell,patchSize + 2*overlap);
154 const int DouleftCellSerialised = peano4::utils::dLinearised(DouleftCell, patchSize + 2*overlap);
155 const int DourightCellSerialised = peano4::utils::dLinearised(DourightCell,patchSize + 2*overlap);
156 for(int i=0; i<59; i++) {
157 gradQ[d*59+i] = ( -1*oldQWithHalo[DourightCellSerialised*(59+n_a_v)+i] + 8*oldQWithHalo[rightCellSerialised*(59+n_a_v)+i] - 8*oldQWithHalo[leftCellSerialised*(59+n_a_v)+i] + 1*oldQWithHalo[DouleftCellSerialised*(59+n_a_v)+i]) / 12.0 / volumeH;
158 }
159 }
160
161 const int cellSerialised = peano4::utils::dLinearised(currentCell, patchSize + 2*overlap);
162
163 double Psi4[n_a_v]={0};
164 double currentPosition[3];
165 for (int d=0; d<3; d++) currentPosition[d]=marker.getOffset()(d)+(cell(d)+0.5)*volumeH;
166 Psi4Calc(Psi4, oldQWithHalo+cellSerialised*(59+n_a_v), gradQ, currentPosition);
167 ::exahype2::enumerator::AoSLexicographicEnumerator enumeratorWithAuxiliaryVariables ( 1, patchSize, 0, 59, n_a_v);
168 fineGridCellCCZ4Q.value[ enumeratorWithAuxiliaryVariables(0,cell,59+0) ] = Psi4[0];
169 fineGridCellCCZ4Q.value[ enumeratorWithAuxiliaryVariables(0,cell,59+1) ] = Psi4[1];
170 }
171 }
172"""
173
174def get_body_of_SommerfeldCondition(scenario, unknowns, auxiliary_variables, restart_time=-1):
175 if restart_time>0:
176 timeStamp="CheckpointTimeStamp"
177 else:
178 timeStamp="-1.0"
179
180 if scenario=="single-puncture":
181 return """
182 ::exahype2::fd::applySommerfeldConditions(
183 [&](
184 const double * __restrict__ Q,
185 const tarch::la::Vector<Dimensions,double>& faceCentre,
186 const tarch::la::Vector<Dimensions,double>& gridCellH,
187 double t,
188 double dt,
189 int normal
190 ) -> double {
191 return repositories::{{SOLVER_INSTANCE}}.maxEigenvalue( Q, faceCentre, gridCellH, t, dt, normal );
192 },
193 [&](
194 double * __restrict__ Q,
195 const tarch::la::Vector<Dimensions,double>& faceCentre,
196 const tarch::la::Vector<Dimensions,double>& gridCellH
197 ) -> void {
198 for (int i=0; i<"""+str(unknowns + auxiliary_variables)+"""; i++) {
199 Q[i] = 0.0;
200 }
201 Q[0] = 1.0; Q[3] = 1.0; Q[5] = 1.0;
202 //const double r=tarch::la::norm2(faceCentre);
203 //Q[16] = 0.5*(1+(1-1.0/2/r)/(1+1.0/2/r));
204 //Q[54] = 1/(1+1.0/2/r)/(1+1.0/2/r);
205 Q[16] = 1.0; Q[54] = 1.0;
206 },
207 marker.x(),
208 marker.h(),
209 {{FACE_METADATA_ACCESSOR}}.getOldTimeStamp(marker.getSelectedFaceNumber()<Dimensions ? 1 : 0),
210 repositories::{{SOLVER_INSTANCE}}.getMinTimeStepSize(),
211 {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}},
212 {{OVERLAP}},
213 {{NUMBER_OF_UNKNOWNS}},
214 {{NUMBER_OF_AUXILIARY_VARIABLES}},
215 marker.getSelectedFaceNumber(),
216 {0.0, 0.0, 0.0},
217 fineGridFace{{UNKNOWN_IDENTIFIER}}Old.value,
218 fineGridFace{{UNKNOWN_IDENTIFIER}}New.value,
219 """+timeStamp+"""
220 );
221"""
222 else:
223 return """
224 ::exahype2::fd::applySommerfeldConditions(
225 [&](
226 const double * __restrict__ Q,
227 const tarch::la::Vector<Dimensions,double>& faceCentre,
228 const tarch::la::Vector<Dimensions,double>& gridCellH,
229 double t,
230 double dt,
231 int normal
232 ) -> double {
233 return repositories::{{SOLVER_INSTANCE}}.maxEigenvalue( Q, faceCentre, gridCellH, t, dt, normal );
234 },
235 [&](
236 double * __restrict__ Q,
237 const tarch::la::Vector<Dimensions,double>& faceCentre,
238 const tarch::la::Vector<Dimensions,double>& gridCellH
239 ) -> void {
240 for (int i=0; i<"""+str(unknowns + auxiliary_variables)+"""; i++) {
241 Q[i] = 0.0;
242 }
243 Q[0] = 1.0; Q[3] = 1.0; Q[5] = 1.0;
244 Q[16] = 1.0; Q[54] = 1.0;
245 },
246 marker.x(),
247 marker.h(),
248 {{FACE_METADATA_ACCESSOR}}.getOldTimeStamp(marker.getSelectedFaceNumber()<Dimensions ? 1 : 0),
249 repositories::{{SOLVER_INSTANCE}}.getMinTimeStepSize(),
250 {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}},
251 {{OVERLAP}},
252 {{NUMBER_OF_UNKNOWNS}},
253 {{NUMBER_OF_AUXILIARY_VARIABLES}},
254 marker.getSelectedFaceNumber(),
255 {0.0, 0.0, 0.0},
256 fineGridFace{{UNKNOWN_IDENTIFIER}}Old.value,
257 fineGridFace{{UNKNOWN_IDENTIFIER}}New.value,
258 """+timeStamp+"""
259 );
260"""
CCZ4Helper.get_body_of_adm_constraints
get_body_of_adm_constraints(patch_size, number_of_output_variable, so_flag)
Definition CCZ4Helper.py:21
CCZ4Helper.get_body_of_SommerfeldCondition
get_body_of_SommerfeldCondition(scenario, unknowns, auxiliary_variables, restart_time=-1)
Definition CCZ4Helper.py:174
CCZ4Helper.get_body_of_enforceCCZ4constraint
get_body_of_enforceCCZ4constraint()
Definition CCZ4Helper.py:1
CCZ4Helper.get_body_of_Psi_Calc
get_body_of_Psi_Calc(patch_size, number_of_output_variable, so_flag)
Definition CCZ4Helper.py:97