domainQ.c
b87f34a0ce7emaster

domainQ.c
View Options

This document is not UTF8. It was detected as Shift JIS and converted to UTF8 for display.

			#include <stdio.h>
			#include <stdlib.h>
			#include <string.h>
			#include <math.h>
			#include <mpi.h>

			#include "allvars.h"
			#include "proto.h"


			#ifdef PY_INTERFACE

			/*! \file domain.c
			* \brief code for domain decomposition
			*
			* This file contains the code for the domain decomposition of the
			* simulation volume. The domains are constructed from disjoint subsets
			* of the leaves of a fiducial top-level tree that covers the full
			* simulation volume. Domain boundaries hence run along tree-node
			* divisions of a fiducial global BH tree. As a result of this method, the
			* tree force are in principle strictly independent of the way the domains
			* are cut. The domain decomposition can be carried out for an arbitrary
			* number of CPUs. Individual domains are not cubical, but spatially
			* coherent since the leaves are traversed in a Peano-Hilbert order and
			* individual domains form segments along this order. This also ensures
			* that each domain has a small surface to volume ratio, which minimizes
			* communication.
			*/

			#define TOPNODEFACTOR 20.0

			#define REDUC_FAC 0.98


			/! toGo[taskNTask + partner] gives the number of particles in task 'task'
			* that have to go to task 'partner'
			*/
			static int toGo, toGoSph;
			static int local_toGo, local_toGoSph;
			static int *list_NumPart;
			static int *list_N_gas;
			static int *list_load;
			static int *list_loadsph;
			static double *list_work;

			static long long maxload, maxloadsph;

			static struct topnode_exchange
			{
			peanokey Startkey;
			int Count;
			}
			toplist, toplist_local;



			/*! This is the main routine for the domain decomposition. It acts as a
			* driver routine that allocates various temporary buffers, maps the
			* particles back onto the periodic box if needed, and then does the
			* domain decomposition, and a final Peano-Hilbert order of all particles
			* as a tuning measure.
			*/
			void domain_DecompositionQ(void)
			{
			double t0, t1;

			#ifdef PMGRID
			//if(All.PM_Ti_endstep == All.Ti_Current)
			{
			All.NumForcesSinceLastDomainDecomp = 1 + All.TotNumPartQ * All.TreeDomainUpdateFrequency;
			/* to make sure that we do a domain decomposition before the PM-force is evaluated.
			this is needed to make sure that the particles are wrapped into the box */
			}
			#endif

			/* Check whether it is really time for a new domain decomposition */
			//if(All.NumForcesSinceLastDomainDecomp > All.TotNumPartQ * All.TreeDomainUpdateFrequency)
			{
			t0 = second();

			#ifdef PERIODIC
			do_box_wrappingQ(); /* map the particles back onto the box */
			#endif
			All.NumForcesSinceLastDomainDecomp = 0;
			TreeReconstructFlag = 1; /* ensures that new tree will be constructed */

			if(ThisTask == 0)
			{
			printf("domain decompositionQ... \n");
			fflush(stdout);
			}

			Key = malloc(sizeof(peanokey) * All.MaxPartQ);
			KeySorted = malloc(sizeof(peanokey) * All.MaxPartQ);

			toGo = malloc(sizeof(int) * NTask * NTask);
			toGoSph = malloc(sizeof(int) * NTask * NTask);
			local_toGo = malloc(sizeof(int) * NTask);
			local_toGoSph = malloc(sizeof(int) * NTask);
			list_NumPart = malloc(sizeof(int) * NTask);
			list_N_gas = malloc(sizeof(int) * NTask);
			list_load = malloc(sizeof(int) * NTask);
			list_loadsph = malloc(sizeof(int) * NTask);
			list_work = malloc(sizeof(double) * NTask);

			MPI_Allgather(&NumPartQ, 1, MPI_INT, list_NumPart, 1, MPI_INT, MPI_COMM_WORLD);
			MPI_Allgather(&N_gasQ, 1, MPI_INT, list_N_gas, 1, MPI_INT, MPI_COMM_WORLD);

			maxload = All.MaxPartQ * REDUC_FAC;
			maxloadsph = All.MaxPartSphQ * REDUC_FAC;

			domain_decomposeQ();

			free(list_work);
			free(list_loadsph);
			free(list_load);
			free(list_N_gas);
			free(list_NumPart);
			free(local_toGoSph);
			free(local_toGo);
			free(toGoSph);
			free(toGo);


			if(ThisTask == 0)
			{
			printf("domain decomposition done. \n");
			fflush(stdout);
			}

			t1 = second();
			All.CPU_Domain += timediff(t0, t1);

			#ifdef PEANOHILBERT
			t0 = second();
			peano_hilbert_orderQ();
			t1 = second();
			All.CPU_Peano += timediff(t0, t1);
			#endif

			free(KeySorted);
			free(Key);
			}

			}



			/*! This function carries out the actual domain decomposition for all
			* particle types. It will try to balance the work-load for each domain,
			* as estimated based on the Q[i]-GravCost values. The decomposition will
			* respect the maximum allowed memory-imbalance given by the value of
			* PartAllocFactor.
			*/
			void domain_decomposeQ(void)
			{
			int i, j, status;
			int ngrp, task, partner, sendcount, recvcount;
			long long sumtogo, sumload;
			int maxload, *temp;
			double sumwork, maxwork;

			for(i = 0; i < 6; i++)
			NtypeLocalQ[i] = 0;

			for(i = 0; i < NumPartQ; i++)
			NtypeLocalQ[Q[i].Type]++;

			/* because Ntype[] is of type `long long', we cannot do a simple
			* MPI_Allreduce() to sum the total particle numbers
			*/
			temp = malloc(NTask * 6 * sizeof(int));
			MPI_Allgather(NtypeLocalQ, 6, MPI_INT, temp, 6, MPI_INT, MPI_COMM_WORLD);
			for(i = 0; i < 6; i++)
			{
			NtypeQ[i] = 0;
			for(j = 0; j < NTask; j++)
			NtypeQ[i] += temp[j * 6 + i];
			}
			free(temp);

			#ifndef UNEQUALSOFTENINGS
			for(i = 0; i < 6; i++)
			if(NtypeQ[i] > 0)
			break;

			for(ngrp = i + 1; ngrp < 6; ngrp++)
			{
			if(NtypeQ[ngrp] > 0)
			if(All.SofteningTableQ[ngrp] != All.SofteningTableQ[i])
			{
			if(ThisTask == 0)
			{
			fprintf(stdout, "Code was not compiled with UNEQUALSOFTENINGS, but some of the\n");
			fprintf(stdout, "softening lengths are unequal nevertheless.\n");
			fprintf(stdout, "This is not allowed.\n");
			}
			endrun(0);
			}
			}
			#endif


			/* determine global dimensions of domain grid */
			domain_findExtentQ();

			domain_determineTopTreeQ();

			/* determine cost distribution in domain grid */
			domain_sumCostQ();

			/* find the split of the domain grid recursively */
			status = domain_findSplitQ(0, NTask, 0, NTopleavesQ - 1);
			if(status != 0)
			{
			if(ThisTask == 0)
			printf("\nNo domain decomposition that stays within memory bounds is possible.\n");
			endrun(0);
			}

			/* now try to improve the work-load balance of the split */
			domain_shiftSplitQ();

			DomainMyStartQ = DomainStartListQ[ThisTask];
			DomainMyLastQ = DomainEndListQ[ThisTask];

			if(ThisTask == 0)
			{
			sumload = maxload = 0;
			sumwork = maxwork = 0;
			for(i = 0; i < NTask; i++)
			{
			sumload += list_load[i];
			sumwork += list_work[i];

			if(list_load[i] > maxload)
			maxload = list_load[i];

			if(list_work[i] > maxwork)
			maxwork = list_work[i];
			}

			printf("work-load balance=%g memory-balance=%g\n",
			maxwork / (sumwork / NTask), maxload / (((double) sumload) / NTask));
			}


			/* determine for each cpu how many particles have to be shifted to other cpus */
			domain_countToGoQ();

			for(i = 0, sumtogo = 0; i < NTask * NTask; i++)
			sumtogo += toGo[i];

			while(sumtogo > 0)
			{
			if(ThisTask == 0)
			{
			printf("exchange of %d%09d particles\n", (int) (sumtogo / 1000000000),
			(int) (sumtogo % 1000000000));
			fflush(stdout);
			}

			for(ngrp = 1; ngrp < (1 << PTask); ngrp++)
			{
			for(task = 0; task < NTask; task++)
			{
			partner = task ^ ngrp;

			if(partner < NTask && task < partner)
			{
			/* treat SPH separately */
			if(All.TotN_gasQ > 0)
			{

			domain_findExchangeNumbersQ(task, partner, 1, &sendcount, &recvcount);

			list_NumPart[task] += recvcount - sendcount;
			list_NumPart[partner] -= recvcount - sendcount;
			list_N_gas[task] += recvcount - sendcount;
			list_N_gas[partner] -= recvcount - sendcount;

			toGo[task * NTask + partner] -= sendcount;
			toGo[partner * NTask + task] -= recvcount;
			toGoSph[task * NTask + partner] -= sendcount;
			toGoSph[partner * NTask + task] -= recvcount;

			if(task == ThisTask) /* actually carry out the exchange */
			domain_exchangeParticlesQ(partner, 1, sendcount, recvcount);
			if(partner == ThisTask)
			domain_exchangeParticlesQ(task, 1, recvcount, sendcount);

			}

			domain_findExchangeNumbersQ(task, partner, 0, &sendcount, &recvcount);

			list_NumPart[task] += recvcount - sendcount;
			list_NumPart[partner] -= recvcount - sendcount;

			toGo[task * NTask + partner] -= sendcount;
			toGo[partner * NTask + task] -= recvcount;

			if(task == ThisTask) /* actually carry out the exchange */
			domain_exchangeParticlesQ(partner, 0, sendcount, recvcount);
			if(partner == ThisTask)
			domain_exchangeParticlesQ(task, 0, recvcount, sendcount);
			}
			}
			}

			for(i = 0, sumtogo = 0; i < NTask * NTask; i++)
			sumtogo += toGo[i];
			}
			}

			/*! This function tries to find a split point in a range of cells in the
			* domain-grid. The range of cells starts at 'first', and ends at 'last'
			* (inclusively). The number of cpus that holds the range is 'ncpu', with
			* the first cpu given by 'cpustart'. If more than 2 cpus are to be split,
			* the function calls itself recursively. The division tries to achieve a
			* best particle-load balance under the constraint that 'maxload' and
			* 'maxloadsph' may not be exceeded, and that each cpu holds at least one
			* cell from the domaingrid. If such a decomposition cannot be achieved, a
			* non-zero error code is returned.
			*
			* After successful completion, DomainMyStartQ[] and DomainMyLastQ[] contain
			* the first and last cell of the domaingrid assigned to the local task
			* for the given type. Also, DomainTaskQ[] contains for each cell the task
			* it was assigned to.
			*/
			int domain_findSplitQ(int cpustart, int ncpu, int first, int last)
			{
			int i, split, ok_left, ok_right;
			long long load, sphload, load_leftOfSplit, sphload_leftOfSplit;
			int ncpu_leftOfSplit;
			double maxAvgLoad_CurrentSplit, maxAvgLoad_NewSplit;


			ncpu_leftOfSplit = ncpu / 2;

			for(i = first, load = 0, sphload = 0; i <= last; i++)
			{
			load += DomainCountQ[i];
			sphload += DomainCountSphQ[i];
			}

			split = first + ncpu_leftOfSplit;

			for(i = first, load_leftOfSplit = sphload_leftOfSplit = 0; i < split; i++)
			{
			load_leftOfSplit += DomainCountQ[i];
			sphload_leftOfSplit += DomainCountSphQ[i];
			}

			/* find the best split point in terms of work-load balance */

			while(split < last - (ncpu - ncpu_leftOfSplit - 1) && split > 0)
			{
			maxAvgLoad_CurrentSplit =
			dmax(load_leftOfSplit / ncpu_leftOfSplit, (load - load_leftOfSplit) / (ncpu - ncpu_leftOfSplit));

			maxAvgLoad_NewSplit =
			dmax((load_leftOfSplit + DomainCountQ[split]) / ncpu_leftOfSplit,
			(load - load_leftOfSplit - DomainCountQ[split]) / (ncpu - ncpu_leftOfSplit));

			if(maxAvgLoad_NewSplit <= maxAvgLoad_CurrentSplit)
			{
			load_leftOfSplit += DomainCountQ[split];
			sphload_leftOfSplit += DomainCountSphQ[split];
			split++;
			}
			else
			break;
			}


			/* we will now have to check whether this solution is possible given the restrictions on the maximum load */

			for(i = first, load_leftOfSplit = 0, sphload_leftOfSplit = 0; i < split; i++)
			{
			load_leftOfSplit += DomainCountQ[i];
			sphload_leftOfSplit += DomainCountSphQ[i];
			}

			if(load_leftOfSplit > maxload * ncpu_leftOfSplit \|\|
			(load - load_leftOfSplit) > maxload * (ncpu - ncpu_leftOfSplit))
			{
			/* we did not find a viable split */
			return -1;
			}

			if(sphload_leftOfSplit > maxloadsph * ncpu_leftOfSplit \|\|
			(sphload - sphload_leftOfSplit) > maxloadsph * (ncpu - ncpu_leftOfSplit))
			{
			/* we did not find a viable split */
			return -1;
			}

			if(ncpu_leftOfSplit >= 2)
			ok_left = domain_findSplitQ(cpustart, ncpu_leftOfSplit, first, split - 1);
			else
			ok_left = 0;

			if((ncpu - ncpu_leftOfSplit) >= 2)
			ok_right = domain_findSplitQ(cpustart + ncpu_leftOfSplit, ncpu - ncpu_leftOfSplit, split, last);
			else
			ok_right = 0;

			if(ok_left == 0 && ok_right == 0)
			{
			/* found a viable split */

			if(ncpu_leftOfSplit == 1)
			{
			for(i = first; i < split; i++)
			DomainTaskQ[i] = cpustart;

			list_load[cpustart] = load_leftOfSplit;
			list_loadsph[cpustart] = sphload_leftOfSplit;
			DomainStartListQ[cpustart] = first;
			DomainEndListQ[cpustart] = split - 1;
			}

			if((ncpu - ncpu_leftOfSplit) == 1)
			{
			for(i = split; i <= last; i++)
			DomainTaskQ[i] = cpustart + ncpu_leftOfSplit;

			list_load[cpustart + ncpu_leftOfSplit] = load - load_leftOfSplit;
			list_loadsph[cpustart + ncpu_leftOfSplit] = sphload - sphload_leftOfSplit;
			DomainStartListQ[cpustart + ncpu_leftOfSplit] = split;
			DomainEndListQ[cpustart + ncpu_leftOfSplit] = last;
			}

			return 0;
			}

			/* we did not find a viable split */
			return -1;
			}



			/*! This function tries to improve the domain decomposition found by
			* domain_findSplitQ() with respect to work-load balance. To this end, the
			* boundaries in the existing domain-split solution (which was found by
			* trying to balance the particle load) are shifted as long as this leads
			* to better work-load while still remaining within the allowed
			* memory-imbalance constraints.
			*/
			void domain_shiftSplitQ(void)
			{
			int i, task, iter = 0, moved;
			double maxw, newmaxw;

			for(task = 0; task < NTask; task++)
			list_work[task] = 0;

			for(i = 0; i < NTopleavesQ; i++)
			list_work[DomainTaskQ[i]] += DomainWorkQ[i];

			do
			{
			for(task = 0, moved = 0; task < NTask - 1; task++)
			{
			maxw = dmax(list_work[task], list_work[task + 1]);

			if(list_work[task] < list_work[task + 1])
			{
			newmaxw = dmax(list_work[task] + DomainWorkQ[DomainStartListQ[task + 1]],
			list_work[task + 1] - DomainWorkQ[DomainStartListQ[task + 1]]);
			if(newmaxw <= maxw)
			{
			if(list_load[task] + DomainCountQ[DomainStartListQ[task + 1]] <= maxload)
			{
			if(list_loadsph[task] + DomainCountSphQ[DomainStartListQ[task + 1]] > maxloadsph)
			continue;

			/* ok, we can move one domain cell from right to left */
			list_work[task] += DomainWorkQ[DomainStartListQ[task + 1]];
			list_load[task] += DomainCountQ[DomainStartListQ[task + 1]];
			list_loadsph[task] += DomainCountSphQ[DomainStartListQ[task + 1]];
			list_work[task + 1] -= DomainWorkQ[DomainStartListQ[task + 1]];
			list_load[task + 1] -= DomainCountQ[DomainStartListQ[task + 1]];
			list_loadsph[task + 1] -= DomainCountSphQ[DomainStartListQ[task + 1]];

			DomainTaskQ[DomainStartListQ[task + 1]] = task;
			DomainStartListQ[task + 1] += 1;
			DomainEndListQ[task] += 1;

			moved++;
			}
			}
			}
			else
			{
			newmaxw = dmax(list_work[task] - DomainWorkQ[DomainEndListQ[task]],
			list_work[task + 1] + DomainWorkQ[DomainEndListQ[task]]);
			if(newmaxw <= maxw)
			{
			if(list_load[task + 1] + DomainCountQ[DomainEndListQ[task]] <= maxload)
			{
			if(list_loadsph[task + 1] + DomainCountSphQ[DomainEndListQ[task]] > maxloadsph)
			continue;

			/* ok, we can move one domain cell from left to right */
			list_work[task] -= DomainWorkQ[DomainEndListQ[task]];
			list_load[task] -= DomainCountQ[DomainEndListQ[task]];
			list_loadsph[task] -= DomainCountSphQ[DomainEndListQ[task]];
			list_work[task + 1] += DomainWorkQ[DomainEndListQ[task]];
			list_load[task + 1] += DomainCountQ[DomainEndListQ[task]];
			list_loadsph[task + 1] += DomainCountSphQ[DomainEndListQ[task]];

			DomainTaskQ[DomainEndListQ[task]] = task + 1;
			DomainEndListQ[task] -= 1;
			DomainStartListQ[task + 1] -= 1;

			moved++;
			}
			}

			}
			}

			iter++;
			}
			while(moved > 0 && iter < 10 * NTopleavesQ);
			}


			/*! This function counts how many particles have to be exchanged between
			* two CPUs according to the domain split. If the CPUs are already quite
			* full and hold data from other CPUs as well, not all the particles may
			* be exchanged at once. In this case the communication phase has to be
			* repeated, until enough of the third-party particles have been moved
			* away such that the decomposition can be completed.
			*/
			void domain_findExchangeNumbersQ(int task, int partner, int sphflag, int send, int recv)
			{
			int numpartA, numpartsphA, ntobesentA, maxsendA, maxsendA_old;
			int numpartB, numpartsphB, ntobesentB, maxsendB, maxsendB_old;

			numpartA = list_NumPart[task];
			numpartsphA = list_N_gas[task];

			numpartB = list_NumPart[partner];
			numpartsphB = list_N_gas[partner];

			if(sphflag == 1)
			{
			ntobesentA = toGoSph[task * NTask + partner];
			ntobesentB = toGoSph[partner * NTask + task];
			}
			else
			{
			ntobesentA = toGo[task * NTask + partner] - toGoSph[task * NTask + partner];
			ntobesentB = toGo[partner * NTask + task] - toGoSph[partner * NTask + task];
			}

			maxsendA = imin(ntobesentA, All.BunchSizeDomain);
			maxsendB = imin(ntobesentB, All.BunchSizeDomain);

			do
			{
			maxsendA_old = maxsendA;
			maxsendB_old = maxsendB;

			maxsendA = imin(All.MaxPartQ - numpartB + maxsendB, maxsendA);
			maxsendB = imin(All.MaxPartQ - numpartA + maxsendA, maxsendB);
			}
			while((maxsendA != maxsendA_old) \|\| (maxsendB != maxsendB_old));


			/* now make also sure that there is enough space for SPH particeles */
			if(sphflag == 1)
			{
			do
			{
			maxsendA_old = maxsendA;
			maxsendB_old = maxsendB;

			maxsendA = imin(All.MaxPartSphQ - numpartsphB + maxsendB, maxsendA);
			maxsendB = imin(All.MaxPartSphQ - numpartsphA + maxsendA, maxsendB);
			}
			while((maxsendA != maxsendA_old) \|\| (maxsendB != maxsendB_old));
			}

			*send = maxsendA;
			*recv = maxsendB;
			}




			/*! This function exchanges particles between two CPUs according to the
			* domain split. In doing this, the memory boundaries which may restrict
			* the exhange process are observed.
			*/
			void domain_exchangeParticlesQ(int partner, int sphflag, int send_count, int recv_count)
			{
			int i, no, n, count, rep;
			MPI_Status status;

			for(n = 0, count = 0; count < send_count && n < NumPartQ; n++)
			{
			if(sphflag)
			{
			if(Q[n].Type != 0)
			continue;
			}
			else
			{
			if(Q[n].Type == 0)
			continue;
			}

			no = 0;

			while(TopNodesQ[no].Daughter >= 0)
			no = TopNodesQ[no].Daughter + (Key[n] - TopNodesQ[no].StartKey) / (TopNodesQ[no].Size / 8);

			no = TopNodesQ[no].Leaf;

			if(DomainTaskQ[no] == partner)
			{
			if(sphflag) /* special reorder routine for SPH particles (need to stay at beginning) */
			{
			DomainPartBufQ[count] = Q[n]; /* copy particle and collect in contiguous memory */
			DomainKeyBufQ[count] = Key[n];
			DomainSphBufQ[count] = SphQ[n];

			Q[n] = Q[N_gasQ - 1];
			Q[N_gasQ - 1] = Q[NumPartQ - 1];

			Key[n] = Key[N_gasQ - 1];
			Key[N_gasQ - 1] = Key[NumPartQ - 1];

			SphQ[n] = SphQ[N_gasQ - 1];

			N_gasQ--;
			}
			else
			{
			DomainPartBufQ[count] = Q[n]; /* copy particle and collect in contiguous memory */
			DomainKeyBufQ[count] = Key[n];
			Q[n] = Q[NumPartQ - 1];
			Key[n] = Key[NumPartQ - 1];
			}

			count++;
			NumPartQ--;
			n--;
			}
			}

			if(count != send_count)
			{
			printf("Houston, we got a problem...\n");
			printf("ThisTask=%d count=%d send_count=%d\n", ThisTask, count, send_count);
			endrun(88);
			}

			/* transmit */

			for(rep = 0; rep < 2; rep++)
			{
			if((rep == 0 && ThisTask < partner) \|\| (rep == 1 && ThisTask > partner))
			{
			if(send_count > 0)
			{
			MPI_Ssend(&DomainPartBufQ[0], send_count * sizeof(struct particle_data), MPI_BYTE, partner,
			TAG_PDATA, MPI_COMM_WORLD);

			MPI_Ssend(&DomainKeyBufQ[0], send_count * sizeof(peanokey), MPI_BYTE, partner, TAG_KEY,
			MPI_COMM_WORLD);

			if(sphflag)
			MPI_Ssend(&DomainSphBufQ[0], send_count * sizeof(struct sph_particle_data), MPI_BYTE, partner,
			TAG_SPHDATA, MPI_COMM_WORLD);
			}
			}

			if((rep == 1 && ThisTask < partner) \|\| (rep == 0 && ThisTask > partner))
			{
			if(recv_count > 0)
			{
			if(sphflag)
			{
			if((NumPartQ - N_gasQ) > recv_count)
			{
			for(i = 0; i < recv_count; i++)
			{
			Q[NumPartQ + i] = Q[N_gasQ + i];
			Key[NumPartQ + i] = Key[N_gasQ + i];
			}
			}
			else
			{
			for(i = NumPartQ - 1; i >= N_gasQ; i--)
			{
			Q[i + recv_count] = Q[i];
			Key[i + recv_count] = Key[i];
			}
			}

			MPI_Recv(&Q[N_gasQ], recv_count * sizeof(struct particle_data), MPI_BYTE, partner, TAG_PDATA,
			MPI_COMM_WORLD, &status);
			MPI_Recv(&Key[N_gasQ], recv_count * sizeof(peanokey), MPI_BYTE, partner, TAG_KEY,
			MPI_COMM_WORLD, &status);
			MPI_Recv(&SphQ[N_gasQ], recv_count * sizeof(struct sph_particle_data), MPI_BYTE, partner,
			TAG_SPHDATA, MPI_COMM_WORLD, &status);

			N_gasQ += recv_count;
			}
			else
			{
			MPI_Recv(&Q[NumPartQ], recv_count * sizeof(struct particle_data), MPI_BYTE, partner,
			TAG_PDATA, MPI_COMM_WORLD, &status);
			MPI_Recv(&Key[NumPartQ], recv_count * sizeof(peanokey), MPI_BYTE, partner,
			TAG_KEY, MPI_COMM_WORLD, &status);
			}

			NumPartQ += recv_count;
			}
			}
			}
			}

			/*! This function determines how many particles that are currently stored
			* on the local CPU have to be moved off according to the domain
			* decomposition.
			*/
			void domain_countToGoQ(void)
			{
			int n, no;

			for(n = 0; n < NTask; n++)
			{
			local_toGo[n] = 0;
			local_toGoSph[n] = 0;
			}

			for(n = 0; n < NumPartQ; n++)
			{
			no = 0;

			while(TopNodesQ[no].Daughter >= 0)
			no = TopNodesQ[no].Daughter + (Key[n] - TopNodesQ[no].StartKey) / (TopNodesQ[no].Size / 8);

			no = TopNodesQ[no].Leaf;

			if(DomainTaskQ[no] != ThisTask)
			{
			local_toGo[DomainTaskQ[no]] += 1;
			if(Q[n].Type == 0)
			local_toGoSph[DomainTaskQ[no]] += 1;
			}
			}

			MPI_Allgather(local_toGo, NTask, MPI_INT, toGo, NTask, MPI_INT, MPI_COMM_WORLD);
			MPI_Allgather(local_toGoSph, NTask, MPI_INT, toGoSph, NTask, MPI_INT, MPI_COMM_WORLD);
			}


			/*! This function walks the global top tree in order to establish the
			* number of leaves it has. These leaves are distributed to different
			* processors.
			*/
			void domain_walktoptreeQ(int no)
			{
			int i;

			if(TopNodesQ[no].Daughter == -1)
			{
			TopNodesQ[no].Leaf = NTopleavesQ;
			NTopleavesQ++;
			}
			else
			{
			for(i = 0; i < 8; i++)
			domain_walktoptreeQ(TopNodesQ[no].Daughter + i);
			}
			}

			/*! This routine bins the particles onto the domain-grid, i.e. it sums up the
			* total number of particles and the total amount of work in each of the
			* domain-cells. This information forms the basis for the actual decision on
			* the adopted domain decomposition.
			*/
			void domain_sumCostQ(void)
			{
			int i, n, no;
			double *local_DomainWorkQ;
			int *local_DomainCountQ;
			int *local_DomainCountSphQ;

			local_DomainWorkQ = malloc(NTopnodesQ * sizeof(double));
			local_DomainCountQ = malloc(NTopnodesQ * sizeof(int));
			local_DomainCountSphQ = malloc(NTopnodesQ * sizeof(int));



			NTopleavesQ = 0;

			domain_walktoptreeQ(0);

			for(i = 0; i < NTopleavesQ; i++)
			{
			local_DomainWorkQ[i] = 0;
			local_DomainCountQ[i] = 0;
			local_DomainCountSphQ[i] = 0;
			}

			if(ThisTask == 0)
			printf("NTopleavesQ= %d\n", NTopleavesQ);

			for(n = 0; n < NumPartQ; n++)
			{
			no = 0;

			while(TopNodesQ[no].Daughter >= 0)
			no = TopNodesQ[no].Daughter + (Key[n] - TopNodesQ[no].StartKey) / (TopNodesQ[no].Size / 8);

			no = TopNodesQ[no].Leaf;

			if(Q[n].Ti_endstep > Q[n].Ti_begstep)
			local_DomainWorkQ[no] += (1.0 + Q[n].GravCost) / (Q[n].Ti_endstep - Q[n].Ti_begstep);
			else
			local_DomainWorkQ[no] += (1.0 + Q[n].GravCost);

			local_DomainCountQ[no] += 1;
			if(Q[n].Type == 0)
			local_DomainCountSphQ[no] += 1;
			}

			MPI_Allreduce(local_DomainWorkQ, DomainWorkQ, NTopleavesQ, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
			MPI_Allreduce(local_DomainCountQ, DomainCountQ, NTopleavesQ, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
			MPI_Allreduce(local_DomainCountSphQ, DomainCountSphQ, NTopleavesQ, MPI_INT, MPI_SUM, MPI_COMM_WORLD);


			free(local_DomainCountSphQ);
			free(local_DomainCountQ);
			free(local_DomainWorkQ);
			}


			/*! This routine finds the extent of the global domain grid.
			*/
			void domain_findExtentQ(void)
			{
			int i, j;
			double len, xmin[3], xmax[3], xmin_glob[3], xmax_glob[3];

			/* determine local extension */
			for(j = 0; j < 3; j++)
			{
			xmin[j] = MAX_REAL_NUMBER;
			xmax[j] = -MAX_REAL_NUMBER;
			}

			for(i = 0; i < NumPartQ; i++)
			{
			for(j = 0; j < 3; j++)
			{
			if(xmin[j] > Q[i].Pos[j])
			xmin[j] = Q[i].Pos[j];

			if(xmax[j] < Q[i].Pos[j])
			xmax[j] = Q[i].Pos[j];
			}
			}

			MPI_Allreduce(xmin, xmin_glob, 3, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
			MPI_Allreduce(xmax, xmax_glob, 3, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);

			len = 0;
			for(j = 0; j < 3; j++)
			if(xmax_glob[j] - xmin_glob[j] > len)
			len = xmax_glob[j] - xmin_glob[j];

			len *= 1.001;

			for(j = 0; j < 3; j++)
			{
			DomainCenterQ[j] = 0.5 * (xmin_glob[j] + xmax_glob[j]);
			DomainCornerQ[j] = 0.5 * (xmin_glob[j] + xmax_glob[j]) - 0.5 * len;
			}

			DomainLenQ = len;
			DomainFacQ = 1.0 / len * (((peanokey) 1) << (BITS_PER_DIMENSION));


			/*
			ici il faudrait plutot tester que tous les Q set P sont dans le meme domaine ...

			il faut faire 軋 ailleurs....

			*/





			if (ThisTask==0)
			if (DomainLenQ!=DomainLen)
			printf("!!!!!!!!!!!!!!!! warning domains are differents !!!!!\n");


			/*
			for(j = 0; j < 3; j++)
			{
			DomainCenterQ[j] = DomainCenter[j];
			DomainCornerQ[j] = DomainCorner[j];
			}

			DomainLenQ = DomainLen;
			DomainFacQ = DomainFac;

			*/


			}


			/*! This function constructs the global top-level tree node that is used
			* for the domain decomposition. This is done by considering the string of
			* Peano-Hilbert keys for all particles, which is recursively chopped off
			* in pieces of eight segments until each segment holds at most a certain
			* number of particles.
			*/
			void domain_determineTopTreeQ(void)
			{
			int i, ntop_local, ntop;
			int ntopnodelist, ntopoffset;

			for(i = 0; i < NumPartQ; i++)
			{
			KeySorted[i] = Key[i] = peano_hilbert_key((Q[i].Pos[0] - DomainCornerQ[0]) * DomainFacQ,
			(Q[i].Pos[1] - DomainCornerQ[1]) * DomainFacQ,
			(Q[i].Pos[2] - DomainCornerQ[2]) * DomainFacQ,
			BITS_PER_DIMENSION);
			}

			qsort(KeySorted, NumPartQ, sizeof(peanokey), domain_compare_key);

			NTopnodesQ = 1;
			TopNodesQ[0].Daughter = -1;
			TopNodesQ[0].Size = PEANOCELLS;
			TopNodesQ[0].StartKey = 0;
			TopNodesQ[0].Count = NumPartQ;
			TopNodesQ[0].Pstart = 0;

			domain_topsplit_localQ(0, 0);

			toplist_local = malloc(NTopnodesQ * sizeof(struct topnode_exchange));

			for(i = 0, ntop_local = 0; i < NTopnodesQ; i++)
			{
			if(TopNodesQ[i].Daughter == -1) /* only use leaves */
			{
			toplist_local[ntop_local].Startkey = TopNodesQ[i].StartKey;
			toplist_local[ntop_local].Count = TopNodesQ[i].Count;
			ntop_local++;
			}
			}

			ntopnodelist = malloc(sizeof(int) * NTask);
			ntopoffset = malloc(sizeof(int) * NTask);

			MPI_Allgather(&ntop_local, 1, MPI_INT, ntopnodelist, 1, MPI_INT, MPI_COMM_WORLD);

			for(i = 0, ntop = 0, ntopoffset[0] = 0; i < NTask; i++)
			{
			ntop += ntopnodelist[i];
			if(i > 0)
			ntopoffset[i] = ntopoffset[i - 1] + ntopnodelist[i - 1];
			}


			toplist = malloc(ntop * sizeof(struct topnode_exchange));

			for(i = 0; i < NTask; i++)
			{
			ntopnodelist[i] *= sizeof(struct topnode_exchange);
			ntopoffset[i] *= sizeof(struct topnode_exchange);
			}

			MPI_Allgatherv(toplist_local, ntop_local * sizeof(struct topnode_exchange), MPI_BYTE,
			toplist, ntopnodelist, ntopoffset, MPI_BYTE, MPI_COMM_WORLD);

			qsort(toplist, ntop, sizeof(struct topnode_exchange), domain_compare_toplist);

			NTopnodesQ = 1;
			TopNodesQ[0].Daughter = -1;
			TopNodesQ[0].Size = PEANOCELLS;
			TopNodesQ[0].StartKey = 0;
			TopNodesQ[0].Count = All.TotNumPartQ;
			TopNodesQ[0].Pstart = 0;
			TopNodesQ[0].Blocks = ntop;

			domain_topsplitQ(0, 0);

			free(toplist);
			free(ntopoffset);
			free(ntopnodelist);
			free(toplist_local);

			}



			/*! This function is responsible for constructing the local top-level
			* Peano-Hilbert segments. A segment is cut into 8 pieces recursively
			* until the number of particles in the segment has fallen below
			* All.TotNumPartQ / (TOPNODEFACTOR * NTask * NTask).
			*/
			void domain_topsplit_localQ(int node, peanokey startkey)
			{
			int i, p, sub, bin;

			if(TopNodesQ[node].Size >= 8)
			{
			TopNodesQ[node].Daughter = NTopnodesQ;

			for(i = 0; i < 8; i++)
			{
			if(NTopnodesQ < MAXTOPNODES)
			{
			sub = TopNodesQ[node].Daughter + i;
			TopNodesQ[sub].Size = TopNodesQ[node].Size / 8;
			TopNodesQ[sub].Count = 0;
			TopNodesQ[sub].Daughter = -1;
			TopNodesQ[sub].StartKey = startkey + i * TopNodesQ[sub].Size;
			TopNodesQ[sub].Pstart = TopNodesQ[node].Pstart;

			NTopnodesQ++;
			}
			else
			{
			printf("task=%d: We are out of Topnodes. Increasing the constant MAXTOPNODES might help.\n",
			ThisTask);
			fflush(stdout);
			endrun(13213);
			}
			}

			for(p = TopNodesQ[node].Pstart; p < TopNodesQ[node].Pstart + TopNodesQ[node].Count; p++)
			{
			bin = (KeySorted[p] - startkey) / (TopNodesQ[node].Size / 8);

			if(bin < 0 \|\| bin > 7)
			{
			printf("task=%d: something odd has happened here. bin=%d\n", ThisTask, bin);
			fflush(stdout);
			endrun(13123123);
			}

			sub = TopNodesQ[node].Daughter + bin;

			if(TopNodesQ[sub].Count == 0)
			TopNodesQ[sub].Pstart = p;

			TopNodesQ[sub].Count++;
			}

			for(i = 0; i < 8; i++)
			{
			sub = TopNodesQ[node].Daughter + i;
			if(TopNodesQ[sub].Count > All.TotNumPartQ / (TOPNODEFACTOR * NTask * NTask))
			domain_topsplit_localQ(sub, TopNodesQ[sub].StartKey);
			}
			}
			}



			/*! This function is responsible for constructing the global top-level tree
			* segments. Starting from a joint list of all local top-level segments,
			* in which mulitple occurences of the same spatial segment have been
			* combined, a segment is subdivided into 8 pieces recursively until the
			* number of particles in each segment has fallen below All.TotNumPartQ /
			* (TOPNODEFACTOR * NTask).
			*/
			void domain_topsplitQ(int node, peanokey startkey)
			{
			int i, p, sub, bin;

			if(TopNodesQ[node].Size >= 8)
			{
			TopNodesQ[node].Daughter = NTopnodesQ;

			for(i = 0; i < 8; i++)
			{
			if(NTopnodesQ < MAXTOPNODES)
			{
			sub = TopNodesQ[node].Daughter + i;
			TopNodesQ[sub].Size = TopNodesQ[node].Size / 8;
			TopNodesQ[sub].Count = 0;
			TopNodesQ[sub].Blocks = 0;
			TopNodesQ[sub].Daughter = -1;
			TopNodesQ[sub].StartKey = startkey + i * TopNodesQ[sub].Size;
			TopNodesQ[sub].Pstart = TopNodesQ[node].Pstart;
			NTopnodesQ++;
			}
			else
			{
			printf("Task=%d: We are out of Topnodes. Increasing the constant MAXTOPNODES might help.\n",
			ThisTask);
			fflush(stdout);
			endrun(137213);
			}
			}

			for(p = TopNodesQ[node].Pstart; p < TopNodesQ[node].Pstart + TopNodesQ[node].Blocks; p++)
			{
			bin = (toplist[p].Startkey - startkey) / (TopNodesQ[node].Size / 8);
			sub = TopNodesQ[node].Daughter + bin;

			if(bin < 0 \|\| bin > 7)
			endrun(77);

			if(TopNodesQ[sub].Blocks == 0)
			TopNodesQ[sub].Pstart = p;

			TopNodesQ[sub].Count += toplist[p].Count;
			TopNodesQ[sub].Blocks++;
			}

			for(i = 0; i < 8; i++)
			{
			sub = TopNodesQ[node].Daughter + i;
			if(TopNodesQ[sub].Count > All.TotNumPartQ / (TOPNODEFACTOR * NTask))
			domain_topsplitQ(sub, TopNodesQ[sub].StartKey);
			}
			}
			}




			#endif

Gear/src/domainQ.cb87f34a0ce7emaster

domainQ.cView Options

Gear/src/domainQ.c
b87f34a0ce7emaster

domainQ.c
View Options