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