/*
MODEL.H file for NFLDLAMBDA model.
Written by Gary Felder (gfelder@physics.stanford.edu)
Last Modified 3/7/01
This file contains the model specific functions and definitions for the model:
V = 1/4 lambda phi^4 + 1/2 g^2 phi^2 chi^2 + 1/2 h_i^2 chi^2 sigma_i^2
Comments: The inflaton potential is taken to be pure lambda phi^4. The inflaton decays (generally resonantly) to the single field chi, which in turn may decay to the fields sigma_i. None of the fields have bare masses and the theory is therefore conformal. In practice this means that after appropriate rescaling of the time, space, and field variables the only effect of the expansion of the universe on the field equations is via the term a''/a, which rapidly becomes negligible after inflation.
We use the default rescalings (see below), which give A=1/f0, B=sqrt(lambda) f0, r=1, s=-1
The potential in program units is given by V_{pr} = 1/4 phi_{pr}^4 + 1/2 gl phi_{pr}^2 chi_{pr}^2 + 1/2 hl_i chi_{pr}^2 sigma_{i,pr}^2
*/
/*
General comments about the model.h file:
This file contains the following functions - all called externally
modelinfo(FILE *info_) outputs information about the model and model-specific parameters to a file.
modelinitialize() performs any model-specific initialization
potential_energy(int term, float *field_values) calculates the average potential energy density term by term. The variable num_potential_terms (just above this function) specifies how many separate potential terms are used in this model.
dvdf(int fld, int i, int j, int k) calculates the potential term in the equation of motion, dV/dfield, for the field fld at the lattice point (i,j,k)
effective_mass(float mass_sq[], float *field_values) calculates the square masses of the fields and puts them in the array mass_sq. The parameter beginning tells the function to use initial field values - if this parameter is zero then the field quantities will be calculated dynamically.
model_output(int flush, char *ext_) allows each model to include its own specialized output function(s). The parameter flush is set to 1 when infrequent calculations are being performed and 0 otherwise. The string ext_ gives the extension for output filenames.
*/
/* - This section should be copied into parameters.h when using this model
// ---Adjustable parameters for NFLDLAMBDA model--- //
const float lambda=9.e-14; // Self-coupling of inflaton
const float gl=200; // Resonance parameter g^2/lambda
// The following array must contain at least as many entries as there are fields. Extra entries will be ignored.
const float hl[]={0.,0.,100.*gl,200.*gl}; // h_i^2/lambda (The first two values are just padding so that the indices match the field indices logically.)
*/
// Rescaling parameters.
// The program variables ("_pr") are defined as
// f_pr = rescale_A a^rescale_r f (f=field value)
// x_pr = rescale_B x (x=distance)
// dt_pr = rescale_B a^rescale_s dt (t=time)
// The constants beta, cpl, and f0 are used to set the variable rescalings rescale_A, rescale_B, rescale_r, and rescale_s.
// These rescaling constants may be reset independently; their settings in terms of beta, cpl, and f0 are suggestions. See the documentation for more details.
// These rescalings are intrinsic to the model and probably shouldn't be changed for individual runs. Adjustable parameters are stored in "parameters.h"
const float beta=4.; // Exponent of the dominant term in the potential
const float cpl=lambda; // Coefficient of the dominant term in the potential (up to numerical factors - see documentation)
const float f0=0.342; // Initial value of phi in Planck units, typically the point at which phi'=0
// By default these are automatically set to A=1/f0, B=sqrt(cpl) f0^(-1+beta/2), R=6/(2+beta), S=3(2-beta)/(2+beta). They may be adjusted to different values, but the relationship S=2R-3 must be maintained for the program equations to remain correct.
const float rescale_A=1./f0;
const float rescale_B=sqrt(cpl)*pow(f0,-1.+beta/2.);
const float rescale_r=6./(2.+beta);
// The value of S in terms of R SHOULD NOT be changed.
const float rescale_s=2.*rescale_r-3.;
// Other global variables
// The array model_vars is intended to hold any model-specific, non-constant global variables.
// These variables should be initialized in modelinitialize() below
// Even if you're not using any, num_model_vars should be at least 1 to keep some compilers happy.
const int num_model_vars=1;
// Model specific variables: None are defined for this model.
extern float model_vars[num_model_vars];
// Macros to make the equations more readable: The values of fld are 0=Phi,1=Chi,1+i=Sigma_i
#define PHI FIELD(0)
#define CHI FIELD(1)
#define SIGMAF FIELD(fld)
#define SIGMAN FIELD(n) // Sometimes sigma needs to be indexed in two different ways in one function
// Model specific details about the run to be output to an information file
inline void modelinfo(FILE *info_)
{
// Name and description of model
fprintf(info_,"N Field Lambda Model\n");
fprintf(info_,"V = 1/4 lambda phi^4 + 1/2 g^2 phi^2 chi^2 + 1/2 h_i^2 chi^2 sigma_i^2\n\n");
// Model specific parameter values
fprintf(info_,"gl (g^2/lambda) = %f\n",gl);
for(int fld=2;fldnumcouplings)
{
printf("Number of sigma fields (%d) exceeds number for which couplings have been defined (%d). Exiting\n",nflds-2,numcouplings-2);
exit(1);
}
}
}
// The constant num_potential_terms must be defined for use by outside functions
// Terms: term=0: 1/4 lambda phi^4 --- term=1: 1/2 g^2 phi^2 chi^2 --- term=2: Sum{i} of 1/2 h_i^2 chi^2 sigma_i^2
const int num_potential_terms=3; // Number of terms that are calculated separately in the potential
// Potential energy terms
// See documentation for normalization of these terms.
// When setting initial conditions field values will be supplied in the array field_values. Otherwise this function will calculate them on the lattice.
inline float potential_energy(int term, float *field_values)
{
DECLARE_INDICES
int fld;
float potential=0.;
float sigma_chi_sq[nflds]; // Average values of sigma_i^2 chi^2 (the first 2 spots in this array are not used)
if(field_values==NULL) // If no values are given calculate averages on the lattice
{
// For calculating term 2 initialize sigma_chi_sq array to 0
if(term==2)
for(fld=2;fld)
hlsigma_sq += hl[fld]*fldsqrd[fld];
mass_sq[1] = gl*fldsqrd[0] + hlsigma_sq;
for(fld=2;fld0) // If there's no expansion don't bother with this.
{
correction = pow(a,2.*rescale_s+2.);
for(fld=0;fld