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