/* Calculates time evolution of initial data (u, u_t)=(A*some function, B*some function) for cubic focusing Klein-Gordon equation u_tt-u_rr=-u+u3/r2 in radial 3d using centered time and centered space discretization. Saves (A,B) in blowup.dat, dispersion.dat, or indecisive.dat depending on outcome */ /***************************************************** ** Parameters **************************************** ******************************************************/ /* minimum amplitude A */ #define MINA 0.0 /* maximum amplitude A */ #define MAXA 7.0 /* stepsize A */ #define INCA 0.05 /* minimum amplitude B */ #define MINB -15.0 /* maximum amplitude B */ #define MAXB 15.0 /* stepsize B */ #define INCB 0.2 /* Grid mesh dr in multiples of DRBASE, see below */ #define DRMULT 40 /* Number of spatial gridpoints */ #define KD 2000 /* Courant constant dt/dr */ #define R 0.9 /* Time steps to calculate */ #define N 6000 /* threshold of L^infty norm for blow up */ #define MAX 100000.0 /* threshold for dispersion */ #define MIN 0.20 /* threshold for dispersion in time, i.e., this is the minimum number of time steps that are calculated until the evolution can be classified as "dispersion" */ #define DISPT 50 /* Maximal number of processes */ #define PROCS 50 /*********************************************************/ #define BLOWUP -1 #define INDEC 0 #define DISP 1 #include #include #include /* define resolution of soliton file */ #define DRBASE 1.0e-4 /* define grid mesh */ #define DR ((double)DRMULT*DRBASE) /****************** initial data ******************/ /* initial position */ double f(long k, double A) { double r=k*DR; /* compute radial coordinate from gridpoint */ return A*r*exp(-sqrt(1+r*r)); } /* initial velocity */ double g(long k, double B) { double r=k*DR; /* compute radial coordinate from gridpoint */ return B*r*sin(6.0*r*r)*exp(-0.1*(r*r-1.0)*(r*r-1.0)); } /***********************************************/ /* discretization of 1d Laplacian */ double D(double u[], long k) { return (u[k+1]-2.0*u[k]+u[k-1])/(DR*DR); } /* definition of the nonlinearity on the right-hand side of the equation */ double F(double u[], long k) { return -u[k]+u[k]*u[k]*u[k]/(k*DR*k*DR); } /* calculates the free energy as in Strauss-Vazquez */ double freeenergy(double u2[], double u1[], long K, double dt) { double en=0.0; /* stores 2 times the energy */ register long k; /* loop variable */ en=0.0; for(k=1; km) m=dum; } return m; } /* calculates time evolution, returns BLOWUP, DISP or INDEC */ int evolution(double A, double B, double u0[], double u1[], double u2[], long K, double dt) { register long k, n; /* loop variables */ double *ut; /* auxiliary pointer for cycling */ double baseenergynorm, basesupnorm; double Egroesse, Supgroesse; long klein; /* copy initial data to u0 */ for(k=1; k1.0) { return BLOWUP; } klein=0; /* calculate further time steps */ for(n=1; n<=N; n++) { /* calculate time step */ for(k=1; k MAX*baseenergynorm) || (Supgroesse*DR>1)) return BLOWUP; /* check for dispersion */ if( Egroesse < MIN*baseenergynorm ) klein=klein+1; if (klein>DISPT) return DISP; /* cycle arrays */ ut=u2; /* 2 -> t */ u2=u0; /* 0 -> 2 */ u0=u1; /* 1 -> 0 */ u1=ut; /* t -> 1 */ } /* if we ever reach this point then the time evolution was indecisive */ return INDEC; } int main(void) { double A, B; /* amplitudes */ int result; /* result of time evolution */ FILE *fdisp, *findec; /* file pointers */ pid_t pid; /* stores return value of fork() */ int status; /* stores status of wait */ long pc; /* blocksize, pointcounter */ double *u0, *u1, *u2; /* Field at 2 different time steps */ int myrank; /* process rank */ long K; /* number of space points actually calculated, equals KD+2N+2 */ double dt; /* time mesh */ /* Set K */ K=KD+N; /* set dt */ dt=R*DR; printf("Calculating %ld points...\n", (long)ceil((MAXA-MINA)/INCA*(MAXB-MINB)/INCB)); fflush(stdout); /* create processes */ for(myrank=0; myrank