Peano
Loading...
Searching...
No Matches
ccz4_archived_24_01_19.py
Go to the documentation of this file.
1import os
2import argparse
3
4import peano4
5import exahype2
6
7import peano4.toolbox.particles
8import dastgen2
9
10from Probe_file_gene import tracer_seeds_generate
11
12from ComputeFirstDerivatives import ComputeFirstDerivativesFD4RK
13
14from peano4.toolbox.blockstructured.DynamicAMR import DynamicAMR
15import numpy as np
16#todo no mpmath here so rkdg not working
17
18def CoordListFromFile(file,scale):
19 f=open(file)
20 data=[]
21 tem=f.readlines()
22 for line in tem:
23 data.append(list(map(float,line.split(' ')[1:])))
24 print(len(data))
25 data=np.array(data)
26 return data*scale
27
28
29modes = {
30 "release": peano4.output.CompileMode.Release,
31 "trace": peano4.output.CompileMode.Trace,
32 "assert": peano4.output.CompileMode.Asserts, "stats": peano4.output.CompileMode.Stats,
33 "debug": peano4.output.CompileMode.Debug,
34}
35
36floatparams = {
37 "GLMc0":1.5, "GLMc":1.2, "GLMd":2.0, "GLMepsA":1.0, "GLMepsP":1.0,
38 "GLMepsD":1.0,
39 "itau":1.0, "k1":0.1, "k2":0.0, "k3":0.5, "eta":1.0,
40 "f":0.75, "g":0.0, "xi":1.0, "e":1.0, "c":1.0, "mu":0.2, "ds":1.0,
41 "sk":1.0, "bs":1.0,
42 "domain_r":0.5, "smoothing":0.0, "KOSigma":4.0 #, \
43 #"itau":1.0, "k1":0.1, "k2":0.0, "k3":0.5, "eta":1.0,
44 #"f":1.0, "g":0.0, "xi":1.0, "e":1.0, "c":1.0, "mu":0.2, "ds":1.0,
45 #"sk":1.0, "bs":1.0#, \
46 #"par_b":666.0, "center_offset_x":-1.0, "center_offset_y":0.0, "center_offset_z":0.0, \
47 #"target_m_plus":1.0, "par_p_plus_x":0.0, "par_p_plus_y":0.0, "par_p_plus_z":0.0, \
48 #"par_s_plus_x":0.0, "par_s_plus_y":0.0, "par_s_plus_z":0.0, \
49 #"target_m_minus":1.0, "par_p_minus_x":0.0, "par_p_minus_y":0.0, "par_p_minus_z":0.0, \
50 #"par_s_minus_x":0.0, "par_s_minus_y":0.0, "par_s_minus_z":0.0, \
51 #"tp_epsilon":1e-6
52}
53
54intparams = {"BBHType":2, "LapseType":1, "tp_grid_setup":0, "swi":99, "ReSwi":0, "source":0}
55#sss
56if __name__ == "__main__":
57 parser = argparse.ArgumentParser(description='ExaHyPE 2 - CCZ4 script')
58 parser.add_argument("-maxh", "--max-h", dest="max_h", type=float, required="True", help="upper limit for refinement. Refers to volume size, i.e. not to patch size" )
59 parser.add_argument("-minh", "--min-h", dest="min_h", type=float, default=0, help="lower limit for refinement (set to 0 to make it equal to max_h - default). Refers to volume size, i.e. not to patch size" )
60 parser.add_argument("-ps", "--patch-size", dest="patch_size", type=int, default=6, help="Patch size, i.e. number of volumes per patch per direction" )
61 parser.add_argument("-plt", "--plot-step-size", dest="plot_step_size", type=float, default=0.04, help="Plot step size (0 to switch it off)" )
62 parser.add_argument("-m", "--mode", dest="mode", default="release", help="|".join(modes.keys()) )
63 parser.add_argument( "--gpu", dest="GPU", default=False, action="store_true", help="Run with accelerator support" )
64 parser.add_argument("-ext", "--extension", dest="extension", choices=["none", "adm", "Psi4","phy-debug", "PT"], default="none", help="Pick extension, i.e. what should be plotted on top. Default is none" )
65 parser.add_argument("-impl", "--implementation", dest="implementation", choices=["fv-global-fixed", "fv-global-adaptive", "fv-global-fixed-enclave", "fv-global-adaptive-enclave","fv-subcyling-adaptive-enclave","fv-local-enclave", "fv-musclhancock", "fd4-rk1-adaptive", "fd4-rk1-adaptive-enclave", "fd4-rk1-fixed", "fd4-rk4-adaptive", "fd4-rk4-fixed", "RKDG"], required="True", help="Pick solver type" )
66 parser.add_argument("-no-pbc", "--no-periodic-boundary-conditions", dest="periodic_bc", action="store_false", default=True, help="switch on or off the periodic BC" )
67 parser.add_argument("-sommerfeld", "--sommerfeld-boundary-conditions", dest="sommerfeld_bc", action="store_true", default=False, help="switch on or off the Sommerfeld radiative BC" )
68 parser.add_argument("-AMRMarker", "--AMR-marker", dest="marker", choices=["none", "poisson"], default="none", help="switch on or off the AMR boundary marker" )
69 parser.add_argument("-et", "--end-time", dest="end_time", type=float, default=1.0, help="End (terminal) time" )
70 parser.add_argument("-pst", "--plot-start-time", dest="plot_start_time", type=float, default=0.0, help="start time for plot" )
71 parser.add_argument("-s", "--scenario", dest="scenario", choices=["gauge", "dia_gauge", "linear", "single-puncture","two-punctures", "flat"], required="True", help="Scenario" )
72 parser.add_argument("-cfl", "--CFL-ratio", dest="cfl", type=float, default=0.1, help="Set CFL ratio" )
73 parser.add_argument("-tracer", "--add-tracer", dest="add_tracer", type=int, default=0, help="Add tracers and specify the seeds. 0-switch off, 1-static point tracer, 2-moving point tracer" )
74 parser.add_argument("-tn", "--tracer-name", dest="tra_name", type=str, default="de", help="name of output tracer file (temporary)" )
75 parser.add_argument("-exn", "--exe-name", dest="exe_name", type=str, default="", help="name of output executable file" )
76 parser.add_argument("-outdir", "--output-directory", dest="path", type=str, default="./", help="specify the output directory, include the patch file and tracer file" )
77 parser.add_argument("-interp", "--interpolation", dest="interpolation", choices=["constant","order-2","linear-constant-extrap","linear-linear-extrap","linear-con-extrap-lin-normal-interp","linear-lin-extrap-lin-normal-interp"], default="linear-lin-extrap-lin-normal-interp", help="interpolation scheme for AMR" )
78 parser.add_argument("-restrict", "--restriction", dest="restriction", choices=["average", "inject" ], default="average", help="restriction scheme for AMR" )
79 parser.add_argument("-dc", "--double-communication", dest="DC_flag", default=False, action="store_true", help="enable double communication per timestep")
80
81
82 for k, v in floatparams.items(): parser.add_argument("--{}".format(k), dest="CCZ4{}".format(k), type=float, default=v, help="default: %(default)s")
83 for k, v in intparams.items():
84 if k=="ReSwi":
85 parser.add_argument("--{}".format(k), dest="CCZ4{}".format(k), type=int, default=v, help="default: %(default)s, choose refinement criterion, 0-no refinement, 1-radius based, 2-SBH phi gradient based, 3-BBH phi gradient based. Notice: 2 and 3 only work with -ext Full")
86 else: parser.add_argument("--{}".format(k), dest="CCZ4{}".format(k), type=int, default=v, help="default: %(default)s")
87
88 args = parser.parse_args()
89
90 SuperClass = None
91
92 print(args.implementation)
93 if args.implementation=="fv-global-fixed":
94 SuperClass = exahype2.solvers.fv.rusanov.GlobalFixedTimeStep
95 if args.implementation=="fv-global-adaptive":
96 SuperClass = exahype2.solvers.fv.rusanov.GlobalAdaptiveTimeStep
97 if args.implementation=="fv-global-fixed-enclave":
98 SuperClass = exahype2.solvers.fv.rusanov.GlobalFixedTimeStepWithEnclaveTasking
99 if args.implementation=="fv-global-adaptive-enclave":
100 SuperClass = exahype2.solvers.fv.rusanov.GlobalAdaptiveTimeStepWithEnclaveTasking
101 if args.implementation=="fv-subcyling-adaptive-enclave":
102 SuperClass = exahype2.solvers.fv.rusanov.SubcyclingAdaptiveTimeStepWithEnclaveTasking
103 if args.implementation=="fv-local-enclave":
104 SuperClass = exahype2.solvers.fv.rusanov.LocalTimeStepWithEnclaveTasking
105 if args.implementation=="fv-musclhancock":
106 SuperClass = exahype2.solvers.fv.musclhancock.GlobalAdaptiveTimeStep
107 if args.implementation=="fd4-rk1-adaptive" or args.implementation=="fd4-rk4-adaptive":
108 SuperClass = exahype2.solvers.rkfd.fd4.GlobalAdaptiveTimeStep
109 if args.implementation=="fd4-rk1-adaptive-enclave":
110 SuperClass = exahype2.solvers.rkfd.fd4.GlobalAdaptiveTimeStepWithEnclaveTasking
111 if args.implementation=="fd4-rk1-fixed" or args.implementation=="fd4-rk4-fixed":
112 SuperClass = exahype2.solvers.rkfd.fd4.GlobalFixedTimeStep
113 if args.implementation=="RKDG":
114 SuperClass = exahype2.solvers.rkdg.rusanov.GlobalAdaptiveTimeStepWithEnclaveTasking
115
116
117 class CCZ4Solver( SuperClass ):
118 def __init__(self, name, patch_size, min_volume_h, max_volume_h, cfl, domain_r, KOSig):
119 unknowns = {
120 "G":6, #0-5
121 "K":6, #6-11
122 "theta":1, #12
123 "Z":3, #13-15
124 "lapse":1, #16
125 "shift":3, #17-19
126 "b":3, #20-22
127 "dLapse":3, #23-25
128 "dxShift":3,#26-28
129 "dyShift":3,#29-31
130 "dzShift":3,#32-34
131 "dxG":6, #35-40
132 "dyG":6, #41-46
133 "dzG":6, #47-52
134 "traceK":1, #53
135 "phi":1, #54
136 "P":3, #55-57
137 "K0":1, #58
138 }
139
140 number_of_unknowns = sum(unknowns.values())
141
142 self._my_user_includes = """
143#include "../CCZ4Kernels.h"
144"""
145
146 if SuperClass==exahype2.solvers.fv.rusanov.GlobalFixedTimeStep:
147 SuperClass.__init__(
148 self,
149 name=name, patch_size=patch_size,
150 unknowns=number_of_unknowns,
151 auxiliary_variables=0,
152 min_volume_h=min_volume_h, max_volume_h=max_volume_h,
153 normalised_time_step_size=1e-2
154 )
155 elif SuperClass==exahype2.solvers.fv.rusanov.GlobalFixedTimeStepWithEnclaveTasking:
156 SuperClass.__init__(
157 self,
158 name=name, patch_size=patch_size,
159 unknowns=number_of_unknowns,
160 auxiliary_variables=0,
161 min_volume_h=min_volume_h, max_volume_h=max_volume_h,
162 time_step_size=1e-2,
163 use_gpu =args.GPU #=="fv-fixed-gpu" else False
164 )
165 elif SuperClass==exahype2.solvers.fv.rusanov.GlobalAdaptiveTimeStep:
166 SuperClass.__init__(
167 self,
168 name=name, patch_size=patch_size,
169 unknowns=number_of_unknowns,
170 auxiliary_variables=0,
171 min_volume_h=min_volume_h, max_volume_h=max_volume_h,
172 time_step_relaxation=cfl
173 )
174 elif SuperClass==exahype2.solvers.rkfd.fd4.GlobalAdaptiveTimeStep:
175 if args.implementation=="fd4-rk1-adaptive" or args.implementation=="fd4-rk1-adaptive-enclave":
176 rk_order = 1
177 else:
178 rk_order = 4
179 SuperClass.__init__(
180 self,
181 name=name, patch_size=patch_size, rk_order=rk_order,
182 unknowns=number_of_unknowns,
183 auxiliary_variables=0,
184 min_meshcell_h=min_volume_h, max_meshcell_h=max_volume_h,
185 time_step_relaxation=cfl,
186 KOSigma=KOSig
187 )
188 elif SuperClass==exahype2.solvers.rkfd.fd4.GlobalAdaptiveTimeStepWithEnclaveTasking:
189 if args.implementation=="fd4-rk1-adaptive" or args.implementation=="fd4-rk1-adaptive-enclave":
190 rk_order = 1
191 else:
192 rk_order = 4
193 SuperClass.__init__(
194 self,
195 name=name, patch_size=patch_size, rk_order=rk_order,
196 unknowns=number_of_unknowns,
197 auxiliary_variables=0,
198 min_meshcell_h=min_volume_h, max_meshcell_h=max_volume_h,
199 time_step_relaxation=cfl,
200 KOSigma=KOSig,
201 pde_terms_without_state=True
202 )
203 elif SuperClass==exahype2.solvers.rkfd.fd4.GlobalFixedTimeStep or SuperClass==exahype2.solvers.rkfd.fd4.GlobalFixedTimeStepWithEnclaveTasking:
204 if args.implementation=="fd4-rk1-fixed":
205 rk_order = 1
206 else:
207 rk_order = 4
208 SuperClass.__init__(
209 self,
210 name=name, patch_size=patch_size, rk_order=rk_order,
211 unknowns=number_of_unknowns,
212 auxiliary_variables=0,
213 min_meshcell_h=min_volume_h, max_meshcell_h=max_volume_h,
214 normalised_time_step_size=0.01,
215 KOSigma=KOSig
216 )
217 elif SuperClass == exahype2.solvers.rkdg.rusanov.GlobalAdaptiveTimeStep:
218 SuperClass.__init__(
219 self,
220 name=name,
221 rk_order = 2,
222 polynomials = exahype2.solvers.GaussLegendreBasis(2),
223 face_projections = exahype2.solvers.rkdg.FaceProjections.Solution,
224 unknowns=number_of_unknowns,
225 auxiliary_variables=0,
226 min_cell_h=min_volume_h, max_cell_h=max_volume_h, #not a good name, but use for now.
227 time_step_relaxation=cfl
228 )
229 else:
230 SuperClass.__init__(
231 self,
232 name=name, patch_size=patch_size,
233 unknowns=number_of_unknowns,
234 auxiliary_variables=0,
235 min_volume_h=min_volume_h, max_volume_h=max_volume_h,
236 time_step_relaxation=cfl
237 #use_gpu =args.GPU #=="fv-fixed-gpu" else False
238# use_gpu = True if args.implementation=="fv-adaptive-gpu" else False
239 )
240
241 self._patch_size = patch_size
242 self._domain_r = domain_r
243
244 self.set_implementation(
245 boundary_conditions=exahype2.solvers.PDETerms.User_Defined_Implementation,
246 ncp=exahype2.solvers.PDETerms.User_Defined_Implementation,
247 flux=exahype2.solvers.PDETerms.None_Implementation,
248 source_term=exahype2.solvers.PDETerms.User_Defined_Implementation,
249 refinement_criterion = exahype2.solvers.PDETerms.User_Defined_Implementation,
250 eigenvalues = exahype2.solvers.PDETerms.User_Defined_Implementation
251 )
252
253 if SuperClass==exahype2.solvers.rkfd.fd4.GlobalAdaptiveTimeStepWithEnclaveTasking:
254 self._fused_compute_kernel_call_cpu = exahype2.solvers.rkfd.fd4.kernels.create_compute_kernel_for_FD4(
255 self._flux_implementation,
256 self._ncp_implementation,
257 self._source_term_implementation,
258 compute_max_eigenvalue_of_next_time_step = True,
259 solver_variant = exahype2.solvers.rkfd.kernels.SolverVariant.Multicore,
260 kernel_variant = exahype2.solvers.rkfd.kernels.KernelVariant.BatchedAoSHeap,
261 KOSigma = self._KO_Sigma
262 )
263
264 self.switch_to_heap_storage(
265 exahype2.solvers.rkfd.Storage.Heap, exahype2.solvers.rkfd.Storage.Heap
266 )
267
268 self.postprocess_updated_patch = """
269{
270"""
271
272 if (
273 SuperClass==exahype2.solvers.rkfd.fd4.GlobalAdaptiveTimeStep
274 or SuperClass==exahype2.solvers.rkfd.fd4.GlobalFixedTimeStep
275 or SuperClass==exahype2.solvers.rkfd.fd4.GlobalAdaptiveTimeStepWithEnclaveTasking
276 #or SuperClass==exahype2.solvers.rkfd.fd4.GlobalFixedTimeStepWithEnclaveTasking no fixed enclave version yet
277 ):
278 self.postprocess_updated_patch += """
279 #if Dimensions==2
280 constexpr int itmax = {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}} * {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}};
281 #endif
282
283 #if Dimensions==3
284 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}};
285 #endif
286"""
287 #elif SuperClass == exahype2.solvers.rkdg.rusanov.GlobalAdaptiveTimeStep: fix it after rkdg is back
288 # self.postprocess_updated_patch = """
289 #constexpr int itmax = ({{ORDER}}+1) * ({{ORDER}}+1) * ({{ORDER}}+1); // only support 3d
290#"""
291 else:
292 self.postprocess_updated_patch += """
293 #if Dimensions==2
294 constexpr int itmax = {{NUMBER_OF_VOLUMES_PER_AXIS}} * {{NUMBER_OF_VOLUMES_PER_AXIS}};
295 #endif
296
297 #if Dimensions==3
298 constexpr int itmax = {{NUMBER_OF_VOLUMES_PER_AXIS}} * {{NUMBER_OF_VOLUMES_PER_AXIS}} * {{NUMBER_OF_VOLUMES_PER_AXIS}};
299 #endif
300"""
301
302 self.postprocess_updated_patch += """
303 int index = 0;
304 for (int i=0;i<itmax;i++)
305 {
306 applications::exahype2::ccz4::enforceCCZ4constraints( newQ+index );
307 index += {{NUMBER_OF_UNKNOWNS}} + {{NUMBER_OF_AUXILIARY_VARIABLES}};
308 }
309 }
310"""
311
312
313 def create_action_sets(self):
314 SuperClass.create_action_sets(self)
315 self._action_set_couple_resolution_transitions_and_handle_dynamic_mesh_refinement.additional_includes += """
316#include "../CCZ4Kernels.h"
317 """
318
319 def get_user_action_set_includes(self):
320 """
321 We take this routine to add a few additional include statements.
322 """
323 return SuperClass.get_user_action_set_includes(self) + self._my_user_includes
324
325 def add_adm_constriants(self):
326 """
327 @add adm constriants (Hamilton and Momentum)
328 """
329 self._auxiliary_variables = 7
330
331 self._my_user_includes += """
332 #include "../CCZ4Kernels.h"
333 """
334
335 self._preprocess_reconstructed_patch = """
336 const int patchSize = """ + str( self._patch.dim[0] ) + """;
337 double volumeH = ::exahype2::fv::getVolumeLength(marker.h(),patchSize);
338
339 const int n_a_v=7;
340 const int overlap=3; //make sure you are using fd4 solver!
341
342 //if (not marker.willBeRefined() and repositories::InstanceOfFiniteVolumeCCZ4.getSolverState()!=FiniteVolumeCCZ4::SolverState::GridConstruction and repositories::InstanceOfFiniteVolumeCCZ4.getSolverState()==FiniteVolumeCCZ4::SolverState::RungeKuttaSubStep0) {
343 if (not marker.willBeRefined() and repositories::InstanceOfFiniteVolumeCCZ4.getSolverState()!=FiniteVolumeCCZ4::SolverState::GridConstruction) {
344 dfor(cell,patchSize) {
345 tarch::la::Vector<Dimensions,int> currentCell = cell + tarch::la::Vector<Dimensions,int>(overlap);
346
347 double gradQ[3*59]={ 0 };
348
349 for (int d=0; d<3; d++) {
350 tarch::la::Vector<Dimensions,int> leftCell = currentCell;
351 tarch::la::Vector<Dimensions,int> rightCell = currentCell;
352 leftCell(d) -= 1;
353 rightCell(d) += 1;
354 const int leftCellSerialised = peano4::utils::dLinearised(leftCell, patchSize + 2*overlap);
355 const int rightCellSerialised = peano4::utils::dLinearised(rightCell,patchSize + 2*overlap);
356 for(int i=0; i<59; i++) {
357 gradQ[d*59+i] = ( oldQWithHalo[rightCellSerialised*(59+n_a_v)+i] - oldQWithHalo[leftCellSerialised*(59+n_a_v)+i] ) / 2.0 / volumeH;
358 }
359 }
360
361 const int cellSerialised = peano4::utils::dLinearised(currentCell, patchSize + 2*overlap);
362
363 double constraints[7]={0,0,0,0,0,0,0};
364 admconstraints(constraints, oldQWithHalo+cellSerialised*(59+n_a_v), gradQ);
365 ::exahype2::enumerator::AoSLexicographicEnumerator enumeratorWithAuxiliaryVariables ( 1, patchSize, 0, 59, n_a_v);
366 for (int i=0;i<7;i++){
367 fineGridCellFiniteVolumeCCZ4Q.value[ enumeratorWithAuxiliaryVariables(0,cell,59+i) ] = constraints[i];
368 }
369 }
370 }
371 """
372
373 self.create_data_structures()
374 self.create_action_sets()
375
376 def add_Psi4W(self):
377 """
378 add psi4 writer
379 """
380 self._auxiliary_variables = 2
381
382 self._my_user_includes += """
383 #include "../libtwopunctures/TP_PunctureTracker.h"
384 #include "../CCZ4Kernels.h"
385 """
386
387 self._preprocess_reconstructed_patch = """
388 const int patchSize = """ + str( self._patch.dim[0] ) + """;
389 double volumeH = ::exahype2::fv::getVolumeLength(marker.h(),patchSize);
390
391 const int n_a_v=2;
392 const int overlap=3; //make sure you are using fd4 solver!
393
394 //if (not marker.willBeRefined() and repositories::InstanceOfFiniteVolumeCCZ4.getSolverState()!=FiniteVolumeCCZ4::SolverState::GridConstruction and repositories::InstanceOfFiniteVolumeCCZ4.getSolverState()==FiniteVolumeCCZ4::SolverState::RungeKuttaSubStep0) {
395 if (not marker.willBeRefined() and repositories::InstanceOfFiniteVolumeCCZ4.getSolverState()!=FiniteVolumeCCZ4::SolverState::GridConstruction) {
396 dfor(cell,patchSize) {
397 tarch::la::Vector<Dimensions,int> currentCell = cell + tarch::la::Vector<Dimensions,int>(overlap);
398
399 double gradQ[3*59]={ 0 };
400
401 for (int d=0; d<3; d++) {
402 tarch::la::Vector<Dimensions,int> leftCell = currentCell;
403 tarch::la::Vector<Dimensions,int> rightCell = currentCell;
404 leftCell(d) -= 1;
405 rightCell(d) += 1;
406 const int leftCellSerialised = peano4::utils::dLinearised(leftCell, patchSize + 2*overlap);
407 const int rightCellSerialised = peano4::utils::dLinearised(rightCell,patchSize + 2*overlap);
408 for(int i=0; i<59; i++) {
409 gradQ[d*59+i] = ( oldQWithHalo[rightCellSerialised*(59+n_a_v)+i] - oldQWithHalo[leftCellSerialised*(59+n_a_v)+i] ) / 2.0 / volumeH;
410 }
411 }
412
413 const int cellSerialised = peano4::utils::dLinearised(currentCell, patchSize + 2*overlap);
414
415 double Psi4[2]={0,0};
416 double currentPosition[3];
417 for (int d=0; d<3; d++) currentPosition[d]=marker.getOffset()(d)+(cell(d)+0.5)*volumeH;
418 Psi4Calc(Psi4, oldQWithHalo+cellSerialised*(59+n_a_v), gradQ, currentPosition);
419 ::exahype2::enumerator::AoSLexicographicEnumerator enumeratorWithAuxiliaryVariables ( 1, patchSize, 0, 59, n_a_v);
420 fineGridCellFiniteVolumeCCZ4Q.value[ enumeratorWithAuxiliaryVariables(0,cell,59+0) ] = Psi4[0];
421 fineGridCellFiniteVolumeCCZ4Q.value[ enumeratorWithAuxiliaryVariables(0,cell,59+1) ] = Psi4[1];
422 }
423 }
424 """
425
426 self.create_data_structures()
427 self.create_action_sets()
428
429 def add_intermediate_quantities(self):
430 """
431 @add adm intermediate quantites
432 """
433 self._auxiliary_variables = 15
434
435 self._my_user_includes += """
436 #include "../CCZ4Kernels.h"
437 """
438
439 self._preprocess_reconstructed_patch = """
440 const int patchSize = """ + str( self._patch.dim[0] ) + """;
441 double volumeH = ::exahype2::fv::getVolumeLength(marker.h(),patchSize);
442
443 const int n_a_v=""" + str( self._auxiliary_variables ) + """;
444 const int overlap=3; //make sure you are using fd4 solver!
445
446 //if (not marker.willBeRefined() and repositories::InstanceOfFiniteVolumeCCZ4.getSolverState()!=FiniteVolumeCCZ4::SolverState::GridConstruction and repositories::InstanceOfFiniteVolumeCCZ4.getSolverState()==FiniteVolumeCCZ4::SolverState::RungeKuttaSubStep0) {
447 if (not marker.willBeRefined() and repositories::InstanceOfFiniteVolumeCCZ4.getSolverState()!=FiniteVolumeCCZ4::SolverState::GridConstruction) {
448 dfor(cell,patchSize) {
449 tarch::la::Vector<Dimensions,int> currentCell = cell + tarch::la::Vector<Dimensions,int>(overlap);
450
451 double gradQ[3*59]={ 0 };
452
453 for (int d=0; d<3; d++) {
454 tarch::la::Vector<Dimensions,int> leftCell = currentCell;
455 tarch::la::Vector<Dimensions,int> leftCell2 = currentCell;
456 tarch::la::Vector<Dimensions,int> rightCell = currentCell;
457 tarch::la::Vector<Dimensions,int> rightCell2 = currentCell;
458 leftCell(d) -= 1; leftCell2(d) -= 2;
459 rightCell(d) += 1; rightCell2(d) += 2;
460 const int leftCellSerialised = peano4::utils::dLinearised(leftCell, patchSize + 2*overlap);
461 const int rightCellSerialised = peano4::utils::dLinearised(rightCell,patchSize + 2*overlap);
462 const int leftCellSerialised2 = peano4::utils::dLinearised(leftCell2, patchSize + 2*overlap);
463 const int rightCellSerialised2 = peano4::utils::dLinearised(rightCell2,patchSize + 2*overlap);
464 for(int i=0; i<59; i++) {
465 gradQ[d*59+i] = ( (2.0/3.0)*oldQWithHalo[rightCellSerialised*(59+n_a_v)+i] - (2.0/3.0)*oldQWithHalo[leftCellSerialised*(59+n_a_v)+i]
466 +(1.0/12.0)*oldQWithHalo[leftCellSerialised2*(59+n_a_v)+i] - (1.0/12.0)*oldQWithHalo[rightCellSerialised2*(59+n_a_v)+i] )/volumeH;
467 }
468 }
469
470 const int cellSerialised = peano4::utils::dLinearised(currentCell, patchSize + 2*overlap);
471
472 double terms[""" + str( self._auxiliary_variables ) + """]={0,0,0,0,0,0,0};
473 TestingOutput(terms, oldQWithHalo+cellSerialised*(59+n_a_v), gradQ);
474 ::exahype2::enumerator::AoSLexicographicEnumerator enumeratorWithAuxiliaryVariables ( 1, patchSize, 0, 59, n_a_v);
475 for (int i=0;i<""" + str( self._auxiliary_variables ) + """;i++){
476 fineGridCellFiniteVolumeCCZ4Q.value[ enumeratorWithAuxiliaryVariables(0,cell,59+i) ] = terms[i];
477 }
478 }
479 }
480 """
481
482 self.create_data_structures()
483 self.create_action_sets()
484
485 def add_puncture_tracker(self):
486 """
487 @add
488 """
489 self._auxiliary_variables = 0
490
491 self._user_solver_includes += """
492 #include "libtwopunctures/TP_PunctureTracker.h"
493 #include "CCZ4Kernels.h"
494 #include <fstream>
495 #include <iomanip>
496 """
497
498# self._user_action_set_includes += """
499# #include "../libtwopunctures/TP_PunctureTracker.h"
500# #include "../CCZ4Kernels.h"
501# #include <fstream>
502# #include <iomanip>
503# """
504
505 self._preprocess_reconstructed_patch = """
506 const int patchSize = """ + str( self._patch.dim[0] ) + """;
507 double volumeH = ::exahype2::fv::getVolumeLength(marker.h(),patchSize);
508
509 std::fstream fin;
510 std::string att="_"""+args.exe_name+""".txt"; std::string p1="puncture1"; std::string p2="puncture2"; std::string tem="ztem";
511 const int n_a_v=""" + str( self._auxiliary_variables ) + """;
512 const double overlap=3;
513
514 if (tarch::la::equals(timeStamp,0.0)){//initialization
515 fin.open((p1+att),std::ios::out|std::ios::trunc);
516 fin << "2.0 0.0 0.0 0.0 0.0" << std::endl;//4.461538
517 fin.close();
518 fin.open((p2+att),std::ios::out|std::ios::trunc);
519 fin << "-2.0 0.0 0.0 0.0 0.0" << std::endl;//-5.538462
520 fin.close();
521 //fin.open((tem+att),std::ios::out|std::ios::trunc);
522 //fin << "tem file" << std::endl;
523 //fin.close();
524 } else {
525 fin.open((p1+att),std::ios::in);
526 std::string pos=getLastLine(fin);
527 fin.close();
528 double coor1[5]={0};//read in previous coordinates
529 CoorReadIn(coor1,pos);
530 fin.open((p2+att),std::ios::in);
531 std::string pos2=getLastLine(fin);
532 fin.close();
533 double coor2[5]={0};
534 CoorReadIn(coor2,pos2);
535 int inter_number=4;
536 if (marker.isContained(coor1)){
537 tarch::la::Vector<Dimensions*2,int> IndexOfCell=FindCellIndex(coor1,marker.getOffset(),volumeH,patchSize); //find where the puncture is
538 tarch::la::Vector<Dimensions,int> IndexForInterpolate[8];
539 FindInterIndex(IndexForInterpolate,IndexOfCell);//find the closest 8 cells
540 double raw[32]={0};
541 for (int i=0;i<8;i++){
542 int Lindex=peano4::utils::dLinearised(IndexForInterpolate[i], patchSize + 2*overlap);
543 for (int j=0;j<Dimensions;j++) {raw[i*inter_number+j]=oldQWithHalo[Lindex*(59+n_a_v)+17+j];} //read in corresponding beta
544 raw[i*inter_number+3]=oldQWithHalo[Lindex*(59+n_a_v)+16];
545 //std::cout << IndexForInterpolate[i](0) << " " << IndexForInterpolate[i](1) << " " << IndexForInterpolate[i](2) << std::endl;
546 //std::cout << raw[i*inter_number+0] << " " << raw[i*inter_number+1] << " " << raw[i*inter_number+2] << " " << raw[i*inter_number+3] << std::endl;
547 //::exahype2::enumerator::AoSLexicographicEnumerator enumeratorWithAuxiliaryVariables ( 1, patchSize, overlap, 59, n_a_v);
548 //std::cout << oldQWithHalo[enumeratorWithAuxiliaryVariables(0,IndexForInterpolate[i],0)] << std::endl;
549 }
550 double result[inter_number];
551 Interpolation(result,IndexForInterpolate,raw,coor1,marker.getOffset(),volumeH,patchSize);
552
553 coor1[0]-=timeStepSize*result[0]; coor1[1]-=timeStepSize*result[1]; coor1[2]-=timeStepSize*result[2];//updates position
554 fin.open((p1+att),std::ios::app);//output
555 fin << std::setprecision (17) << coor1[0] << " " << coor1[1];
556 fin << " " << coor1[2] << " " << timeStamp << " " << result[3] << std::endl;
557 fin.close();
558 //fin.open((tem+att),std::ios::app);
559 //fin << "for cellindex" << IndexOfCell(0) << " " << IndexOfCell(1) << " " << IndexOfCell(2) << " " << IndexOfCell(3) << " " << IndexOfCell(4) << " " << IndexOfCell(5) << std::endl;
560 //for (int i=0;i<8;i++){
561 //fin << IndexForInterpolate[i](0) << " " << IndexForInterpolate[i](1) << " " << IndexForInterpolate[i](2) << std::endl;
562 //fin << raw[i*inter_number+0] << " " << raw[i*inter_number+1] << " " << raw[i*inter_number+2] << " " << raw[i*inter_number+3] << std::endl;
563 //}
564
565 //fin << "after inter" << result[0] << " " << result[1] << " " << result[2] << " " << result[3] << " " << timeStamp << std::endl;
566 //fin.close();
567 }
568 if (marker.isContained(coor2)){//do the same for the second puncutre
569 tarch::la::Vector<Dimensions*2,int> IndexOfCell=FindCellIndex(coor2,marker.getOffset(),volumeH,patchSize);
570 tarch::la::Vector<Dimensions,int> IndexForInterpolate[8];
571 FindInterIndex(IndexForInterpolate,IndexOfCell);
572 double raw[8*inter_number];
573 for (int i=0;i<8;i++){
574 int Lindex=peano4::utils::dLinearised(IndexForInterpolate[i], patchSize + 2*overlap);
575 for (int j=0;j<Dimensions;j++) {raw[i*inter_number+j]=oldQWithHalo[Lindex*(59+n_a_v)+17+j];}
576 raw[i*inter_number+3]=oldQWithHalo[Lindex*(59+n_a_v)+16];
577 }
578 double result[inter_number];
579 Interpolation(result,IndexForInterpolate,raw,coor2,marker.getOffset(),volumeH,patchSize);
580
581 coor2[0]-=timeStepSize*result[0]; coor2[1]-=timeStepSize*result[1]; coor2[2]-=timeStepSize*result[2];
582 fin.open((p2+att),std::ios::app);
583 fin << std::setprecision (17) << coor2[0] << " " << coor2[1];
584 fin << " " << coor2[2] << " " << timeStamp << " " << result[3] << std::endl;
585 fin.close();
586 }
587 }
588 """
589
590 self.create_data_structures()
591 self.create_action_sets()
592
593
594
597 userinfo = []
598 exe="peano4"
599
600 if args.exe_name!="":
601 exe += "_"
602 exe += args.exe_name
603 if not args.tra_name=="de":
604 exe += "_" + args.tra_name
605 project = exahype2.Project( ["applications", "exahype2", "ccz4"], "ccz4", executable=exe)
606
607
610 is_aderdg = False
611 is_rkdg = False
612 solver_name = "CCZ4"
613 try:
614 if SuperClass==exahype2.solvers.aderdg.NonFusedGenericRusanovFixedTimeStepSize:
615 is_aderdg = True
616 order = 3
617 unknowns = 59
618 time_step_size = 0.001
619 except Exception as e:
620 pass
621 #msg = "Warning: ADER-DG no supported on this machine"
622 #print(msg)
623 #userinfo.append((msg,e))
624
625 try:
626 if SuperClass==exahype2.solvers.rkdg.rusanov.GlobalAdaptiveTimeStep:
627 is_rkdg = True
628 except Exception as e:
629 pass
630 msg = "Warning: RKDG not supported on this machine"
631 print(msg)
632 userinfo.append((msg,e))
633
634 if is_aderdg:
635 solver_name = "ADERDG" + solver_name
636 elif is_rkdg:
637 solver_name = "RKDG" + solver_name
638 else:
639 solver_name = "FiniteVolume" + solver_name
640
641 if args.DC_flag==True:
642 args.cfl=args.cfl/2
643
644 min_h = args.min_h
645 if min_h <=0.0:
646 print( "No minimal mesh size chosen. Set it to max mesh size (regular grid)" )
647 min_h = args.max_h
648
649 if is_aderdg:
650 my_solver = exahype2.solvers.aderdg.NonFusedGenericRusanovFixedTimeStepSize(
651 solver_name, order, unknowns, 0, #auxiliary_variables
652 exahype2.solvers.aderdg.Polynomials.Gauss_Legendre,
653 min_h, args.max_h, time_step_size,
654 flux = None,
655 ncp = exahype2.solvers.PDETerms.User_Defined_Implementation,
656 sources = exahype2.solvers.PDETerms.User_Defined_Implementation
657 )
658 else:
659 my_solver = CCZ4Solver(solver_name, args.patch_size, min_h, args.max_h, args.cfl, args.CCZ4domain_r,args.CCZ4KOSigma)
660 userinfo.append(("CFL ratio set as "+str(args.cfl), None))
661
662 userinfo.append(("The solver is "+args.implementation, None))
663
664
667 print(args.marker)
668 if args.marker=="poisson":
669 amr_label = exahype2.solvers.elliptic.AMRMarker( name="AMRMarker" )
670 project.add_solver( amr_label )
671 my_solver._auxiliary_variables+=1
672 print( "Add Poisson AMR marker for variable {}".format(my_solver._auxiliary_variables-1) )
673 amr_label.couple_with_FV_solver( my_solver, my_solver._auxiliary_variables-1 )
674
675
678
679 if args.interpolation=="order-2":
680 my_solver.overlap=2
681
682 if args.interpolation=="constant":
683 my_solver.interpolation = "piecewise_constant"
684 print( "Interpolation rule: piecewise_constant" )
685 if args.interpolation=="linear-constant-extrap":
686 my_solver.interpolation = "linear_with_constant_extrapolation"
687 print( "Interpolation rule: linear constant extrapolation" )
688 if args.interpolation=="linear-linear-extrap":
689 my_solver.interpolation = "linear_with_linear_extrapolation"
690 print( "Interpolation rule: linear extrapolation" )
691 if args.interpolation=="linear-con-extrap-lin-normal-interp":
692 my_solver.interpolation = "linear_with_constant_extrapolation_and_linear_normal_interpolation"
693 print( "Interpolation rule: linear+constant extrapolation and linear normal interpolation" )
694
695 if args.interpolation=="order-2":
696 my_solver.interpolation = "linear"
697
698 #Below are the I&R set functions for FD4
699 #These functions here do two things:
700 #1. it set the solver.interpolation parameters similar to what we do for fv, so the code know what template shoule be named.
701 #2. it computes and write the corresponding matrices for named schemes into the solver head files.
702 #notice point 2 is not needed as for build-in function (for fv) the matrices are already hard-coded.
703 tem_interp=["TP_constant","TP_linear_with_linear_extrap_normal_interp"]
704 tem_restrict=["TP_inject_normal_extrap","TP_average_normal_extrap"]
705 if ("fd4-rk" in args.implementation):
706 exahype2.solvers.rkfd.fd4.switch_to_FD4_tensor_product_interpolation( my_solver, tem_interp[1] )
707 exahype2.solvers.rkfd.fd4.switch_to_FD4_tensor_product_restriction( my_solver, tem_restrict[1] )
708 userinfo.append(("FD4 Interpolation: " + tem_interp[1] + " & Restriction: " + tem_restrict[1], None))
709 #if args.interpolation=="linear+enforce" or args.interpolation=="linear-slow+enforce":
710 # self.point_wise_postprocessing_of_interpolation = "examples::exahype2::ccz4::enforceCCZ4constraints(targetVolume)" )
711 #if args.restriction=="average+enforce" or args.restriction=="inject+enforce":
712 # self.point_wise_postprocessing_of_restriction = "examples::exahype2::ccz4::enforceCCZ4constraints(targetVolume)"
713
714
715
716
719 if args.extension=="adm":
720 my_solver.add_adm_constriants()
721 if args.extension=="Psi4":
722 my_solver.add_Psi4W()
723 if args.extension=="phy-debug":
724 my_solver.add_intermediate_quantities()
725 if args.extension=="PT":
726 my_solver.add_puncture_tracker()
727
728 #my_solver._auxiliary_variables = 59*3
729
730
733 for k, v in intparams.items():
734 intparams.update({k:eval("args.CCZ4{}".format(k))})
735 for k, v in floatparams.items():
736 floatparams.update({k:eval("args.CCZ4{}".format(k))})
737
738 if args.scenario=="two-punctures":
739 msg = "Periodic BC deactivated because you pick Puncture scenario\nInitialize binary black holes"
740 print(msg)
741 periodic_boundary_conditions = [False,False,False]
742 intparams.update({"swi":2}) #notice it may change, see according function in FiniteVolumeCCZ4.cpp
743 print(intparams)
744 userinfo.append((msg,None))
745 elif args.scenario=="single-puncture":
746 msg = "Periodic BC deactivated because you pick Puncture scenario\nInitialize single black hole"
747 print(msg)
748 periodic_boundary_conditions = [False,False,False]
749 intparams.update({"swi":0})
750 userinfo.append((msg,None))
751 elif args.periodic_bc==True:
752 msg = "Periodic BC set"
753 print(msg)
754 periodic_boundary_conditions = [True,True,True] # Periodic BC
755 userinfo.append((msg,None))
756 else:
757 msg = "WARNING: Periodic BC deactivated by hand"
758 print(msg)
759 periodic_boundary_conditions = [False,False,False]
760 userinfo.append((msg,None))
761
762 if args.sommerfeld_bc==True:
763 msg = "set Sommerfeld boundary condition"
764 userinfo.append((msg,None))
765 periodic_boundary_conditions = [False,False,False]
766 if args.scenario=="single-puncture": #please notice the default black hole mass (1.0) is hard-coded. need to be modified accordingly if necessary.
767 my_solver._action_set_handle_boundary.TemplateHandleBoundary_KernelCalls = """
768 ::exahype2::fd::applySommerfeldConditions(
769 [&](
770 const double * __restrict__ Q,
771 const tarch::la::Vector<Dimensions,double>& faceCentre,
772 const tarch::la::Vector<Dimensions,double>& gridCellH,
773 double t,
774 double dt,
775 int normal
776 ) -> double {
777 return repositories::{{SOLVER_INSTANCE}}.maxEigenvalue( Q, faceCentre, gridCellH, t, dt, normal );
778 },
779 [&](
780 double * __restrict__ Q,
781 const tarch::la::Vector<Dimensions,double>& faceCentre,
782 const tarch::la::Vector<Dimensions,double>& gridCellH
783 ) -> void {
784 for (int i=0; i<"""+str(my_solver._unknowns+my_solver._auxiliary_variables)+"""; i++) {
785 Q[i] = 0.0;
786 }
787 Q[0] = 1.0; Q[3] = 1.0; Q[5] = 1.0;
788 //const double r=tarch::la::norm2(faceCentre);
789 //Q[16] = 0.5*(1+(1-1.0/2/r)/(1+1.0/2/r));
790 //Q[54] = 1/(1+1.0/2/r)/(1+1.0/2/r);
791 Q[16] = 1.0; Q[54] = 1.0;
792 },
793 marker.x(),
794 marker.h(),
795 {{FACE_METADATA_ACCESSOR}}.getOldTimeStamp(marker.getSelectedFaceNumber()<Dimensions ? 1 : 0),
796 repositories::{{SOLVER_INSTANCE}}.getMinTimeStepSize(),
797 {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}},
798 {{OVERLAP}},
799 {{NUMBER_OF_UNKNOWNS}},
800 {{NUMBER_OF_AUXILIARY_VARIABLES}},
801 marker.getSelectedFaceNumber(),
802 {0.0, 0.0, 0.0},
803 fineGridFace{{UNKNOWN_IDENTIFIER}}Old.value,
804 fineGridFace{{UNKNOWN_IDENTIFIER}}New.value
805 );
806 """
807 else:
808 my_solver._action_set_handle_boundary.TemplateHandleBoundary_KernelCalls = """
809 ::exahype2::fd::applySommerfeldConditions(
810 [&](
811 const double * __restrict__ Q,
812 const tarch::la::Vector<Dimensions,double>& faceCentre,
813 const tarch::la::Vector<Dimensions,double>& gridCellH,
814 double t,
815 double dt,
816 int normal
817 ) -> double {
818 return repositories::{{SOLVER_INSTANCE}}.maxEigenvalue( Q, faceCentre, gridCellH, t, dt, normal );
819 },
820 [&](
821 double * __restrict__ Q,
822 const tarch::la::Vector<Dimensions,double>& faceCentre,
823 const tarch::la::Vector<Dimensions,double>& gridCellH
824 ) -> void {
825 for (int i=0; i<"""+str(my_solver._unknowns+my_solver._auxiliary_variables)+"""; i++) {
826 Q[i] = 0.0;
827 }
828 Q[0] = 1.0; Q[3] = 1.0; Q[5] = 1.0;
829 Q[16] = 1.0; Q[54] = 1.0;
830 },
831 marker.x(),
832 marker.h(),
833 {{FACE_METADATA_ACCESSOR}}.getOldTimeStamp(marker.getSelectedFaceNumber()<Dimensions ? 1 : 0),
834 repositories::{{SOLVER_INSTANCE}}.getMinTimeStepSize(),
835 {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}},
836 {{OVERLAP}},
837 {{NUMBER_OF_UNKNOWNS}},
838 {{NUMBER_OF_AUXILIARY_VARIABLES}},
839 marker.getSelectedFaceNumber(),
840 {0.0, 0.0, 0.0},
841 fineGridFace{{UNKNOWN_IDENTIFIER}}Old.value,
842 fineGridFace{{UNKNOWN_IDENTIFIER}}New.value
843 );
844 """
845
846
847
848 solverconstants=""
849
850 if args.scenario=="gauge":
851 solverconstants+= "static constexpr int Scenario=0; /* Gauge wave */ \n "
852 userinfo.append(("picking gauge wave scenario",None))
853 floatparams.update({"sk":0.0}); floatparams.update({"bs":0.0})
854 intparams.update({"LapseType":0})
855 elif args.scenario=="dia_gauge":
856 solverconstants+= "static constexpr int Scenario=3; /* Diagonal Gauge wave */ \n "
857 userinfo.append(("picking diagonal gauge wave scenario",None))
858 floatparams.update({"sk":0.0}); floatparams.update({"bs":0.0})
859 intparams.update({"LapseType":0})
860 elif args.scenario=="linear":
861 solverconstants+= "static constexpr int Scenario=1; /* Linear wave */ \n "
862 userinfo.append(("picking linear wave scenario",None))
863 floatparams.update({"sk":0.0}); floatparams.update({"bs":0.0})
864 intparams.update({"LapseType":0})
865 elif args.scenario=="flat":
866 solverconstants+= "static constexpr int Scenario=4; /* flat spacetime */ \n "
867 userinfo.append(("picking flat scenario",None))
868 elif (args.scenario=="two-punctures") or (args.scenario=="single-puncture"):
869 solverconstants+= "static constexpr int Scenario=2; /* Two-puncture */ \n"
870 userinfo.append(("picking black hole scenario",None))
871 else:
872 raise Exception( "Scenario " + args.scenario + " is now unknown")
873
874 for k, v in floatparams.items(): solverconstants+= "static constexpr double {} = {};\n".format("CCZ4{}".format(k), v)
875 for k, v in intparams.items(): solverconstants+= "static constexpr int {} = {};\n".format("CCZ4{}".format(k), v)
876 my_solver.add_solver_constants(solverconstants)
877
878 project.add_solver(my_solver)
879
880 build_mode = modes[args.mode]
881
882 dimensions = 3
883
884
887 floatparams.update({"domain_r":args.CCZ4domain_r})
888 dr=floatparams["domain_r"]
889 offset=[-dr, -dr, -dr]; domain_size=[2*dr, 2*dr, 2*dr]
890 msg = "domain set as "+str(offset)+" and "+str(domain_size)
891 print(msg)
892 userinfo.append((msg,None))
893
894 project.set_global_simulation_parameters(
895 dimensions, # dimensions
896 offset, domain_size,
897 args.end_time, # end time
898 args.plot_start_time, args.plot_step_size, # snapshots
899 periodic_boundary_conditions,
900 8 # plotter precision
901 )
902 userinfo.append(("plot start time: "+str(args.plot_start_time)+", plot step size: "+str(args.plot_step_size),None))
903 userinfo.append(("Terminal time: "+str(args.end_time),None))
904
905 project.set_Peano4_installation("../../..", build_mode)
906
907
910 path="./"
911 if not args.path=="./":
912 path=args.path
913 #path="/cosma5/data/durham/dc-zhan3/bbh-c5-1"
914 #path="/cosma6/data/dp004/dc-zhan3/exahype2/sbh-fv3"
915 project.set_output_path(path)
916 probe_point = [-12,-12,0.0]
917 project.add_plot_filter( probe_point,[24.0,24.0,0.01],1 )
918 if args.extension=="Psi4":
919 my_solver.select_dofs_to_print = [0,12,16,17,18,53,54,59,60]
920 elif args.extension=="adm" or args.extension=="phy-debug":
921 pass
922 else:
923 pass#my_solver.select_dofs_to_print = [0,12,16,17,18,53,54]
924
925
926 #project.set_load_balancing("toolbox::loadbalancing::strategies::RecursiveSubdivision","(new ::exahype2::LoadBalancingConfiguration(0.95,100,-1,32))" )
927 #project.set_load_balancing("toolbox::loadbalancing::strategies::RecursiveBipartition","(new ::exahype2::LoadBalancingConfiguration(0.9,1000,16))" )
928 #project.set_load_balancing("toolbox::loadbalancing::strategies::RecursiveBipartition","(new ::exahype2::LoadBalancingConfiguration(0.98))" )
929 #project.set_load_balancing("toolbox::loadbalancing::strategies::SpreadOut","(new ::exahype2::LoadBalancingConfiguration(0.95))" )
930 project.set_load_balancing("toolbox::loadbalancing::strategies::SpreadOutHierarchically","(new ::exahype2::LoadBalancingConfiguration(0.98))" )
931 #project.set_load_balancing("toolbox::loadbalancing::strategies::cascade::SpreadOut_RecursiveBipartition","(new ::exahype2::LoadBalancingConfiguration(0.98))" )
932 #project.set_load_balancing("toolbox::loadbalancing::strategies::cascade::SpreadOut_SplitOversizedTree",
933 # "new ::exahype2::LoadBalancingConfiguration(0.98, 1, 8,::exahype2::LoadBalancingConfiguration::UseNumberOfThreads), new toolbox::loadbalancing::metrics::CellCount()",)
934
935
938 if args.add_tracer==1:
939 tracer_particles = project.add_tracer( name="Tracer",attribute_count=61 )
940 init_action_set = exahype2.tracer.InsertParticlesByCoordinates(
941 particle_set=tracer_particles, coordinates=[[0,0,8.0],[0,0,-8.0],[8.0,8.0,0],[-8.0,-8.0,0]])
942 # particle_set=tracer_particles, coordinates=[[-0.1,0,0],[0.1,0,0],[0,0.1,0],[0,-0.1,0]])
943 #init_action_set = exahype2.tracer.InsertParticlesFromFile(
944 # particle_set=tracer_particles, filename="t-design.dat", scale_factor=abs(offset[0])*0.8)
945 #"Gauss_Legendre_quadrature.dat" #"t-design.dat"
946 init_action_set.priority = 0
947 project.add_action_set_to_initialisation( init_action_set )
948
949
950 # Not required anymore, as will be done by tracer routine itself. See
951 # exahype2.Project.add_tracer()
952 # -------------------------------------------------------------------
953 #project.add_action_set_to_timestepping(
954 # peano4.toolbox.particles.UpdateParallelState(
955 # particle_set=tracer_particles
956 # )
957 #)
958
959 tracer_action_set=exahype2.tracer.FiniteVolumesTracing(tracer_particles,
960 my_solver,
961 project_on_tracer_properties_kernel="::exahype2::fv::projectAllValuesOntoParticle_piecewiseLinear",
962 )
963 tracer_action_set.priority = my_solver._action_set_update_cell.priority+1
964 project.add_action_set_to_timestepping(tracer_action_set)
965 project.add_action_set_to_initialisation(tracer_action_set)
966
967 dump_action_set=exahype2.tracer.DumpTracerIntoDatabase(
968 particle_set=tracer_particles,
969 solver=my_solver,
970 filename=args.path+"Tracer-test",
971 number_of_entries_between_two_db_flushes=1000,
972 output_precision=8,
973 data_delta_between_two_snapsots=1e16,
974 position_delta_between_two_snapsots=1e16,
975 time_delta_between_two_snapsots=0.02
976 )
977 dump_action_set.priority = my_solver._action_set_update_cell.priority+2
978 project.add_action_set_to_timestepping(dump_action_set)
979
980 if args.add_tracer==2:
981 tracer_particles = project.add_tracer( name="Tracer",attribute_count=59 )
982 init_action_set = exahype2.tracer.InsertParticlesByCoordinates(
983 particle_set=tracer_particles, coordinates=[[2.0,0,0],[-2.0,0,0]])
984 init_action_set.priority = 0
985 project.add_action_set_to_initialisation( init_action_set )
986
987 # Not required anymore, as will be done by tracer routine itself. See
988 # exahype2.Project.add_tracer()
989 # -------------------------------------------------------------------
990 #project.add_action_set_to_timestepping(
991 # peano4.toolbox.particles.UpdateParallelState(
992 # particle_set=tracer_particles
993 # )
994 #)
995
996 tracer_action_set=exahype2.tracer.FiniteVolumesTracing(tracer_particles,
997 my_solver,
998 project_on_tracer_properties_kernel="::exahype2::fv::projectAllValuesOntoParticle_piecewiseLinear_explicit_Euler",
999 projection_kernel_arguments="""
1000 marker,
1001 {{PATCH_SIZE}},
1002 {{NUMBER_OF_UNKNOWNS}}+{{NUMBER_OF_AUXILIARY_VARIABLES}},
1003 fineGridCell{{SOLVER_NAME}}Q.value,
1004 fineGridCell{{SOLVER_NAME}}CellLabel.getTimeStepSize(),
1005 {17,18,19},
1006 {-1,-1,-1},
1007 *p
1008 """
1009 )
1010 tracer_action_set.priority = my_solver._action_set_compute_final_linear_combination.priority+1
1011 project.add_action_set_to_timestepping(tracer_action_set)
1012 project.add_action_set_to_initialisation(tracer_action_set)
1013
1014 dump_action_set=exahype2.tracer.DumpTracerIntoDatabase(
1015 particle_set=tracer_particles,
1016 solver=my_solver,
1017 filename=args.path+"Tracer-test",
1018 number_of_entries_between_two_db_flushes=500,
1019 output_precision=8,
1020 data_delta_between_two_snapsots=1e16,
1021 position_delta_between_two_snapsots=1e16,
1022 time_delta_between_two_snapsots=0.02
1023 )
1024 dump_action_set.priority = my_solver._action_set_compute_final_linear_combination.priority+2
1025 project.add_action_set_to_timestepping(dump_action_set)
1026
1027
1030 peano4_project = project.generate_Peano4_project(verbose=True)
1031
1032 if args.scenario=="gauge" or args.scenario=="linear" or args.scenario=="dia_gauge" or args.scenario=="flat":
1033 pass
1034 elif args.scenario=="two-punctures" or args.scenario=="single-puncture":
1035 #
1036 # There are two different things to do when we pick a scneario: We have
1037 # to configure the solver accordingly (and keep in mind that every solver
1038 # needs its own config), and then we might have to adopt the build
1039 # environment.
1040 #
1041 peano4_project.output.makefile.add_linker_flag( "-lm -lgsl -lgslcblas" )
1042 peano4_project.output.makefile.add_cpp_file( "libtwopunctures/TP_Utilities.cpp" )
1043 peano4_project.output.makefile.add_cpp_file( "libtwopunctures/TP_Parameters.cpp" )
1044 peano4_project.output.makefile.add_cpp_file( "libtwopunctures/TP_Logging.cpp" )
1045 peano4_project.output.makefile.add_cpp_file( "libtwopunctures/TwoPunctures.cpp" )
1046 peano4_project.output.makefile.add_cpp_file( "libtwopunctures/CoordTransf.cpp" )
1047 peano4_project.output.makefile.add_cpp_file( "libtwopunctures/Equations.cpp" )
1048 peano4_project.output.makefile.add_cpp_file( "libtwopunctures/FuncAndJacobian.cpp" )
1049 peano4_project.output.makefile.add_cpp_file( "libtwopunctures/Newton.cpp" )
1050 peano4_project.output.makefile.add_CXX_flag( "-DIncludeTwoPunctures" )
1051 else:
1052 raise Exception( "Scenario " + args.scenario + " is now unknown")
1053
1054 peano4_project.output.makefile.add_CXX_flag( "-DCCZ4EINSTEIN" )
1055 peano4_project.output.makefile.add_cpp_file( "InitialValues.cpp" )
1056 peano4_project.output.makefile.add_cpp_file( "CCZ4Kernels.cpp" )
1057
1058 # NOTE these lines are required to build with the fortran routines --- this will also require to uncomment some
1059 # includes etc
1060 #
1061 # peano4_project.output.makefile.add_Fortran_flag( "-DCCZ4EINSTEIN -DDim3" )
1062 # peano4_project.output.makefile.add_Fortran_flag( "-lstdc++ -fdefault-real-8 -fdefault-double-8 -cpp -std=legacy -ffree-line-length-512 -fPIC" )
1063 # peano4_project.output.makefile.add_CXX_flag( "-fPIE -DCCZ4EINSTEIN" )
1064 # peano4_project.output.makefile.add_linker_flag( "-lstdc++ -fPIC -lgfortran" )
1065 # peano4_project.output.makefile.add_linker_flag( "-lstdc++ -fPIC -L/usr/lib/x86_64-linux-gnu -lgfortran" )
1066
1067 # This might work for Intel (not tested)
1068 #peano4_project.output.makefile.add_Fortran_flag( "-r8 -cpp -auto -qopenmp-simd -O2" )
1069 #peano4_project.output.makefile.add_linker_flag( "-lstdc++ -fPIC" )
1070 # you might need -lifcore
1071
1072 # peano4_project.output.makefile.add_Fortran_module( "MainVariables.f90" )
1073
1074 # peano4_project.output.makefile.add_Fortran_files(
1075 # ["PDE.f90 ", "EinsteinConstraints.f90 ", "Properties.f90","ADMConstraints.f90"]
1076 # )
1077
1078
1079 # peano4_project.constants.export_string( "Scenario", "gaugewave-c++" )
1080 # peano4_project.constants.add_target_begin()
1081 # for k, v in floatparams.items(): peano4_project.constants.export_constexpr_with_type("CCZ4{}".format(k), eval('args.CCZ4{}'.format(k)), "double")
1082 # peano4_project.constants.add_target_end()
1083 # peano4_project.constants.add_target_begin()
1084 # for k, v in intparams.items(): peano4_project.constants.export_constexpr_with_type("CCZ4{}".format(k), eval('args.CCZ4{}'.format(k)), "int")
1085 # peano4_project.constants.add_target_end()
1086 # project.export_const_with_type("NumberOfUnknowns", 59, "int")
1087 #peano4_project.constants.export_string( "Scenario", "linearwave-c++" )
1088
1089
1090 if args.DC_flag==True:
1091 userinfo.append(("Enable double communication, make sure you are using the Second formulism.", None))
1092 additional_mesh_traversal = peano4.solversteps.Step( name = "AdditionalMeshTraversal",
1093 add_user_defined_actions=False,
1094 )
1095 project_patch_onto_faces = exahype2.solvers.rkfd.actionsets.ProjectPatchOntoFaces(my_solver)
1096 project_patch_onto_faces.guards = [ my_solver._store_cell_data_default_guard() for x in range(0,my_solver.number_of_Runge_Kutta_steps()+1) ]
1097
1098 roll_over_projected_faces = exahype2.solvers.rkfd.actionsets.RollOverUpdatedFace(my_solver,
1099 my_solver._store_face_data_default_guard(),
1100 )
1101
1102 dynamic_AMR = exahype2.solvers.rkfd.actionsets.DynamicAMR(
1103 solver=my_solver,
1104 interpolation=my_solver.interpolation,
1105 restriction=my_solver.restriction,
1106 )
1107
1108 additional_mesh_traversal.add_action_set( ComputeFirstDerivativesFD4RK(solver=my_solver,
1109 is_enclave_solver = False,
1110 ))
1111 additional_mesh_traversal.add_action_set( project_patch_onto_faces )
1112 additional_mesh_traversal.add_action_set( roll_over_projected_faces )
1113 additional_mesh_traversal.add_action_set( dynamic_AMR )
1114
1115 project.register_new_user_defined_algorithmic_step( additional_mesh_traversal )
1116 peano4_project.solversteps.add_step( additional_mesh_traversal )
1117
1118
1119 peano4_project.generate( throw_away_data_after_generation=False )
1120 peano4_project.build( make_clean_first = True )
1121
1122 #Report the application information
1123 userinfo.append(("the executable file name: "+exe, None))
1124 userinfo.append(("output directory: "+path, None))
1125 print("=========================================================")
1126 if not args.add_tracer==0:
1127 userinfo.append(("tracer output file: "+args.path+"Tracer-test", None))
1128 if len(userinfo) >0:
1129 print("The building information:")
1130 for msg, exception in userinfo:
1131 if exception is None:
1132 print(msg)
1133 else: print(msg, "Exception: {}".format(exception))
1134 print(intparams)
1135 print(floatparams)
1136 print("=========================================================")
ComputeFirstDerivatives.ComputeFirstDerivativesFD4RK
Explicit reconstruction of first derivatives for FD4 discretisation.
Definition ComputeFirstDerivatives.py:55
ccz4_archived_24_01_19.CCZ4Solver
Definition ccz4_archived_24_01_19.py:117
ccz4_archived_24_01_19.CCZ4Solver._source_term_implementation
_source_term_implementation
Definition ccz4_archived_24_01_19.py:257
ccz4_archived_24_01_19.CCZ4Solver._domain_r
_domain_r
Definition ccz4_archived_24_01_19.py:242
ccz4_archived_24_01_19.CCZ4Solver.postprocess_updated_patch
postprocess_updated_patch
Definition ccz4_archived_24_01_19.py:268
ccz4_archived_24_01_19.CCZ4Solver._patch_size
_patch_size
Definition ccz4_archived_24_01_19.py:241
ccz4_archived_24_01_19.CCZ4Solver._flux_implementation
_flux_implementation
Definition ccz4_archived_24_01_19.py:255
ccz4_archived_24_01_19.CCZ4Solver.__init__
__init__(self, name, patch_size, min_volume_h, max_volume_h, cfl, domain_r, KOSig)
Definition ccz4_archived_24_01_19.py:118
ccz4_archived_24_01_19.CCZ4Solver._my_user_includes
_my_user_includes
Definition ccz4_archived_24_01_19.py:142
ccz4_archived_24_01_19.CCZ4Solver._fused_compute_kernel_call_cpu
_fused_compute_kernel_call_cpu
Definition ccz4_archived_24_01_19.py:254
ccz4_archived_24_01_19.CCZ4Solver._ncp_implementation
_ncp_implementation
Definition ccz4_archived_24_01_19.py:256
exahype2.Project.Project
ExaHyPE 2 project.
Definition Project.py:14
exahype2.solvers.LagrangeBasisWithDiagonalMassMatrix.GaussLegendreBasis
The Gauss-Legendre Basis is by construction the only basis which yields diagonal mass matrices.
Definition LagrangeBasisWithDiagonalMassMatrix.py:609
exahype2.solvers.elliptic.AMRMarker.AMRMarker
A very simple Poisson equation system solver to identify AMR transitions.
Definition AMRMarker.py:24
exahype2.solvers.rkdg.rusanov.GlobalAdaptiveTimeStep.GlobalAdaptiveTimeStep
RKDG solver with Rusanov Riemann solver employing global adaptive time stepping.
Definition GlobalAdaptiveTimeStep.py:23
exahype2.solvers.rkfd.actionsets.DynamicAMR.DynamicAMR
The handling of (dynamically) adaptive meshes for finite differences.
Definition DynamicAMR.py:11
exahype2.solvers.rkfd.actionsets.ProjectPatchOntoFaces.ProjectPatchOntoFaces
Project patch data onto faces, so the faces hold valid data which can feed into the subsequent Runge-...
Definition ProjectPatchOntoFaces.py:13
exahype2.tracer.FiniteVolumesTracing.FiniteVolumesTracing
Particle tracing over the Finite Volumes solver.
Definition FiniteVolumesTracing.py:12
peano4.toolbox.blockstructured.DynamicAMR.DynamicAMR
This class assumes that you have an 2MxNxN patch on your faces.
Definition DynamicAMR.py:9
ccz4_archived_24_01_19.add_Psi4W
add_Psi4W(self)
add psi4 writer
Definition ccz4_archived_24_01_19.py:376
ccz4_archived_24_01_19.CoordListFromFile
CoordListFromFile(file, scale)
Definition ccz4_archived_24_01_19.py:18
ccz4_archived_24_01_19.add_adm_constriants
add_adm_constriants(self)
@add adm constriants (Hamilton and Momentum)
Definition ccz4_archived_24_01_19.py:325
ccz4_archived_24_01_19.add_puncture_tracker
add_puncture_tracker(self)
@add
Definition ccz4_archived_24_01_19.py:485
ccz4_archived_24_01_19.get_user_action_set_includes
get_user_action_set_includes(self)
We take this routine to add a few additional include statements.
Definition ccz4_archived_24_01_19.py:319
ccz4_archived_24_01_19.add_intermediate_quantities
add_intermediate_quantities(self)
@add adm intermediate quantites
Definition ccz4_archived_24_01_19.py:429
ccz4_archived_24_01_19.create_action_sets
create_action_sets(self)
Definition ccz4_archived_24_01_19.py:313
exahype2.solvers.rkfd.fd4.kernels.create_compute_kernel_for_FD4
create_compute_kernel_for_FD4(flux_implementation, ncp_implementation, source_implementation, compute_max_eigenvalue_of_next_time_step, solver_variant, kernel_variant, KOSigma)
I return only the unqualified function call, i.e.
Definition kernels.py:20