testing/icgen/random_verschillende_resoluties_N-GenIC/gadget3_64/conduction.c
void conduction(void)
{
  int i, iter;
  double delta0, delta1, alpha, beta, a3inv, dt, rel_change, loc_max_rel_change, glob_max_rel_change;
  double sumnew, sumold, sumtransfer, sumnew_tot, sumold_tot, sumtransfer_tot;

  if(ThisTask == 0)
    {
      printf("Start thermal conduction...\n");
      fflush(stdout);
    }

  Energy = (double *) mymalloc(N_gas * sizeof(double));
  EnergyOld = (double *) mymalloc(N_gas * sizeof(double));
  Residual = (double *) mymalloc(N_gas * sizeof(double));
  DVec = (double *) mymalloc(N_gas * sizeof(double));
  QVec = (double *) mymalloc(N_gas * sizeof(double));

  Kappa = (double *) mymalloc(N_gas * sizeof(double));


  if(All.ComovingIntegrationOn)
    a3inv = 1 / (All.Time * All.Time * All.Time);
  else
    a3inv = 1.0;


  dt = (All.Conduction_Ti_endstep - All.Conduction_Ti_begstep) * All.Timebase_interval;

  if(All.ComovingIntegrationOn)
    dt *= All.Time / hubble_function(All.Time);

  if(ThisTask == 0)
    {
      printf("dt=%g\n", dt);
    }

  /* First, let's compute the thermal energies per unit mass and 
     conductivities for all particles */

  for(i = 0; i < N_gas; i++)
    {
      if(P[i].Type == 0)
	{
	  /* this gives the thermal energy per unit mass for particle i */
	  Energy[i] = EnergyOld[i] = SphP[i].Entropy *
	    pow(SphP[i].d.Density * a3inv, GAMMA_MINUS1) / GAMMA_MINUS1;

#ifdef CONDUCTION_CONSTANT
	  Kappa[i] = All.ConductionCoeff;
#else
	  Kappa[i] = All.ConductionCoeff * pow(EnergyOld[i], 2.5);
#ifdef CONDUCTION_SATURATION
	  electron_free_path = All.ElectronFreePathFactor * Energy[i] * Energy[i] / (SphP[i].Density * a3inv);
	  temp_scale_length =
	    atime * fabs(smoothentr) / sqrt(gradentr[0] * gradentr[0] +
					    gradentr[1] * gradentr[1] + gradentr[2] * gradentr[2]);
	  Kappa[i] /= (1 + 4.2 * electron_free_path / temp_scale_length);
#endif
#endif

#ifdef SFR
	  if(SphP[i].Density * a3inv >= All.PhysDensThresh)
	    Kappa[i] = 0;
#endif

	  /* we'll factor the timestep into the conductivities, for simplicity */
	  Kappa[i] *= dt;
	}
    }



  /* Let's start the Conjugate Gradient Algorithm */

  /* Initialization */

  conduction_matrix_multiply(EnergyOld, Residual);

  for(i = 0; i < N_gas; i++)
    {
      if(P[i].Type == 0)
	{
	  Residual[i] = EnergyOld[i] - Residual[i];
	  DVec[i] = Residual[i];
	}
    }

  delta1 = conduction_vector_multiply(Residual, Residual);
  delta0 = delta1;

  iter = 0;			/* iteration counter */
  glob_max_rel_change = 1 + COND_ITER_ACCURACY;	/* to make sure that we enter the iteration */

  while(iter < MAX_COND_ITER && glob_max_rel_change > COND_ITER_ACCURACY && delta1 > 0)
    {
      conduction_matrix_multiply(DVec, QVec);

      alpha = delta1 / conduction_vector_multiply(DVec, QVec);

      for(i = 0, loc_max_rel_change = 0; i < N_gas; i++)
	{
	  Energy[i] += alpha * DVec[i];
	  Residual[i] -= alpha * QVec[i];

	  rel_change = alpha * DVec[i] / Energy[i];
	  if(loc_max_rel_change < rel_change)
	    loc_max_rel_change = rel_change;
	}

      delta0 = delta1;
      delta1 = conduction_vector_multiply(Residual, Residual);

      beta = delta1 / delta0;

      for(i = 0; i < N_gas; i++)
	DVec[i] = Residual[i] + beta * DVec[i];

      iter++;


      MPI_Allreduce(&loc_max_rel_change, &glob_max_rel_change, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);

      if(ThisTask == 0)
	{
	  printf("conduction iter=%d  delta1=%g delta1/delta0=%g  max-rel-change=%g\n",
		 iter, delta1, delta1 / delta0, glob_max_rel_change);
	  fflush(stdout);
	}
    }


  /* Now we have the solution vector in Energy[] */
  /* assign it to the entropies, and update the pressure */

  for(i = 0, sumnew = sumold = sumtransfer = 0; i < N_gas; i++)
    {
      if(P[i].Type == 0)
	{
	  sumnew += P[i].Mass * Energy[i];
	  sumold += P[i].Mass * EnergyOld[i];
	  sumtransfer += P[i].Mass * fabs(Energy[i] - EnergyOld[i]);

	  SphP[i].Entropy = Energy[i] / pow(SphP[i].d.Density * a3inv, GAMMA_MINUS1) * GAMMA_MINUS1;
	}
    }

  MPI_Allreduce(&sumnew, &sumnew_tot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
  MPI_Allreduce(&sumold, &sumold_tot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
  MPI_Allreduce(&sumtransfer, &sumtransfer_tot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);

  if(ThisTask == 0)
    {
      printf("\nconduction finished. energy_before=%g energy_after=%g  rel-change=%g rel-transfer=%g\n\n",
	     sumold_tot, sumnew_tot, (sumnew_tot - sumold_tot) / sumold_tot, sumtransfer_tot / sumold_tot);
      fflush(stdout);
    }

  myfree(Kappa);
  myfree(QVec);
  myfree(DVec);
  myfree(Residual);
  myfree(EnergyOld);
  myfree(Energy);
}
/* Note: the conduction equation we solve is really formulated with u instead of T, i.e.
   the factor (gamma-1)*mu*mp/k_B that converts from T to u is implicitely absorbed in a 
   redefinition of kappa */
the conduction equation we solve is really formulated with u instead of T, i.e.
the factor (gamma-1)*mu*mp/k_B that converts from T to u is implicitely absorbed in a
redefinition of kappa
the conduction equation we solve is really formulated with u instead of T i . e . the factor ( gamma - 1 ) * mu * mp / k_B that converts from T to u is implicitely absorbed in a redefinition of kappa
void conduction(void)
{
  int i, iter;
  double delta0, delta1, alpha, beta, a3inv, dt, rel_change, loc_max_rel_change, glob_max_rel_change;
  double sumnew, sumold, sumtransfer, sumnew_tot, sumold_tot, sumtransfer_tot;
  if(ThisTask == 0)
    {
      printf("Start thermal conduction...\n");
      fflush(stdout);
    }
  Energy = (double *) mymalloc(N_gas * sizeof(double));
  EnergyOld = (double *) mymalloc(N_gas * sizeof(double));
  Residual = (double *) mymalloc(N_gas * sizeof(double));
  DVec = (double *) mymalloc(N_gas * sizeof(double));
  QVec = (double *) mymalloc(N_gas * sizeof(double));
  Kappa = (double *) mymalloc(N_gas * sizeof(double));
  if(All.ComovingIntegrationOn)
    a3inv = 1 / (All.Time * All.Time * All.Time);
  else
    a3inv = 1.0;
  dt = (All.Conduction_Ti_endstep - All.Conduction_Ti_begstep) * All.Timebase_interval;
  if(All.ComovingIntegrationOn)
    dt *= All.Time / hubble_function(All.Time);
  if(ThisTask == 0)
    {
      printf("dt=%g\n", dt);
    }
  for(i = 0; i < N_gas; i++)
    {
      if(P[i].Type == 0)
	{
	  Energy[i] = EnergyOld[i] = SphP[i].Entropy *
	    pow(SphP[i].d.Density * a3inv, GAMMA_MINUS1) / GAMMA_MINUS1;
#ifdef CONDUCTION_CONSTANT
	  Kappa[i] = All.ConductionCoeff;
#else
	  Kappa[i] = All.ConductionCoeff * pow(EnergyOld[i], 2.5);
#ifdef CONDUCTION_SATURATION
	  electron_free_path = All.ElectronFreePathFactor * Energy[i] * Energy[i] / (SphP[i].Density * a3inv);
	  temp_scale_length =
	    atime * fabs(smoothentr) / sqrt(gradentr[0] * gradentr[0] +
					    gradentr[1] * gradentr[1] + gradentr[2] * gradentr[2]);
	  Kappa[i] /= (1 + 4.2 * electron_free_path / temp_scale_length);
#endif
#endif
#ifdef SFR
	  if(SphP[i].Density * a3inv >= All.PhysDensThresh)
	    Kappa[i] = 0;
#endif
	  Kappa[i] *= dt;
	}
    }
  conduction_matrix_multiply(EnergyOld, Residual);
  for(i = 0; i < N_gas; i++)
    {
      if(P[i].Type == 0)
	{
	  Residual[i] = EnergyOld[i] - Residual[i];
	  DVec[i] = Residual[i];
	}
    }
  delta1 = conduction_vector_multiply(Residual, Residual);
  delta0 = delta1;
  iter = 0;			
  glob_max_rel_change = 1 + COND_ITER_ACCURACY;	
  while(iter < MAX_COND_ITER && glob_max_rel_change > COND_ITER_ACCURACY && delta1 > 0)
    {
      conduction_matrix_multiply(DVec, QVec);
      alpha = delta1 / conduction_vector_multiply(DVec, QVec);
      for(i = 0, loc_max_rel_change = 0; i < N_gas; i++)
	{
	  Energy[i] += alpha * DVec[i];
	  Residual[i] -= alpha * QVec[i];
	  rel_change = alpha * DVec[i] / Energy[i];
	  if(loc_max_rel_change < rel_change)
	    loc_max_rel_change = rel_change;
	}
      delta0 = delta1;
      delta1 = conduction_vector_multiply(Residual, Residual);
      beta = delta1 / delta0;
      for(i = 0; i < N_gas; i++)
	DVec[i] = Residual[i] + beta * DVec[i];
      iter++;
      MPI_Allreduce(&loc_max_rel_change, &glob_max_rel_change, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
      if(ThisTask == 0)
	{
	  printf("conduction iter=%d  delta1=%g delta1/delta0=%g  max-rel-change=%g\n",
		 iter, delta1, delta1 / delta0, glob_max_rel_change);
	  fflush(stdout);
	}
    }
  for(i = 0, sumnew = sumold = sumtransfer = 0; i < N_gas; i++)
    {
      if(P[i].Type == 0)
	{
	  sumnew += P[i].Mass * Energy[i];
	  sumold += P[i].Mass * EnergyOld[i];
	  sumtransfer += P[i].Mass * fabs(Energy[i] - EnergyOld[i]);
	  SphP[i].Entropy = Energy[i] / pow(SphP[i].d.Density * a3inv, GAMMA_MINUS1) * GAMMA_MINUS1;
	}
    }
  MPI_Allreduce(&sumnew, &sumnew_tot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
  MPI_Allreduce(&sumold, &sumold_tot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
  MPI_Allreduce(&sumtransfer, &sumtransfer_tot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
  if(ThisTask == 0)
    {
      printf("\nconduction finished. energy_before=%g energy_after=%g  rel-change=%g rel-transfer=%g\n\n",
	     sumold_tot, sumnew_tot, (sumnew_tot - sumold_tot) / sumold_tot, sumtransfer_tot / sumold_tot);
      fflush(stdout);
    }
  myfree(Kappa);
  myfree(QVec);
  myfree(DVec);
  myfree(Residual);
  myfree(EnergyOld);
  myfree(Energy);
}
void conduction ( void ) { int i , iter ; double delta0 , delta1 , alpha , beta , a3inv , dt , rel_change , loc_max_rel_change , glob_max_rel_change ; double sumnew , sumold , sumtransfer , sumnew_tot , sumold_tot , sumtransfer_tot ; if ( ThisTask == 0 ) { printf ( " \n " ) ; fflush ( stdout ) ; } Energy = ( double * ) mymalloc ( N_gas * sizeof ( double ) ) ; EnergyOld = ( double * ) mymalloc ( N_gas * sizeof ( double ) ) ; Residual = ( double * ) mymalloc ( N_gas * sizeof ( double ) ) ; DVec = ( double * ) mymalloc ( N_gas * sizeof ( double ) ) ; QVec = ( double * ) mymalloc ( N_gas * sizeof ( double ) ) ; Kappa = ( double * ) mymalloc ( N_gas * sizeof ( double ) ) ; if ( All . ComovingIntegrationOn ) a3inv = 1 / ( All . Time * All . Time * All . Time ) ; else a3inv = 1.0 ; dt = ( All . Conduction_Ti_endstep - All . Conduction_Ti_begstep ) * All . Timebase_interval ; if ( All . ComovingIntegrationOn ) dt *= All . Time / hubble_function ( All . Time ) ; if ( ThisTask == 0 ) { printf ( " \n " , dt ) ; } for ( i = 0 ; i < N_gas ; i ++ ) { if ( P [ i ] . Type == 0 ) { Energy [ i ] = EnergyOld [ i ] = SphP [ i ] . Entropy * pow ( SphP [ i ] . d . Density * a3inv , GAMMA_MINUS1 ) / GAMMA_MINUS1 ; #ifdef CONDUCTION_CONSTANT Kappa [ i ] = All . ConductionCoeff ; #else Kappa [ i ] = All . ConductionCoeff * pow ( EnergyOld [ i ] , 2.5 ) ; #ifdef CONDUCTION_SATURATION electron_free_path = All . ElectronFreePathFactor * Energy [ i ] * Energy [ i ] / ( SphP [ i ] . Density * a3inv ) ; temp_scale_length = atime * fabs ( smoothentr ) / sqrt ( gradentr [ 0 ] * gradentr [ 0 ] + gradentr [ 1 ] * gradentr [ 1 ] + gradentr [ 2 ] * gradentr [ 2 ] ) ; Kappa [ i ] /= ( 1 + 4.2 * electron_free_path / temp_scale_length ) ; #endif #endif #ifdef SFR if ( SphP [ i ] . Density * a3inv >= All . PhysDensThresh ) Kappa [ i ] = 0 ; #endif Kappa [ i ] *= dt ; } } conduction_matrix_multiply ( EnergyOld , Residual ) ; for ( i = 0 ; i < N_gas ; i ++ ) { if ( P [ i ] . Type == 0 ) { Residual [ i ] = EnergyOld [ i ] - Residual [ i ] ; DVec [ i ] = Residual [ i ] ; } } delta1 = conduction_vector_multiply ( Residual , Residual ) ; delta0 = delta1 ; iter = 0 ; glob_max_rel_change = 1 + COND_ITER_ACCURACY ; while ( iter < MAX_COND_ITER && glob_max_rel_change > COND_ITER_ACCURACY && delta1 > 0 ) { conduction_matrix_multiply ( DVec , QVec ) ; alpha = delta1 / conduction_vector_multiply ( DVec , QVec ) ; for ( i = 0 , loc_max_rel_change = 0 ; i < N_gas ; i ++ ) { Energy [ i ] += alpha * DVec [ i ] ; Residual [ i ] -= alpha * QVec [ i ] ; rel_change = alpha * DVec [ i ] / Energy [ i ] ; if ( loc_max_rel_change < rel_change ) loc_max_rel_change = rel_change ; } delta0 = delta1 ; delta1 = conduction_vector_multiply ( Residual , Residual ) ; beta = delta1 / delta0 ; for ( i = 0 ; i < N_gas ; i ++ ) DVec [ i ] = Residual [ i ] + beta * DVec [ i ] ; iter ++ ; MPI_Allreduce ( & loc_max_rel_change , & glob_max_rel_change , 1 , MPI_DOUBLE , MPI_MAX , MPI_COMM_WORLD ) ; if ( ThisTask == 0 ) { printf ( " \n " , iter , delta1 , delta1 / delta0 , glob_max_rel_change ) ; fflush ( stdout ) ; } } for ( i = 0 , sumnew = sumold = sumtransfer = 0 ; i < N_gas ; i ++ ) { if ( P [ i ] . Type == 0 ) { sumnew += P [ i ] . Mass * Energy [ i ] ; sumold += P [ i ] . Mass * EnergyOld [ i ] ; sumtransfer += P [ i ] . Mass * fabs ( Energy [ i ] - EnergyOld [ i ] ) ; SphP [ i ] . Entropy = Energy [ i ] / pow ( SphP [ i ] . d . Density * a3inv , GAMMA_MINUS1 ) * GAMMA_MINUS1 ; } } MPI_Allreduce ( & sumnew , & sumnew_tot , 1 , MPI_DOUBLE , MPI_SUM , MPI_COMM_WORLD ) ; MPI_Allreduce ( & sumold , & sumold_tot , 1 , MPI_DOUBLE , MPI_SUM , MPI_COMM_WORLD ) ; MPI_Allreduce ( & sumtransfer , & sumtransfer_tot , 1 , MPI_DOUBLE , MPI_SUM , MPI_COMM_WORLD ) ; if ( ThisTask == 0 ) { printf ( " \n \n \n " , sumold_tot , sumnew_tot , ( sumnew_tot - sumold_tot ) / sumold_tot , sumtransfer_tot / sumold_tot ) ; fflush ( stdout ) ; } myfree ( Kappa ) ; myfree ( QVec ) ; myfree ( DVec ) ; myfree ( Residual ) ; myfree ( EnergyOld ) ; myfree ( Energy ) ; }
Note: the conduction equation we solve is really formulated with u instead of T, i.e.
Note : the conduction equation we solve is really formulated with u instead of T i . e .
/* First, let's compute the thermal energies per unit mass and 
     conductivities for all particles */ /* this gives the thermal energy per unit mass for particle i */ /* we'll factor the timestep into the conductivities, for simplicity */ /* Let's start the Conjugate Gradient Algorithm */ /* Initialization */ /* iteration counter */ /* to make sure that we enter the iteration */ /* Now we have the solution vector in Energy[] */ /* assign it to the entropies, and update the pressure */