From e4ef95d2227e74735c6c6e60df1e9cc07ab1d480 Mon Sep 17 00:00:00 2001 From: Vinpasso Date: Sat, 7 Dec 2019 16:17:05 +0100 Subject: [PATCH] Import changes to lulesh from main repository --- CMakeLists.txt | 29 + laik-lulesh-repartition.cc | 345 +- laik_partitioners.cc | 4 +- laik_vector.cc | 92 +- laik_vector.h | 18 +- laik_vector_comm_exclusive_halo.cc | 9 +- laik_vector_comm_exclusive_halo.h | 3 +- laik_vector_comm_overlapping_overlapping.cc | 10 +- laik_vector_comm_overlapping_overlapping.h | 4 +- laik_vector_repart_exclusive.cc | 138 +- laik_vector_repart_exclusive.h | 36 +- laik_vector_repart_overlapping.cc | 52 +- laik_vector_repart_overlapping.h | 9 +- lulesh-util.cc | 35 +- lulesh.cc | 4633 ++++++++++--------- lulesh.h | 23 +- lulesh_tuple.h | 5 +- 17 files changed, 2998 insertions(+), 2447 deletions(-) create mode 100755 CMakeLists.txt diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100755 index 0000000..c0992ff --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,29 @@ +# Build the C examples +FILE(GLOB SRCS *.cc) +FILE(GLOB HEADERS *.h) + +add_executable ("lulesh" + ${SRCS} + ${HEADERS} + "../util/fault-tolerance-options.h" + "../util/fault-tolerance-options.c" + ) + +target_compile_options("lulesh" + PRIVATE + "-Wall" + "-Wextra" + ) + +add_definitions("-DUSE_MPI=1" "-DREPARTITIONING" "-DFAULT_TOLERANCE") + +foreach(SOURCE ${SRCS} ${HEADERS}) + set_source_files_properties( ${SOURCE} PROPERTIES COMPILE_FLAGS "-Wno-error=implicit-fallthrough -Wno-error=unused-parameter -Wno-error=unused-but-set-parameter -Wno-error=unknown-pragmas" ) +endforeach() + +target_link_libraries ("lulesh" + PRIVATE "laik" + PRIVATE "m" + PRIVATE "mpi" + ) + diff --git a/laik-lulesh-repartition.cc b/laik-lulesh-repartition.cc index ab72dce..fd5b48b 100644 --- a/laik-lulesh-repartition.cc +++ b/laik-lulesh-repartition.cc @@ -1,106 +1,219 @@ #include "lulesh.h" -#include +#include "laik_partitioners.h" #include -void Domain::re_distribute_data_structures(Laik_Group* new_group, Laik_Partitioning* p_exclusive, Laik_Partitioning* p_halo, Laik_Partitioning* p_overlapping, Laik_Transition *t_to_exclusive, Laik_Transition *t_to_halo, Laik_Transition *t_to_overlapping_init, Laik_Transition *t_to_overlapping_reduce){ +void +Domain::re_distribute_data_structures(Laik_Group *new_group, Laik_Partitioning *p_exclusive, Laik_Partitioning *p_halo, + Laik_Partitioning *p_overlapping, Laik_Transition *t_to_exclusive, + Laik_Transition *t_to_halo, Laik_Transition *t_to_overlapping_init, + Laik_Transition *t_to_overlapping_reduce, bool suppressSwitchToP1) { #ifdef REPARTITIONING - m_x.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce); - m_y.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce); - m_z.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce); - m_xd.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce); - m_yd.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce); - m_zd.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce); - m_xdd.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce); - m_ydd.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce); - m_zdd.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce); + m_x.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce, suppressSwitchToP1); + m_y.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce, suppressSwitchToP1); + m_z.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce, suppressSwitchToP1); + m_xd.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce, suppressSwitchToP1); + m_yd.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce, suppressSwitchToP1); + m_zd.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce, suppressSwitchToP1); + m_xdd.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce, suppressSwitchToP1); + m_ydd.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce, suppressSwitchToP1); + m_zdd.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce, suppressSwitchToP1); #endif - m_fx.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce); - m_fy.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce); - m_fz.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce); - m_nodalMass.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce); + m_fx.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce, suppressSwitchToP1); + m_fy.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce, suppressSwitchToP1); + m_fz.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce, suppressSwitchToP1); + m_nodalMass.migrate(new_group, p_overlapping, p_overlapping, t_to_overlapping_init, t_to_overlapping_reduce, suppressSwitchToP1); #ifdef REPARTITIONING - m_dxx.migrate(new_group, p_exclusive,nullptr, nullptr, nullptr); - m_dyy.migrate(new_group, p_exclusive,nullptr, nullptr, nullptr); - m_dzz.migrate(new_group, p_exclusive,nullptr, nullptr, nullptr); + m_dxx.migrate(new_group, p_exclusive, nullptr, nullptr, nullptr, suppressSwitchToP1); + m_dyy.migrate(new_group, p_exclusive, nullptr, nullptr, nullptr, suppressSwitchToP1); + m_dzz.migrate(new_group, p_exclusive, nullptr, nullptr, nullptr, suppressSwitchToP1); #endif - m_delv_xi.migrate(new_group, p_exclusive, p_halo, t_to_exclusive, t_to_halo); - m_delv_eta.migrate(new_group, p_exclusive, p_halo, t_to_exclusive, t_to_halo); - m_delv_zeta.migrate(new_group, p_exclusive, p_halo, t_to_exclusive, t_to_halo); + m_delv_xi.migrate(new_group, p_exclusive, p_halo, t_to_exclusive, t_to_halo, suppressSwitchToP1); + m_delv_eta.migrate(new_group, p_exclusive, p_halo, t_to_exclusive, t_to_halo, suppressSwitchToP1); + m_delv_zeta.migrate(new_group, p_exclusive, p_halo, t_to_exclusive, t_to_halo, suppressSwitchToP1); #ifdef REPARTITIONING - m_delx_xi.migrate(new_group, p_exclusive,nullptr, nullptr, nullptr); - m_delx_eta.migrate(new_group, p_exclusive,nullptr, nullptr, nullptr); - m_delx_zeta.migrate(new_group, p_exclusive,nullptr, nullptr, nullptr); - m_e.migrate(new_group, p_exclusive,nullptr, nullptr, nullptr); - m_p.migrate(new_group, p_exclusive,nullptr, nullptr, nullptr); - m_q.migrate(new_group, p_exclusive,nullptr, nullptr, nullptr); - m_ql.migrate(new_group, p_exclusive,nullptr, nullptr, nullptr); - m_qq.migrate(new_group, p_exclusive,nullptr, nullptr, nullptr); - m_v.migrate(new_group, p_exclusive,nullptr, nullptr, nullptr); - m_volo.migrate(new_group, p_exclusive,nullptr, nullptr, nullptr); + m_delx_xi.migrate(new_group, p_exclusive, nullptr, nullptr, nullptr, suppressSwitchToP1); + m_delx_eta.migrate(new_group, p_exclusive, nullptr, nullptr, nullptr, suppressSwitchToP1); + m_delx_zeta.migrate(new_group, p_exclusive, nullptr, nullptr, nullptr, suppressSwitchToP1); + m_e.migrate(new_group, p_exclusive, nullptr, nullptr, nullptr, suppressSwitchToP1); + m_p.migrate(new_group, p_exclusive, nullptr, nullptr, nullptr, suppressSwitchToP1); + m_q.migrate(new_group, p_exclusive, nullptr, nullptr, nullptr, suppressSwitchToP1); + m_ql.migrate(new_group, p_exclusive, nullptr, nullptr, nullptr, suppressSwitchToP1); + m_qq.migrate(new_group, p_exclusive, nullptr, nullptr, nullptr, suppressSwitchToP1); + m_v.migrate(new_group, p_exclusive, nullptr, nullptr, nullptr, suppressSwitchToP1); + m_volo.migrate(new_group, p_exclusive, nullptr, nullptr, nullptr, suppressSwitchToP1); //m_vnew.migrate(new_group, p_exclusive,nullptr, nullptr, nullptr); - m_delv.migrate(new_group, p_exclusive,nullptr, nullptr, nullptr); - m_vdov.migrate(new_group, p_exclusive,nullptr, nullptr, nullptr); - m_arealg.migrate(new_group, p_exclusive,nullptr, nullptr, nullptr); - m_ss.migrate(new_group, p_exclusive,nullptr, nullptr, nullptr); - m_elemMass.migrate(new_group, p_exclusive,nullptr, nullptr, nullptr); + m_delv.migrate(new_group, p_exclusive, nullptr, nullptr, nullptr, suppressSwitchToP1); + m_vdov.migrate(new_group, p_exclusive, nullptr, nullptr, nullptr, suppressSwitchToP1); + m_arealg.migrate(new_group, p_exclusive, nullptr, nullptr, nullptr, suppressSwitchToP1); + m_ss.migrate(new_group, p_exclusive, nullptr, nullptr, nullptr, suppressSwitchToP1); + m_elemMass.migrate(new_group, p_exclusive, nullptr, nullptr, nullptr, suppressSwitchToP1); #endif - this -> world = new_group; + this->world = new_group; } -void init_config_params(Laik_Group* group, int& b,int& f, int& d, int& u, int& l, int& r){ +Real_t checkpointDeltaTimeMultLB, checkpointDeltaTimeMultUB, checkpointDTCourant, checkpointDTHydro, checkpointDTMax; +Real_t checkpointTime, checkpointDeltaTime; +Int_t checkpointCycle; + +void Domain::createCheckpoints(std::vector &checkpoints, int redundancyCount, int rotationDistance) { +#ifdef REPARTITIONING + checkpoints.push_back(m_x.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_y.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_z.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_xd.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_yd.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_zd.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_xdd.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_ydd.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_zdd.checkpoint(redundancyCount, rotationDistance)); +#endif + checkpoints.push_back(m_fx.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_fy.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_fz.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_nodalMass.checkpoint(redundancyCount, rotationDistance)); +#ifdef REPARTITIONING + checkpoints.push_back(m_dxx.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_dyy.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_dzz.checkpoint(redundancyCount, rotationDistance)); +#endif + checkpoints.push_back(m_delv_xi.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_delv_eta.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_delv_zeta.checkpoint(redundancyCount, rotationDistance)); +#ifdef REPARTITIONING + checkpoints.push_back(m_delx_xi.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_delx_eta.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_delx_zeta.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_e.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_p.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_q.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_ql.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_qq.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_v.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_volo.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_delv.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_vdov.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_arealg.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_ss.checkpoint(redundancyCount, rotationDistance)); + checkpoints.push_back(m_elemMass.checkpoint(redundancyCount, rotationDistance)); +#endif + + checkpointDeltaTimeMultLB = deltatimemultlb(); + checkpointDeltaTimeMultUB = deltatimemultub(); + checkpointDTCourant = dtcourant(); + checkpointDTHydro = dthydro(); + checkpointDTMax = dtmax(); + checkpointTime = time(); + checkpointCycle = cycle(); + checkpointDeltaTime = deltatime(); +} + +int Domain::restore(std::vector &checkpoints, Laik_Group *newGroup) { + deltatimemultlb() = checkpointDeltaTimeMultLB; + deltatimemultub() = checkpointDeltaTimeMultUB; + dtcourant() = checkpointDTCourant; + dthydro() = checkpointDTHydro; + dtmax() = checkpointDTMax; + time() = checkpointTime; + cycle() = checkpointCycle; + deltatime() = checkpointDeltaTime; + + int index = 0; +#ifdef REPARTITIONING + m_x.restore(checkpoints[index++], newGroup); + m_y.restore(checkpoints[index++], newGroup); + m_z.restore(checkpoints[index++], newGroup); + m_xd.restore(checkpoints[index++], newGroup); + m_yd.restore(checkpoints[index++], newGroup); + m_zd.restore(checkpoints[index++], newGroup); + m_xdd.restore(checkpoints[index++], newGroup); + m_ydd.restore(checkpoints[index++], newGroup); + m_zdd.restore(checkpoints[index++], newGroup); +#endif + m_fx.restore(checkpoints[index++], newGroup); + m_fy.restore(checkpoints[index++], newGroup); + m_fz.restore(checkpoints[index++], newGroup); + m_nodalMass.restore(checkpoints[index++], newGroup); +#ifdef REPARTITIONING + m_dxx.restore(checkpoints[index++], newGroup); + m_dyy.restore(checkpoints[index++], newGroup); + m_dzz.restore(checkpoints[index++], newGroup); +#endif + m_delv_xi.restore(checkpoints[index++], newGroup); + m_delv_eta.restore(checkpoints[index++], newGroup); + m_delv_zeta.restore(checkpoints[index++], newGroup); +#ifdef REPARTITIONING + m_delx_xi.restore(checkpoints[index++], newGroup); + m_delx_eta.restore(checkpoints[index++], newGroup); + m_delx_zeta.restore(checkpoints[index++], newGroup); + m_e.restore(checkpoints[index++], newGroup); + m_p.restore(checkpoints[index++], newGroup); + m_q.restore(checkpoints[index++], newGroup); + m_ql.restore(checkpoints[index++], newGroup); + m_qq.restore(checkpoints[index++], newGroup); + m_v.restore(checkpoints[index++], newGroup); + m_volo.restore(checkpoints[index++], newGroup); + m_delv.restore(checkpoints[index++], newGroup); + m_vdov.restore(checkpoints[index++], newGroup); + m_arealg.restore(checkpoints[index++], newGroup); + m_ss.restore(checkpoints[index++], newGroup); + m_elemMass.restore(checkpoints[index++], newGroup); +#endif + return index; +} + +void init_config_params(Laik_Group *group, int &b, int &f, int &d, int &u, int &l, int &r) { int col, row, plane, side; InitMeshDecomp(laik_size(group), laik_myid(group), &col, &row, &plane, &side); - b=1; - f=1; - d=1; - u=1; - l=1; - r=1; + b = 1; + f = 1; + d = 1; + u = 1; + l = 1; + r = 1; - if (col==0) { - l=0; + if (col == 0) { + l = 0; } - if (col==side-1) { - r=0; + if (col == side - 1) { + r = 0; } - if (row==0) { - d=0; + if (row == 0) { + d = 0; } - if (row==side-1) { - u=0; + if (row == side - 1) { + u = 0; } - if (plane==0) { - b=0; + if (plane == 0) { + b = 0; } - if (plane==side-1) { - f=0; + if (plane == side - 1) { + f = 0; } //state = 0; } -void create_partitionings_and_transitions( Laik_Group *&world, - Laik_Space *&indexSpaceElements, - Laik_Space *&indexSpaceNodes, - Laik_Space *&indexSapceDt, - Laik_Partitioning *&exclusivePartitioning, - Laik_Partitioning *&haloPartitioning, - Laik_Partitioning *&overlapingPartitioning, - Laik_Partitioning *&allPartitioning, - Laik_Transition *&transitionToExclusive, - Laik_Transition *&transitionToHalo, - Laik_Transition *&transitionToOverlappingInit, - Laik_Transition *&transitionToOverlappingReduce) -{ +void create_partitionings_and_transitions(Laik_Group *&world, + Laik_Space *&indexSpaceElements, + Laik_Space *&indexSpaceNodes, + Laik_Space *&indexSapceDt, + Laik_Partitioning *&exclusivePartitioning, + Laik_Partitioning *&haloPartitioning, + Laik_Partitioning *&overlapingPartitioning, + Laik_Partitioning *&allPartitioning, + Laik_Transition *&transitionToExclusive, + Laik_Transition *&transitionToHalo, + Laik_Transition *&transitionToOverlappingInit, + Laik_Transition *&transitionToOverlappingReduce) { int halo_depth = 1; // depth of halos used in partitioner algorithms // run partitioner algorithms to @@ -115,9 +228,10 @@ void create_partitionings_and_transitions( Laik_Group *&world, // the iterations. Only for re-partitioning // they have to be recalculated exclusivePartitioning = laik_new_partitioning(exclusive_partitioner(), world, indexSpaceElements, 0); - haloPartitioning = laik_new_partitioning(overlaping_partitioner(halo_depth), world, indexSpaceElements, exclusivePartitioning); - overlapingPartitioning =laik_new_partitioning(overlaping_reduction_partitioner(halo_depth), - world, indexSpaceNodes, 0); + haloPartitioning = laik_new_partitioning(overlaping_partitioner(halo_depth), world, indexSpaceElements, + exclusivePartitioning); + overlapingPartitioning = laik_new_partitioning(overlaping_reduction_partitioner(halo_depth), + world, indexSpaceNodes, 0); // create all partitioning for dt to perform reductions allPartitioning = laik_new_partitioning(laik_All, world, indexSapceDt, 0); @@ -151,38 +265,87 @@ void remove_partitionings_and_transitions(Laik_Partitioning *&exclusivePartition Laik_Transition *&transitionToExclusive, Laik_Transition *&transitionToHalo, Laik_Transition *&transitionToOverlappingInit, - Laik_Transition *&transitionToOverlappingReduce) -{ - laik_free_partitioning (exclusivePartitioning); - laik_free_partitioning (haloPartitioning); - laik_free_partitioning (overlapingPartitioning); - laik_free_partitioning (allPartitioning); - free (transitionToExclusive); - free (transitionToHalo); - free (transitionToOverlappingInit); - free (transitionToOverlappingReduce); + Laik_Transition *&transitionToOverlappingReduce) { + laik_free_partitioning(exclusivePartitioning); + laik_free_partitioning(haloPartitioning); + laik_free_partitioning(overlapingPartitioning); + laik_free_partitioning(allPartitioning); + free(transitionToExclusive); + free(transitionToHalo); + free(transitionToOverlappingInit); + free(transitionToOverlappingReduce); +} + +void calculate_removing_list(Laik_Group *world, cmdLineOpts &opts, double side, double &newside, int &diffsize, + int *&removeList) { + int cursize = laik_size(world); + if(opts.repart != 0) { + // Planned repartitioning (not fault tolerant) + newside = cbrt(opts.repart); + if (newside - ((int) floor(newside + 0.1)) != 0) { + MPI_Abort(MPI_COMM_WORLD, -1); + } + diffsize = cursize - opts.repart; + } else { + //Fault tolerance + //TODO + printf("Not implemented\n"); + abort(); + } + removeList = (int *) malloc(diffsize * sizeof(int)); + for (int i = 0; i < diffsize; i++) { + removeList[i] = i + opts.repart; + } + + // check if the repartitioning scenario is valid (the target process group is a cubic int and + // the total size of elements in domain still valid) + double new_nx = (double) opts.nx * (double) side / (double) newside; + double verifier; + + if (modf(new_nx, &verifier) != 0.0) { + std::cout << "Repartitioning is not allowed for inbalanced domains after repartitioning. \n"; + MPI_Abort(MPI_COMM_WORLD, -1); + } } -void calculate_removing_list(Laik_Group* world, cmdLineOpts& opts, double side, double& newside, int& diffsize, int *&removeList) -{ +void calculate_removing_list_ft(Laik_Group *world, cmdLineOpts &opts, double side, double &newside, int &diffsize, + int *&removeList, int *nodeStatuses) { int cursize = laik_size(world); newside = cbrt(opts.repart); - if (newside - ((int) floor(newside+0.1)) != 0){ - MPI_Abort(MPI_COMM_WORLD,-1); + if (newside - ((int) floor(newside + 0.1)) != 0) { + MPI_Abort(MPI_COMM_WORLD, -1); } diffsize = cursize - opts.repart; - removeList = (int*) malloc (diffsize * sizeof(int)); - for (int i=0; i -#include +#include "laik_partitioners.h" +#include "lulesh.h" /** * @brief Exclusiv Partitioner diff --git a/laik_vector.cc b/laik_vector.cc index eefb26f..8678ea0 100644 --- a/laik_vector.cc +++ b/laik_vector.cc @@ -1,14 +1,17 @@ -#include +#include "laik_vector.h" #include "laik_vector_comm_exclusive_halo.h" #include "laik_vector_comm_overlapping_overlapping.h" #include "laik_vector_repart_exclusive.h" #include "laik_vector_repart_overlapping.h" -#include -#include +#include "laik_partitioners.h" +#include "lulesh.h" #include #include #include +#include +#include +#include template laik_vector::laik_vector(Laik_Instance* inst, Laik_Group* world, Laik_Space* indexSpace, Laik_Partitioning *p1, Laik_Partitioning *p2, Laik_Transition* t1, Laik_Transition* t2, Laik_ReductionOperation operation):reduction_operation(operation){ @@ -44,6 +47,89 @@ void laik_vector::test_print(){ template void laik_vector::clear(){} +template +void laik_vector::copyLaikDataToVector(std::vector &data_vector) { + uint64_t cnt; + T *base; + // copy the data back into the stl vectors +// int nSlices = laik_my_slicecount(this->p1); + int nSlices = laik_my_slicecount(laik_data_get_partitioning(this->data)); + for (int n = 0; n < nSlices; n++) { + assert(this->data != NULL); + assert(laik_data_get_partitioning(this->data) != nullptr); + assert(laik_my_slicecount(laik_data_get_partitioning(this->data)) == nSlices); + laik_get_map_1d(this->data, n, (void **) &base, &cnt); + uint64_t elemOffset = n * cnt; + laik_log(LAIK_LL_Debug, "Copy LAIK data to vector: vector (capacity) %zu data %" PRIu64 + " offset %" PRIu64 " length %" PRIu64, data_vector.capacity(), cnt, + elemOffset, cnt); + assert(elemOffset >= 0 && elemOffset + cnt <= data_vector.capacity()); + memcpy(&data_vector[0] + elemOffset, base, cnt * sizeof(T)); + //std::copy(data_vector.begin() + n*count ,data_vector.begin() + (n+1)*count-1 , base); + } +} + +template +void laik_vector::copyVectorToLaikData(std::vector &data_vector) { + uint64_t cnt; + T *base; + // copy the data from stl vector into the laik container + int nSlices = laik_my_slicecount(this->p1); + for (int n = 0; n < nSlices; n++) { + laik_get_map_1d(this->data, n, (void **) &base, &cnt); +// laik_log(LAIK_LL_Info, "Copy vector to LAIK data: vector (size) %lu data %lu", data_vector.size(), cnt); + assert(n * cnt >= 0 && n * cnt + cnt <= data_vector.capacity()); + memcpy(base, &data_vector[0] + n * cnt, cnt * sizeof(T)); + //std::copy( base, base + cnt, data_vector.begin() + n*count ); + } +} + +template +void laik_vector::resizeVector(std::vector &data_vector) {// resize vector + uint64_t cnt; + T* base; + assert(laik_my_mapcount(laik_data_get_partitioning(this->data)) == 1); + laik_get_map_1d(this->data, 0, (void **)&base, &cnt); + int s = cnt*cnt*cnt; + data_vector.resize(s); +} + +template +void laik_vector::resizeVectorToLaikData(std::vector &data_vector) {// resize vector + uint64_t cnt = 0; + for (int i = 0; i < laik_my_slicecount(laik_data_get_partitioning(data)); ++i) { + cnt += laik_slice_size(laik_taskslice_get_slice(laik_my_slice(laik_data_get_partitioning(data), i))); + } + laik_log(LAIK_LL_Info, "Resizing vector from %zu to %" PRIu64, data_vector.capacity(), cnt); + data_vector.resize(cnt); +} + +template +void laik_vector::prepareMigration(bool suppressDataSwitchToP1) { + if(!suppressDataSwitchToP1) { + laik_switchto_partitioning(this->data, this->p1, LAIK_DF_None, LAIK_RO_None); + } +} + +#ifdef FAULT_TOLERANCE +template +Laik_Checkpoint * laik_vector::checkpoint(int redundancyCount, int rotationDistance) { +// std::cout << "Creating checkpoint of " << laik_my_slicecount(laik_data_get_partitioning(data)) << " slices." << std::endl; + return laik_checkpoint_create(data, laik_Master, redundancyCount, rotationDistance, + laik_data_get_group(data), LAIK_RO_Min); +} + +template +void laik_vector::restore(Laik_Checkpoint *checkpoint, Laik_Group *newGroup) { + // Set partitioning to backup partitioning so that it can be migrated later + assert(checkpoint->data != nullptr && laik_data_get_partitioning(checkpoint->data) != nullptr); + Laik_Partitioning* newPartitioning = laik_new_partitioning(laik_Master, newGroup, indexSpace, nullptr); + laik_switchto_partitioning(data, newPartitioning, LAIK_DF_None, LAIK_RO_None); +// laik_partitioning_migrate(laik_data_get_partitioning(checkpoint->data), newGroup); + laik_checkpoint_restore(checkpoint, data); +} + +#endif template class laik_vector; diff --git a/laik_vector.h b/laik_vector.h index 243de96..270685c 100644 --- a/laik_vector.h +++ b/laik_vector.h @@ -20,6 +20,7 @@ extern "C"{ template class laik_vector { + public: /** * @brief laik_vector constructor @@ -67,13 +68,26 @@ class laik_vector * @param t_new_1 transition to p_new_1 * @param t_new_2 transition to p_new_2 */ - virtual void migrate(Laik_Group* new_group, Laik_Partitioning* p_new_1, Laik_Partitioning* p_new_2, Laik_Transition* t_new_1, Laik_Transition* t_new_2) = 0; + virtual void + migrate(Laik_Group *new_group, Laik_Partitioning *p_new_1, Laik_Partitioning *p_new_2, Laik_Transition *t_new_1, + Laik_Transition *t_new_2, bool suppressSwitchToP1) = 0; /** * @brief clearing laik_vectors */ void clear(); + void copyLaikDataToVector(std::vector &data_vector); + void copyVectorToLaikData(std::vector &data_vector); + + void resizeVector(std::vector&); + void resizeVectorToLaikData(std::vector&); + +#ifdef FAULT_TOLERANCE + virtual Laik_Checkpoint *checkpoint(int redundancyCount, int rotationDistance); + virtual void restore(Laik_Checkpoint *checkpoint, Laik_Group *newGroup); +#endif + protected: // members from laik Laik_Instance* inst; // laik context @@ -111,6 +125,8 @@ class laik_vector * @brief test_print printing laik_vector for debug */ void test_print(); + + void prepareMigration(bool suppressDataSwitchToP1); }; #endif // LAIK_VECTOR diff --git a/laik_vector_comm_exclusive_halo.cc b/laik_vector_comm_exclusive_halo.cc index 22c582d..9ac154e 100644 --- a/laik_vector_comm_exclusive_halo.cc +++ b/laik_vector_comm_exclusive_halo.cc @@ -1,6 +1,6 @@ #include"laik_vector_comm_exclusive_halo.h" -#include -#include +#include "laik_partitioners.h" +#include "lulesh.h" #include #include #include @@ -248,12 +248,13 @@ void laik_vector_comm_exclusive_halo::switch_to_p2(){ } template -void laik_vector_comm_exclusive_halo::migrate(Laik_Group* new_group, Laik_Partitioning* p_new_1, Laik_Partitioning* p_new_2, Laik_Transition* t_new_1, Laik_Transition* t_new_2){ +void laik_vector_comm_exclusive_halo::migrate(Laik_Group* new_group, Laik_Partitioning* p_new_1, Laik_Partitioning* p_new_2, Laik_Transition* t_new_1, Laik_Transition* t_new_2, + bool suppressSwitchToP1){ uint64_t cnt; int* base; //int slice = 0; - laik_switchto_partitioning(this->data, this->p1, LAIK_DF_None, LAIK_RO_None); + this->prepareMigration(suppressSwitchToP1); // use the reservation API to precalculate the pointers Laik_Reservation* reservation = laik_reservation_new(this->data); diff --git a/laik_vector_comm_exclusive_halo.h b/laik_vector_comm_exclusive_halo.h index 0a43c6c..2d4927e 100644 --- a/laik_vector_comm_exclusive_halo.h +++ b/laik_vector_comm_exclusive_halo.h @@ -14,7 +14,8 @@ class laik_vector_comm_exclusive_halo:public laik_vector void resize(int count) override; void switch_to_p1() override; void switch_to_p2() override; - void migrate(Laik_Group* new_group, Laik_Partitioning* p_new_1, Laik_Partitioning* p_new_2, Laik_Transition* t_new_1, Laik_Transition* t_new_2) override; + void migrate(Laik_Group* new_group, Laik_Partitioning* p_new_1, Laik_Partitioning* p_new_2, Laik_Transition* t_new_1, Laik_Transition* t_new_2, + bool suppressSwitchToP1) override; }; template diff --git a/laik_vector_comm_overlapping_overlapping.cc b/laik_vector_comm_overlapping_overlapping.cc index 5f281ec..e3f596e 100644 --- a/laik_vector_comm_overlapping_overlapping.cc +++ b/laik_vector_comm_overlapping_overlapping.cc @@ -1,6 +1,6 @@ #include"laik_vector_comm_overlapping_overlapping.h" -#include -#include +#include "laik_partitioners.h" +#include "lulesh.h" #include #include #include @@ -99,12 +99,14 @@ void laik_vector_comm_overlapping_overlapping::switch_to_p2(){ } template -void laik_vector_comm_overlapping_overlapping::migrate(Laik_Group* new_group, Laik_Partitioning* p_new_1, Laik_Partitioning* p_new_2, Laik_Transition* t_new_1, Laik_Transition* t_new_2){ +void laik_vector_comm_overlapping_overlapping::migrate(Laik_Group *new_group, Laik_Partitioning *p_new_1, + Laik_Partitioning *p_new_2, Laik_Transition *t_new_1, + Laik_Transition *t_new_2, bool suppressSwitchToP1) { uint64_t cnt; int* base; //int slice = 0; - laik_switchto_partitioning(this->data, this->p1, LAIK_DF_None, LAIK_RO_Min); + this->prepareMigration(suppressSwitchToP1); Laik_Reservation* reservation = laik_reservation_new(this->data); laik_reservation_add(reservation, p_new_1); diff --git a/laik_vector_comm_overlapping_overlapping.h b/laik_vector_comm_overlapping_overlapping.h index 28a1ead..890c23b 100644 --- a/laik_vector_comm_overlapping_overlapping.h +++ b/laik_vector_comm_overlapping_overlapping.h @@ -14,7 +14,9 @@ class laik_vector_comm_overlapping_overlapping:public laik_vector void resize(int count) override; void switch_to_p1() override; void switch_to_p2() override; - void migrate(Laik_Group* new_group, Laik_Partitioning* p_new_1, Laik_Partitioning* p_new_2, Laik_Transition* t_new_1, Laik_Transition* t_new_2) override; + void migrate(Laik_Group *new_group, Laik_Partitioning *p_new_1, + Laik_Partitioning *p_new_2, Laik_Transition *t_new_1, + Laik_Transition *t_new_2, bool suppressSwitchToP1) override; }; template diff --git a/laik_vector_repart_exclusive.cc b/laik_vector_repart_exclusive.cc index 51b83c8..f9271c0 100644 --- a/laik_vector_repart_exclusive.cc +++ b/laik_vector_repart_exclusive.cc @@ -1,6 +1,6 @@ #include"laik_vector_repart_exclusive.h" -#include -#include +#include "laik_partitioners.h" +#include "lulesh.h" #include #include #include @@ -9,102 +9,106 @@ // implementation of laik_vector with exclusive partitioning (elem partitioning) // for repartitioning of exclusive data structs // //////////////////////////////////////////////////////////////////////// -template +template laik_vector_repart_exclusive::laik_vector_repart_exclusive(Laik_Instance *inst, - Laik_Group *world, - Laik_Space* indexSpace, - Laik_Partitioning *p1, Laik_Partitioning *p2, Laik_Transition* t1, Laik_Transition* t2, Laik_ReductionOperation operation):laik_vector(inst,world, indexSpace, p1, p2, t1, t2, operation){} -template -void laik_vector_repart_exclusive::resize(int count){ - - int s = count / laik_size(this->world); + Laik_Group *world, + Laik_Space *indexSpace, + Laik_Partitioning *p1, Laik_Partitioning *p2, + Laik_Transition *t1, Laik_Transition *t2, + Laik_ReductionOperation operation):laik_vector(inst, + world, + indexSpace, + p1, p2, + t1, t2, + operation) {} + +template +void laik_vector_repart_exclusive::resize(int count) { + + int s = count / laik_size(this->world); data_vector.resize(s); - this -> size = count; + this->size = count; - if (std::is_same ::value) { - this->data = laik_new_data(this->indexSpace, laik_Double ); + if (std::is_same::value) { + this->data = laik_new_data(this->indexSpace, laik_Double); - } - else if (std::is_same ::value){ - this->data = laik_new_data(this->indexSpace, laik_Int64 ); + } else if (std::is_same::value) { + this->data = laik_new_data(this->indexSpace, laik_Int64); } laik_switchto_partitioning(this->data, this->p1, LAIK_DF_None, this->reduction_operation); uint64_t cnt; - T* base; + T *base; int nSlices = laik_my_slicecount(this->p1); - for (int n = 0; n < nSlices; ++n) - { - laik_get_map_1d(this->data, n, (void **)&base, &cnt); + for (int n = 0; n < nSlices; ++n) { + laik_get_map_1d(this->data, n, (void **) &base, &cnt); } - this -> count = cnt; + this->count = cnt; } -template -T* laik_vector_repart_exclusive::calc_pointer(int idx, int state){ +template +T *laik_vector_repart_exclusive::calc_pointer(int idx, int state) { return nullptr; } -template -void laik_vector_repart_exclusive::precalculate_base_pointers(){ +template +void laik_vector_repart_exclusive::precalculate_base_pointers() { } -template -void laik_vector_repart_exclusive::switch_to_p1(){ +template +void laik_vector_repart_exclusive::switch_to_p1() { } -template -void laik_vector_repart_exclusive::switch_to_p2(){ +template +void laik_vector_repart_exclusive::switch_to_p2() { } -template -void laik_vector_repart_exclusive::migrate(Laik_Group* new_group, Laik_Partitioning* p_new_1, Laik_Partitioning* p_new_2, Laik_Transition* t_new_1, Laik_Transition* t_new_2){ - uint64_t cnt; - T* base; - int nSlices; - - this -> state = 0; - - laik_switchto_partitioning(this->data, this->p1, LAIK_DF_Preserve, LAIK_RO_None); - // copy the data from stl vector into the laik container - nSlices = laik_my_slicecount(this->p1); - for (int n = 0; n < nSlices; n++) - { - laik_get_map_1d(this->data, n, (void **)&base, &cnt); - memcpy(base, &data_vector[0] + n*cnt, cnt*sizeof(T)); - //std::copy( base, base + cnt, data_vector.begin() + n*count ); - } +template +void +laik_vector_repart_exclusive::migrate(Laik_Group *new_group, Laik_Partitioning *p_new_1, Laik_Partitioning *p_new_2, + Laik_Transition *t_new_1, Laik_Transition *t_new_2, bool suppressSwitchToP1) { - // perform switches for communication + this->state = 0; + + this->prepareMigration(suppressSwitchToP1); + + this->copyVectorToLaikData(data_vector); + // perform switches for communication laik_switchto_partitioning(this->data, p_new_1, LAIK_DF_Preserve, LAIK_RO_None); - this -> world = new_group; - if (laik_myid(this->world)<0) - return ; + this->world = new_group; + if (laik_myid(this->world) < 0) + return; - this -> p1=p_new_1; - this -> p2=p_new_2; - this -> t1=t_new_1; - this -> t2=t_new_2; + this->p1 = p_new_1; + this->p2 = p_new_2; + this->t1 = t_new_1; + this->t2 = t_new_2; + this->resizeVector(data_vector); - // resize vector - laik_get_map_1d(this->data, 0, (void **)&base, &cnt); - int s = cnt*cnt*cnt; - data_vector.resize(s); + this->copyLaikDataToVector(data_vector); +} - // copy the data back into the stl vecotrs - nSlices = laik_my_slicecount(this->p1); - for (int n = 0; n < nSlices; n++) - { - laik_get_map_1d(this->data, n, (void **)&base, &cnt); - memcpy(&data_vector[0] + n*cnt, base, cnt*sizeof(T)); - //std::copy(data_vector.begin() + n*count ,data_vector.begin() + (n+1)*count-1 , base); - } + +#ifdef FAULT_TOLERANCE +template +Laik_Checkpoint * laik_vector_repart_exclusive::checkpoint(int redundancyCount, int rotationDistance) { + this->copyVectorToLaikData(data_vector); + return laik_vector::checkpoint(redundancyCount, rotationDistance); +} + +template +void laik_vector_repart_exclusive::restore(Laik_Checkpoint *checkpoint, Laik_Group *newGroup) { + laik_vector::restore(checkpoint, newGroup); + this->resizeVectorToLaikData(data_vector); + this->copyLaikDataToVector(data_vector); } +#endif -template class laik_vector_repart_exclusive; +template +class laik_vector_repart_exclusive; diff --git a/laik_vector_repart_exclusive.h b/laik_vector_repart_exclusive.h index 230c31f..bf66019 100644 --- a/laik_vector_repart_exclusive.h +++ b/laik_vector_repart_exclusive.h @@ -3,27 +3,43 @@ #include"laik_vector.h" -template -class laik_vector_repart_exclusive:public laik_vector -{ +template +class laik_vector_repart_exclusive : public laik_vector { public: - laik_vector_repart_exclusive(Laik_Instance* inst, Laik_Group* world, Laik_Space* indexSpace, Laik_Partitioning *p1, Laik_Partitioning *p2, Laik_Transition* t1, Laik_Transition* t2, Laik_ReductionOperation operation = LAIK_RO_None); - inline T& operator [](int idx) override; - T* calc_pointer(int idx, int state); + laik_vector_repart_exclusive(Laik_Instance *inst, Laik_Group *world, Laik_Space *indexSpace, Laik_Partitioning *p1, + Laik_Partitioning *p2, Laik_Transition *t1, Laik_Transition *t2, + Laik_ReductionOperation operation = LAIK_RO_None); + + inline T &operator[](int idx) override; + + T *calc_pointer(int idx, int state); + void precalculate_base_pointers() override; + void resize(int count) override; + void switch_to_p1() override; + void switch_to_p2() override; - void migrate(Laik_Group* new_group, Laik_Partitioning* p_new_1, Laik_Partitioning* p_new_2, Laik_Transition* t_new_1, Laik_Transition* t_new_2) override; + + void + migrate(Laik_Group *new_group, Laik_Partitioning *p_new_1, Laik_Partitioning *p_new_2, + Laik_Transition *t_new_1, Laik_Transition *t_new_2, bool suppressSwitchToP1) override; + +#ifdef FAULT_TOLERANCE + Laik_Checkpoint *checkpoint(int redundancyCount, int rotationDistance) override; + void restore(Laik_Checkpoint *checkpoint, Laik_Group *newGroup) override; +#endif private: std::vector data_vector; + }; -template +template inline -T& laik_vector_repart_exclusive::operator [](int idx){ - return this -> data_vector[idx]; +T &laik_vector_repart_exclusive::operator[](int idx) { + return this->data_vector[idx]; } #endif // LAIK_VECTOR_REPART_EXCLUSIVE_H diff --git a/laik_vector_repart_overlapping.cc b/laik_vector_repart_overlapping.cc index 9ac3f32..6eb53b0 100644 --- a/laik_vector_repart_overlapping.cc +++ b/laik_vector_repart_overlapping.cc @@ -1,6 +1,7 @@ #include"laik_vector_repart_overlapping.h" -#include -#include +#include "laik_partitioners.h" +#include "lulesh.h" +#include "laik_vector.h" #include #include #include @@ -59,22 +60,15 @@ void laik_vector_repart_overlapping::switch_to_p2(){ } template -void laik_vector_repart_overlapping::migrate(Laik_Group* new_group, Laik_Partitioning* p_new_1, Laik_Partitioning* p_new_2, Laik_Transition* t_new_1, Laik_Transition* t_new_2){ - uint64_t cnt; - T* base; - int nSlices; +void laik_vector_repart_overlapping::migrate(Laik_Group *new_group, Laik_Partitioning *p_new_1, + Laik_Partitioning *p_new_2, Laik_Transition *t_new_1, + Laik_Transition *t_new_2, bool suppressSwitchToP1) { this -> state = 0; - laik_switchto_partitioning(this->data, this->p1, LAIK_DF_None, LAIK_RO_Min); - // copy the data from stl vector into the laik container - nSlices = laik_my_slicecount(this->p1); - for (int n = 0; n < nSlices; n++) - { - laik_get_map_1d(this->data, n, (void **)&base, &cnt); - memcpy(base, &data_vector[0] + n*cnt, cnt*sizeof(T)); - //std::copy( base, base + cnt, data_vector.begin() + n*count ); - } + this->prepareMigration(suppressSwitchToP1); + + this->copyVectorToLaikData(data_vector); // perform switches for communication laik_switchto_partitioning(this->data, p_new_1, LAIK_DF_Preserve, LAIK_RO_Min); @@ -87,21 +81,25 @@ void laik_vector_repart_overlapping::migrate(Laik_Group* new_group, Laik_Part this -> p2=p_new_2; this -> t1=t_new_1; this -> t2=t_new_2; + this->resizeVector(data_vector); - // resize vector - laik_get_map_1d(this->data, 0, (void **)&base, &cnt); - int s = cnt*cnt*cnt; - data_vector.resize(s); + this->copyLaikDataToVector(data_vector); - // copy the data back into the stl vecotrs - nSlices = laik_my_slicecount(this->p1); - for (int n = 0; n < nSlices; n++) - { - laik_get_map_1d(this->data, n, (void **)&base, &cnt); - memcpy(&data_vector[0] + n*cnt, base, cnt*sizeof(T)); - //std::copy(data_vector.begin() + n*count ,data_vector.begin() + (n+1)*count-1 , base); - } +} + +#ifdef FAULT_TOLERANCE +template +Laik_Checkpoint * laik_vector_repart_overlapping::checkpoint(int redundancyCount, int rotationDistance) { + this->copyVectorToLaikData(data_vector); + return laik_vector::checkpoint(redundancyCount, rotationDistance); +} +template +void laik_vector_repart_overlapping::restore(Laik_Checkpoint *checkpoint, Laik_Group *newGroup) { + laik_vector::restore(checkpoint, newGroup); + this->resizeVectorToLaikData(data_vector); + this->copyLaikDataToVector(data_vector); } +#endif template class laik_vector_repart_overlapping; diff --git a/laik_vector_repart_overlapping.h b/laik_vector_repart_overlapping.h index 7c5e3e6..808cf20 100644 --- a/laik_vector_repart_overlapping.h +++ b/laik_vector_repart_overlapping.h @@ -14,7 +14,14 @@ class laik_vector_repart_overlapping:public laik_vector void resize(int count) override; void switch_to_p1() override; void switch_to_p2() override; - void migrate(Laik_Group* new_group, Laik_Partitioning* p_new_1, Laik_Partitioning* p_new_2, Laik_Transition* t_new_1, Laik_Transition* t_new_2) override; + void migrate(Laik_Group *new_group, Laik_Partitioning *p_new_1, + Laik_Partitioning *p_new_2, Laik_Transition *t_new_1, + Laik_Transition *t_new_2, bool suppressSwitchToP1) override; + +#ifdef FAULT_TOLERANCE + Laik_Checkpoint *checkpoint(int redundancyCount, int rotationDistance) override; + void restore(Laik_Checkpoint *checkpoint, Laik_Group *newGroup) override; +#endif private: std::vector data_vector; diff --git a/lulesh-util.cc b/lulesh-util.cc index a7ad101..70361f1 100644 --- a/lulesh-util.cc +++ b/lulesh-util.cc @@ -7,6 +7,7 @@ #endif #include "lulesh.h" + /* Helper function for converting strings to ints, with error checking */ int StrToInt(const char *token, int *retVal) { @@ -16,7 +17,7 @@ int StrToInt(const char *token, int *retVal) if (token == NULL) return 0 ; - + c = token ; *retVal = (int)strtol(c, &endptr, decimal_base) ; if((endptr != c) && ((*endptr == ' ') || (*endptr == '\0'))) @@ -42,8 +43,10 @@ static void PrintCommandLineOptions(char *execname, int myRank) printf(" -v : Output viz file (requires compiling with -DVIZ_MESH\n"); printf(" -repart : enable repartitioning by defining the number of target group\n"); printf(" -repart_cycle : cycle at which repartitioning happens\n"); + printf(" -no-timers : disable MPI timing functionality to make output stable\n"); printf(" -h : This message\n"); - printf("\n\n"); + printf(FAULT_TOLERANCE_OPTIONS_HELP); + printf("\n\n"); } } @@ -51,7 +54,7 @@ static void ParseError(const char *message, int myRank) { if (myRank == 0) { printf("%s\n", message); -#if USE_MPI +#if USE_MPI MPI_Abort(MPI_COMM_WORLD, -1); #else exit(-1); @@ -59,8 +62,12 @@ static void ParseError(const char *message, int myRank) } } -void ParseCommandLineOptions(int argc, char *argv[], - int myRank, struct cmdLineOpts *opts) +void printOptions(cmdLineOpts* opts) { + laik_log(LAIK_LL_Info, "LAIK LULESH Options:\nQuiet: %i, %i iterations, size %i, %i regions, balance %i, cost %i, %i files, %i progress, %i viz output, %i repartitioning (at cycle %i).\n", opts->quiet, opts->its, opts->nx, opts->numReg, opts->balance, opts->cost, opts->numFiles, opts->showProg, opts->viz, opts->repart, opts->cycle); +} + +void +ParseCommandLineOptions(int argc, char *argv[], int myRank, struct cmdLineOpts *opts, FaultToleranceOptions *ftOptions) { if(argc > 1) { int i = 1; @@ -143,7 +150,7 @@ void ParseCommandLineOptions(int argc, char *argv[], } /* -v */ else if (strcmp(argv[i], "-v") == 0) { -#if VIZ_MESH +#if VIZ_MESH opts->viz = 1; #else ParseError("Use of -v requires compiling with -DVIZ_MESH\n", myRank); @@ -170,10 +177,18 @@ void ParseCommandLineOptions(int argc, char *argv[], } i+=2; } + else if (parseFaultToleranceOptionsProxy(argc, argv, &i, myRank, ftOptions)) { + // Successfully parsed arg, need to shift by one more arg since lulesh doesn't do any shifting on its own + i++; + } + else if(strcmp(argv[i], "-no-timers") == 0) { + opts->noTimer = true; + i++; + } /* -h */ else if (strcmp(argv[i], "-h") == 0) { PrintCommandLineOptions(argv[0], myRank); -#if USE_MPI +#if USE_MPI MPI_Abort(MPI_COMM_WORLD, 0); #else exit(0); @@ -182,11 +197,13 @@ void ParseCommandLineOptions(int argc, char *argv[], else { char msg[80]; PrintCommandLineOptions(argv[0], myRank); - sprintf(msg, "ERROR: Unknown command line argument: %s\n", argv[i]); + sprintf(msg, "ERROR: Unknown command line argument at position %i: %s\n", i, argv[i]); ParseError(msg, myRank); } } } + + printOptions(opts); } ///////////////////////////////////////////////////////////////////// @@ -198,7 +215,7 @@ void VerifyAndWriteFinalOutput(Real_t elapsed_time, { // GrindTime1 only takes a single domain into account, and is thus a good way to measure // processor speed indepdendent of MPI parallelism. - // GrindTime2 takes into account speedups from MPI parallelism + // GrindTime2 takes into account speedups from MPI parallelism Real_t grindTime1 = ((elapsed_time*1e6)/locDom.cycle())/(nx*nx*nx); Real_t grindTime2 = ((elapsed_time*1e6)/locDom.cycle())/(nx*nx*nx*numRanks); diff --git a/lulesh.cc b/lulesh.cc index 2bf0908..1a281db 100644 --- a/lulesh.cc +++ b/lulesh.cc @@ -155,43 +155,45 @@ Additional BSD Notice #include #include #include +#include #if _OPENMP # include #endif /*laik headers*/ -extern "C"{ +extern "C" { #include "laik.h" #include "laik-backend-mpi.h" +#include "laik-backend-tcp.h" +#include "../util/fault-tolerance-options.h" } //porting to laik //#include "laik_vector.h" // includes lulesh.h -#include -#include -#include +#include "lulesh.h" +#include "laik_partitioners.h" +#include "laik_vector.h" /*********************************/ /* Data structure implementation */ /*********************************/ /* might want to add access methods so that memory can be */ + /* better managed, as in luleshFT */ -template -T *Allocate(size_t size) -{ - return static_cast(malloc(sizeof(T)*size)) ; +template +T *Allocate(size_t size) { + return static_cast(malloc(sizeof(T) * size)); } -template -void Release(T **ptr) -{ - if (*ptr != NULL) { - free(*ptr) ; - *ptr = NULL ; - } +template +void Release(T **ptr) { + if (*ptr != NULL) { + free(*ptr); + *ptr = NULL; + } } @@ -200,114 +202,115 @@ void Release(T **ptr) /* Work Routines */ +double laik_reduce(double input, Int_t myRank, Laik_Data *laikTimer, Laik_Partitioning *all, Laik_Partitioning *master); + +double getMPI_Wtime(cmdLineOpts* opts); + static inline -void TimeIncrement(Domain& domain, Laik_Data* laikDt, Laik_Partitioning* allPartitioning, Real_t &gnewdt) -{ - Real_t targetdt = domain.stoptime() - domain.time() ; - - if ((domain.dtfixed() <= Real_t(0.0)) && (domain.cycle() != Int_t(0))) { - Real_t ratio ; - Real_t olddt = domain.deltatime() ; - - /* This will require a reduction in parallel */ - gnewdt = Real_t(1.0e+20) ; - Real_t newdt ; - if (domain.dtcourant() < gnewdt) { - gnewdt = domain.dtcourant() / Real_t(2.0) ; - } - if (domain.dthydro() < gnewdt) { - gnewdt = domain.dthydro() * Real_t(2.0) / Real_t(3.0) ; - } +void TimeIncrement(Domain &domain, Laik_Data *laikDt, Laik_Partitioning *allPartitioning, Real_t &gnewdt) { + Real_t targetdt = domain.stoptime() - domain.time(); + + if ((domain.dtfixed() <= Real_t(0.0)) && (domain.cycle() != Int_t(0))) { + Real_t ratio; + Real_t olddt = domain.deltatime(); + + /* This will require a reduction in parallel */ + gnewdt = Real_t(1.0e+20); + Real_t newdt; + if (domain.dtcourant() < gnewdt) { + gnewdt = domain.dtcourant() / Real_t(2.0); + } + if (domain.dthydro() < gnewdt) { + gnewdt = domain.dthydro() * Real_t(2.0) / Real_t(3.0); + } -#if USE_MPI - /* +#if USE_MPI + /* MPI_Allreduce(&gnewdt, &newdt, 1, ((sizeof(Real_t) == 4) ? MPI_FLOAT : MPI_DOUBLE), MPI_MIN, MPI_COMM_WORLD) ; */ - laik_switchto_partitioning(laikDt, allPartitioning, LAIK_DF_Preserve, LAIK_RO_Min); - newdt = gnewdt; + laik_switchto_partitioning(laikDt, allPartitioning, LAIK_DF_Preserve, LAIK_RO_Min); + newdt = gnewdt; #else - newdt = gnewdt; + newdt = gnewdt; #endif - - ratio = newdt / olddt ; - if (ratio >= Real_t(1.0)) { - if (ratio < domain.deltatimemultlb()) { - newdt = olddt ; - } - else if (ratio > domain.deltatimemultub()) { - newdt = olddt*domain.deltatimemultub() ; - } - } - if (newdt > domain.dtmax()) { - newdt = domain.dtmax() ; - } - domain.deltatime() = newdt ; - } + ratio = newdt / olddt; + if (ratio >= Real_t(1.0)) { + if (ratio < domain.deltatimemultlb()) { + newdt = olddt; + } else if (ratio > domain.deltatimemultub()) { + newdt = olddt * domain.deltatimemultub(); + } + } - /* TRY TO PREVENT VERY SMALL SCALING ON THE NEXT CYCLE */ - if ((targetdt > domain.deltatime()) && - (targetdt < (Real_t(4.0) * domain.deltatime() / Real_t(3.0))) ) { - targetdt = Real_t(2.0) * domain.deltatime() / Real_t(3.0) ; - } + if (newdt > domain.dtmax()) { + newdt = domain.dtmax(); + } + domain.deltatime() = newdt; + } - if (targetdt < domain.deltatime()) { - domain.deltatime() = targetdt ; - } + /* TRY TO PREVENT VERY SMALL SCALING ON THE NEXT CYCLE */ + if ((targetdt > domain.deltatime()) && + (targetdt < (Real_t(4.0) * domain.deltatime() / Real_t(3.0)))) { + targetdt = Real_t(2.0) * domain.deltatime() / Real_t(3.0); + } + + if (targetdt < domain.deltatime()) { + domain.deltatime() = targetdt; + } - domain.time() += domain.deltatime() ; + domain.time() += domain.deltatime(); - ++domain.cycle() ; + ++domain.cycle(); } /******************************************/ static inline void CollectDomainNodesToElemNodes(Domain &domain, - const Index_t* elemToNode, + const Index_t *elemToNode, Real_t elemX[8], Real_t elemY[8], - Real_t elemZ[8]) -{ - Index_t nd0i = elemToNode[0] ; - Index_t nd1i = elemToNode[1] ; - Index_t nd2i = elemToNode[2] ; - Index_t nd3i = elemToNode[3] ; - Index_t nd4i = elemToNode[4] ; - Index_t nd5i = elemToNode[5] ; - Index_t nd6i = elemToNode[6] ; - Index_t nd7i = elemToNode[7] ; - - elemX[0] = domain.x(nd0i); - elemX[1] = domain.x(nd1i); - elemX[2] = domain.x(nd2i); - elemX[3] = domain.x(nd3i); - elemX[4] = domain.x(nd4i); - elemX[5] = domain.x(nd5i); - elemX[6] = domain.x(nd6i); - elemX[7] = domain.x(nd7i); - - elemY[0] = domain.y(nd0i); - elemY[1] = domain.y(nd1i); - elemY[2] = domain.y(nd2i); - elemY[3] = domain.y(nd3i); - elemY[4] = domain.y(nd4i); - elemY[5] = domain.y(nd5i); - elemY[6] = domain.y(nd6i); - elemY[7] = domain.y(nd7i); - - elemZ[0] = domain.z(nd0i); - elemZ[1] = domain.z(nd1i); - elemZ[2] = domain.z(nd2i); - elemZ[3] = domain.z(nd3i); - elemZ[4] = domain.z(nd4i); - elemZ[5] = domain.z(nd5i); - elemZ[6] = domain.z(nd6i); - elemZ[7] = domain.z(nd7i); + Real_t elemZ[8]) { + Index_t nd0i = elemToNode[0]; + Index_t nd1i = elemToNode[1]; + Index_t nd2i = elemToNode[2]; + Index_t nd3i = elemToNode[3]; + Index_t nd4i = elemToNode[4]; + Index_t nd5i = elemToNode[5]; + Index_t nd6i = elemToNode[6]; + Index_t nd7i = elemToNode[7]; + + elemX[0] = domain.x(nd0i); + elemX[1] = domain.x(nd1i); + elemX[2] = domain.x(nd2i); + elemX[3] = domain.x(nd3i); + elemX[4] = domain.x(nd4i); + elemX[5] = domain.x(nd5i); + elemX[6] = domain.x(nd6i); + elemX[7] = domain.x(nd7i); + + elemY[0] = domain.y(nd0i); + elemY[1] = domain.y(nd1i); + elemY[2] = domain.y(nd2i); + elemY[3] = domain.y(nd3i); + elemY[4] = domain.y(nd4i); + elemY[5] = domain.y(nd5i); + elemY[6] = domain.y(nd6i); + elemY[7] = domain.y(nd7i); + + elemZ[0] = domain.z(nd0i); + elemZ[1] = domain.z(nd1i); + elemZ[2] = domain.z(nd2i); + elemZ[3] = domain.z(nd3i); + elemZ[4] = domain.z(nd4i); + elemZ[5] = domain.z(nd5i); + elemZ[6] = domain.z(nd6i); + elemZ[7] = domain.z(nd7i); } @@ -316,109 +319,119 @@ void CollectDomainNodesToElemNodes(Domain &domain, static inline void InitStressTermsForElems(Domain &domain, Real_t *sigxx, Real_t *sigyy, Real_t *sigzz, - Index_t numElem) -{ - // - // pull in the stresses appropriate to the hydro integration - // + Index_t numElem) { + // + // pull in the stresses appropriate to the hydro integration + // #pragma omp parallel for firstprivate(numElem) - for (Index_t i = 0 ; i < numElem ; ++i){ - sigxx[i] = sigyy[i] = sigzz[i] = - domain.p(i) - domain.q(i) ; - } + for (Index_t i = 0; i < numElem; ++i) { + sigxx[i] = sigyy[i] = sigzz[i] = -domain.p(i) - domain.q(i); + } } /******************************************/ static inline -void CalcElemShapeFunctionDerivatives( Real_t const x[], - Real_t const y[], - Real_t const z[], - Real_t b[][8], - Real_t* const volume ) -{ - const Real_t x0 = x[0] ; const Real_t x1 = x[1] ; - const Real_t x2 = x[2] ; const Real_t x3 = x[3] ; - const Real_t x4 = x[4] ; const Real_t x5 = x[5] ; - const Real_t x6 = x[6] ; const Real_t x7 = x[7] ; - - const Real_t y0 = y[0] ; const Real_t y1 = y[1] ; - const Real_t y2 = y[2] ; const Real_t y3 = y[3] ; - const Real_t y4 = y[4] ; const Real_t y5 = y[5] ; - const Real_t y6 = y[6] ; const Real_t y7 = y[7] ; - - const Real_t z0 = z[0] ; const Real_t z1 = z[1] ; - const Real_t z2 = z[2] ; const Real_t z3 = z[3] ; - const Real_t z4 = z[4] ; const Real_t z5 = z[5] ; - const Real_t z6 = z[6] ; const Real_t z7 = z[7] ; - - Real_t fjxxi, fjxet, fjxze; - Real_t fjyxi, fjyet, fjyze; - Real_t fjzxi, fjzet, fjzze; - Real_t cjxxi, cjxet, cjxze; - Real_t cjyxi, cjyet, cjyze; - Real_t cjzxi, cjzet, cjzze; - - fjxxi = Real_t(.125) * ( (x6-x0) + (x5-x3) - (x7-x1) - (x4-x2) ); - fjxet = Real_t(.125) * ( (x6-x0) - (x5-x3) + (x7-x1) - (x4-x2) ); - fjxze = Real_t(.125) * ( (x6-x0) + (x5-x3) + (x7-x1) + (x4-x2) ); - - fjyxi = Real_t(.125) * ( (y6-y0) + (y5-y3) - (y7-y1) - (y4-y2) ); - fjyet = Real_t(.125) * ( (y6-y0) - (y5-y3) + (y7-y1) - (y4-y2) ); - fjyze = Real_t(.125) * ( (y6-y0) + (y5-y3) + (y7-y1) + (y4-y2) ); - - fjzxi = Real_t(.125) * ( (z6-z0) + (z5-z3) - (z7-z1) - (z4-z2) ); - fjzet = Real_t(.125) * ( (z6-z0) - (z5-z3) + (z7-z1) - (z4-z2) ); - fjzze = Real_t(.125) * ( (z6-z0) + (z5-z3) + (z7-z1) + (z4-z2) ); - - /* compute cofactors */ - cjxxi = (fjyet * fjzze) - (fjzet * fjyze); - cjxet = - (fjyxi * fjzze) + (fjzxi * fjyze); - cjxze = (fjyxi * fjzet) - (fjzxi * fjyet); - - cjyxi = - (fjxet * fjzze) + (fjzet * fjxze); - cjyet = (fjxxi * fjzze) - (fjzxi * fjxze); - cjyze = - (fjxxi * fjzet) + (fjzxi * fjxet); - - cjzxi = (fjxet * fjyze) - (fjyet * fjxze); - cjzet = - (fjxxi * fjyze) + (fjyxi * fjxze); - cjzze = (fjxxi * fjyet) - (fjyxi * fjxet); - - /* calculate partials : +void CalcElemShapeFunctionDerivatives(Real_t const x[], + Real_t const y[], + Real_t const z[], + Real_t b[][8], + Real_t *const volume) { + const Real_t x0 = x[0]; + const Real_t x1 = x[1]; + const Real_t x2 = x[2]; + const Real_t x3 = x[3]; + const Real_t x4 = x[4]; + const Real_t x5 = x[5]; + const Real_t x6 = x[6]; + const Real_t x7 = x[7]; + + const Real_t y0 = y[0]; + const Real_t y1 = y[1]; + const Real_t y2 = y[2]; + const Real_t y3 = y[3]; + const Real_t y4 = y[4]; + const Real_t y5 = y[5]; + const Real_t y6 = y[6]; + const Real_t y7 = y[7]; + + const Real_t z0 = z[0]; + const Real_t z1 = z[1]; + const Real_t z2 = z[2]; + const Real_t z3 = z[3]; + const Real_t z4 = z[4]; + const Real_t z5 = z[5]; + const Real_t z6 = z[6]; + const Real_t z7 = z[7]; + + Real_t fjxxi, fjxet, fjxze; + Real_t fjyxi, fjyet, fjyze; + Real_t fjzxi, fjzet, fjzze; + Real_t cjxxi, cjxet, cjxze; + Real_t cjyxi, cjyet, cjyze; + Real_t cjzxi, cjzet, cjzze; + + fjxxi = Real_t(.125) * ((x6 - x0) + (x5 - x3) - (x7 - x1) - (x4 - x2)); + fjxet = Real_t(.125) * ((x6 - x0) - (x5 - x3) + (x7 - x1) - (x4 - x2)); + fjxze = Real_t(.125) * ((x6 - x0) + (x5 - x3) + (x7 - x1) + (x4 - x2)); + + fjyxi = Real_t(.125) * ((y6 - y0) + (y5 - y3) - (y7 - y1) - (y4 - y2)); + fjyet = Real_t(.125) * ((y6 - y0) - (y5 - y3) + (y7 - y1) - (y4 - y2)); + fjyze = Real_t(.125) * ((y6 - y0) + (y5 - y3) + (y7 - y1) + (y4 - y2)); + + fjzxi = Real_t(.125) * ((z6 - z0) + (z5 - z3) - (z7 - z1) - (z4 - z2)); + fjzet = Real_t(.125) * ((z6 - z0) - (z5 - z3) + (z7 - z1) - (z4 - z2)); + fjzze = Real_t(.125) * ((z6 - z0) + (z5 - z3) + (z7 - z1) + (z4 - z2)); + + /* compute cofactors */ + cjxxi = (fjyet * fjzze) - (fjzet * fjyze); + cjxet = -(fjyxi * fjzze) + (fjzxi * fjyze); + cjxze = (fjyxi * fjzet) - (fjzxi * fjyet); + + cjyxi = -(fjxet * fjzze) + (fjzet * fjxze); + cjyet = (fjxxi * fjzze) - (fjzxi * fjxze); + cjyze = -(fjxxi * fjzet) + (fjzxi * fjxet); + + cjzxi = (fjxet * fjyze) - (fjyet * fjxze); + cjzet = -(fjxxi * fjyze) + (fjyxi * fjxze); + cjzze = (fjxxi * fjyet) - (fjyxi * fjxet); + + /* calculate partials : this need only be done for l = 0,1,2,3 since , by symmetry , (6,7,4,5) = - (0,1,2,3) . */ - b[0][0] = - cjxxi - cjxet - cjxze; - b[0][1] = cjxxi - cjxet - cjxze; - b[0][2] = cjxxi + cjxet - cjxze; - b[0][3] = - cjxxi + cjxet - cjxze; - b[0][4] = -b[0][2]; - b[0][5] = -b[0][3]; - b[0][6] = -b[0][0]; - b[0][7] = -b[0][1]; - - b[1][0] = - cjyxi - cjyet - cjyze; - b[1][1] = cjyxi - cjyet - cjyze; - b[1][2] = cjyxi + cjyet - cjyze; - b[1][3] = - cjyxi + cjyet - cjyze; - b[1][4] = -b[1][2]; - b[1][5] = -b[1][3]; - b[1][6] = -b[1][0]; - b[1][7] = -b[1][1]; - - b[2][0] = - cjzxi - cjzet - cjzze; - b[2][1] = cjzxi - cjzet - cjzze; - b[2][2] = cjzxi + cjzet - cjzze; - b[2][3] = - cjzxi + cjzet - cjzze; - b[2][4] = -b[2][2]; - b[2][5] = -b[2][3]; - b[2][6] = -b[2][0]; - b[2][7] = -b[2][1]; - - /* calculate jacobian determinant (volume) */ - //laik_log((Laik_LogLevel)2,"Debug %f\n", *volume); - - *volume = Real_t(8.) * ( fjxet * cjxet + fjyet * cjyet + fjzet * cjzet); + b[0][0] = -cjxxi - cjxet - cjxze; + b[0][1] = cjxxi - cjxet - cjxze; + b[0][2] = cjxxi + cjxet - cjxze; + b[0][3] = -cjxxi + cjxet - cjxze; + b[0][4] = -b[0][2]; + b[0][5] = -b[0][3]; + b[0][6] = -b[0][0]; + b[0][7] = -b[0][1]; + + b[1][0] = -cjyxi - cjyet - cjyze; + b[1][1] = cjyxi - cjyet - cjyze; + b[1][2] = cjyxi + cjyet - cjyze; + b[1][3] = -cjyxi + cjyet - cjyze; + b[1][4] = -b[1][2]; + b[1][5] = -b[1][3]; + b[1][6] = -b[1][0]; + b[1][7] = -b[1][1]; + + b[2][0] = -cjzxi - cjzet - cjzze; + b[2][1] = cjzxi - cjzet - cjzze; + b[2][2] = cjzxi + cjzet - cjzze; + b[2][3] = -cjzxi + cjzet - cjzze; + b[2][4] = -b[2][2]; + b[2][5] = -b[2][3]; + b[2][6] = -b[2][0]; + b[2][7] = -b[2][1]; + + /* calculate jacobian determinant (volume) */ + //laik_log((Laik_LogLevel)2,"Debug %f\n", *volume); + + *volume = Real_t(8.) * (fjxet * cjxet + fjyet * cjyet + fjzet * cjzet); } /******************************************/ @@ -431,32 +444,31 @@ void SumElemFaceNormal(Real_t *normalX0, Real_t *normalY0, Real_t *normalZ0, const Real_t x0, const Real_t y0, const Real_t z0, const Real_t x1, const Real_t y1, const Real_t z1, const Real_t x2, const Real_t y2, const Real_t z2, - const Real_t x3, const Real_t y3, const Real_t z3) -{ - Real_t bisectX0 = Real_t(0.5) * (x3 + x2 - x1 - x0); - Real_t bisectY0 = Real_t(0.5) * (y3 + y2 - y1 - y0); - Real_t bisectZ0 = Real_t(0.5) * (z3 + z2 - z1 - z0); - Real_t bisectX1 = Real_t(0.5) * (x2 + x1 - x3 - x0); - Real_t bisectY1 = Real_t(0.5) * (y2 + y1 - y3 - y0); - Real_t bisectZ1 = Real_t(0.5) * (z2 + z1 - z3 - z0); - Real_t areaX = Real_t(0.25) * (bisectY0 * bisectZ1 - bisectZ0 * bisectY1); - Real_t areaY = Real_t(0.25) * (bisectZ0 * bisectX1 - bisectX0 * bisectZ1); - Real_t areaZ = Real_t(0.25) * (bisectX0 * bisectY1 - bisectY0 * bisectX1); - - *normalX0 += areaX; - *normalX1 += areaX; - *normalX2 += areaX; - *normalX3 += areaX; - - *normalY0 += areaY; - *normalY1 += areaY; - *normalY2 += areaY; - *normalY3 += areaY; - - *normalZ0 += areaZ; - *normalZ1 += areaZ; - *normalZ2 += areaZ; - *normalZ3 += areaZ; + const Real_t x3, const Real_t y3, const Real_t z3) { + Real_t bisectX0 = Real_t(0.5) * (x3 + x2 - x1 - x0); + Real_t bisectY0 = Real_t(0.5) * (y3 + y2 - y1 - y0); + Real_t bisectZ0 = Real_t(0.5) * (z3 + z2 - z1 - z0); + Real_t bisectX1 = Real_t(0.5) * (x2 + x1 - x3 - x0); + Real_t bisectY1 = Real_t(0.5) * (y2 + y1 - y3 - y0); + Real_t bisectZ1 = Real_t(0.5) * (z2 + z1 - z3 - z0); + Real_t areaX = Real_t(0.25) * (bisectY0 * bisectZ1 - bisectZ0 * bisectY1); + Real_t areaY = Real_t(0.25) * (bisectZ0 * bisectX1 - bisectX0 * bisectZ1); + Real_t areaZ = Real_t(0.25) * (bisectX0 * bisectY1 - bisectY0 * bisectX1); + + *normalX0 += areaX; + *normalX1 += areaX; + *normalX2 += areaX; + *normalX3 += areaX; + + *normalY0 += areaY; + *normalY1 += areaY; + *normalY2 += areaY; + *normalY3 += areaY; + + *normalZ0 += areaZ; + *normalZ1 += areaZ; + *normalZ2 += areaZ; + *normalZ3 += areaZ; } /******************************************/ @@ -467,168 +479,162 @@ void CalcElemNodeNormals(Real_t pfx[8], Real_t pfz[8], const Real_t x[8], const Real_t y[8], - const Real_t z[8]) -{ - for (Index_t i = 0 ; i < 8 ; ++i) { - pfx[i] = Real_t(0.0); - pfy[i] = Real_t(0.0); - pfz[i] = Real_t(0.0); - } - /* evaluate face one: nodes 0, 1, 2, 3 */ - SumElemFaceNormal(&pfx[0], &pfy[0], &pfz[0], - &pfx[1], &pfy[1], &pfz[1], - &pfx[2], &pfy[2], &pfz[2], - &pfx[3], &pfy[3], &pfz[3], - x[0], y[0], z[0], x[1], y[1], z[1], - x[2], y[2], z[2], x[3], y[3], z[3]); - /* evaluate face two: nodes 0, 4, 5, 1 */ - SumElemFaceNormal(&pfx[0], &pfy[0], &pfz[0], - &pfx[4], &pfy[4], &pfz[4], - &pfx[5], &pfy[5], &pfz[5], - &pfx[1], &pfy[1], &pfz[1], - x[0], y[0], z[0], x[4], y[4], z[4], - x[5], y[5], z[5], x[1], y[1], z[1]); - /* evaluate face three: nodes 1, 5, 6, 2 */ - SumElemFaceNormal(&pfx[1], &pfy[1], &pfz[1], - &pfx[5], &pfy[5], &pfz[5], - &pfx[6], &pfy[6], &pfz[6], - &pfx[2], &pfy[2], &pfz[2], - x[1], y[1], z[1], x[5], y[5], z[5], - x[6], y[6], z[6], x[2], y[2], z[2]); - /* evaluate face four: nodes 2, 6, 7, 3 */ - SumElemFaceNormal(&pfx[2], &pfy[2], &pfz[2], - &pfx[6], &pfy[6], &pfz[6], - &pfx[7], &pfy[7], &pfz[7], - &pfx[3], &pfy[3], &pfz[3], - x[2], y[2], z[2], x[6], y[6], z[6], - x[7], y[7], z[7], x[3], y[3], z[3]); - /* evaluate face five: nodes 3, 7, 4, 0 */ - SumElemFaceNormal(&pfx[3], &pfy[3], &pfz[3], - &pfx[7], &pfy[7], &pfz[7], - &pfx[4], &pfy[4], &pfz[4], - &pfx[0], &pfy[0], &pfz[0], - x[3], y[3], z[3], x[7], y[7], z[7], - x[4], y[4], z[4], x[0], y[0], z[0]); - /* evaluate face six: nodes 4, 7, 6, 5 */ - SumElemFaceNormal(&pfx[4], &pfy[4], &pfz[4], - &pfx[7], &pfy[7], &pfz[7], - &pfx[6], &pfy[6], &pfz[6], - &pfx[5], &pfy[5], &pfz[5], - x[4], y[4], z[4], x[7], y[7], z[7], - x[6], y[6], z[6], x[5], y[5], z[5]); + const Real_t z[8]) { + for (Index_t i = 0; i < 8; ++i) { + pfx[i] = Real_t(0.0); + pfy[i] = Real_t(0.0); + pfz[i] = Real_t(0.0); + } + /* evaluate face one: nodes 0, 1, 2, 3 */ + SumElemFaceNormal(&pfx[0], &pfy[0], &pfz[0], + &pfx[1], &pfy[1], &pfz[1], + &pfx[2], &pfy[2], &pfz[2], + &pfx[3], &pfy[3], &pfz[3], + x[0], y[0], z[0], x[1], y[1], z[1], + x[2], y[2], z[2], x[3], y[3], z[3]); + /* evaluate face two: nodes 0, 4, 5, 1 */ + SumElemFaceNormal(&pfx[0], &pfy[0], &pfz[0], + &pfx[4], &pfy[4], &pfz[4], + &pfx[5], &pfy[5], &pfz[5], + &pfx[1], &pfy[1], &pfz[1], + x[0], y[0], z[0], x[4], y[4], z[4], + x[5], y[5], z[5], x[1], y[1], z[1]); + /* evaluate face three: nodes 1, 5, 6, 2 */ + SumElemFaceNormal(&pfx[1], &pfy[1], &pfz[1], + &pfx[5], &pfy[5], &pfz[5], + &pfx[6], &pfy[6], &pfz[6], + &pfx[2], &pfy[2], &pfz[2], + x[1], y[1], z[1], x[5], y[5], z[5], + x[6], y[6], z[6], x[2], y[2], z[2]); + /* evaluate face four: nodes 2, 6, 7, 3 */ + SumElemFaceNormal(&pfx[2], &pfy[2], &pfz[2], + &pfx[6], &pfy[6], &pfz[6], + &pfx[7], &pfy[7], &pfz[7], + &pfx[3], &pfy[3], &pfz[3], + x[2], y[2], z[2], x[6], y[6], z[6], + x[7], y[7], z[7], x[3], y[3], z[3]); + /* evaluate face five: nodes 3, 7, 4, 0 */ + SumElemFaceNormal(&pfx[3], &pfy[3], &pfz[3], + &pfx[7], &pfy[7], &pfz[7], + &pfx[4], &pfy[4], &pfz[4], + &pfx[0], &pfy[0], &pfz[0], + x[3], y[3], z[3], x[7], y[7], z[7], + x[4], y[4], z[4], x[0], y[0], z[0]); + /* evaluate face six: nodes 4, 7, 6, 5 */ + SumElemFaceNormal(&pfx[4], &pfy[4], &pfz[4], + &pfx[7], &pfy[7], &pfz[7], + &pfx[6], &pfy[6], &pfz[6], + &pfx[5], &pfy[5], &pfz[5], + x[4], y[4], z[4], x[7], y[7], z[7], + x[6], y[6], z[6], x[5], y[5], z[5]); } /******************************************/ static inline -void SumElemStressesToNodeForces( const Real_t B[][8], - const Real_t stress_xx, - const Real_t stress_yy, - const Real_t stress_zz, - Real_t fx[], Real_t fy[], Real_t fz[] ) -{ - for(Index_t i = 0; i < 8; i++) { - fx[i] = -( stress_xx * B[0][i] ); - fy[i] = -( stress_yy * B[1][i] ); - fz[i] = -( stress_zz * B[2][i] ); - } +void SumElemStressesToNodeForces(const Real_t B[][8], + const Real_t stress_xx, + const Real_t stress_yy, + const Real_t stress_zz, + Real_t fx[], Real_t fy[], Real_t fz[]) { + for (Index_t i = 0; i < 8; i++) { + fx[i] = -(stress_xx * B[0][i]); + fy[i] = -(stress_yy * B[1][i]); + fz[i] = -(stress_zz * B[2][i]); + } } /******************************************/ static inline -void IntegrateStressForElems( Domain &domain, - Real_t *sigxx, Real_t *sigyy, Real_t *sigzz, - Real_t *determ, Index_t numElem, Index_t numNode) -{ +void IntegrateStressForElems(Domain &domain, + Real_t *sigxx, Real_t *sigyy, Real_t *sigzz, + Real_t *determ, Index_t numElem, Index_t numNode) { #if _OPENMP - Index_t numthreads = omp_get_max_threads(); + Index_t numthreads = omp_get_max_threads(); #else - Index_t numthreads = 1; + Index_t numthreads = 1; #endif - Index_t numElem8 = numElem * 8 ; - Real_t *fx_elem; - Real_t *fy_elem; - Real_t *fz_elem; - Real_t fx_local[8] ; - Real_t fy_local[8] ; - Real_t fz_local[8] ; + Index_t numElem8 = numElem * 8; + Real_t *fx_elem; + Real_t *fy_elem; + Real_t *fz_elem; + Real_t fx_local[8]; + Real_t fy_local[8]; + Real_t fz_local[8]; - if (numthreads > 1) { - fx_elem = Allocate(numElem8) ; - fy_elem = Allocate(numElem8) ; - fz_elem = Allocate(numElem8) ; - } - // loop over all elements + if (numthreads > 1) { + fx_elem = Allocate(numElem8); + fy_elem = Allocate(numElem8); + fz_elem = Allocate(numElem8); + } + // loop over all elements #pragma omp parallel for firstprivate(numElem) - for( Index_t k=0 ; k 1) { - // Eliminate thread writing conflicts at the nodes by giving - // each element its own copy to write to - SumElemStressesToNodeForces( B, sigxx[k], sigyy[k], sigzz[k], - &fx_elem[k*8], - &fy_elem[k*8], - &fz_elem[k*8] ) ; - } - else { - SumElemStressesToNodeForces( B, sigxx[k], sigyy[k], sigzz[k], - fx_local, fy_local, fz_local ) ; - - // copy nodal force contributions to global force arrray. - for( Index_t lnode=0 ; lnode<8 ; ++lnode ) { - Index_t gnode = elemToNode[lnode]; - domain.fx(gnode) += fx_local[lnode]; - domain.fy(gnode) += fy_local[lnode]; - domain.fz(gnode) += fz_local[lnode]; - } + CalcElemNodeNormals(B[0], B[1], B[2], + x_local, y_local, z_local); + + if (numthreads > 1) { + // Eliminate thread writing conflicts at the nodes by giving + // each element its own copy to write to + SumElemStressesToNodeForces(B, sigxx[k], sigyy[k], sigzz[k], + &fx_elem[k * 8], + &fy_elem[k * 8], + &fz_elem[k * 8]); + } else { + SumElemStressesToNodeForces(B, sigxx[k], sigyy[k], sigzz[k], + fx_local, fy_local, fz_local); + + // copy nodal force contributions to global force arrray. + for (Index_t lnode = 0; lnode < 8; ++lnode) { + Index_t gnode = elemToNode[lnode]; + domain.fx(gnode) += fx_local[lnode]; + domain.fy(gnode) += fy_local[lnode]; + domain.fz(gnode) += fz_local[lnode]; + } + } } - } - if (numthreads > 1) { - // If threaded, then we need to copy the data out of the temporary - // arrays used above into the final forces field + if (numthreads > 1) { + // If threaded, then we need to copy the data out of the temporary + // arrays used above into the final forces field #pragma omp parallel for firstprivate(numNode) - for( Index_t gnode=0 ; gnode 1) { - fx_elem = Allocate(numElem8) ; - fy_elem = Allocate(numElem8) ; - fz_elem = Allocate(numElem8) ; - } + Real_t *fx_elem; + Real_t *fy_elem; + Real_t *fz_elem; - Real_t gamma[4][8]; - - gamma[0][0] = Real_t( 1.); - gamma[0][1] = Real_t( 1.); - gamma[0][2] = Real_t(-1.); - gamma[0][3] = Real_t(-1.); - gamma[0][4] = Real_t(-1.); - gamma[0][5] = Real_t(-1.); - gamma[0][6] = Real_t( 1.); - gamma[0][7] = Real_t( 1.); - gamma[1][0] = Real_t( 1.); - gamma[1][1] = Real_t(-1.); - gamma[1][2] = Real_t(-1.); - gamma[1][3] = Real_t( 1.); - gamma[1][4] = Real_t(-1.); - gamma[1][5] = Real_t( 1.); - gamma[1][6] = Real_t( 1.); - gamma[1][7] = Real_t(-1.); - gamma[2][0] = Real_t( 1.); - gamma[2][1] = Real_t(-1.); - gamma[2][2] = Real_t( 1.); - gamma[2][3] = Real_t(-1.); - gamma[2][4] = Real_t( 1.); - gamma[2][5] = Real_t(-1.); - gamma[2][6] = Real_t( 1.); - gamma[2][7] = Real_t(-1.); - gamma[3][0] = Real_t(-1.); - gamma[3][1] = Real_t( 1.); - gamma[3][2] = Real_t(-1.); - gamma[3][3] = Real_t( 1.); - gamma[3][4] = Real_t( 1.); - gamma[3][5] = Real_t(-1.); - gamma[3][6] = Real_t( 1.); - gamma[3][7] = Real_t(-1.); + if (numthreads > 1) { + fx_elem = Allocate(numElem8); + fy_elem = Allocate(numElem8); + fz_elem = Allocate(numElem8); + } + + Real_t gamma[4][8]; + + gamma[0][0] = Real_t(1.); + gamma[0][1] = Real_t(1.); + gamma[0][2] = Real_t(-1.); + gamma[0][3] = Real_t(-1.); + gamma[0][4] = Real_t(-1.); + gamma[0][5] = Real_t(-1.); + gamma[0][6] = Real_t(1.); + gamma[0][7] = Real_t(1.); + gamma[1][0] = Real_t(1.); + gamma[1][1] = Real_t(-1.); + gamma[1][2] = Real_t(-1.); + gamma[1][3] = Real_t(1.); + gamma[1][4] = Real_t(-1.); + gamma[1][5] = Real_t(1.); + gamma[1][6] = Real_t(1.); + gamma[1][7] = Real_t(-1.); + gamma[2][0] = Real_t(1.); + gamma[2][1] = Real_t(-1.); + gamma[2][2] = Real_t(1.); + gamma[2][3] = Real_t(-1.); + gamma[2][4] = Real_t(1.); + gamma[2][5] = Real_t(-1.); + gamma[2][6] = Real_t(1.); + gamma[2][7] = Real_t(-1.); + gamma[3][0] = Real_t(-1.); + gamma[3][1] = Real_t(1.); + gamma[3][2] = Real_t(-1.); + gamma[3][3] = Real_t(1.); + gamma[3][4] = Real_t(1.); + gamma[3][5] = Real_t(-1.); + gamma[3][6] = Real_t(1.); + gamma[3][7] = Real_t(-1.); /*************************************************/ /* compute the hourglass modes */ #pragma omp parallel for firstprivate(numElem, hourg) - for(Index_t i2=0;i2 1) { - fx_local = &fx_elem[i3] ; - fx_local[0] = hgfx[0]; - fx_local[1] = hgfx[1]; - fx_local[2] = hgfx[2]; - fx_local[3] = hgfx[3]; - fx_local[4] = hgfx[4]; - fx_local[5] = hgfx[5]; - fx_local[6] = hgfx[6]; - fx_local[7] = hgfx[7]; - - fy_local = &fy_elem[i3] ; - fy_local[0] = hgfy[0]; - fy_local[1] = hgfy[1]; - fy_local[2] = hgfy[2]; - fy_local[3] = hgfy[3]; - fy_local[4] = hgfy[4]; - fy_local[5] = hgfy[5]; - fy_local[6] = hgfy[6]; - fy_local[7] = hgfy[7]; - - fz_local = &fz_elem[i3] ; - fz_local[0] = hgfz[0]; - fz_local[1] = hgfz[1]; - fz_local[2] = hgfz[2]; - fz_local[3] = hgfz[3]; - fz_local[4] = hgfz[4]; - fz_local[5] = hgfz[5]; - fz_local[6] = hgfz[6]; - fz_local[7] = hgfz[7]; - } - else { - domain.fx(n0si2) += hgfx[0]; - domain.fy(n0si2) += hgfy[0]; - domain.fz(n0si2) += hgfz[0]; - - domain.fx(n1si2) += hgfx[1]; - domain.fy(n1si2) += hgfy[1]; - domain.fz(n1si2) += hgfz[1]; - - domain.fx(n2si2) += hgfx[2]; - domain.fy(n2si2) += hgfy[2]; - domain.fz(n2si2) += hgfz[2]; - - domain.fx(n3si2) += hgfx[3]; - domain.fy(n3si2) += hgfy[3]; - domain.fz(n3si2) += hgfz[3]; - - domain.fx(n4si2) += hgfx[4]; - domain.fy(n4si2) += hgfy[4]; - domain.fz(n4si2) += hgfz[4]; - - domain.fx(n5si2) += hgfx[5]; - domain.fy(n5si2) += hgfy[5]; - domain.fz(n5si2) += hgfz[5]; - - domain.fx(n6si2) += hgfx[6]; - domain.fy(n6si2) += hgfy[6]; - domain.fz(n6si2) += hgfz[6]; - - domain.fx(n7si2) += hgfx[7]; - domain.fy(n7si2) += hgfy[7]; - domain.fz(n7si2) += hgfz[7]; - } - } + /* compute forces */ + /* store forces into h arrays (force arrays) */ + + ss1 = domain.ss(i2); + mass1 = domain.elemMass(i2); + volume13 = CBRT(determ[i2]); + + Index_t n0si2 = elemToNode[0]; + Index_t n1si2 = elemToNode[1]; + Index_t n2si2 = elemToNode[2]; + Index_t n3si2 = elemToNode[3]; + Index_t n4si2 = elemToNode[4]; + Index_t n5si2 = elemToNode[5]; + Index_t n6si2 = elemToNode[6]; + Index_t n7si2 = elemToNode[7]; + + xd1[0] = domain.xd(n0si2); + xd1[1] = domain.xd(n1si2); + xd1[2] = domain.xd(n2si2); + xd1[3] = domain.xd(n3si2); + xd1[4] = domain.xd(n4si2); + xd1[5] = domain.xd(n5si2); + xd1[6] = domain.xd(n6si2); + xd1[7] = domain.xd(n7si2); + + yd1[0] = domain.yd(n0si2); + yd1[1] = domain.yd(n1si2); + yd1[2] = domain.yd(n2si2); + yd1[3] = domain.yd(n3si2); + yd1[4] = domain.yd(n4si2); + yd1[5] = domain.yd(n5si2); + yd1[6] = domain.yd(n6si2); + yd1[7] = domain.yd(n7si2); + + zd1[0] = domain.zd(n0si2); + zd1[1] = domain.zd(n1si2); + zd1[2] = domain.zd(n2si2); + zd1[3] = domain.zd(n3si2); + zd1[4] = domain.zd(n4si2); + zd1[5] = domain.zd(n5si2); + zd1[6] = domain.zd(n6si2); + zd1[7] = domain.zd(n7si2); + + coefficient = -hourg * Real_t(0.01) * ss1 * mass1 / volume13; + + CalcElemFBHourglassForce(xd1, yd1, zd1, + hourgam, + coefficient, hgfx, hgfy, hgfz); + + // With the threaded version, we write into local arrays per elem + // so we don't have to worry about race conditions + if (numthreads > 1) { + fx_local = &fx_elem[i3]; + fx_local[0] = hgfx[0]; + fx_local[1] = hgfx[1]; + fx_local[2] = hgfx[2]; + fx_local[3] = hgfx[3]; + fx_local[4] = hgfx[4]; + fx_local[5] = hgfx[5]; + fx_local[6] = hgfx[6]; + fx_local[7] = hgfx[7]; + + fy_local = &fy_elem[i3]; + fy_local[0] = hgfy[0]; + fy_local[1] = hgfy[1]; + fy_local[2] = hgfy[2]; + fy_local[3] = hgfy[3]; + fy_local[4] = hgfy[4]; + fy_local[5] = hgfy[5]; + fy_local[6] = hgfy[6]; + fy_local[7] = hgfy[7]; + + fz_local = &fz_elem[i3]; + fz_local[0] = hgfz[0]; + fz_local[1] = hgfz[1]; + fz_local[2] = hgfz[2]; + fz_local[3] = hgfz[3]; + fz_local[4] = hgfz[4]; + fz_local[5] = hgfz[5]; + fz_local[6] = hgfz[6]; + fz_local[7] = hgfz[7]; + } else { + domain.fx(n0si2) += hgfx[0]; + domain.fy(n0si2) += hgfy[0]; + domain.fz(n0si2) += hgfz[0]; + + domain.fx(n1si2) += hgfx[1]; + domain.fy(n1si2) += hgfy[1]; + domain.fz(n1si2) += hgfz[1]; + + domain.fx(n2si2) += hgfx[2]; + domain.fy(n2si2) += hgfy[2]; + domain.fz(n2si2) += hgfz[2]; + + domain.fx(n3si2) += hgfx[3]; + domain.fy(n3si2) += hgfy[3]; + domain.fz(n3si2) += hgfz[3]; + + domain.fx(n4si2) += hgfx[4]; + domain.fy(n4si2) += hgfy[4]; + domain.fz(n4si2) += hgfz[4]; + + domain.fx(n5si2) += hgfx[5]; + domain.fy(n5si2) += hgfy[5]; + domain.fz(n5si2) += hgfz[5]; + + domain.fx(n6si2) += hgfx[6]; + domain.fy(n6si2) += hgfy[6]; + domain.fz(n6si2) += hgfz[6]; + + domain.fx(n7si2) += hgfx[7]; + domain.fy(n7si2) += hgfy[7]; + domain.fz(n7si2) += hgfz[7]; + } + } - if (numthreads > 1) { - // Collect the data from the local arrays into the final force arrays + if (numthreads > 1) { + // Collect the data from the local arrays into the final force arrays #pragma omp parallel for firstprivate(numNode) - for( Index_t gnode=0 ; gnode(numElem8) ; - Real_t *dvdy = Allocate(numElem8) ; - Real_t *dvdz = Allocate(numElem8) ; - Real_t *x8n = Allocate(numElem8) ; - Real_t *y8n = Allocate(numElem8) ; - Real_t *z8n = Allocate(numElem8) ; - - /* start loop over elements */ +void CalcHourglassControlForElems(Domain &domain, + Real_t determ[], Real_t hgcoef) { + Index_t numElem = domain.numElem(); + Index_t numElem8 = numElem * 8; + Real_t *dvdx = Allocate(numElem8); + Real_t *dvdy = Allocate(numElem8); + Real_t *dvdz = Allocate(numElem8); + Real_t *x8n = Allocate(numElem8); + Real_t *y8n = Allocate(numElem8); + Real_t *z8n = Allocate(numElem8); + + /* start loop over elements */ #pragma omp parallel for firstprivate(numElem) - for (Index_t i=0 ; i Real_t(0.) ) { - CalcFBHourglassForceForElems( domain, - determ, x8n, y8n, z8n, dvdx, dvdy, dvdz, - hgcoef, numElem, domain.numNode()) ; - } + if (hgcoef > Real_t(0.)) { + CalcFBHourglassForceForElems(domain, + determ, x8n, y8n, z8n, dvdx, dvdy, dvdz, + hgcoef, numElem, domain.numNode()); + } - Release(&z8n) ; - Release(&y8n) ; - Release(&x8n) ; - Release(&dvdz) ; - Release(&dvdy) ; - Release(&dvdx) ; + Release(&z8n); + Release(&y8n); + Release(&x8n); + Release(&dvdz); + Release(&dvdy); + Release(&dvdx); - return ; + return; } /******************************************/ static inline -void CalcVolumeForceForElems(Domain& domain) -{ - Index_t numElem = domain.numElem() ; - if (numElem != 0) { - Real_t hgcoef = domain.hgcoef() ; - Real_t *sigxx = Allocate(numElem) ; - Real_t *sigyy = Allocate(numElem) ; - Real_t *sigzz = Allocate(numElem) ; - Real_t *determ = Allocate(numElem) ; - - /* Sum contributions to total stress tensor */ - InitStressTermsForElems(domain, sigxx, sigyy, sigzz, numElem); - - // call elemlib stress integration loop to produce nodal forces from - // material stresses. - IntegrateStressForElems( domain, - sigxx, sigyy, sigzz, determ, numElem, - domain.numNode()) ; - - /* +void CalcVolumeForceForElems(Domain &domain) { + Index_t numElem = domain.numElem(); + if (numElem != 0) { + Real_t hgcoef = domain.hgcoef(); + Real_t *sigxx = Allocate(numElem); + Real_t *sigyy = Allocate(numElem); + Real_t *sigzz = Allocate(numElem); + Real_t *determ = Allocate(numElem); + + /* Sum contributions to total stress tensor */ + InitStressTermsForElems(domain, sigxx, sigyy, sigzz, numElem); + + // call elemlib stress integration loop to produce nodal forces from + // material stresses. + IntegrateStressForElems(domain, + sigxx, sigyy, sigzz, determ, numElem, + domain.numNode()); + + /* if (domain.cycle() == 6){ laik_log((Laik_LogLevel)2,"Debug: numNodes: %d\n",domain.numNode()); for (int k = 0; k < numElem; ++k) { @@ -1134,64 +1132,63 @@ void CalcVolumeForceForElems(Domain& domain) } */ - // check for negative element volume + // check for negative element volume #pragma omp parallel for firstprivate(numElem) - for ( Index_t k=0 ; k 0) { - const Real_t deltatime = domain.deltatime() ; + Index_t numElem = domain.numElem(); + if (numElem > 0) { + const Real_t deltatime = domain.deltatime(); - //domain.AllocateStrains(numElem); + //domain.AllocateStrains(numElem); - CalcKinematicsForElems(domain, vnew, deltatime, numElem) ; + CalcKinematicsForElems(domain, vnew, deltatime, numElem); - // element loop to do some stuff not included in the elemlib function. + // element loop to do some stuff not included in the elemlib function. #pragma omp parallel for firstprivate(numElem) - for ( Index_t k=0 ; k monoq_max_slope) phixi = monoq_max_slope; - - - /* phieta */ - norm = Real_t(1.) / ( domain.delv_eta(i) + ptiny ) ; - - switch (bcMask & ETA_M) { - case ETA_M_COMM: /* needs comm data */ - case 0: delvm = domain.delv_eta(domain.letam(i)) ; break ; - case ETA_M_SYMM: delvm = domain.delv_eta(i) ; break ; - case ETA_M_FREE: delvm = Real_t(0.0) ; break ; - default: fprintf(stderr, "Error in switch at %s line %d\n", - __FILE__, __LINE__); - delvm = 0; /* ERROR - but quiets the compiler */ - break; - } - switch (bcMask & ETA_P) { - case ETA_P_COMM: /* needs comm data */ - case 0: delvp = domain.delv_eta(domain.letap(i)) ; break ; - case ETA_P_SYMM: delvp = domain.delv_eta(i) ; break ; - case ETA_P_FREE: delvp = Real_t(0.0) ; break ; - default: fprintf(stderr, "Error in switch at %s line %d\n", - __FILE__, __LINE__); - delvp = 0; /* ERROR - but quiets the compiler */ - break; - } - - delvm = delvm * norm ; - delvp = delvp * norm ; + delvm = delvm * norm; + delvp = delvp * norm; + + phixi = Real_t(.5) * (delvm + delvp); + + delvm *= monoq_limiter_mult; + delvp *= monoq_limiter_mult; + + if (delvm < phixi) phixi = delvm; + if (delvp < phixi) phixi = delvp; + if (phixi < Real_t(0.)) phixi = Real_t(0.); + if (phixi > monoq_max_slope) phixi = monoq_max_slope; + + + /* phieta */ + norm = Real_t(1.) / (domain.delv_eta(i) + ptiny); + + switch (bcMask & ETA_M) { + case ETA_M_COMM: /* needs comm data */ + case 0: + delvm = domain.delv_eta(domain.letam(i)); + break; + case ETA_M_SYMM: + delvm = domain.delv_eta(i); + break; + case ETA_M_FREE: + delvm = Real_t(0.0); + break; + default: + fprintf(stderr, "Error in switch at %s line %d\n", + __FILE__, __LINE__); + delvm = 0; /* ERROR - but quiets the compiler */ + break; + } + switch (bcMask & ETA_P) { + case ETA_P_COMM: /* needs comm data */ + case 0: + delvp = domain.delv_eta(domain.letap(i)); + break; + case ETA_P_SYMM: + delvp = domain.delv_eta(i); + break; + case ETA_P_FREE: + delvp = Real_t(0.0); + break; + default: + fprintf(stderr, "Error in switch at %s line %d\n", + __FILE__, __LINE__); + delvp = 0; /* ERROR - but quiets the compiler */ + break; + } - phieta = Real_t(.5) * ( delvm + delvp ) ; + delvm = delvm * norm; + delvp = delvp * norm; + + phieta = Real_t(.5) * (delvm + delvp); + + delvm *= monoq_limiter_mult; + delvp *= monoq_limiter_mult; + + if (delvm < phieta) phieta = delvm; + if (delvp < phieta) phieta = delvp; + if (phieta < Real_t(0.)) phieta = Real_t(0.); + if (phieta > monoq_max_slope) phieta = monoq_max_slope; + + /* phizeta */ + norm = Real_t(1.) / (domain.delv_zeta(i) + ptiny); + + switch (bcMask & ZETA_M) { + case ZETA_M_COMM: /* needs comm data */ + case 0: + delvm = domain.delv_zeta(domain.lzetam(i)); + break; + case ZETA_M_SYMM: + delvm = domain.delv_zeta(i); + break; + case ZETA_M_FREE: + delvm = Real_t(0.0); + break; + default: + fprintf(stderr, "Error in switch at %s line %d\n", + __FILE__, __LINE__); + delvm = 0; /* ERROR - but quiets the compiler */ + break; + } + switch (bcMask & ZETA_P) { + case ZETA_P_COMM: /* needs comm data */ + case 0: + delvp = domain.delv_zeta(domain.lzetap(i)); + break; + case ZETA_P_SYMM: + delvp = domain.delv_zeta(i); + break; + case ZETA_P_FREE: + delvp = Real_t(0.0); + break; + default: + fprintf(stderr, "Error in switch at %s line %d\n", + __FILE__, __LINE__); + delvp = 0; /* ERROR - but quiets the compiler */ + break; + } - delvm *= monoq_limiter_mult ; - delvp *= monoq_limiter_mult ; + delvm = delvm * norm; + delvp = delvp * norm; - if ( delvm < phieta ) phieta = delvm ; - if ( delvp < phieta ) phieta = delvp ; - if ( phieta < Real_t(0.)) phieta = Real_t(0.) ; - if ( phieta > monoq_max_slope) phieta = monoq_max_slope; + phizeta = Real_t(.5) * (delvm + delvp); - /* phizeta */ - norm = Real_t(1.) / ( domain.delv_zeta(i) + ptiny ) ; + delvm *= monoq_limiter_mult; + delvp *= monoq_limiter_mult; - switch (bcMask & ZETA_M) { - case ZETA_M_COMM: /* needs comm data */ - case 0: delvm = domain.delv_zeta(domain.lzetam(i)) ; break ; - case ZETA_M_SYMM: delvm = domain.delv_zeta(i) ; break ; - case ZETA_M_FREE: delvm = Real_t(0.0) ; break ; - default: fprintf(stderr, "Error in switch at %s line %d\n", - __FILE__, __LINE__); - delvm = 0; /* ERROR - but quiets the compiler */ - break; - } - switch (bcMask & ZETA_P) { - case ZETA_P_COMM: /* needs comm data */ - case 0: delvp = domain.delv_zeta(domain.lzetap(i)) ; break ; - case ZETA_P_SYMM: delvp = domain.delv_zeta(i) ; break ; - case ZETA_P_FREE: delvp = Real_t(0.0) ; break ; - default: fprintf(stderr, "Error in switch at %s line %d\n", - __FILE__, __LINE__); - delvp = 0; /* ERROR - but quiets the compiler */ - break; - } + if (delvm < phizeta) phizeta = delvm; + if (delvp < phizeta) phizeta = delvp; + if (phizeta < Real_t(0.)) phizeta = Real_t(0.); + if (phizeta > monoq_max_slope) phizeta = monoq_max_slope; - delvm = delvm * norm ; - delvp = delvp * norm ; + /* Remove length scale */ - phizeta = Real_t(.5) * ( delvm + delvp ) ; + if (domain.vdov(i) > Real_t(0.)) { + qlin = Real_t(0.); + qquad = Real_t(0.); + } else { + Real_t delvxxi = domain.delv_xi(i) * domain.delx_xi(i); + Real_t delvxeta = domain.delv_eta(i) * domain.delx_eta(i); + Real_t delvxzeta = domain.delv_zeta(i) * domain.delx_zeta(i); - delvm *= monoq_limiter_mult ; - delvp *= monoq_limiter_mult ; + if (delvxxi > Real_t(0.)) delvxxi = Real_t(0.); + if (delvxeta > Real_t(0.)) delvxeta = Real_t(0.); + if (delvxzeta > Real_t(0.)) delvxzeta = Real_t(0.); - if ( delvm < phizeta ) phizeta = delvm ; - if ( delvp < phizeta ) phizeta = delvp ; - if ( phizeta < Real_t(0.)) phizeta = Real_t(0.); - if ( phizeta > monoq_max_slope ) phizeta = monoq_max_slope; + Real_t rho = domain.elemMass(i) / (domain.volo(i) * vnew[i]); - /* Remove length scale */ + qlin = -qlc_monoq * rho * + (delvxxi * (Real_t(1.) - phixi) + + delvxeta * (Real_t(1.) - phieta) + + delvxzeta * (Real_t(1.) - phizeta)); - if ( domain.vdov(i) > Real_t(0.) ) { - qlin = Real_t(0.) ; - qquad = Real_t(0.) ; - } - else { - Real_t delvxxi = domain.delv_xi(i) * domain.delx_xi(i) ; - Real_t delvxeta = domain.delv_eta(i) * domain.delx_eta(i) ; - Real_t delvxzeta = domain.delv_zeta(i) * domain.delx_zeta(i) ; - - if ( delvxxi > Real_t(0.) ) delvxxi = Real_t(0.) ; - if ( delvxeta > Real_t(0.) ) delvxeta = Real_t(0.) ; - if ( delvxzeta > Real_t(0.) ) delvxzeta = Real_t(0.) ; - - Real_t rho = domain.elemMass(i) / (domain.volo(i) * vnew[i]) ; - - qlin = -qlc_monoq * rho * - ( delvxxi * (Real_t(1.) - phixi) + - delvxeta * (Real_t(1.) - phieta) + - delvxzeta * (Real_t(1.) - phizeta) ) ; - - qquad = qqc_monoq * rho * - ( delvxxi*delvxxi * (Real_t(1.) - phixi*phixi) + - delvxeta*delvxeta * (Real_t(1.) - phieta*phieta) + - delvxzeta*delvxzeta * (Real_t(1.) - phizeta*phizeta) ) ; - } + qquad = qqc_monoq * rho * + (delvxxi * delvxxi * (Real_t(1.) - phixi * phixi) + + delvxeta * delvxeta * (Real_t(1.) - phieta * phieta) + + delvxzeta * delvxzeta * (Real_t(1.) - phizeta * phizeta)); + } - domain.qq(i) = qquad ; - domain.ql(i) = qlin ; - } + domain.qq(i) = qquad; + domain.ql(i) = qlin; + } } /******************************************/ static inline -void CalcMonotonicQForElems(Domain& domain, Real_t vnew[]) -{ - // - // initialize parameters - // - const Real_t ptiny = Real_t(1.e-36) ; - - // - // calculate the monotonic q for all regions - // - for (Index_t r=0 ; r 0) { - CalcMonotonicQRegionForElems(domain, r, vnew, ptiny) ; - } - } +void CalcMonotonicQForElems(Domain &domain, Real_t vnew[]) { + // + // initialize parameters + // + const Real_t ptiny = Real_t(1.e-36); + + // + // calculate the monotonic q for all regions + // + for (Index_t r = 0; r < domain.numReg(); ++r) { + + if (domain.regElemSize(r) > 0) { + CalcMonotonicQRegionForElems(domain, r, vnew, ptiny); + } + } } /******************************************/ static inline -void CalcQForElems(Domain& domain, Real_t vnew[]) -{ - // - // MONOTONIC Q option - // +void CalcQForElems(Domain &domain, Real_t vnew[]) { + // + // MONOTONIC Q option + // + + Index_t numElem = domain.numElem(); - Index_t numElem = domain.numElem() ; + if (numElem != 0) { - if (numElem != 0) { + //Int_t allElem = numElem + /* local elem */ + // 2*domain.sizeX()*domain.sizeY() + /* plane ghosts */ + // 2*domain.sizeX()*domain.sizeZ() + /* row ghosts */ + // 2*domain.sizeY()*domain.sizeZ() ; /* col ghosts */ + //domain.AllocateGradients(numElem, numRanks, allElem); - //Int_t allElem = numElem + /* local elem */ - // 2*domain.sizeX()*domain.sizeY() + /* plane ghosts */ - // 2*domain.sizeX()*domain.sizeZ() + /* row ghosts */ - // 2*domain.sizeY()*domain.sizeZ() ; /* col ghosts */ - //domain.AllocateGradients(numElem, numRanks, allElem); +#if USE_MPI + //CommRecv(domain, MSG_MONOQ, 3, + // domain.sizeX(), domain.sizeY(), domain.sizeZ(), + // true, true) ; + domain.get_delv_xi().switch_to_p1(); + domain.get_delv_eta().switch_to_p1(); + domain.get_delv_zeta().switch_to_p1(); +#endif -#if USE_MPI - //CommRecv(domain, MSG_MONOQ, 3, - // domain.sizeX(), domain.sizeY(), domain.sizeZ(), - // true, true) ; - domain.get_delv_xi().switch_to_p1(); - domain.get_delv_eta().switch_to_p1(); - domain.get_delv_zeta().switch_to_p1(); -#endif + /* Calculate velocity gradients */ + CalcMonotonicQGradientsForElems(domain, vnew); - /* Calculate velocity gradients */ - CalcMonotonicQGradientsForElems(domain, vnew); +#if USE_MPI -#if USE_MPI + domain.get_delv_xi().switch_to_p2(); + domain.get_delv_eta().switch_to_p2(); + domain.get_delv_zeta().switch_to_p2(); - domain.get_delv_xi().switch_to_p2(); - domain.get_delv_eta().switch_to_p2(); - domain.get_delv_zeta().switch_to_p2(); + //Domain_member fieldData[3] ; - //Domain_member fieldData[3] ; - - /* Transfer veloctiy gradients in the first order elements */ - /* problem->commElements->Transfer(CommElements::monoQ) ; */ + /* Transfer veloctiy gradients in the first order elements */ + /* problem->commElements->Transfer(CommElements::monoQ) ; */ - //fieldData[0] = &Domain::delv_xi ; - //fieldData[1] = &Domain::delv_eta ; - //fieldData[2] = &Domain::delv_zeta ; + //fieldData[0] = &Domain::delv_xi ; + //fieldData[1] = &Domain::delv_eta ; + //fieldData[2] = &Domain::delv_zeta ; - //CommSend(domain, MSG_MONOQ, 3, fieldData, - // domain.sizeX(), domain.sizeY(), domain.sizeZ(), - // true, true) ; + //CommSend(domain, MSG_MONOQ, 3, fieldData, + // domain.sizeX(), domain.sizeY(), domain.sizeZ(), + // true, true) ; - //CommMonoQ(domain) ; -#endif + //CommMonoQ(domain) ; +#endif - CalcMonotonicQForElems(domain, vnew) ; + CalcMonotonicQForElems(domain, vnew); - // Free up memory - //domain.DeallocateGradients(); + // Free up memory + //domain.DeallocateGradients(); - /* + /* for (Index_t i=0; i<20; ++i) { laik_log((Laik_LogLevel)2, "i: %d, value: %f", i, domain.delv_xi(i)); } */ - /* Don't allow excessive artificial viscosity */ - Index_t idx = -1; - for (Index_t i=0; i domain.qstop() ) { - idx = i ; - break ; - } - } + /* Don't allow excessive artificial viscosity */ + Index_t idx = -1; + for (Index_t i = 0; i < numElem; ++i) { + if (domain.q(i) > domain.qstop()) { + idx = i; + break; + } + } - if(idx >= 0) { + if (idx >= 0) { #if USE_MPI - laik_log((Laik_LogLevel)2,"Debug 4\n"); - MPI_Barrier(MPI_COMM_WORLD); - MPI_Abort(MPI_COMM_WORLD, QStopError) ; + laik_log(LAIK_LL_Warning, "Debug 4\n"); + MPI_Barrier(MPI_COMM_WORLD); + MPI_Abort(MPI_COMM_WORLD, QStopError); #else - exit(QStopError); + exit(QStopError); #endif - } - } + } + } } /******************************************/ static inline -void CalcPressureForElems(Real_t* p_new, Real_t* bvc, - Real_t* pbvc, Real_t* e_old, - Real_t* compression, Real_t *vnewc, +void CalcPressureForElems(Real_t *p_new, Real_t *bvc, + Real_t *pbvc, Real_t *e_old, + Real_t *compression, Real_t *vnewc, Real_t pmin, Real_t p_cut, Real_t eosvmax, - Index_t length, std::vector regElemList) -{ + Index_t length, std::vector regElemList) { #pragma omp parallel for firstprivate(length) - for (Index_t i = 0; i < length ; ++i) { - Real_t c1s = Real_t(2.0)/Real_t(3.0) ; - bvc[i] = c1s * (compression[i] + Real_t(1.)); - pbvc[i] = c1s; - } + for (Index_t i = 0; i < length; ++i) { + Real_t c1s = Real_t(2.0) / Real_t(3.0); + bvc[i] = c1s * (compression[i] + Real_t(1.)); + pbvc[i] = c1s; + } #pragma omp parallel for firstprivate(length, pmin, p_cut, eosvmax) - for (Index_t i = 0 ; i < length ; ++i){ - Index_t elem = regElemList[i]; - - p_new[i] = bvc[i] * e_old[i] ; + for (Index_t i = 0; i < length; ++i) { + Index_t elem = regElemList[i]; - if (FABS(p_new[i]) < p_cut ) - p_new[i] = Real_t(0.0) ; + p_new[i] = bvc[i] * e_old[i]; - if ( vnewc[elem] >= eosvmax ) /* impossible condition here? */ - p_new[i] = Real_t(0.0) ; + if (FABS(p_new[i]) < p_cut) + p_new[i] = Real_t(0.0); - if (p_new[i] < pmin) - p_new[i] = pmin ; - } + if (vnewc[elem] >= eosvmax) /* impossible condition here? */ + p_new[i] = Real_t(0.0); + + if (p_new[i] < pmin) + p_new[i] = pmin; + } } /******************************************/ static inline -void CalcEnergyForElems(Real_t* p_new, Real_t* e_new, Real_t* q_new, - Real_t* bvc, Real_t* pbvc, - Real_t* p_old, Real_t* e_old, Real_t* q_old, - Real_t* compression, Real_t* compHalfStep, - Real_t* vnewc, Real_t* work, Real_t* delvc, Real_t pmin, - Real_t p_cut, Real_t e_cut, Real_t q_cut, Real_t emin, - Real_t* qq_old, Real_t* ql_old, +void CalcEnergyForElems(Real_t *p_new, Real_t *e_new, Real_t *q_new, + Real_t *bvc, Real_t *pbvc, + Real_t *p_old, Real_t *e_old, Real_t *q_old, + Real_t *compression, Real_t *compHalfStep, + Real_t *vnewc, Real_t *work, Real_t *delvc, Real_t pmin, + Real_t p_cut, Real_t e_cut, Real_t q_cut, Real_t emin, + Real_t *qq_old, Real_t *ql_old, Real_t rho0, Real_t eosvmax, - Index_t length, std::vector regElemList) -{ - Real_t *pHalfStep = Allocate(length) ; + Index_t length, std::vector regElemList) { + Real_t *pHalfStep = Allocate(length); #pragma omp parallel for firstprivate(length, emin) - for (Index_t i = 0 ; i < length ; ++i) { - e_new[i] = e_old[i] - Real_t(0.5) * delvc[i] * (p_old[i] + q_old[i]) - + Real_t(0.5) * work[i]; + for (Index_t i = 0; i < length; ++i) { + e_new[i] = e_old[i] - Real_t(0.5) * delvc[i] * (p_old[i] + q_old[i]) + + Real_t(0.5) * work[i]; - if (e_new[i] < emin ) { - e_new[i] = emin ; - } - } + if (e_new[i] < emin) { + e_new[i] = emin; + } + } - CalcPressureForElems(pHalfStep, bvc, pbvc, e_new, compHalfStep, vnewc, - pmin, p_cut, eosvmax, length, regElemList); + CalcPressureForElems(pHalfStep, bvc, pbvc, e_new, compHalfStep, vnewc, + pmin, p_cut, eosvmax, length, regElemList); #pragma omp parallel for firstprivate(length, rho0) - for (Index_t i = 0 ; i < length ; ++i) { - Real_t vhalf = Real_t(1.) / (Real_t(1.) + compHalfStep[i]) ; - - if ( delvc[i] > Real_t(0.) ) { - q_new[i] /* = qq_old[i] = ql_old[i] */ = Real_t(0.) ; - } - else { - Real_t ssc = ( pbvc[i] * e_new[i] - + vhalf * vhalf * bvc[i] * pHalfStep[i] ) / rho0 ; - - if ( ssc <= Real_t(.1111111e-36) ) { - ssc = Real_t(.3333333e-18) ; - } else { - ssc = SQRT(ssc) ; - } + for (Index_t i = 0; i < length; ++i) { + Real_t vhalf = Real_t(1.) / (Real_t(1.) + compHalfStep[i]); + + if (delvc[i] > Real_t(0.)) { + q_new[i] /* = qq_old[i] = ql_old[i] */ = Real_t(0.); + } else { + Real_t ssc = (pbvc[i] * e_new[i] + + vhalf * vhalf * bvc[i] * pHalfStep[i]) / rho0; + + if (ssc <= Real_t(.1111111e-36)) { + ssc = Real_t(.3333333e-18); + } else { + ssc = SQRT(ssc); + } - q_new[i] = (ssc*ql_old[i] + qq_old[i]) ; - } + q_new[i] = (ssc * ql_old[i] + qq_old[i]); + } - e_new[i] = e_new[i] + Real_t(0.5) * delvc[i] - * ( Real_t(3.0)*(p_old[i] + q_old[i]) - - Real_t(4.0)*(pHalfStep[i] + q_new[i])) ; - } + e_new[i] = e_new[i] + Real_t(0.5) * delvc[i] + * (Real_t(3.0) * (p_old[i] + q_old[i]) + - Real_t(4.0) * (pHalfStep[i] + q_new[i])); + } #pragma omp parallel for firstprivate(length, emin, e_cut) - for (Index_t i = 0 ; i < length ; ++i) { + for (Index_t i = 0; i < length; ++i) { - e_new[i] += Real_t(0.5) * work[i]; + e_new[i] += Real_t(0.5) * work[i]; - if (FABS(e_new[i]) < e_cut) { - e_new[i] = Real_t(0.) ; - } - if ( e_new[i] < emin ) { - e_new[i] = emin ; - } - } + if (FABS(e_new[i]) < e_cut) { + e_new[i] = Real_t(0.); + } + if (e_new[i] < emin) { + e_new[i] = emin; + } + } - CalcPressureForElems(p_new, bvc, pbvc, e_new, compression, vnewc, - pmin, p_cut, eosvmax, length, regElemList); + CalcPressureForElems(p_new, bvc, pbvc, e_new, compression, vnewc, + pmin, p_cut, eosvmax, length, regElemList); #pragma omp parallel for firstprivate(length, rho0, emin, e_cut) - for (Index_t i = 0 ; i < length ; ++i){ - const Real_t sixth = Real_t(1.0) / Real_t(6.0) ; - Index_t elem = regElemList[i]; - Real_t q_tilde ; - - if (delvc[i] > Real_t(0.)) { - q_tilde = Real_t(0.) ; - } - else { - Real_t ssc = ( pbvc[i] * e_new[i] - + vnewc[elem] * vnewc[elem] * bvc[i] * p_new[i] ) / rho0 ; - - if ( ssc <= Real_t(.1111111e-36) ) { - ssc = Real_t(.3333333e-18) ; - } else { - ssc = SQRT(ssc) ; - } + for (Index_t i = 0; i < length; ++i) { + const Real_t sixth = Real_t(1.0) / Real_t(6.0); + Index_t elem = regElemList[i]; + Real_t q_tilde; + + if (delvc[i] > Real_t(0.)) { + q_tilde = Real_t(0.); + } else { + Real_t ssc = (pbvc[i] * e_new[i] + + vnewc[elem] * vnewc[elem] * bvc[i] * p_new[i]) / rho0; + + if (ssc <= Real_t(.1111111e-36)) { + ssc = Real_t(.3333333e-18); + } else { + ssc = SQRT(ssc); + } - q_tilde = (ssc*ql_old[i] + qq_old[i]) ; - } + q_tilde = (ssc * ql_old[i] + qq_old[i]); + } - e_new[i] = e_new[i] - ( Real_t(7.0)*(p_old[i] + q_old[i]) - - Real_t(8.0)*(pHalfStep[i] + q_new[i]) - + (p_new[i] + q_tilde)) * delvc[i]*sixth ; + e_new[i] = e_new[i] - (Real_t(7.0) * (p_old[i] + q_old[i]) + - Real_t(8.0) * (pHalfStep[i] + q_new[i]) + + (p_new[i] + q_tilde)) * delvc[i] * sixth; - if (FABS(e_new[i]) < e_cut) { - e_new[i] = Real_t(0.) ; - } - if ( e_new[i] < emin ) { - e_new[i] = emin ; - } - } + if (FABS(e_new[i]) < e_cut) { + e_new[i] = Real_t(0.); + } + if (e_new[i] < emin) { + e_new[i] = emin; + } + } - CalcPressureForElems(p_new, bvc, pbvc, e_new, compression, vnewc, - pmin, p_cut, eosvmax, length, regElemList); + CalcPressureForElems(p_new, bvc, pbvc, e_new, compression, vnewc, + pmin, p_cut, eosvmax, length, regElemList); #pragma omp parallel for firstprivate(length, rho0, q_cut) - for (Index_t i = 0 ; i < length ; ++i){ - Index_t elem = regElemList[i]; + for (Index_t i = 0; i < length; ++i) { + Index_t elem = regElemList[i]; - if ( delvc[i] <= Real_t(0.) ) { - Real_t ssc = ( pbvc[i] * e_new[i] - + vnewc[elem] * vnewc[elem] * bvc[i] * p_new[i] ) / rho0 ; + if (delvc[i] <= Real_t(0.)) { + Real_t ssc = (pbvc[i] * e_new[i] + + vnewc[elem] * vnewc[elem] * bvc[i] * p_new[i]) / rho0; - if ( ssc <= Real_t(.1111111e-36) ) { - ssc = Real_t(.3333333e-18) ; - } else { - ssc = SQRT(ssc) ; - } + if (ssc <= Real_t(.1111111e-36)) { + ssc = Real_t(.3333333e-18); + } else { + ssc = SQRT(ssc); + } - q_new[i] = (ssc*ql_old[i] + qq_old[i]) ; + q_new[i] = (ssc * ql_old[i] + qq_old[i]); - if (FABS(q_new[i]) < q_cut) q_new[i] = Real_t(0.) ; - } - } + if (FABS(q_new[i]) < q_cut) q_new[i] = Real_t(0.); + } + } - Release(&pHalfStep) ; + Release(&pHalfStep); - return ; + return; } /******************************************/ @@ -2290,263 +2301,257 @@ void CalcSoundSpeedForElems(Domain &domain, Real_t *vnewc, Real_t rho0, Real_t *enewc, Real_t *pnewc, Real_t *pbvc, Real_t *bvc, Real_t ss4o3, - Index_t len, std::vector regElemList) -{ + Index_t len, std::vector regElemList) { #pragma omp parallel for firstprivate(rho0, ss4o3) - for (Index_t i = 0; i < len ; ++i) { - Index_t elem = regElemList[i]; - Real_t ssTmp = (pbvc[i] * enewc[i] + vnewc[elem] * vnewc[elem] * - bvc[i] * pnewc[i]) / rho0; - if (ssTmp <= Real_t(.1111111e-36)) { - ssTmp = Real_t(.3333333e-18); - } - else { - ssTmp = SQRT(ssTmp); - } - domain.ss(elem) = ssTmp ; - } + for (Index_t i = 0; i < len; ++i) { + Index_t elem = regElemList[i]; + Real_t ssTmp = (pbvc[i] * enewc[i] + vnewc[elem] * vnewc[elem] * + bvc[i] * pnewc[i]) / rho0; + if (ssTmp <= Real_t(.1111111e-36)) { + ssTmp = Real_t(.3333333e-18); + } else { + ssTmp = SQRT(ssTmp); + } + domain.ss(elem) = ssTmp; + } } /******************************************/ static inline -void EvalEOSForElems(Domain& domain, Real_t *vnewc, - Int_t numElemReg, std::vector regElemList, Int_t rep) -{ - Real_t e_cut = domain.e_cut() ; - Real_t p_cut = domain.p_cut() ; - Real_t ss4o3 = domain.ss4o3() ; - Real_t q_cut = domain.q_cut() ; - - Real_t eosvmax = domain.eosvmax() ; - Real_t eosvmin = domain.eosvmin() ; - Real_t pmin = domain.pmin() ; - Real_t emin = domain.emin() ; - Real_t rho0 = domain.refdens() ; - - // These temporaries will be of different size for - // each call (due to different sized region element - // lists) - Real_t *e_old = Allocate(numElemReg) ; - Real_t *delvc = Allocate(numElemReg) ; - Real_t *p_old = Allocate(numElemReg) ; - Real_t *q_old = Allocate(numElemReg) ; - Real_t *compression = Allocate(numElemReg) ; - Real_t *compHalfStep = Allocate(numElemReg) ; - Real_t *qq_old = Allocate(numElemReg) ; - Real_t *ql_old = Allocate(numElemReg) ; - Real_t *work = Allocate(numElemReg) ; - Real_t *p_new = Allocate(numElemReg) ; - Real_t *e_new = Allocate(numElemReg) ; - Real_t *q_new = Allocate(numElemReg) ; - Real_t *bvc = Allocate(numElemReg) ; - Real_t *pbvc = Allocate(numElemReg) ; - - //loop to add load imbalance based on region number - for(Int_t j = 0; j < rep; j++) { - /* compress data, minimal set */ +void EvalEOSForElems(Domain &domain, Real_t *vnewc, + Int_t numElemReg, std::vector regElemList, Int_t rep) { + Real_t e_cut = domain.e_cut(); + Real_t p_cut = domain.p_cut(); + Real_t ss4o3 = domain.ss4o3(); + Real_t q_cut = domain.q_cut(); + + Real_t eosvmax = domain.eosvmax(); + Real_t eosvmin = domain.eosvmin(); + Real_t pmin = domain.pmin(); + Real_t emin = domain.emin(); + Real_t rho0 = domain.refdens(); + + // These temporaries will be of different size for + // each call (due to different sized region element + // lists) + Real_t *e_old = Allocate(numElemReg); + Real_t *delvc = Allocate(numElemReg); + Real_t *p_old = Allocate(numElemReg); + Real_t *q_old = Allocate(numElemReg); + Real_t *compression = Allocate(numElemReg); + Real_t *compHalfStep = Allocate(numElemReg); + Real_t *qq_old = Allocate(numElemReg); + Real_t *ql_old = Allocate(numElemReg); + Real_t *work = Allocate(numElemReg); + Real_t *p_new = Allocate(numElemReg); + Real_t *e_new = Allocate(numElemReg); + Real_t *q_new = Allocate(numElemReg); + Real_t *bvc = Allocate(numElemReg); + Real_t *pbvc = Allocate(numElemReg); + + //loop to add load imbalance based on region number + for (Int_t j = 0; j < rep; j++) { + /* compress data, minimal set */ #pragma omp parallel - { + { #pragma omp for nowait firstprivate(numElemReg) - for (Index_t i=0; i eosvmax or v < eosvmin */ - if ( eosvmin != Real_t(0.) ) { + for (Index_t i = 0; i < numElemReg; ++i) { + Index_t elem = regElemList[i]; + Real_t vchalf; + compression[i] = Real_t(1.) / vnewc[elem] - Real_t(1.); + vchalf = vnewc[elem] - delvc[i] * Real_t(.5); + compHalfStep[i] = Real_t(1.) / vchalf - Real_t(1.); + } + + /* Check for v > eosvmax or v < eosvmin */ + if (eosvmin != Real_t(0.)) { #pragma omp for nowait firstprivate(numElemReg, eosvmin) - for(Index_t i=0 ; i= eosvmax) { /* impossible due to calling func? */ - p_old[i] = Real_t(0.) ; - compression[i] = Real_t(0.) ; - compHalfStep[i] = Real_t(0.) ; - } + for (Index_t i = 0; i < numElemReg; ++i) { + Index_t elem = regElemList[i]; + if (vnewc[elem] >= eosvmax) { /* impossible due to calling func? */ + p_old[i] = Real_t(0.); + compression[i] = Real_t(0.); + compHalfStep[i] = Real_t(0.); + } + } } - } #pragma omp for nowait firstprivate(numElemReg) - for (Index_t i = 0 ; i < numElemReg ; ++i) { - work[i] = Real_t(0.) ; - } - } - CalcEnergyForElems(p_new, e_new, q_new, bvc, pbvc, - p_old, e_old, q_old, compression, compHalfStep, - vnewc, work, delvc, pmin, - p_cut, e_cut, q_cut, emin, - qq_old, ql_old, rho0, eosvmax, - numElemReg, regElemList); - } + for (Index_t i = 0; i < numElemReg; ++i) { + work[i] = Real_t(0.); + } + } + CalcEnergyForElems(p_new, e_new, q_new, bvc, pbvc, + p_old, e_old, q_old, compression, compHalfStep, + vnewc, work, delvc, pmin, + p_cut, e_cut, q_cut, emin, + qq_old, ql_old, rho0, eosvmax, + numElemReg, regElemList); + } #pragma omp parallel for firstprivate(numElemReg) - for (Index_t i=0; i eosvmax) - vnew[i] = eosvmax ; - } - } + for (Index_t i = 0; i < numElem; ++i) { + if (vnew[i] > eosvmax) + vnew[i] = eosvmax; + } + } - // This check may not make perfect sense in LULESH, but - // it's representative of something in the full code - - // just leave it in, please + // This check may not make perfect sense in LULESH, but + // it's representative of something in the full code - + // just leave it in, please #pragma omp for nowait firstprivate(numElem) - for (Index_t i=0; i eosvmax) - vc = eosvmax ; - } - if (vc <= 0.) { -#if USE_MPI - laik_log((Laik_LogLevel)2,"Debug 5\n"); - MPI_Barrier(MPI_COMM_WORLD); - MPI_Abort(MPI_COMM_WORLD, VolumeError) ; + for (Index_t i = 0; i < numElem; ++i) { + Real_t vc = domain.v(i); + if (eosvmin != Real_t(0.)) { + if (vc < eosvmin) + vc = eosvmin; + } + if (eosvmax != Real_t(0.)) { + if (vc > eosvmax) + vc = eosvmax; + } + if (vc <= 0.) { +#if USE_MPI + laik_log(LAIK_LL_Warning, "Debug 5\n"); + MPI_Barrier(MPI_COMM_WORLD); + MPI_Abort(MPI_COMM_WORLD, VolumeError); #else - exit(VolumeError); + exit(VolumeError); #endif - } - } - } + } + } + } - for (Int_t r=0 ; r regElemList = domain.regElemlist(r); - Int_t rep; - //Determine load imbalance for this region - //round down the number with lowest cost - if(r < domain.numReg()/2) - rep = 1; - //you don't get an expensive region unless you at least have 5 regions - else if(r < (domain.numReg() - (domain.numReg()+15)/20)) - rep = 1 + domain.cost(); - //very expensive regions - else - rep = 10 * (1+ domain.cost()); - EvalEOSForElems(domain, vnew, numElemReg, regElemList, rep); - } + for (Int_t r = 0; r < domain.numReg(); r++) { + Index_t numElemReg = domain.regElemSize(r); + std::vector regElemList = domain.regElemlist(r); + Int_t rep; + //Determine load imbalance for this region + //round down the number with lowest cost + if (r < domain.numReg() / 2) + rep = 1; + //you don't get an expensive region unless you at least have 5 regions + else if (r < (domain.numReg() - (domain.numReg() + 15) / 20)) + rep = 1 + domain.cost(); + //very expensive regions + else + rep = 10 * (1 + domain.cost()); + EvalEOSForElems(domain, vnew, numElemReg, regElemList, rep); + } - } + } } /******************************************/ static inline void UpdateVolumesForElems(Domain &domain, Real_t *vnew, - Real_t v_cut, Index_t length) -{ - if (length != 0) { + Real_t v_cut, Index_t length) { + if (length != 0) { #pragma omp parallel for firstprivate(length, v_cut) - for(Index_t i=0 ; i(numElem) ; /* new relative vol -- temp */ +void LagrangeElements(Domain &domain, Index_t numElem) { + Real_t *vnew = Allocate(numElem); /* new relative vol -- temp */ - CalcLagrangeElements(domain, vnew) ; + CalcLagrangeElements(domain, vnew); - /* Calculate Q. (Monotonic q option requires communication) */ - CalcQForElems(domain, vnew) ; + /* Calculate Q. (Monotonic q option requires communication) */ + CalcQForElems(domain, vnew); - ApplyMaterialPropertiesForElems(domain, vnew) ; + ApplyMaterialPropertiesForElems(domain, vnew); - UpdateVolumesForElems(domain, vnew, - domain.v_cut(), numElem) ; + UpdateVolumesForElems(domain, vnew, + domain.v_cut(), numElem); - Release(&vnew); + Release(&vnew); } /******************************************/ @@ -2554,10 +2559,9 @@ void LagrangeElements(Domain& domain, Index_t numElem) static inline void CalcCourantConstraintForElems(Domain &domain, Index_t length, std::vector regElemlist, - Real_t qqc, Real_t& dtcourant) -{ -#if _OPENMP - Index_t threads = omp_get_max_threads(); + Real_t qqc, Real_t &dtcourant) { +#if _OPENMP + Index_t threads = omp_get_max_threads(); static Index_t *courant_elem_per_thread; static Real_t *dtcourant_per_thread; static bool first = true; @@ -2567,62 +2571,62 @@ void CalcCourantConstraintForElems(Domain &domain, Index_t length, first = false; } #else - Index_t threads = 1; - Index_t courant_elem_per_thread[1]; - Real_t dtcourant_per_thread[1]; + Index_t threads = 1; + Index_t courant_elem_per_thread[1]; + Real_t dtcourant_per_thread[1]; #endif #pragma omp parallel firstprivate(length, qqc) - { - Real_t qqc2 = Real_t(64.0) * qqc * qqc ; - Real_t dtcourant_tmp = dtcourant; - Index_t courant_elem = -1 ; + { + Real_t qqc2 = Real_t(64.0) * qqc * qqc; + Real_t dtcourant_tmp = dtcourant; + Index_t courant_elem = -1; #if _OPENMP - Index_t thread_num = omp_get_thread_num(); + Index_t thread_num = omp_get_thread_num(); #else - Index_t thread_num = 0; -#endif - -#pragma omp for - for (Index_t i = 0 ; i < length ; ++i) { - Index_t indx = regElemlist[i] ; - Real_t dtf = domain.ss(indx) * domain.ss(indx) ; - - if ( domain.vdov(indx) < Real_t(0.) ) { - dtf = dtf - + qqc2 * domain.arealg(indx) * domain.arealg(indx) - * domain.vdov(indx) * domain.vdov(indx) ; - } - - dtf = SQRT(dtf) ; - dtf = domain.arealg(indx) / dtf ; - - if (domain.vdov(indx) != Real_t(0.)) { - if ( dtf < dtcourant_tmp ) { - dtcourant_tmp = dtf ; - courant_elem = indx ; + Index_t thread_num = 0; +#endif + +#pragma omp for + for (Index_t i = 0; i < length; ++i) { + Index_t indx = regElemlist[i]; + Real_t dtf = domain.ss(indx) * domain.ss(indx); + + if (domain.vdov(indx) < Real_t(0.)) { + dtf = dtf + + qqc2 * domain.arealg(indx) * domain.arealg(indx) + * domain.vdov(indx) * domain.vdov(indx); } - } - } - dtcourant_per_thread[thread_num] = dtcourant_tmp ; - courant_elem_per_thread[thread_num] = courant_elem ; - } + dtf = SQRT(dtf); + dtf = domain.arealg(indx) / dtf; - for (Index_t i = 1; i < threads; ++i) { - if (dtcourant_per_thread[i] < dtcourant_per_thread[0] ) { - dtcourant_per_thread[0] = dtcourant_per_thread[i]; - courant_elem_per_thread[0] = courant_elem_per_thread[i]; - } - } + if (domain.vdov(indx) != Real_t(0.)) { + if (dtf < dtcourant_tmp) { + dtcourant_tmp = dtf; + courant_elem = indx; + } + } + } - if (courant_elem_per_thread[0] != -1) { - dtcourant = dtcourant_per_thread[0] ; - } + dtcourant_per_thread[thread_num] = dtcourant_tmp; + courant_elem_per_thread[thread_num] = courant_elem; + } - return ; + for (Index_t i = 1; i < threads; ++i) { + if (dtcourant_per_thread[i] < dtcourant_per_thread[0]) { + dtcourant_per_thread[0] = dtcourant_per_thread[i]; + courant_elem_per_thread[0] = courant_elem_per_thread[i]; + } + } + + if (courant_elem_per_thread[0] != -1) { + dtcourant = dtcourant_per_thread[0]; + } + + return; } @@ -2630,10 +2634,9 @@ void CalcCourantConstraintForElems(Domain &domain, Index_t length, static inline void CalcHydroConstraintForElems(Domain &domain, Index_t length, - std::vector regElemlist, Real_t dvovmax, Real_t& dthydro) -{ -#if _OPENMP - Index_t threads = omp_get_max_threads(); + std::vector regElemlist, Real_t dvovmax, Real_t &dthydro) { +#if _OPENMP + Index_t threads = omp_get_max_threads(); static Index_t *hydro_elem_per_thread; static Real_t *dthydro_per_thread; static bool first = true; @@ -2643,100 +2646,99 @@ void CalcHydroConstraintForElems(Domain &domain, Index_t length, first = false; } #else - Index_t threads = 1; - Index_t hydro_elem_per_thread[1]; - Real_t dthydro_per_thread[1]; + Index_t threads = 1; + Index_t hydro_elem_per_thread[1]; + Real_t dthydro_per_thread[1]; #endif #pragma omp parallel firstprivate(length, dvovmax) - { - Real_t dthydro_tmp = dthydro ; - Index_t hydro_elem = -1 ; + { + Real_t dthydro_tmp = dthydro; + Index_t hydro_elem = -1; -#if _OPENMP - Index_t thread_num = omp_get_thread_num(); -#else - Index_t thread_num = 0; -#endif +#if _OPENMP + Index_t thread_num = omp_get_thread_num(); +#else + Index_t thread_num = 0; +#endif #pragma omp for - for (Index_t i = 0 ; i < length ; ++i) { - Index_t indx = regElemlist[i] ; + for (Index_t i = 0; i < length; ++i) { + Index_t indx = regElemlist[i]; - if (domain.vdov(indx) != Real_t(0.)) { - Real_t dtdvov = dvovmax / (FABS(domain.vdov(indx))+Real_t(1.e-20)) ; + if (domain.vdov(indx) != Real_t(0.)) { + Real_t dtdvov = dvovmax / (FABS(domain.vdov(indx)) + Real_t(1.e-20)); - if ( dthydro_tmp > dtdvov ) { - dthydro_tmp = dtdvov ; - hydro_elem = indx ; + if (dthydro_tmp > dtdvov) { + dthydro_tmp = dtdvov; + hydro_elem = indx; + } } - } - } + } - dthydro_per_thread[thread_num] = dthydro_tmp ; - hydro_elem_per_thread[thread_num] = hydro_elem ; - } + dthydro_per_thread[thread_num] = dthydro_tmp; + hydro_elem_per_thread[thread_num] = hydro_elem; + } - for (Index_t i = 1; i < threads; ++i) { - if(dthydro_per_thread[i] < dthydro_per_thread[0]) { - dthydro_per_thread[0] = dthydro_per_thread[i]; - hydro_elem_per_thread[0] = hydro_elem_per_thread[i]; - } - } + for (Index_t i = 1; i < threads; ++i) { + if (dthydro_per_thread[i] < dthydro_per_thread[0]) { + dthydro_per_thread[0] = dthydro_per_thread[i]; + hydro_elem_per_thread[0] = hydro_elem_per_thread[i]; + } + } - if (hydro_elem_per_thread[0] != -1) { - dthydro = dthydro_per_thread[0] ; - } + if (hydro_elem_per_thread[0] != -1) { + dthydro = dthydro_per_thread[0]; + } - return ; + return; } /******************************************/ static inline -void CalcTimeConstraintsForElems(Domain& domain) { +void CalcTimeConstraintsForElems(Domain &domain) { - // Initialize conditions to a very large value - domain.dtcourant() = 1.0e+20; - domain.dthydro() = 1.0e+20; + // Initialize conditions to a very large value + domain.dtcourant() = 1.0e+20; + domain.dthydro() = 1.0e+20; - for (Index_t r=0 ; r < domain.numReg() ; ++r) { - /* evaluate time constraint */ - CalcCourantConstraintForElems(domain, domain.regElemSize(r), + for (Index_t r = 0; r < domain.numReg(); ++r) { + /* evaluate time constraint */ + CalcCourantConstraintForElems(domain, domain.regElemSize(r), + domain.regElemlist(r), + domain.qqc(), + domain.dtcourant()); + + /* check hydro constraint */ + CalcHydroConstraintForElems(domain, domain.regElemSize(r), domain.regElemlist(r), - domain.qqc(), - domain.dtcourant()) ; - - /* check hydro constraint */ - CalcHydroConstraintForElems(domain, domain.regElemSize(r), - domain.regElemlist(r), - domain.dvovmax(), - domain.dthydro()) ; - } + domain.dvovmax(), + domain.dthydro()); + } } /******************************************/ static inline -void LagrangeLeapFrog(Domain& domain) -{ +void LagrangeLeapFrog(Domain &domain) { #ifdef SEDOV_SYNC_POS_VEL_LATE - Domain_member fieldData[6] ; + Domain_member fieldData[6] ; #endif - /* calculate nodal forces, accelerations, velocities, positions, with + /* calculate nodal forces, accelerations, velocities, positions, with * applied boundary conditions and slide surface considerations */ - LagrangeNodal(domain); + LagrangeNodal(domain); #ifdef SEDOV_SYNC_POS_VEL_LATE #endif - /* calculate element quantities (i.e. velocity gradient & q), and update + /* calculate element quantities (i.e. velocity gradient & q), and update * material states */ - LagrangeElements(domain, domain.numElem()); + LagrangeElements(domain, domain.numElem()); -#if USE_MPI +#if USE_MPI #ifdef SEDOV_SYNC_POS_VEL_LATE - CommRecv(domain, MSG_SYNC_POS_VEL, 6, + CommRecv(domain, MSG_SYNC_POS_VEL, 6, domain.sizeX() + 1, domain.sizeY() + 1, domain.sizeZ() + 1, false, false) ; @@ -2751,153 +2753,175 @@ void LagrangeLeapFrog(Domain& domain) domain.sizeX() + 1, domain.sizeY() + 1, domain.sizeZ() + 1, false, false) ; #endif -#endif +#endif - CalcTimeConstraintsForElems(domain); + CalcTimeConstraintsForElems(domain); -#if USE_MPI +#if USE_MPI #ifdef SEDOV_SYNC_POS_VEL_LATE - CommSyncPosVel(domain) ; + CommSyncPosVel(domain) ; +#endif #endif -#endif } /******************************************/ -int main(int argc, char *argv[]) -{ - Domain *locDom ; - Int_t numRanks ; - Int_t myRank ; - struct cmdLineOpts opts; +int main(int argc, char *argv[]) { + Domain *locDom; + Int_t numRanks; + Int_t myRank; + struct cmdLineOpts opts; + + +#if USE_MPI + //Domain_member fieldData ; -#if USE_MPI - //Domain_member fieldData ; + Laik_Instance *inst = laik_init_mpi(&argc, &argv); + Laik_Group *world = laik_world(inst); - Laik_Instance* inst = laik_init_mpi(&argc, &argv); - Laik_Group* world = laik_world(inst); + TRACE_INIT(laik_myid(world)); + TRACE_EVENT_START("INIT", ""); - numRanks = laik_size(world); - myRank = laik_myid(world); + numRanks = laik_size(world); + myRank = laik_myid(world); - //MPI_Comm_size(MPI_COMM_WORLD, &numRanks) ; - //MPI_Comm_rank(MPI_COMM_WORLD, &myRank) ; + //MPI_Comm_size(MPI_COMM_WORLD, &numRanks) ; + //MPI_Comm_rank(MPI_COMM_WORLD, &myRank) ; #else - numRanks = 1; + numRanks = 1; myRank = 0; -#endif - - /* Set defaults that can be overridden by command line opts */ - opts.its = 9999999; - opts.nx = 30; - opts.numReg = 11; - opts.numFiles = (int)(numRanks+10)/9; - opts.showProg = 0; - opts.quiet = 0; - opts.viz = 0; - opts.balance = 1; - opts.cost = 1; - opts.repart = 0; - opts.cycle = 5; - - ParseCommandLineOptions(argc, argv, myRank, &opts); - - if ((myRank == 0) && (opts.quiet == 0)) { - printf("Running problem size %d^3 per domain until completion\n", opts.nx); - printf("Num processors: %d\n", numRanks); +#endif + +#ifdef FAULT_TOLERANCE + std::vector checkpoints; +#endif + + /* Set defaults that can be overridden by command line opts */ + opts.its = 9999999; + opts.nx = 30; + opts.numReg = 11; + opts.numFiles = (int) (numRanks + 10) / 9; + opts.showProg = 0; + opts.quiet = 0; + opts.viz = 0; + opts.balance = 1; + opts.cost = 1; + opts.repart = 0; + opts.cycle = 5; + opts.noTimer = false; + + FaultToleranceOptions ftOptions = FaultToleranceOptionsDefault; + + ParseCommandLineOptions(argc, argv, myRank, &opts, &ftOptions); + + if ((myRank == 0) && (opts.quiet == 0)) { + printf("Running problem size %d^3 per domain until completion\n", opts.nx); + printf("Num processors: %d\n", numRanks); #if _OPENMP - printf("Num threads: %d\n", omp_get_max_threads()); + printf("Num threads: %d\n", omp_get_max_threads()); #endif - printf("Total number of elements: %lld\n\n", (long long int)(numRanks*opts.nx*opts.nx*opts.nx)); - printf("To run other sizes, use -s .\n"); - printf("To run a fixed number of iterations, use -i .\n"); - printf("To run a more or less balanced region set, use -b .\n"); - printf("To change the relative costs of regions, use -c .\n"); - printf("To print out progress, use -p\n"); - printf("To write an output file for VisIt, use -v\n"); - printf("See help (-h) for more options\n\n"); - + printf("Total number of elements: %lld\n\n", (long long int) (numRanks * opts.nx * opts.nx * opts.nx)); + printf("To run other sizes, use -s .\n"); + printf("To run a fixed number of iterations, use -i .\n"); + printf("To run a more or less balanced region set, use -b .\n"); + printf("To change the relative costs of regions, use -c .\n"); + printf("To print out progress, use -p\n"); + printf("To write an output file for VisIt, use -v\n"); + printf("See help (-h) for more options\n\n"); + #ifdef PERFORMANCE - printf(" \n=== TESTS ON PERFORMANCE MODE === \n"); + printf(" \n=== TESTS ON PERFORMANCE MODE === \n"); #endif #ifdef REPARTITIONING - printf(" \n=== TESTS ON REPARTITIONING MODE === \n"); + printf(" \n=== TESTS ON REPARTITIONING MODE === \n"); #endif - } + } - // Set up the mesh and decompose. Assumes regular cubes for now - Int_t col, row, plane, side; - InitMeshDecomp(numRanks, myRank, &col, &row, &plane, &side); - - // laik related code - int numElem = opts.nx * opts.nx * opts.nx; // local number of elements - int NumNodes=(opts.nx*side+1)*(opts.nx*side+1)*(opts.nx*side+1); // global number of nodes - - // index spaces used for partitioning of Element and Node data structures - Laik_Space *indexSpaceElements = nullptr; // Elements - Laik_Space *indexSpaceNodes = nullptr; // Nodes - Laik_Space* indexSapceDt = nullptr; // for reduction of dt - - // partitioning of Element and Node data structures - Laik_Partitioning *exclusivePartitioning = nullptr; - Laik_Partitioning *haloPartitioning = nullptr; - Laik_Partitioning *overlapingPartitioning = nullptr; - Laik_Partitioning *allPartitioning = nullptr; - // partitioning objects used after re-partitioning - Laik_Partitioning *exclusivePartitioning2 = nullptr; - Laik_Partitioning *haloPartitioning2 = nullptr; - Laik_Partitioning *overlapingPartitioning2 = nullptr; - Laik_Partitioning *allPartitioning2 = nullptr; - - // transition objects among partitionings. - Laik_Transition *transitionToExclusive = nullptr; - Laik_Transition *transitionToHalo = nullptr; - Laik_Transition *transitionToOverlappingInit = nullptr; - Laik_Transition *transitionToOverlappingReduce = nullptr; - // transition objects used after repartitioning - Laik_Transition *transitionToExclusive2 = nullptr; - Laik_Transition *transitionToHalo2 = nullptr; - Laik_Transition *transitionToOverlappingInit2 = nullptr; - Laik_Transition *transitionToOverlappingReduce2 = nullptr; - - // Laik data container for dt - Laik_Data* laikDt = nullptr; - - // create index spaces for elements, nodes and dt - indexSpaceElements = laik_new_space_1d(inst, numRanks*numElem); - indexSpaceNodes = laik_new_space_1d(inst, NumNodes); - indexSapceDt = laik_new_space_1d(inst, 1); - - // create partitionings and transition objects - create_partitionings_and_transitions(world, - indexSpaceElements, indexSpaceNodes, indexSapceDt, - exclusivePartitioning, haloPartitioning, overlapingPartitioning, allPartitioning, - transitionToExclusive, transitionToHalo, transitionToOverlappingInit, transitionToOverlappingReduce); - - // create laik_data for dt - // Note the rest of containers are - // created inside inside Domain by - // initializing laik_vectors - laikDt = laik_new_data(indexSapceDt, laik_Double); - - // for dt we need to get the base pointer to and - // send it to TimeIncrement function to - // perform the local updates and reduction - laik_switchto_partitioning(laikDt, allPartitioning, LAIK_DF_Init, LAIK_RO_Min); - uint64_t dt_count; double* dt_base; - laik_get_map_1d(laikDt, 0, (void**) &dt_base, &dt_count); - - // Build the main data structure and initialize it - // pass the laik_inst and laik_world and partitionings - // and transitions to Domain since Domain contains - // all the data structures and inititlizes them. - locDom = new Domain(numRanks, col, row, plane, opts.nx, - side, opts.numReg, opts.balance, opts.cost, - inst, world, - indexSpaceElements, indexSpaceNodes, - exclusivePartitioning, haloPartitioning, overlapingPartitioning, - transitionToExclusive, transitionToHalo, transitionToOverlappingInit, - transitionToOverlappingReduce) ; + // Set up the mesh and decompose. Assumes regular cubes for now + Int_t col, row, plane, side; + InitMeshDecomp(numRanks, myRank, &col, &row, &plane, &side); + + // laik related code + int numElem = opts.nx * opts.nx * opts.nx; // local number of elements + int NumNodes = (opts.nx * side + 1) * (opts.nx * side + 1) * (opts.nx * side + 1); // global number of nodes + + // index spaces used for partitioning of Element and Node data structures + Laik_Space *indexSpaceElements = nullptr; // Elements + Laik_Space *indexSpaceNodes = nullptr; // Nodes + Laik_Space *indexSapceDt = nullptr; // for reduction of dt + + // partitioning of Element and Node data structures + Laik_Partitioning *exclusivePartitioning = nullptr; + Laik_Partitioning *haloPartitioning = nullptr; + Laik_Partitioning *overlapingPartitioning = nullptr; + Laik_Partitioning *allPartitioning = nullptr; + // partitioning objects used after re-partitioning + Laik_Partitioning *exclusivePartitioning2 = nullptr; + Laik_Partitioning *haloPartitioning2 = nullptr; + Laik_Partitioning *overlapingPartitioning2 = nullptr; + Laik_Partitioning *allPartitioning2 = nullptr; + + // transition objects among partitionings. + Laik_Transition *transitionToExclusive = nullptr; + Laik_Transition *transitionToHalo = nullptr; + Laik_Transition *transitionToOverlappingInit = nullptr; + Laik_Transition *transitionToOverlappingReduce = nullptr; + // transition objects used after repartitioning + Laik_Transition *transitionToExclusive2 = nullptr; + Laik_Transition *transitionToHalo2 = nullptr; + Laik_Transition *transitionToOverlappingInit2 = nullptr; + Laik_Transition *transitionToOverlappingReduce2 = nullptr; + + // Laik data container for dt + Laik_Data *laikDt = nullptr; + + Laik_Data *laikReductionAuxiliaryContainer = nullptr; + + // create index spaces for elements, nodes and dt + indexSpaceElements = laik_new_space_1d(inst, numRanks * numElem); + indexSpaceNodes = laik_new_space_1d(inst, NumNodes); + indexSapceDt = laik_new_space_1d(inst, 1); + + // create partitionings and transition objects + create_partitionings_and_transitions(world, + indexSpaceElements, indexSpaceNodes, indexSapceDt, + exclusivePartitioning, haloPartitioning, overlapingPartitioning, + allPartitioning, + transitionToExclusive, transitionToHalo, transitionToOverlappingInit, + transitionToOverlappingReduce); + + // create laik_data for dt + // Note the rest of containers are + // created inside inside Domain by + // initializing laik_vectors + laikDt = laik_new_data(indexSapceDt, laik_Double); + + laikReductionAuxiliaryContainer = laik_new_data_1d(inst, laik_Double, 1); + Laik_Partitioning *all, *master; + all = laik_new_partitioning(laik_All, world, laik_data_get_space(laikReductionAuxiliaryContainer), nullptr); + master = laik_new_partitioning(laik_Master, world, laik_data_get_space(laikReductionAuxiliaryContainer), nullptr); + laik_switchto_partitioning(laikReductionAuxiliaryContainer, all, LAIK_DF_None, LAIK_RO_None); + + // for dt we need to get the base pointer to and + // send it to TimeIncrement function to + // perform the local updates and reduction + laik_switchto_partitioning(laikDt, allPartitioning, LAIK_DF_Init, LAIK_RO_Min); + uint64_t dt_count; + double *dt_base; + assert(laik_my_slicecount(laik_data_get_partitioning(laikDt)) == 1); + laik_get_map_1d(laikDt, 0, (void **) &dt_base, &dt_count); + + // Build the main data structure and initialize it + // pass the laik_inst and laik_world and partitionings + // and transitions to Domain since Domain contains + // all the data structures and inititlizes them. + locDom = new Domain(numRanks, col, row, plane, opts.nx, + side, opts.numReg, opts.balance, opts.cost, + inst, world, + indexSpaceElements, indexSpaceNodes, + exclusivePartitioning, haloPartitioning, overlapingPartitioning, + transitionToExclusive, transitionToHalo, transitionToOverlappingInit, + transitionToOverlappingReduce); #if USE_MPI /* @@ -2918,105 +2942,228 @@ int main(int argc, char *argv[]) MPI_Barrier(MPI_COMM_WORLD); */ -#endif - - // BEGIN timestep to solution */ -#if USE_MPI - double start = MPI_Wtime(); - double start2 = MPI_Wtime(); - int repart_iter = 0; +#endif + + // BEGIN timestep to solution */ +#if USE_MPI + double start = getMPI_Wtime(&opts); + double start2 = getMPI_Wtime(&opts); + int repart_iter = 0; #else - timeval start; + timeval start; gettimeofday(&start, NULL) ; #endif - while((locDom->time() < locDom->stoptime()) && (locDom->cycle() < opts.its)) { + TRACE_EVENT_END("INIT", ""); + + //TODO do not terminate early in case of fault + while ( /*(locDom->time() < locDom->stoptime()) && */ (locDom->cycle() < opts.its)) { // TODO remove the following def #define USE_MPI 1 #if USE_MPI - - // check if repartitioning has to be done in this iteration - // if so, do repartitoing before the actual iteration - // after repartitioning continue from the current iteration - if (opts.repart>0 && locDom->cycle() == opts.cycle) - { - double intermediate_timer = MPI_Wtime() - start2; - double itG; - MPI_Reduce(&intermediate_timer, &itG, 1, MPI_DOUBLE, - MPI_MAX, 0, MPI_COMM_WORLD); - - if ((myRank == 0) && (opts.quiet == 0)) { - printf("Starting Repartitioning, current runtime = %f s\n", itG); - printf("Previous Time per Iteration: %f s \n", itG/locDom->cycle()); - } - - laik_log((Laik_LogLevel)2,"Before repart\n" ); - laik_log((Laik_LogLevel)2,"My ID in main world: %d\n", laik_myid(world) ); - - double repart_start = MPI_Wtime(); - - // calculate number of removing tasks and a list of them - double newside; int diffsize; int *removeList=nullptr; - calculate_removing_list(world, opts, side, newside, diffsize, removeList); - - // update number of local number of elements per edge in each task - opts.nx = opts.nx * side / (int) newside; - - // create the shrinked group by removeing tasks in removeList - Laik_Group* shrinked_group = laik_new_shrinked_group(world, diffsize, removeList); - - // create partitionings and transition objects - // on the target (shrinked) process group - // we switch to these partitionings and - // the old objects are not used anymore. - create_partitionings_and_transitions(shrinked_group, - indexSpaceElements, indexSpaceNodes, indexSapceDt, - exclusivePartitioning2, haloPartitioning2, overlapingPartitioning2, allPartitioning2, - transitionToExclusive2, transitionToHalo2, transitionToOverlappingInit2, transitionToOverlappingReduce2); - - // migrate data for all the data structures - locDom->re_distribute_data_structures(shrinked_group, exclusivePartitioning2, haloPartitioning2, overlapingPartitioning2, transitionToExclusive2, transitionToHalo2, transitionToOverlappingInit2, transitionToOverlappingReduce2); - - // processes that are not part of the new (shrinked) - // process group have to exit the main loop - if (laik_myid(shrinked_group)==-1) { - break ; - } - - // repartition and geting the base pointer for only for dt - laik_switchto_partitioning(laikDt, allPartitioning2, LAIK_DF_None, LAIK_RO_Min); - laik_switchto_partitioning(laikDt, allPartitioning2, LAIK_DF_Preserve, LAIK_RO_Min); - laik_get_map_1d(laikDt, 0, (void**) &dt_base, &dt_count); - - // update the working process group in codes - world = shrinked_group; - - // update required values according to the new process group - InitMeshDecomp(laik_size(world), laik_myid(world), &col, &row, &plane, &side); - // update Domain according to the new process group - locDom -> re_init_domain(laik_size(world), col, row, plane, opts.nx, - side, opts.numReg, opts.balance, opts.cost); - - double duration = MPI_Wtime() - repart_start; - - laik_log((Laik_LogLevel)2,"After repart\n"); - if (opts.quiet == 0) - printf("Repartition Done in %f s on Rank %d\n", duration, laik_myid(world)); - - start2 = MPI_Wtime(); - repart_iter = locDom->cycle(); - - // remove the old partitioning and transition objects - remove_partitionings_and_transitions(exclusivePartitioning, haloPartitioning, overlapingPartitioning, allPartitioning, - transitionToExclusive, transitionToHalo, transitionToOverlappingInit, transitionToOverlappingReduce); - - // update the working allPartitioning used later in - // reduction for dt - allPartitioning = allPartitioning2; - - free (removeList); - } + + laik_set_iteration(inst, locDom->cycle()); + + // Execute any pre planned failures + exitIfFailureIteration(locDom->cycle(), &ftOptions, inst); + + // check if repartitioning has to be done in this iteration + // if so, do repartitioning before the actual iteration + // after repartitioning continue from the current iteration + if ((opts.repart > 0 && locDom->cycle() == opts.cycle) || (isFaultToleranceActive(&ftOptions) && locDom->cycle() % ftOptions.failureCheckFrequency == 0)) { + std::vector nodeStatuses; + int failedCount = -1; + // Check if a node has failed, then do restore. Else, do a checkpoint. + if(isFaultToleranceActive(&ftOptions)){ + if(myRank == 0) { + laik_log(LAIK_LL_Info, "Checking for failed nodes."); + } + TRACE_EVENT_S("FAILURE-DETECT", ""); + nodeStatuses.reserve(laik_size(world)); + failedCount = laik_failure_check_nodes(inst, world, &nodeStatuses[0]); + laik_log(LAIK_LL_Info, "Rank %i detected %i node failures.", myRank, failedCount); + } + if (!isFaultToleranceActive(&ftOptions) || failedCount > 0) { + TRACE_EVENT_START("RESTORE", ""); + double intermediate_timer = getMPI_Wtime(&opts) - start2; + double itG = -1; + +// MPI_Reduce(&intermediate_timer, &itG, 1, MPI_DOUBLE, +// MPI_MAX, 0, MPI_COMM_WORLD); + itG = laik_reduce(intermediate_timer, myRank, laikReductionAuxiliaryContainer, all, master); + + if ((myRank == 0) && (opts.quiet == 0)) { + printf("Starting Repartitioning, current runtime = %f s\n", itG); + printf("Previous Time per Iteration: %f s \n", itG / locDom->cycle()); + } + + laik_log((Laik_LogLevel) 2, "Before repart\n"); + laik_log((Laik_LogLevel) 2, "My ID in main world: %d\n", laik_myid(world)); + + double repart_start = getMPI_Wtime(&opts); + + // calculate number of removing tasks and a list of them + double newside; + int diffsize; + int *removeList = nullptr; + Laik_Group *shrinked_group; + + //Todo this was switched logic. Based on previous working version, this was changed back. + if (isFaultToleranceActive(&ftOptions)) { + laik_log(LAIK_LL_Info, "Fault tolerance recovery repartitioning pre-step initiated"); +// calculate_removing_list_ft(world, opts, side, newside, diffsize, removeList, &nodeStatuses[0]); + + // Only remove the actually failed before restoring + laik_failure_eliminate_nodes(inst, failedCount, &nodeStatuses[0]); + shrinked_group = laik_world_fault_tolerant(inst); + + + for(auto& checkpoint : checkpoints) { + laik_checkpoint_remove_failed_slices(checkpoint, world, &nodeStatuses[0]); + } + + laik_log(LAIK_LL_Info, "Restoring checkpoints."); + int restored = locDom->restore(checkpoints, shrinked_group); + laik_log(LAIK_LL_Info, "Restored %i checkpoints.", restored); + + // Make sure the dt container is also switched to the temporary world before doing the second switch. + laik_log(LAIK_LL_Debug, "Switching the dt container to temporary group."); + laik_switchto_new_partitioning(laikDt, shrinked_group, laik_All, LAIK_DF_None, LAIK_RO_None); + laik_log(LAIK_LL_Debug, "Switched the dt container to temporary group."); + + laik_log(LAIK_LL_Debug, "Switching the auxiliary container to temporary group."); + laik_switchto_new_partitioning(laikReductionAuxiliaryContainer, shrinked_group, laik_All, LAIK_DF_None, LAIK_RO_None); + laik_log(LAIK_LL_Debug, "Switched the auxiliary container to temporary group."); + + // Fake new temporary world + world = shrinked_group; + laik_log(LAIK_LL_Info, "Installed fake world."); + } + + calculate_removing_list(world, opts, side, newside, diffsize, removeList); + + // create the shrinked group by removing tasks in removeList + shrinked_group = laik_new_shrinked_group(world, diffsize, removeList); + + if(laik_size(world) - diffsize == 1) { + laik_log(LAIK_LL_Info, "Only 1 process left in group, disabling redundancy."); + ftOptions.redundancyCount = 0; + } + + // update number of local number of elements per edge in each task + opts.nx = opts.nx * side / (int) newside; + + + // create partitionings and transition objects + // on the target (shrinked) process group + // we switch to these partitionings and + // the old objects are not used anymore. + laik_log(LAIK_LL_Info, "Creating partitionings and transitions."); + create_partitionings_and_transitions(shrinked_group, + indexSpaceElements, indexSpaceNodes, indexSapceDt, + exclusivePartitioning2, haloPartitioning2, overlapingPartitioning2, + allPartitioning2, + transitionToExclusive2, transitionToHalo2, + transitionToOverlappingInit2, transitionToOverlappingReduce2); + master = laik_new_partitioning(laik_Master, world, laik_data_get_space(laikReductionAuxiliaryContainer), nullptr); + all = laik_new_partitioning(laik_All, world, laik_data_get_space(laikReductionAuxiliaryContainer), nullptr); + + // migrate data for all the data structures + laik_log(LAIK_LL_Info, "Redistributing data structures."); + //Todo: This was also switched sign, see other TODO + locDom->re_distribute_data_structures(shrinked_group, exclusivePartitioning2, haloPartitioning2, + overlapingPartitioning2, transitionToExclusive2, + transitionToHalo2, transitionToOverlappingInit2, + transitionToOverlappingReduce2, + isFaultToleranceActive(&ftOptions)); + laik_switchto_partitioning(laikReductionAuxiliaryContainer, all, LAIK_DF_None, LAIK_RO_None); + + // processes that are not part of the new (shrinked) + // process group have to exit the main loop + if (laik_myid(shrinked_group) == -1) { + laik_log(LAIK_LL_Info, "Rank %i eliminated. Good bye.", myRank); + break; + } + + // repartition and geting the base pointer for only for dt + laik_switchto_partitioning(laikDt, allPartitioning2, LAIK_DF_None, LAIK_RO_Min); + laik_switchto_partitioning(laikDt, allPartitioning2, LAIK_DF_Preserve, LAIK_RO_Min); + assert(laik_my_slicecount(laik_data_get_partitioning(laikDt)) == 1); + laik_get_map_1d(laikDt, 0, (void **) &dt_base, &dt_count); + + // update the working process group in codes + world = shrinked_group; + std::cout << "New world size: " << laik_size(world) << std::endl; + + // update required values according to the new process group + InitMeshDecomp(laik_size(world), laik_myid(world), &col, &row, &plane, &side); + // update Domain according to the new process group + locDom->re_init_domain(laik_size(world), col, row, plane, opts.nx, + side, opts.numReg, opts.balance, opts.cost); + + double duration = getMPI_Wtime(&opts) - repart_start; + + laik_log((Laik_LogLevel) 2, "After repart\n"); + if (opts.quiet == 0) + printf("Repartition Done in %f s on Rank %d\n", duration, laik_myid(world)); + + start2 = getMPI_Wtime(&opts); + repart_iter = locDom->cycle(); + + // remove the old partitioning and transition objects + remove_partitionings_and_transitions(exclusivePartitioning, haloPartitioning, overlapingPartitioning, + allPartitioning, + transitionToExclusive, transitionToHalo, + transitionToOverlappingInit, transitionToOverlappingReduce); + + // update the working allPartitioning used later in + // reduction for dt + allPartitioning = allPartitioning2; + + free(removeList); + + TRACE_EVENT_END("RESTORE", ""); + } else if(failedCount == 0 && locDom->cycle() % ftOptions.checkpointFrequency == 0) { + if(!ftOptions.delayCheckpointRelease) { + if(myRank == 0) { + laik_log(LAIK_LL_Info, "Freeing %zu checkpoints.\n", checkpoints.size()); + } + + for (auto & checkpoint : checkpoints) { + laik_checkpoint_free(checkpoint); + } + checkpoints.clear(); + } + + if(myRank == 0) { + laik_log(LAIK_LL_Info, "Creating checkpoints."); + } + TRACE_EVENT_START("CHECKPOINT", ""); + std::vector newCheckpoints; + locDom->createCheckpoints(newCheckpoints, ftOptions.redundancyCount, ftOptions.rotationDistance); + TRACE_EVENT_END("CHECKPOINT", ""); + if(myRank == 0) { + laik_log(LAIK_LL_Info, "Finished creating checkpoints."); + } + + if(ftOptions.delayCheckpointRelease) { + if(myRank == 0) { + laik_log(LAIK_LL_Info, "Freeing %zu checkpoints.\n", checkpoints.size()); + } + + for (auto & checkpoint : checkpoints) { + laik_checkpoint_free(checkpoint); + } + checkpoints.clear(); + } + + checkpoints = newCheckpoints; + + // Make sure error handler is installed so we can make use of the checkpoints + laik_error_handler_set(inst, laik_failure_default_error_handler); + } + } #endif /* @@ -3034,61 +3181,109 @@ int main(int argc, char *argv[]) */ - TimeIncrement(*locDom, laikDt, allPartitioning, *dt_base) ; - LagrangeLeapFrog(*locDom) ; + TimeIncrement(*locDom, laikDt, allPartitioning, *dt_base); + LagrangeLeapFrog(*locDom); - if ((opts.showProg != 0) && (opts.quiet == 0) && (myRank == 0)) { - printf("cycle = %d, time = %e, dt=%e\n", - locDom->cycle(), double(locDom->time()), double(locDom->deltatime()) ) ; - } - } + if ((opts.showProg != 0) && (opts.quiet == 0) && (myRank == 0)) { + printf("cycle = %d, time = %e, dt=%e\n", + locDom->cycle(), double(locDom->time()), double(locDom->deltatime())); + } + } - // Use reduced max elapsed time - double elapsed_time; -#if USE_MPI - elapsed_time = MPI_Wtime() - start; - double alternate_time = MPI_Wtime() - start2; + TRACE_EVENT_START("FINALIZE", ""); + + // Use reduced max elapsed time + double elapsed_time; +#if USE_MPI + elapsed_time = getMPI_Wtime(&opts) - start; + double alternate_time = getMPI_Wtime(&opts) - start2; #else - timeval end; + timeval end; gettimeofday(&end, NULL) ; elapsed_time = (double)(end.tv_sec - start.tv_sec) + ((double)(end.tv_usec - start.tv_usec))/1000000 ; #endif - double elapsed_timeG; -#if USE_MPI + double elapsed_timeG; +#if USE_MPI double elapsed_timeG2; - MPI_Reduce(&elapsed_time, &elapsed_timeG, 1, MPI_DOUBLE, - MPI_MAX, 0, MPI_COMM_WORLD); - MPI_Reduce(&alternate_time, &elapsed_timeG2, 1, MPI_DOUBLE, - MPI_MAX, 0, MPI_COMM_WORLD); +// MPI_Reduce(&elapsed_time, &elapsed_timeG, 1, MPI_DOUBLE, +// MPI_MAX, 0, MPI_COMM_WORLD); +// MPI_Reduce(&alternate_time, &elapsed_timeG2, 1, MPI_DOUBLE, +// MPI_MAX, 0, MPI_COMM_WORLD); + elapsed_timeG = laik_reduce(elapsed_time, myRank, laikReductionAuxiliaryContainer, all, master); + elapsed_timeG2 = laik_reduce(alternate_time, myRank, laikReductionAuxiliaryContainer, all, master); #else - elapsed_timeG = elapsed_time; + elapsed_timeG = elapsed_time; #endif - // Write out final viz file */ - if (opts.viz) { - DumpToVisit(*locDom, opts.numFiles, myRank, numRanks) ; - } - - if ((myRank == 0) && (opts.quiet == 0)) { - VerifyAndWriteFinalOutput(elapsed_timeG, *locDom, opts.nx, numRanks); + // Write out final viz file */ + if (opts.viz) { + DumpToVisit(*locDom, opts.numFiles, myRank, numRanks); + } + + if ((myRank == 0) && (opts.quiet == 0)) { + VerifyAndWriteFinalOutput(elapsed_timeG, *locDom, opts.nx, numRanks); #ifdef USE_MPI - if (opts.repart>0){ - printf("Time After Repartitioing: %f, per Iteration %f\n", elapsed_timeG2, elapsed_timeG2/(locDom->cycle() - repart_iter)); - } + if (opts.repart > 0) { + printf("Time After Repartitioning: %f, per Iteration %f\n", elapsed_timeG2, + elapsed_timeG2 / (locDom->cycle() - repart_iter)); + } #endif - } + } #if USE_MPI - // free laik_objects - if (opts.repart>0) - remove_partitionings_and_transitions(exclusivePartitioning2, haloPartitioning2, overlapingPartitioning2, allPartitioning2, - transitionToExclusive2, transitionToHalo2, transitionToOverlappingInit2, transitionToOverlappingReduce2); - else - remove_partitionings_and_transitions(exclusivePartitioning, haloPartitioning, overlapingPartitioning, allPartitioning, - transitionToExclusive, transitionToHalo, transitionToOverlappingInit, transitionToOverlappingReduce); - - laik_finalize(inst); + // free laik_objects + if (exclusivePartitioning2 != nullptr) + remove_partitionings_and_transitions(exclusivePartitioning2, haloPartitioning2, overlapingPartitioning2, + allPartitioning2, + transitionToExclusive2, transitionToHalo2, transitionToOverlappingInit2, + transitionToOverlappingReduce2); + else + remove_partitionings_and_transitions(exclusivePartitioning, haloPartitioning, overlapingPartitioning, + allPartitioning, + transitionToExclusive, transitionToHalo, transitionToOverlappingInit, + transitionToOverlappingReduce); + + laik_finalize(inst); #endif - return 0 ; + TRACE_EVENT_END("FINALIZE", ""); + return 0; +} + +double getMPI_Wtime(cmdLineOpts* opts) { + if(opts->noTimer) { + return 0; + } + return MPI_Wtime(); +} + +double laik_reduce(double input, Int_t myRank, Laik_Data *laikTimer, Laik_Partitioning *all, Laik_Partitioning *master) { + double output; + double* buffer; + uint64_t count; + laik_switchto_partitioning(laikTimer, all, LAIK_DF_None, LAIK_RO_None); + assert(laik_my_slicecount(laik_data_get_partitioning(laikTimer)) == 1); + laik_get_map_1d(laikTimer, 0, reinterpret_cast(&buffer), &count); + if(count != 1) { + std::cerr << "Laik reduce received incorrect partitioning size\n" << std::endl; + MPI_Abort(MPI_COMM_WORLD, -1); + } + buffer[0] = input; + laik_switchto_partitioning(laikTimer, master, LAIK_DF_Preserve, LAIK_RO_Max); + if(myRank == 0) { + assert(laik_my_slicecount(laik_data_get_partitioning(laikTimer)) == 1); + laik_get_map_1d(laikTimer, 0, reinterpret_cast(buffer), &count); + if(count != 1) { + std::cerr << "Laik reduce received incorrect partitioning size\n" << std::endl; + MPI_Abort(MPI_COMM_WORLD, -1); + } + output = buffer[0]; + } else { + output = -1; + } + return output; } + +bool parseFaultToleranceOptionsProxy(int argc, char **argv, int *arg, int rank, FaultToleranceOptions *ftOptions) { + return parseFaultToleranceOptions(argc, argv, arg, rank, ftOptions); +} \ No newline at end of file diff --git a/lulesh.h b/lulesh.h index 8b660b4..f14c764 100644 --- a/lulesh.h +++ b/lulesh.h @@ -26,11 +26,12 @@ #include #include -#include +#include "laik_vector.h" #include "laik_vector_comm_exclusive_halo.h" #include "laik_vector_comm_overlapping_overlapping.h" #include "laik_vector_repart_exclusive.h" #include "laik_vector_repart_overlapping.h" +#include "../util/fault-tolerance-options.h" //************************************************** // Allow flexibility for arithmetic representations @@ -571,9 +572,17 @@ m_delx_zeta.resize(numElem) ; * @param t_to_overlapping_init * @param t_to_overlapping_reduce */ - void re_distribute_data_structures(Laik_Group* new_group, Laik_Partitioning* p_exclusive, Laik_Partitioning* p_halo, Laik_Partitioning* p_overlapping, Laik_Transition *t_to_exclusive, Laik_Transition *t_to_halo, Laik_Transition *t_to_overlapping_init, Laik_Transition *t_to_overlapping_reduce); + void re_distribute_data_structures(Laik_Group *new_group, Laik_Partitioning *p_exclusive, Laik_Partitioning *p_halo, + Laik_Partitioning *p_overlapping, Laik_Transition *t_to_exclusive, + Laik_Transition *t_to_halo, Laik_Transition *t_to_overlapping_init, + Laik_Transition *t_to_overlapping_reduce, bool suppressSwitchToP1); + +#ifdef FAULT_TOLERANCE + void createCheckpoints(std::vector &checkpoints, int redundancyCount, int rotationDistance); + int restore(std::vector &checkpoints, Laik_Group *newGroup); +#endif - private: +private: void BuildMesh(Int_t nx, Int_t edgeNodes, Int_t edgeElems); void SetupThreadSupportStructures(); @@ -779,6 +788,7 @@ m_delx_zeta.resize(numElem) ; Index_t m_rowMin, m_rowMax; Index_t m_colMin, m_colMax; Index_t m_planeMin, m_planeMax ; + } ; typedef Real_t &(Domain::* Domain_member )(Index_t) ; @@ -795,6 +805,7 @@ struct cmdLineOpts { Int_t balance; // -b Int_t repart; // -repart Int_t cycle; // -repart_cycle + bool noTimer; }; @@ -807,8 +818,8 @@ Real_t CalcElemVolume( const Real_t x[8], const Real_t z[8]); // lulesh-util -void ParseCommandLineOptions(int argc, char *argv[], - Int_t myRank, struct cmdLineOpts *opts); +void +ParseCommandLineOptions(int argc, char *argv[], int myRank, struct cmdLineOpts *opts, FaultToleranceOptions *ftOptions); void VerifyAndWriteFinalOutput(Real_t elapsed_time, Domain& locDom, Int_t nx, @@ -898,5 +909,7 @@ void remove_partitionings_and_transitions(Laik_Partitioning *&exclusivePartition * @param removeList */ void calculate_removing_list(Laik_Group* world, cmdLineOpts& opts, double side, double& newside, int& diffsize, int *&removeList); +void calculate_removing_list_ft(Laik_Group* world, cmdLineOpts& opts, double side, double& newside, int& diffsize, int *&removeList, int* nodeStatuses); +bool parseFaultToleranceOptionsProxy(int argc, char **argv, int *arg, int rank, FaultToleranceOptions *ftOptions); #endif diff --git a/lulesh_tuple.h b/lulesh_tuple.h index 9ad7412..734eb0b 100644 --- a/lulesh_tuple.h +++ b/lulesh_tuple.h @@ -25,6 +25,7 @@ #include #include +#include "../util/fault-tolerance-options.h" //************************************************** // Allow flexibility for arithmetic representations @@ -583,8 +584,8 @@ Real_t CalcElemVolume( const Real_t x[8], const Real_t z[8]); // lulesh-util -void ParseCommandLineOptions(int argc, char *argv[], - Int_t myRank, struct cmdLineOpts *opts); +void +ParseCommandLineOptions(int argc, char *argv[], int myRank, struct cmdLineOpts *opts, FaultToleranceOptions *ftOptions); void VerifyAndWriteFinalOutput(Real_t elapsed_time, Domain& locDom, Int_t nx,