Peano
Loading...
Searching...
No Matches
SBH.py
Go to the documentation of this file.
1# This file is part of the Peano project. For conditions of distribution and
2# use, please see the copyright notice at www.peano-framework.org
3import os
4import argparse
5
6import peano4
7import exahype2
8import dastgen2
9
11import numpy as np
12
13import jinja2
14
15
17
18from enum import Enum
19
20# See comments in README.dox
21# export PYTHONPATH=../../../../python
22# export PYTHONPATH=$PYTHONPATH:../../../../applications/exahype2/ccz4
23from CCZ4Solver import CCZ4Solver_FV_GlobalAdaptiveTimeStepWithEnclaveTasking
24from CCZ4Solver import CCZ4Solver_FD4_GlobalAdaptiveTimeStepWithEnclaveTasking
25from CCZ4Solver import (
26 CCZ4Solver_FD4_SecondOrderFormulation_GlobalAdaptiveTimeStepWithEnclaveTasking,
27)
28
29from CCZ4Solver import CCZ4Solver_FV_GlobalAdaptiveTimeStep
30from CCZ4Solver import CCZ4Solver_FD4_GlobalAdaptiveTimeStep
31
32
33# for the coupling
34from exahype2.solvers.fv.actionsets.AbstractFVActionSet import AbstractFVActionSet
35
36BlackHoleRegion = 0.1
37
38
40 NONE = "::peano4::utils::LoopPlacement::Serial"
41 PARALLEL_FOR = "::peano4::utils::LoopPlacement::Nested"
42 SUBTASKS = "::peano4::utils::LoopPlacement::SpreadOut"
43
44
46 """!
47
48 Construct the Finite Volume (limiter) scheme
49
50 We assume that the underlying Finite Differences scheme has a patch
51 size of 6x6x6. To make the Finite Volume scheme's time stepping (and
52 accuracy) match this patch size, we have to employ a 16 times finer
53 mesh.
54
55 It is interesting to see that the limiter does not really have a min
56 and max mesh size. The point is that the higher order solver dictates
57 the mesh structure, and we then follow this structure with the
58 Finite Volume scheme.
59
60 """
61
63 self,
64 name,
65 patch_size,
66 amend_priorities: bool,
67 parallelisation_of_kernels: KernelParallelisation,
68 ):
69 """!
70
71 Construct the limiter
72
73 patch_size: Integer
74 Pass in the patch size of the FD4 scheme or, if you are using RKDG,
75 hand in the number 1. The Finite Volume patch then will be 16 times
76 finer.
77
78 """
79 SubPatchSize = int(patch_size * 4 * 4)
80
81 super(Limiter, self).__init__(
82 name=name + "_FV",
83 patch_size=SubPatchSize,
84 min_volume_h=1.0 / 65536,
85 max_volume_h=65536,
86 pde_terms_without_state=True,
87 )
88 self.double_constants["BlackHoleFVRegion"] = BlackHoleRegion
89
91#include "toolbox/blockstructured/Restriction.h"
92#include "toolbox/blockstructured/Interpolation.h"
93"""
94
95 if amend_priorities:
96 self.enclave_task_priorityenclave_task_priority = "tarch::multicore::Task::DefaultPriority+1"
97
98 # self._fused_compute_kernel_call_gpu = exahype2.solvers.fv.rusanov.kernels.create_compute_Riemann_kernel_for_Rusanov(
99 # self._flux_implementation, self._ncp_implementation, self._source_term_implementation,
100 # compute_max_eigenvalue_of_next_time_step=True,
101 # solver_variant = exahype2.solvers.fv.rusanov.kernels.SolverVariant.Accelerator,
102 # kernel_variant = exahype2.solvers.fv.rusanov.kernels.KernelVariant.GenericAoS
103 # )
104
106 initial_conditions=exahype2.solvers.PDETerms.User_Defined_Implementation,
107 refinement_criterion=exahype2.solvers.PDETerms.Empty_Implementation,
108 boundary_conditions=exahype2.solvers.PDETerms.Empty_Implementation,
109 )
110
111 if parallelisation_of_kernels == KernelParallelisation.SUBTASKS:
112 # self._compute_kernel_call = exahype2.solvers.fv.rusanov.kernels.create_compute_Riemann_kernel_for_Rusanov(
113 # self._flux_implementation, self._ncp_implementation, self._source_term_implementation,
114 # compute_max_eigenvalue_of_next_time_step=True,
115 # solver_variant = exahype2.solvers.fv.rusanov.kernels.SolverVariant.WithVirtualFunctions,
116 # kernel_variant = exahype2.solvers.fv.rusanov.kernels.KernelVariant.PatchWiseAoSHeap
117 # )
118
123 compute_max_eigenvalue_of_next_time_step=True,
124 solver_variant=exahype2.solvers.fv.rusanov.kernels.SolverVariant.Multicore,
125 kernel_variant=exahype2.solvers.fv.rusanov.kernels.KernelVariant.PatchWiseAoSHeap,
126 )
127
130
132 """!
133
134 Mask out exterior cells
135
136 """
137 return (
138 """("""
139 + super(Limiter, self)._store_cell_data_default_guard()
140 + " and repositories::instanceOf"
141 + self._name_name_name
142 + ".isCellOverlappingWithBHImpactArea(marker) )"
143 )
144
146 return (
147 """("""
148 + super(Limiter, self)._load_cell_data_default_guard()
149 + " and repositories::instanceOf"
150 + self._name_name_name
151 + ".isCellOverlappingWithBHImpactArea(marker) )"
152 )
153
155 return (
156 """("""
158 + " and repositories::instanceOf"
159 + self._name_name_name
160 + ".isCellOverlappingWithBHImpactArea(marker) )"
161 )
162
164 return "{} and repositories::instanceOf{}.isOneAdjacentCellOverlappingWithBHImpactArea(marker)".format(
167 )
168
170 return "({} and repositories::instanceOf{}.areBothAdjacentCellsOverlappingWithBHImpactArea(marker))".format(
171 super(Limiter, self)._store_face_data_default_guard(), self._name_name_name
172 )
173
175 return "({} and repositories::instanceOf{}.areBothAdjacentCellsOverlappingWithBHImpactArea(marker))".format(
176 super(Limiter, self)._load_face_data_default_guard(), self._name_name_name
177 )
178
180 """!
181
182 Not really a lot of things to do here. The only exception that is
183 really important is that we have to ensure that we only solve stuff
184 inside the local domain of the FV. By default, ExaHyPE 2 solves the
185 PDE everywhere. If data is not stored persistently or loaded from
186 the persistent stacks, it still solves, as it then would assume that
187 such data arises from dynamic AMR. In this particular case, we have
188 to really mask out certain subdomains.
189
190 It is not just a nice optimisation to do so. It is absolutely key,
191 as the application of the compute kernel on garbage would mean that
192 we end up with invalid eigenvalues.
193
194 """
195 super(Limiter, self).create_action_sets()
196
197 self._action_set_update_cell_action_set_update_cell.guard = "({} and repositories::instanceOf{}.isCellOverlappingWithBHImpactArea(marker))".format(
199 )
200 self._action_set_merge_enclave_task_outcome.guard = "({} and repositories::instanceOf{}.isCellOverlappingWithBHImpactArea(marker))".format(
202 )
203
204
206 """!
207
208 4th order Finite Differences solver with a limiter
209
210 The FD scheme is our primary solver, i.e. the one we wanna use (almost)
211 everywhere. Hence, we have to be sure that we use the right boundary
212 conditions. In our case, that should be Sommerfeld ones. As we work
213 with a limited FD4 solver, we have to ensure that we get the
214 @ref benchmarks_exahype2_ccz4_single_black_hole "solver coupling" right.
215
216 The primary solver is a plain CCZ4 FD4 solver with enclave tasking. There
217 are some modification though:
218
219 1. We ensure that we use Sommerfeld boundary conditions. This means that
220 we have to replace the whole generic boundary treatment with a bespoke
221 action set.
222 2. The coupling has to be injected. We model the coupling as separate
223 (enclave) tasks to avoid that we totally serialise the solver steps
224 whenever we encounter a coupling.
225
226 """
227
229 self,
230 name,
231 patch_size,
232 min_cell_size,
233 max_cell_size,
234 parallelisation_of_interpolation: KernelParallelisation,
235 parallelisation_of_kernels: KernelParallelisation,
236 interpolation_method,
237 ):
239 self._parallelisation_of_interpolation = parallelisation_of_interpolation
240 self.interpolation_method = interpolation_method
241
242 CCZ4Solver_FD4_GlobalAdaptiveTimeStepWithEnclaveTasking.__init__(
243 self,
244 name=name + "_FD4",
245 patch_size=patch_size,
246 rk_order=1,
247 min_meshcell_h=min_cell_size / patch_size,
248 max_meshcell_h=max_cell_size / patch_size,
249 pde_terms_without_state=False,
250 )
251
252 self._user_action_set_includes += """
253#include "toolbox/blockstructured/Projection.h"
254#include "toolbox/blockstructured/Restriction.h"
255#include "toolbox/blockstructured/Interpolation.h"
256#include "toolbox/blockstructured/Copy.h"
257"""
258 if self.interpolation_method == "matrix":
260 """
261#include "matrixdata/Interpolator_"""
262 + str(patch_size)
263 + """_"""
264 + str(patch_size * 16)
265 + """.h"
266"""
267 )
268 self.double_constants["BlackHoleFVRegion"] = BlackHoleRegion
269
271 initial_conditions=exahype2.solvers.PDETerms.User_Defined_Implementation,
272 refinement_criterion=exahype2.solvers.PDETerms.User_Defined_Implementation,
273 boundary_conditions=exahype2.solvers.PDETerms.None_Implementation,
274 )
276
277 if parallelisation_of_kernels == KernelParallelisation.SUBTASKS:
278 print(
279 "WARNING: Unfortunately, we have not yet written parallelised kernels for FD4"
280 )
281 pass
282
284
285 if self.interpolation_method == "matrix":
286 solver_matrix_interpolator = exahype2.solvers.SolverMatrixInterpolator(
287 patch_size, patch_size * 16, 3, 1
288 )
289 solver_matrix_interpolator.generateData()
290
292 """!
293
294 Tailor action set behaviour
295
296 We first make a few additional cells skeleton cells. The rationale
297 behind additional skeletons is given in the @ref benchmarks_exahype2_ccz4_single_black_hole "generic overview".
298 Given the first remark there on FD4-FV coupling, one would be tempted
299 to use the predicate
300
301 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
302 self._action_set_update_cell.additional_skeleton_guard = " " "(
303 repositories::instanceOf" " " + self._name_without_FD4_extension + " " "_FV.isCellOverlappingWithBHImpactArea(marker)
304 and
305 not repositories::instanceOf" " " + self._name_without_FD4_extension + " " "_FV.areAllFaceConnectedCellsOverlappingWithBHImpactArea(marker)
306 )
307 " " "
308 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
309
310 Once we study the other items (notably the fourth), we see that it is
311 reasonable to make all the overlap region a skeleton within the FD4
312 solver.
313
314 """
315 super(FD4SolverWithLimiter, self).create_action_sets()
316
317 self._action_set_handle_boundary.TemplateHandleBoundary_KernelCalls = (
318 """
319 double Qinf[59]={1.0, 0.0, 0.0, 1.0, 0.0, 1.0, //q0-5
320 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, //q6-11
321 0.0, 0.0, 0.0, 0.0, //q12-15
322 1.0, 0.0, 0.0, 0.0, //q16-19
323 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, //q20-25
324 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, //q26-34
325 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
326 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, //q35-52
327 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 //q53-58
328 }; //approximate background solution at infinity
329 ::exahype2::fd::applySommerfeldConditions(
330 [&](
331 const double * __restrict__ Q,
332 const tarch::la::Vector<Dimensions,double>& faceCentre,
333 const tarch::la::Vector<Dimensions,double>& gridCellH,
334 double t,
335 double dt,
336 int normal
337 ) -> double {
338 return repositories::{{SOLVER_INSTANCE}}.maxEigenvalue( Q, faceCentre, gridCellH, t, dt, normal );
339 },
340 [&](
341 double * __restrict__ Q,
342 const tarch::la::Vector<Dimensions,double>& faceCentre,
343 const tarch::la::Vector<Dimensions,double>& gridCellH
344 ) -> void {
345 for (int i=0; i<"""
347 + """; i++) {
348 Q[i] = 0.0;
349 }
350 Q[0] = 1.0; Q[3] = 1.0; Q[5] = 1.0;
351 Q[16] = 0.95; Q[54] = 0.95;
352 },
353 marker.x(),
354 marker.h(),
355 {{FACE_METADATA_ACCESSOR}}.getOldTimeStamp(marker.getSelectedFaceNumber()<Dimensions ? 1 : 0),
356 repositories::{{SOLVER_INSTANCE}}.getMinTimeStepSize(),
357 {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}},
358 {{OVERLAP}},
359 {{NUMBER_OF_UNKNOWNS}},
360 {{NUMBER_OF_AUXILIARY_VARIABLES}},
361 marker.getSelectedFaceNumber(),
362 {0.0, 0.0, 0.0},
363 fineGridFace{{UNKNOWN_IDENTIFIER}}Old.value,
364 fineGridFace{{UNKNOWN_IDENTIFIER}}New.value
365 );
366 """
367 )
368
370 solver=self,
371 enclave_task_cell_label="fineGridCell"
373 + "_FVCellLabel",
374 compute_kernel_implementation="""
375const int sizeOfPatch = (repositories::instanceOf"""
377 + """_FV.NumberOfFiniteVolumesPerAxisPerPatch+2)
378 * (repositories::instanceOf"""
380 + """_FV.NumberOfFiniteVolumesPerAxisPerPatch+2)
381 * (repositories::instanceOf"""
383 + """_FV.NumberOfFiniteVolumesPerAxisPerPatch+2)
384 * (repositories::instanceOf"""
386 + """_FV.NumberOfUnknowns + repositories::instanceOf"""
388 + """_FV.NumberOfAuxiliaryVariables);
389const int sizeOfFace = 2
390 * (repositories::instanceOf"""
392 + """_FV.NumberOfFiniteVolumesPerAxisPerPatch+2)
393 * (repositories::instanceOf"""
395 + """_FV.NumberOfFiniteVolumesPerAxisPerPatch+2)
396 * (repositories::instanceOf"""
398 + """_FV.NumberOfUnknowns + repositories::instanceOf"""
400 + """_FV.NumberOfAuxiliaryVariables);
401
402double* interpolatedFVDataWithHalo = ::tarch::allocateMemory<double>(sizeOfPatch, ::tarch::MemoryLocation::Heap);
403
404bool faceIsReal0 = repositories::instanceOf"""
406 + """_FV.areBothAdjacentCellsOverlappingWithBHImpactArea(marker,0);
407bool faceIsReal1 = repositories::instanceOf"""
409 + """_FV.areBothAdjacentCellsOverlappingWithBHImpactArea(marker,1);
410bool faceIsReal2 = repositories::instanceOf"""
412 + """_FV.areBothAdjacentCellsOverlappingWithBHImpactArea(marker,2);
413bool faceIsReal3 = repositories::instanceOf"""
415 + """_FV.areBothAdjacentCellsOverlappingWithBHImpactArea(marker,3);
416bool faceIsReal4 = repositories::instanceOf"""
418 + """_FV.areBothAdjacentCellsOverlappingWithBHImpactArea(marker,4);
419bool faceIsReal5 = repositories::instanceOf"""
421 + """_FV.areBothAdjacentCellsOverlappingWithBHImpactArea(marker,5);
422
423double* realFace0 = faceIsReal0 ? fineGridFaces"""
425 + """_FVQNew(0).value : nullptr;
426double* realFace1 = faceIsReal1 ? fineGridFaces"""
428 + """_FVQNew(1).value : nullptr;
429double* realFace2 = faceIsReal2 ? fineGridFaces"""
431 + """_FVQNew(2).value : nullptr;
432double* realFace3 = faceIsReal3 ? fineGridFaces"""
434 + """_FVQNew(3).value : nullptr;
435double* realFace4 = faceIsReal4 ? fineGridFaces"""
437 + """_FVQNew(4).value : nullptr;
438double* realFace5 = faceIsReal5 ? fineGridFaces"""
440 + """_FVQNew(5).value : nullptr;
441
442double* dummyFace0 = faceIsReal0 ? nullptr : ::tarch::allocateMemory<double>(sizeOfFace, ::tarch::MemoryLocation::Heap);
443double* dummyFace1 = faceIsReal1 ? nullptr : ::tarch::allocateMemory<double>(sizeOfFace, ::tarch::MemoryLocation::Heap);
444double* dummyFace2 = faceIsReal2 ? nullptr : ::tarch::allocateMemory<double>(sizeOfFace, ::tarch::MemoryLocation::Heap);
445double* dummyFace3 = faceIsReal3 ? nullptr : ::tarch::allocateMemory<double>(sizeOfFace, ::tarch::MemoryLocation::Heap);
446double* dummyFace4 = faceIsReal4 ? nullptr : ::tarch::allocateMemory<double>(sizeOfFace, ::tarch::MemoryLocation::Heap);
447double* dummyFace5 = faceIsReal5 ? nullptr : ::tarch::allocateMemory<double>(sizeOfFace, ::tarch::MemoryLocation::Heap);
448
449::toolbox::blockstructured::interpolateCellDataAssociatedToVolumesIntoOverlappingCell_"""
451 + """(
452 repositories::instanceOf"""
454 + """_FD4.NumberOfGridCellsPerPatchPerAxis,
455 repositories::instanceOf"""
457 + """_FV.NumberOfFiniteVolumesPerAxisPerPatch,
458 3, // halo in FD4
459 1, // halo in FV
460 repositories::instanceOf"""
462 + """_FV.NumberOfUnknowns + repositories::instanceOf"""
464 + """_FV.NumberOfAuxiliaryVariables,
465 """
466 + (
467 """InterpolationMatrixData::Data,
468 InterpolationMatrixData::Indices,
469 InterpolationMatrixData::Indptr,"""
470 if self.interpolation_method == "matrix"
471 else """"""
472 )
473 + """
474 oldQWithHalo,
475 interpolatedFVDataWithHalo,
476 """
477 + str(self._parallelisation_of_interpolation.value)
478 + """
479);
480
481::toolbox::blockstructured::projectPatchHaloOntoFaces(
482 repositories::instanceOf"""
484 + """_FV.NumberOfFiniteVolumesPerAxisPerPatch,
485 1, // int haloSize,
486 repositories::instanceOf"""
488 + """_FV.NumberOfUnknowns,
489 repositories::instanceOf"""
491 + """_FV.NumberOfAuxiliaryVariables,
492 interpolatedFVDataWithHalo,
493 faceIsReal0 ? realFace0 : dummyFace0,
494 faceIsReal1 ? realFace1 : dummyFace1,
495 faceIsReal2 ? realFace2 : dummyFace2,
496 faceIsReal3 ? realFace3 : dummyFace3,
497 faceIsReal4 ? realFace4 : dummyFace4,
498 faceIsReal5 ? realFace5 : dummyFace5
499);
500
501::tarch::freeMemory(interpolatedFVDataWithHalo, tarch::MemoryLocation::Heap );
502if (dummyFace0!=nullptr) ::tarch::freeMemory(dummyFace0, tarch::MemoryLocation::Heap );
503if (dummyFace1!=nullptr) ::tarch::freeMemory(dummyFace1, tarch::MemoryLocation::Heap );
504if (dummyFace2!=nullptr) ::tarch::freeMemory(dummyFace2, tarch::MemoryLocation::Heap );
505if (dummyFace3!=nullptr) ::tarch::freeMemory(dummyFace3, tarch::MemoryLocation::Heap );
506if (dummyFace4!=nullptr) ::tarch::freeMemory(dummyFace4, tarch::MemoryLocation::Heap );
507if (dummyFace5!=nullptr) ::tarch::freeMemory(dummyFace5, tarch::MemoryLocation::Heap );
508""",
509 )
510
512 "("
514 + """
515 and
516 repositories::instanceOf"""
518 + """_FV.isCellOverlappingWithBHImpactArea(marker)
519 and
520 not repositories::instanceOf"""
522 + """_FV.areAllFaceConnectedCellsOverlappingWithBHImpactArea(marker)
523 and
524 not marker.willBeRefined()
525 )"""
526 )
527
529 self,
530 """
531if (
532 repositories::instanceOf{0}_FV.isCellOverlappingWithBHImpactArea(marker)
533 and
534 not marker.willBeRefined()
535) {{
536 logTraceIn( "touchCellFirstTime(...)-inject" );
537
538 repositories::instanceOf{0}_FV.incNumberOfPatches();
539
540 //if ( repositories::instanceOf{0}_FV.areAllFaceConnectedCellsOverlappingWithBHImpactArea(marker) ) {{
541 ::toolbox::blockstructured::restrictCellIntoOverlappingCell_inject(
542 repositories::instanceOf{0}_FV.NumberOfFiniteVolumesPerAxisPerPatch,
543 repositories::instanceOf{0}_FD4.NumberOfGridCellsPerPatchPerAxis,
544 repositories::instanceOf{0}_FV.NumberOfUnknowns + repositories::instanceOf{0}_FV.NumberOfAuxiliaryVariables,
545 fineGridCell{0}_FVQ.value,
546 fineGridCell{0}_FD4Q.value
547 );
548 /*}}
549 else {{
550 ::toolbox::blockstructured::restrictCellIntoOverlappingCell_inject_and_average(
551 repositories::instanceOf{0}_FV.NumberOfFiniteVolumesPerAxisPerPatch,
552 repositories::instanceOf{0}_FD4.NumberOfGridCellsPerPatchPerAxis,
553 repositories::instanceOf{0}_FV.NumberOfUnknowns + repositories::instanceOf{0}_FV.NumberOfAuxiliaryVariables,
554 fineGridCell{0}_FVQ.value,
555 fineGridCell{0}_FD4Q.value
556 );
557 }}*/
558
559 logTraceOut( "touchCellFirstTime(...)-inject" );
560}}
561""".format(
563 ),
564 )
565
566
568 """!
569
570 Construct 4th order Finite Differences solver without a limiter
571
572 """
573
575 self,
576 name,
577 patch_size,
578 min_cell_size,
579 max_cell_size,
580 ):
582
583 CCZ4Solver_FD4_GlobalAdaptiveTimeStepWithEnclaveTasking.__init__(
584 self,
585 name=name + "_FD4",
586 patch_size=patch_size,
587 rk_order=1,
588 min_meshcell_h=min_cell_size / patch_size,
589 max_meshcell_h=max_cell_size / patch_size,
590 pde_terms_without_state=False,
591 )
592
593 self._user_action_set_includes += """
594#include "toolbox/blockstructured/Projection.h"
595#include "toolbox/blockstructured/Restriction.h"
596#include "toolbox/blockstructured/Interpolation.h"
597#include "toolbox/blockstructured/Copy.h"
598"""
599
601 self.double_constants["BlackHoleFVRegion"] = BlackHoleRegion
602
604 initial_conditions=exahype2.solvers.PDETerms.User_Defined_Implementation,
605 refinement_criterion=exahype2.solvers.PDETerms.User_Defined_Implementation,
606 boundary_conditions=exahype2.solvers.PDETerms.None_Implementation,
607 )
608
611
613 """!
614
615 Tailor action set behaviour
616
617 We first make a few additional cells skeleton cells. The rationale
618 behind additional skeletons is given in the @ref benchmarks_exahype2_ccz4_single_black_hole "generic overview".
619 Given the first remark there on FD4-FV coupling, one would be tempted
620 to use the predicate
621
622 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
623 self._action_set_update_cell.additional_skeleton_guard = " " "(
624 repositories::instanceOf" " " + self._name_without_FD4_extension + " " "_FV.isCellOverlappingWithBHImpactArea(marker)
625 and
626 not repositories::instanceOf" " " + self._name_without_FD4_extension + " " "_FV.areAllFaceConnectedCellsOverlappingWithBHImpactArea(marker)
627 )
628 " " "
629 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
630
631 Once we study the other items (notably the fourth), we see that it is
632 reasonable to make all the overlap region a skeleton within the FD4
633 solver.
634
635 """
636 super(FD4SolverWithoutLimiter, self).create_action_sets()
637
638 self._action_set_handle_boundary.TemplateHandleBoundary_KernelCalls = (
639 """
640 double Qinf[59]={1.0, 0.0, 0.0, 1.0, 0.0, 1.0, //q0-5
641 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, //q6-11
642 0.0, 0.0, 0.0, 0.0, //q12-15
643 1.0, 0.0, 0.0, 0.0, //q16-19
644 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, //q20-25
645 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, //q26-34
646 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
647 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, //q35-52
648 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 //q53-58
649 }; //approximate background solution at infinity
650 ::exahype2::fd::applySommerfeldConditions(
651 [&](
652 const double * __restrict__ Q,
653 const tarch::la::Vector<Dimensions,double>& faceCentre,
654 const tarch::la::Vector<Dimensions,double>& gridCellH,
655 double t,
656 double dt,
657 int normal
658 ) -> double {
659 return repositories::{{SOLVER_INSTANCE}}.maxEigenvalue( Q, faceCentre, gridCellH, t, dt, normal );
660 },
661 [&](
662 double * __restrict__ Q,
663 const tarch::la::Vector<Dimensions,double>& faceCentre,
664 const tarch::la::Vector<Dimensions,double>& gridCellH
665 ) -> void {
666 for (int i=0; i<"""
668 + """; i++) {
669 Q[i] = 0.0;
670 }
671 Q[0] = 1.0; Q[3] = 1.0; Q[5] = 1.0;
672 Q[16] = 0.95; Q[54] = 0.95;
673 },
674 marker.x(),
675 marker.h(),
676 {{FACE_METADATA_ACCESSOR}}.getOldTimeStamp(marker.getSelectedFaceNumber()<Dimensions ? 1 : 0),
677 repositories::{{SOLVER_INSTANCE}}.getMinTimeStepSize(),
678 {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}},
679 {{OVERLAP}},
680 {{NUMBER_OF_UNKNOWNS}},
681 {{NUMBER_OF_AUXILIARY_VARIABLES}},
682 marker.getSelectedFaceNumber(),
683 {0.0, 0.0, 0.0},
684 fineGridFace{{UNKNOWN_IDENTIFIER}}Old.value,
685 fineGridFace{{UNKNOWN_IDENTIFIER}}New.value
686 );
687 """
688 )
689
690
692 """!
693
694 A finite volume solver
695
696 This solver is not appropriate to simulate black holes as a stand-alone
697 solver, as it is way too diffusive. If you use it without another scheme,
698 you typically see the black hole disappear after a brief period. So we
699 have it in here merely for performance tests.
700
701 """
702
704 self,
705 name,
706 patch_size,
707 min_cell_size,
708 max_cell_size,
709 ):
710 """!
711
712 Construct the Finite Volume solver
713
714 @param patch_size: Integer
715 Defines how big the individual patches are. If you pass in 10, each
716 Finite Volume patch will have the dimensions 10x10x10.
717 @param min_cell_size: Float
718 This parameter refers to the cell size, i.e. the size of a whole
719 patch. We use this one here, to make the signature the same as for
720 the FD and DG solver variants. The superclass constructor argues
721 over finite volume sizes, and we hence have to recalibrate this
722 parameter with patch_size.
723
724 """
725
726 super(FVSolver, self).__init__(
727 name=name + "_FV",
728 patch_size=patch_size,
729 min_volume_h=min_cell_size / patch_size,
730 max_volume_h=max_cell_size / patch_size,
731 pde_terms_without_state=True,
732 )
733
734 self._user_action_set_includes += """
735#include "toolbox/blockstructured/Projection.h"
736#include "toolbox/blockstructured/Restriction.h"
737#include "toolbox/blockstructured/Interpolation.h"
738#include "toolbox/blockstructured/Copy.h"
739"""
740
742 self.double_constants["BlackHoleFVRegion"] = BlackHoleRegion
743
745 initial_conditions=exahype2.solvers.PDETerms.User_Defined_Implementation,
746 refinement_criterion=exahype2.solvers.PDETerms.Empty_Implementation,
747 boundary_conditions=exahype2.solvers.PDETerms.Empty_Implementation,
748 )
749
752
754 """!
755
756 Tailor action set behaviour
757
758 We first make a few additional cells skeleton cells. The rationale
759 behind additional skeletons is given in the @ref benchmarks_exahype2_ccz4_single_black_hole "generic overview".
760 Given the first remark there on FD4-FV coupling, one would be tempted
761 to use the predicate
762
763 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
764 self._action_set_update_cell.additional_skeleton_guard = " " "(
765 repositories::instanceOf" " " + self._name_without_FD4_extension + " " "_FV.isCellOverlappingWithBHImpactArea(marker)
766 and
767 not repositories::instanceOf" " " + self._name_without_FD4_extension + " " "_FV.areAllFaceConnectedCellsOverlappingWithBHImpactArea(marker)
768 )
769 " " "
770 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
771
772 Once we study the other items (notably the fourth), we see that it is
773 reasonable to make all the overlap region a skeleton within the FD4
774 solver.
775
776 """
777 super(FVSolver, self).create_action_sets()
778
779
781 """!
782
783 Update preconfigured solver parameters
784
785 The default parameters of CCZ4 are tailored towards gauge waves or similar.
786 For the single black hole, two parameters have to be changed. That's bs and
787 sk which both have to be 1.0.
788
789 """
790 solver.double_constants["CCZ4bs"] = 1.0
791 solver.double_constants["CCZ4sk"] = 1.0
add_all_solver_constants(self)
Add domain-specific constants.
CCZ4 solver using fourth-order finite differences and global adaptive time stepping incl enclave task...
CCZ4 solver using finite volumes and global adaptive time stepping incl enclave tasking.
4th order Finite Differences solver with a limiter
Definition SBH.py:205
__init__(self, name, patch_size, min_cell_size, max_cell_size, KernelParallelisation parallelisation_of_interpolation, KernelParallelisation parallelisation_of_kernels, interpolation_method)
Constructor.
Definition SBH.py:237
create_action_sets(self)
Tailor action set behaviour.
Definition SBH.py:291
Construct 4th order Finite Differences solver without a limiter.
Definition SBH.py:567
create_action_sets(self)
Tailor action set behaviour.
Definition SBH.py:612
__init__(self, name, patch_size, min_cell_size, max_cell_size)
Constructor.
Definition SBH.py:580
A finite volume solver.
Definition SBH.py:691
__init__(self, name, patch_size, min_cell_size, max_cell_size)
Construct the Finite Volume solver.
Definition SBH.py:709
create_action_sets(self)
Tailor action set behaviour.
Definition SBH.py:753
Construct the Finite Volume (limiter) scheme.
Definition SBH.py:45
_flux_implementation
Definition SBH.py:120
create_action_sets(self)
Not really a lot of things to do here.
Definition SBH.py:179
_source_term_implementation
Definition SBH.py:122
_fused_compute_kernel_call_cpu
Definition SBH.py:119
_provide_face_data_to_compute_kernels_default_guard(self)
Definition SBH.py:163
_store_cell_data_default_guard(self)
Mask out exterior cells.
Definition SBH.py:131
_ncp_implementation
Definition SBH.py:121
_provide_cell_data_to_compute_kernels_default_guard(self)
Default logic when to create cell data or not.
Definition SBH.py:154
_store_face_data_default_guard(self)
Extend the guard via ands only.
Definition SBH.py:169
_load_cell_data_default_guard(self)
Extend the guard via ands only.
Definition SBH.py:145
__init__(self, name, patch_size, bool amend_priorities, KernelParallelisation parallelisation_of_kernels)
Construct the limiter.
Definition SBH.py:68
enclave_task_priority
Definition SBH.py:96
_load_face_data_default_guard(self)
Extend the guard via ands only.
Definition SBH.py:174
set_implementation(self, boundary_conditions, refinement_criterion, initial_conditions, memory_location, use_split_loop, additional_action_set_includes, additional_user_includes)
If you pass in User_Defined, then the generator will create C++ stubs that you have to befill manuall...
set_implementation(self, flux=None, ncp=None, eigenvalues=None, boundary_conditions=None, refinement_criterion=None, initial_conditions=None, source_term=None, memory_location=None, use_split_loop=False, additional_action_set_includes="", additional_user_includes="")
If you pass in User_Defined, then the generator will create C++ stubs that you have to befill manuall...
set_implementation(self, flux, ncp, source_term, eigenvalues, boundary_conditions, refinement_criterion, initial_conditions, memory_location, additional_action_set_includes, additional_user_includes)
If you pass in User_Defined, then the generator will create C++ stubs that you have to befill manuall...
set_implementation(self, flux=None, ncp=None, source_term=None, eigenvalues=None, boundary_conditions=None, refinement_criterion=None, initial_conditions=None, memory_location=None, additional_action_set_includes="", additional_user_includes="", KOSigma=None)
If you pass in User_Defined, then the generator will create C++ stubs that you have to befill manuall...
update_solver_parameters_for_single_black_hole(solver)
Update preconfigured solver parameters.
Definition SBH.py:780
create_compute_Riemann_kernel_for_Rusanov(flux_implementation, ncp_implementation, source_implementation, compute_max_eigenvalue_of_next_time_step, SolverVariant solver_variant, KernelVariant kernel_variant)
Return only the unqualified function call, i.e., without any namespaces.
Definition kernels.py:34