init_entropy.cc 6.7 KB
Newer Older
Volker Springel's avatar
Volker Springel committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
/*******************************************************************************
 * \copyright   This file is part of the GADGET4 N-body/SPH code developed
 * \copyright   by Volker Springel. Copyright (C) 2014-2020 by Volker Springel
 * \copyright   (vspringel@mpa-garching.mpg.de) and all contributing authors.
 *******************************************************************************/

/*! \file  init_entropy.cc
 *
 *  \brief initialization code for the entropy variable of the particles
 */

#include "gadgetconfig.h"

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

#include "../data/allvars.h"
#include "../data/dtypes.h"
#include "../data/mymalloc.h"
#include "../domain/domain.h"
#include "../io/io.h"
#include "../logs/timer.h"
#include "../main/simulation.h"
#include "../mpi_utils/mpi_utils.h"
#include "../ngbtree/ngbtree.h"
#include "../sort/peano.h"
#include "../sph/kernel.h"
#include "../system/system.h"
#include "../time_integration/timestep.h"

#ifdef PRESSURE_ENTROPY_SPH

#define MAX_ITER_ENTROPY 100
#define ENTROPY_TOLERANCE 1.0e-5

/*! \file init_entropy.c
 *  \brief SPH entropy computation from internal energies for pressure-entropy formulation of SPH
 *
 */
void sph::init_entropy(void)
{
  TIMER_STORE;
  TIMER_START(CPU_DENSITY);

  D->mpi_printf("SPH-INIT-ENTROPY: Begin entropy calculation.  (presently allocated=%g MB)\n", Mem.getAllocatedBytesInMB());

  /* Create list of targets. We do this here to simplify the treatment later on */
  int *targetList = (int *)Mem.mymalloc("targetlist", Tp->NumGas * sizeof(int));

  int ndensities = 0;

  for(int i = 0; i < Tp->TimeBinsHydro.NActiveParticles; i++)
    {
      int target = Tp->TimeBinsHydro.ActiveParticleList[i];
      if(target < 0 || target >= Tp->NumGas)
        Terminate("target=%d i=%d\n", target, i);
      targetList[ndensities++] = target;
    }

  int iter = 0;

  // let's grab at most half the still available memory for imported points and nodes
  int nspace = (0.33 * Mem.FreeBytes) / (sizeof(ngbnode) + 8 * sizeof(foreign_sphpoint_data));

  MaxForeignNodes  = nspace;
  MaxForeignPoints = 8 * nspace;
  NumForeignNodes  = 0;
  NumForeignPoints = 0;

  sum_NumForeignNodes  = 0;
  sum_NumForeignPoints = 0;

  /* the following two arrays hold imported tree nodes and imported points to augment the local tree */
  Foreign_Nodes  = (ngbnode *)Mem.mymalloc_movable(&Foreign_Nodes, "Foreign_Nodes", MaxForeignNodes * sizeof(ngbnode));
  Foreign_Points = (foreign_sphpoint_data *)Mem.mymalloc_movable(&Foreign_Points, "Foreign_Points",
                                                                 MaxForeignPoints * sizeof(foreign_sphpoint_data));



  double tstart = Logs.second();

85
86
87
88
  int global_left_particles = 0;

  MPI_Allreduce(&ndensities, &global_left_particles, 1, MPI_INT, MPI_SUM, D->Communicator);

Volker Springel's avatar
Volker Springel committed
89
90
91
92
  do
    {
      double t0 = Logs.second();

93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
 /*  Since EntropyToInvGammaPred of remote particles can change, we have to import the particles in every iteration */

      MaxForeignNodes  = nspace;
      MaxForeignPoints = 8 * nspace;
      NumForeignNodes  = 0;
      NumForeignPoints = 0;

      sum_NumForeignNodes  = 0;
      sum_NumForeignPoints = 0;

      tree_initialize_leaf_node_access_info();

      max_ncycles = 0;

      prepare_shared_memory_access();

Volker Springel's avatar
Volker Springel committed
109
110
111
      /* now do the primary work with this call */
      densities_determine(ndensities, targetList);

112
113
114
115
      MPI_Allreduce(MPI_IN_PLACE, &max_ncycles, 1, MPI_INT, MPI_MAX, D->Communicator);

      cleanup_shared_memory_access();

Volker Springel's avatar
Volker Springel committed
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
      /* do final operations on results */

      double entropy_old;

      int npleft = 0;

      for(int i = 0; i < ndensities; i++)
        {
          int target = targetList[i];
          if(target >= 0)
            {
              if(Tp->P[target].getType() != 0)
                Terminate("P[target].getType() != 0");

              sph_particle_data *SphP = Tp->SphP;

              if(SphP[target].EntropyToInvGammaPred > 0 && SphP[target].Density > 0)
                {
                  entropy_old = SphP[target].Entropy;
                  SphP[target].PressureSphDensity /= SphP[target].EntropyToInvGammaPred;
                  SphP[target].Entropy =
                      GAMMA_MINUS1 * SphP[target].EntropyPred / pow(SphP[target].PressureSphDensity * All.cf_a3inv, GAMMA_MINUS1);
                  SphP[target].EntropyToInvGammaPred = pow(SphP[target].Entropy, 1.0 / GAMMA);
                }
              else
                {
                  entropy_old                        = SphP[target].Entropy;
                  SphP[target].PressureSphDensity    = 0;
                  SphP[target].Entropy               = 0;
                  SphP[target].EntropyToInvGammaPred = 0;
                }
              /* entropy has not converged yet */
              if(fabs(entropy_old - SphP[target].Entropy) > ENTROPY_TOLERANCE * entropy_old)
                targetList[npleft++] = target;
            }
        }

      ndensities = npleft;

155
156
      MPI_Allreduce(&ndensities, &global_left_particles, 1, MPI_INT, MPI_SUM, D->Communicator);

Volker Springel's avatar
Volker Springel committed
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
      double t1 = Logs.second();

      if(npleft > 0)
        {
          iter++;

          D->mpi_printf("SPH-INIT-ENTROPY: ngb iteration %4d: took %8.3f  , need to repeat for %012lld local particles.\n", iter,
                        Logs.timediff(t0, t1), npleft);

          if(iter > MAXITER)
            Terminate("failed to converge in neighbour iteration in density()\n");
        }
      else
        D->mpi_printf("SPH-INIT-ENTROPY: ngb iteration %4d: took %8.3f\n", ++iter, Logs.timediff(t0, t1));
    }
172
  while(global_left_particles > 0);
Volker Springel's avatar
Volker Springel committed
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209


  TIMER_STOP(CPU_DENSITY);


  /* free temporary buffers */

  Mem.myfree(Foreign_Points);
  Mem.myfree(Foreign_Nodes);

  Mem.myfree(targetList);

  double tb = Logs.second();

  D->mpi_printf("SPH-INIT-ENTROPY: entropy calculation is done. took: %8.3f\n", Logs.timediff(tstart, tb));
}

/*! \brief This function is used to find the initial entropy^invgamma for each SPH
 *  particle in the pressure-entropy formulation of SPH in case the ICs
 *  file contains internal energies.
 */
void sph::setup_entropy_to_invgamma(void)
{
  All.set_cosmo_factors_for_current_time();

  /* Initialize entropy and entropy^invgamma with a fist guess coming from standard SPH density estimate. */
  /* EntropyPred is untouched since it contains the internal energies needed for the iterative process; it */
  /* will be set in init.c to the correct value.                                                          */
  for(int i = 0; i < Tp->NumGas; i++)
    {
      Tp->SphP[i].Entropy = GAMMA_MINUS1 * Tp->SphP[i].EntropyPred / pow(Tp->SphP[i].Density * All.cf_a3inv, GAMMA_MINUS1);
      Tp->SphP[i].EntropyToInvGammaPred = pow(Tp->SphP[i].Entropy, 1.0 / GAMMA);
    }

  init_entropy();
}
#endif