119 def __init__(self, name, patch_size, min_volume_h, max_volume_h, cfl, domain_r, KOSig):
141 number_of_unknowns = sum(unknowns.values())
144#include "../CCZ4Kernels.h"
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
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,
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
176 if args.implementation==
"fd4-rk1-adaptive" or args.implementation==
"fd4-rk1-adaptive-enclave":
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,
190 if args.implementation==
"fd4-rk1-adaptive" or args.implementation==
"fd4-rk1-adaptive-enclave":
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,
202 pde_terms_without_state=
True
205 if args.implementation==
"fd4-rk1-fixed":
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,
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,
228 time_step_relaxation=cfl
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
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
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
265 self.switch_to_heap_storage(
True,
True)
279 constexpr int itmax = {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}} * {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}};
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}};
293 constexpr int itmax = {{NUMBER_OF_VOLUMES_PER_AXIS}} * {{NUMBER_OF_VOLUMES_PER_AXIS}};
297 constexpr int itmax = {{NUMBER_OF_VOLUMES_PER_AXIS}} * {{NUMBER_OF_VOLUMES_PER_AXIS}} * {{NUMBER_OF_VOLUMES_PER_AXIS}};
303 for (int i=0;i<itmax;i++)
305 applications::exahype2::ccz4::enforceCCZ4constraints( newQ+index );
306 index += {{NUMBER_OF_UNKNOWNS}} + {{NUMBER_OF_AUXILIARY_VARIABLES}};
313 SuperClass.create_action_sets(self)
314 self._action_set_couple_resolution_transitions_and_handle_dynamic_mesh_refinement.additional_includes +=
"""
315#include "../CCZ4Kernels.h"
320 We take this routine to add a few additional include statements.
322 return SuperClass.get_user_action_set_includes(self) + self._my_user_includes
329 Add the constraint verification code
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
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
341 self._auxiliary_variables = 6+3+9
343 self._my_user_includes +=
"""
344 #include "../CCZ4Kernels.h"
345 #include "../AbstractFiniteVolumeCCZ4.h"
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];
357 for (int d=0; d<3; d++) currentPosition[d]=marker.getOffset()(d)+(cell(d)+0.5)*volumeH;
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.
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 };
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;
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]);
389 double constraints[6]={ 0 };
392 const int cellSerialised = peano4::utils::dLinearised(currentCell, patchSize + 2*1);
394 admconstraints(constraints,oldQWithHalo+cellSerialised*(59+n_a_v),gradQ);
396 for(int i=0;i<6;i++){
397 oldQWithHalo[cellSerialised*(59+n_a_v)+59+i] = constraints[i];
400 double ncpsum=0, ncp[1]={0};
401 double output[2+9]={ 0 };
402 for (int d=0; d<3; d++)
404 ThetaOutputNCP(ncp,faceQ+2*d*59,faceDeltaQ+2*d*59,d);
406 ThetaOutputNCP(ncp,faceQ+2*d*59+59,faceDeltaQ+2*d*59+59,d);
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];
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])+"]");}
422 ///////////////////////////////////////////////////
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;
431 oldQWithHalo[cellSerialised*(59+n_a_v)+59+6]=pow((P1*P1+P2*P2+P3*P3),0.5);//gradient of phi
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
441 self.create_data_structures()
442 self.create_action_sets()
449 Add the constraint verification code
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
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
461 self._auxiliary_variables = 59*3
463 self._preprocess_reconstructed_patch =
"""
464 const int patchSize = """ + str( self._patch.dim[0] ) +
""";
465 int aux_var=""" + str( self._auxiliary_variables ) +
""";
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);
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;
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;
491 self.create_data_structures()
492 self.create_action_sets()
499 add puncutrestracker, constraints writer, psi4 writer, AMRflag(currently using gradient of phi)
501 self._auxiliary_variables = 9
503 self._my_user_includes +=
"""
504 #include "../libtwopunctures/TP_PunctureTracker.h"
505 #include "../CCZ4Kernels.h"
510 self._preprocess_reconstructed_patch =
"""
511 const int patchSize = """ + str( self._patch.dim[0] ) +
""";
512 double volumeH = ::exahype2::fv::getVolumeLength(marker.h(),patchSize);
515 std::string att="_"""+args.tra_name+
""".txt"; std::string p1="puncture1"; std::string p2="puncture2"; std::string tem="ztem";
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
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
525 //fin.open((tem+att),std::ios::out|std::ios::trunc);
526 //fin << "tem file" << std::endl;
529 fin.open((p1+att),std::ios::in);
530 std::string pos=getLastLine(fin);
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);
538 CoorReadIn(coor2,pos2);
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];
550 double result[inter_number];
551 Interpolation(result,IndexForInterpolate,raw,coor1,marker.getOffset(),volumeH,patchSize);
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;
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;
565 fin << "after inter" << result[0] << " " << result[1] << " " << result[2] << " " << t << std::endl;
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];
578 double result[inter_number];
579 Interpolation(result,IndexForInterpolate,raw,coor2,marker.getOffset(),volumeH,patchSize);
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;
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);
594 double gradQ[3*59]={ 0 };
596 for (int d=0; d<3; d++) {
597 tarch::la::Vector<Dimensions,int> leftCell = currentCell;
598 tarch::la::Vector<Dimensions,int> rightCell = currentCell;
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;
607 double constraints[n_a_v-1]={ 0 };
608 const int cellSerialised = peano4::utils::dLinearised(currentCell, patchSize + 2*1);
610 admconstraints(constraints,oldQWithHalo+cellSerialised*(59+n_a_v),gradQ);
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];
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;
622 oldQWithHalo[cellSerialised*(59+n_a_v)+59+6]=pow((P1*P1+P2*P2+P3*P3),0.5);
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
633 self.create_data_structures()
634 self.create_action_sets()
640 self._auxiliary_variables = 2
642 self._my_user_includes +=
"""
643 #include "../libtwopunctures/TP_PunctureTracker.h"
644 #include "../CCZ4Kernels.h"
647 self._preprocess_reconstructed_patch =
"""
648 const int patchSize = """ + str( self._patch.dim[0] ) +
""";
649 double volumeH = ::exahype2::fv::getVolumeLength(marker.h(),patchSize);
652 const int overlap=3; //make sure you are using fd4 solver!
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);
659 double gradQ[3*59]={ 0 };
661 for (int d=0; d<3; d++) {
662 tarch::la::Vector<Dimensions,int> leftCell = currentCell;
663 tarch::la::Vector<Dimensions,int> rightCell = currentCell;
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;
673 const int cellSerialised = peano4::utils::dLinearised(currentCell, patchSize + 2*overlap);
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];
686 self.create_data_structures()
687 self.create_action_sets()
695 if args.exe_name!=
"":
698 if not args.tra_name==
"de":
699 exe +=
"_" + args.tra_name
709 if SuperClass==exahype2.solvers.aderdg.NonFusedGenericRusanovFixedTimeStepSize:
713 time_step_size = 0.001
714 except Exception
as e:
723 except Exception
as e:
725 msg =
"Warning: RKDG not supported on this machine"
727 userinfo.append((msg,e))
730 solver_name =
"ADERDG" + solver_name
732 solver_name =
"RKDG" + solver_name
734 solver_name =
"FiniteVolume" + solver_name
739 print(
"No minimal mesh size chosen. Set it to max mesh size (regular grid)" )
743 my_solver = exahype2.solvers.aderdg.NonFusedGenericRusanovFixedTimeStepSize(
744 solver_name, order, unknowns, 0,
745 exahype2.solvers.aderdg.Polynomials.Gauss_Legendre,
746 min_h, args.max_h, time_step_size,
748 ncp = exahype2.solvers.PDETerms.User_Defined_Implementation,
749 sources = exahype2.solvers.PDETerms.User_Defined_Implementation
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))
755 userinfo.append((
"The solver is "+args.implementation,
None))
761 if args.marker==
"poisson":
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 )
772 if args.interpolation==
"order-2":
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" )
788 if args.interpolation==
"order-2":
789 my_solver.interpolation =
"linear"
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))
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()
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))})
831 if args.scenario==
"two-punctures":
832 msg =
"Periodic BC deactivated because you pick Puncture scenario\nInitialize binary black holes"
834 periodic_boundary_conditions = [
False,
False,
False]
835 intparams.update({
"swi":2})
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"
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"
847 periodic_boundary_conditions = [
True,
True,
True]
848 userinfo.append((msg,
None))
850 msg =
"WARNING: Periodic BC deactivated by hand"
852 periodic_boundary_conditions = [
False,
False,
False]
853 userinfo.append((msg,
None))
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(
872 const double * __restrict__ Q,
873 const tarch::la::Vector<Dimensions,double>& faceCentre,
874 const tarch::la::Vector<Dimensions,double>& gridCellH,
879 return repositories::{{SOLVER_INSTANCE}}.maxEigenvalue( Q, faceCentre, gridCellH, t, dt, normal );
882 double * __restrict__ Q,
883 const tarch::la::Vector<Dimensions,double>& faceCentre,
884 const tarch::la::Vector<Dimensions,double>& gridCellH
886 for (int i=0; i<"""+str(my_solver._unknowns+my_solver._auxiliary_variables)+
"""; i++) {
889 Q[0] = 1.0; Q[3] = 1.0; Q[5] = 1.0;
890 Q[16] = 1.0; Q[54] = 1.0;
894 {{FACE_METADATA_ACCESSOR}}.getOldTimeStamp(marker.getSelectedFaceNumber()<Dimensions ? 1 : 0),
895 repositories::{{SOLVER_INSTANCE}}.getMinTimeStepSize(),
896 {{NUMBER_OF_GRID_CELLS_PER_PATCH_PER_AXIS}},
898 {{NUMBER_OF_UNKNOWNS}},
899 {{NUMBER_OF_AUXILIARY_VARIABLES}},
900 marker.getSelectedFaceNumber(),
902 fineGridFace{{UNKNOWN_IDENTIFIER}}Old.value,
903 fineGridFace{{UNKNOWN_IDENTIFIER}}New.value
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))
932 raise Exception(
"Scenario " + args.scenario +
" is now unknown")
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)
938 project.add_solver(my_solver)
940 build_mode = modes[args.mode]
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)
952 userinfo.append((msg,
None))
954 project.set_global_simulation_parameters(
958 args.plot_start_time, args.plot_step_size,
959 periodic_boundary_conditions,
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))
965 project.set_Peano4_installation(
"../../..", build_mode)
971 if not args.path==
"./":
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]
981 my_solver.select_dofs_to_print = [0,12,16,17,18,53,54]
988 project.set_load_balancing(
"toolbox::loadbalancing::strategies::SpreadOutHierarchically",
"(new ::exahype2::LoadBalancingConfiguration(0.98))" )
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
1008 time_stepping_kernel=
"toolbox::particles::staticPosition",
1009 project_on_tracer_properties_kernel=
"::exahype2::fv::projectAllValuesOntoParticle_piecewiseLinear"))
1012 project.add_action_set_to_timestepping(exahype2.tracer.DumpTrajectoryIntoDatabase(
1013 particle_set=tracer_particles,
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
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 )
1028 project.add_action_set_to_timestepping(
1030 tracer_particles,my_solver,
1031 [17,18,19],range(65),-1,
1033 time_stepping_kernel=
"toolbox::particles::StaticPosition"
1037 project.add_action_set_to_timestepping(
1038 peano4.toolbox.particles.UpdateParallelState(
1039 particle_set=tracer_particles
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)
1045 if args.add_tracer==4:
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))
1052 if path==
"./": path1=
"."
1054 project.add_action_set_to_timestepping(exahype2.tracer.DumpTrajectoryIntoDatabase(
1055 particle_set=tracer_particles,
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
1074 peano4_project = project.generate_Peano4_project(verbose=
True)
1076 if args.scenario==
"gauge" or args.scenario==
"linear" or args.scenario==
"dia_gauge" or args.scenario==
"flat":
1078 elif args.scenario==
"two-punctures" or args.scenario==
"single-puncture":
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" )
1096 raise Exception(
"Scenario " + args.scenario +
" is now unknown")
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" )
1133 peano4_project.generate( throw_away_data_after_generation=
False )
1135 peano4_project.build( make_clean_first =
True )
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:
1148 else: print(msg,
"Exception: {}".format(exception))
1151 print(
"=========================================================")