I am doing my homework to realize a C programe and during my work there is one step that I need to get the interger part of the double type numbers. So I choose ceil() or floor() in to realize this. But the output is unpredictable and far from expected.
The whole program is the following :
/*
************************************************************************
Includes
************************************************************************
*/
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>
#include <time.h>
#include <math.h>
/* Solvent Particle Propersities*/
typedef struct
{
double vx,vy,rx,ry ; /* velocity AND position */
double cell; /* index of grid */
}solvent_p;
/* Cell Propersities*/
typedef struct
{
double vcm_x,vcm_y ; /* center of mass velocity */
int number; /* number of solvent in the cell */
double roation_angle; /* roation_angle of the cell */
}cell_p;
/* periodic boundary condition judging and adjusting fuction PBC */
void PBC(solvent_p *sol)
{
double L = 20.0; // box size 20
if(sol->rx >20) sol->rx=sol->rx-L;
if(sol->rx < 0) sol->rx=sol->rx+L;
if(sol->ry >20) sol->ry=sol->ry-L;
if(sol->ry < 0) sol->ry=sol->ry+L;
}
int main(int argc, char **argv)
{
// Randome setup generates random numbers from GSL functions
const gsl_rng_type * T;
gsl_rng * r;
T = gsl_rng_default;
gsl_rng_default_seed = ((unsigned long)(time(NULL))); //设seed值为当前时间
r = gsl_rng_alloc (T);
solvent_p solvent[4000];
int i,j,k,index=0;
cell_p grid[400];
double alpha=90.0; //roation angle
/* Iniinitializing solvent
* All vx=0
* half vy = sqrt(3) and half vy=-sqrt(3) to make momentum zero and another requirement is the overall energy is 6000 equals energy of temperature=1 with 4000 solvent 3NkT/2 ,assuming mass of solvent = 1 ,which is all a test quantity
* rx and ry are random number generated by GSL library
* cell=20*(ry+rx) is an index of which cell the solvent is in
*/
for(i=0;i<=3999;i++)
{
solvent[i].vx=0.0;
if(i<=1999)
solvent[i].vy=sqrt(3);
else
solvent[i].vy=-sqrt(3);
solvent[i].rx =20.0 * gsl_rng_uniform_pos(r);
solvent[i].ry =20.0 * gsl_rng_uniform_pos(r);
//printf("%f \t %f \n",solvent[i].rx,solvent[i].ry);
solvent[i].cell=20*(floor(solvent[i].ry))+floor(solvent[i].rx)+1;
}
// grid setting up
for (i=0;i<=399;i++)
{
grid[i].vcm_x=0;
grid[i].vcm_y=0;
grid[i].number=0;
grid[i].roation_angle=0.0;
}
/* Begining Simulation Work
*
* Fist process is preparing the system to equilibrium for 10000 processes
*
* the whole process involving two steps steaming and collision and the two steps are conducted one by one in our simulation
* time duration for steaming is 0.1 which is assigned for Molecular Dynamics and time duration for collision is ignored
*
*
*/
for(i=0;i<=9999;i++)
{
double temp;
double delta_t_MD=0.1; //time step
temp=gsl_rng_uniform_pos(r);
double rand_rx = (temp < 0.5) ? temp : ((-1) * temp ); // randomshift rx;
temp=gsl_rng_uniform_pos(r);
double rand_ry = (temp < 0.5) ? temp : ((-1) * temp ); // randomshift ry
//steaming
for(j=0;j<=3999;j++)
{
//printf("%d \n",j);
printf("1 :%d \t %f \t %f \n",j,solvent[j].rx,solvent[j].ry);
solvent[j].rx=solvent[j].rx+solvent[j].vx * delta_t_MD;
solvent[j].ry=solvent[j].ry+solvent[j].vy * delta_t_MD;
printf("2: %d \t %f \t %f \n",j,solvent[j].rx,solvent[j].ry);
//randomshift
solvent[j].rx=solvent[j].rx+rand_rx;
solvent[j].ry=solvent[j].ry+rand_ry;
printf("3: %d \t %f \t %f \n",j,solvent[j].rx,solvent[j].ry);
// periodic boundary condition
PBC(&solvent[j]);
printf("4: %d \t %f \t %f \n",j,solvent[j].rx,solvent[j].ry);
solvent[j].cell=20*(ceil(solvent[j].ry)-1)+ceil(solvent[j].rx);
printf("5: %d \t %f \t %f \t %f \t %f \n",j,solvent[j].rx,solvent[j].ry,ceil(solvent[j].rx),ceil(solvent[j].ry));
index = (int)(solvent[j].cell);
//printf("%d \t %d \t %f \t %f \t %f \n",j,index,solvent[j].cell,ceil(solvent[j].rx),ceil(solvent[j].ry));
if ((index>=0) &&(index<=400))
{
grid[index].vcm_x=grid[index].vcm_x+solvent[i].vx;
grid[index].vcm_y=grid[index].vcm_y+solvent[i].vy;
grid[index].number=grid[index].number+1;
}
}
// caculating vcm
for (k=0;k<=399;k++)
{
if (grid[k].number >= 1)
{
grid[k].vcm_x=grid[k].vcm_x/grid[k].number;
grid[k].vcm_y=grid[k].vcm_y/grid[k].number;
}
double temp;
temp=gsl_rng_uniform_pos(r);
grid[k].roation_angle = (temp < 0.5) ? alpha : ((-1) * alpha );
}
//collsion
}
gsl_rng_free (r); // free RNG
return 0;
}
Sorry it is some extent long so I did not put in first.And something did not finished but the program framework is set up.
my programe is like this:
printf("4: %d \t %f \t %f \n",j,solvent[j].rx,solvent[i].ry);
solvent[j].cell=20*(floor(solvent[j].ry))+floor(solvent[j].rx)+1;
printf("5: %d \t %f \t %f \t %f \t %f \n",j,solvent[j].rx,solvent[i].ry,floor(solvent[j].rx),floor(solvent[j].ry));
And although something as I wanted something more is wrong and below is I choose some parts of the output:
4: 3993 3.851240 17.047031
5: 3993 3.851240 17.047031 3.000000 9.000000
although floor(solvent[j].rx) is right but floor(solvent[j].ry) is totally wrong.
And the final result of my program is
Segmentation fault (core dumped)
------------------
(program exited with code: 139)
How to fix this one? Is there anything I use was wrong?
And to further test the problem I tried ceil() function in and the program and result is like this
program:
printf("4: %d \t %f \t %f \n",j,solvent[j].rx,solvent[i].ry);
solvent[j].cell=20*(ceil(solvent[j].ry)-1)+ceil(solvent[j].rx);
printf("5: %d \t %f \t %f \t %f \t %f \n",j,solvent[j].rx,solvent[i].ry,ceil(solvent[j].rx),ceil(solvent[j].ry));
result:
2: 3999 14.571205 2.837654
4: 3999 14.571205 2.837654
5: 3999 14.571205 2.837654 15.000000 15.000000
So use the nearest output as an example to illustrate my desire result is:
a= 14.571205, ceil(a)=15.00 b= 2.837654 , ceil(b)=3.00 not the 15.000 in the output
And the problem worsening is that when I just use a and b to test ceil() everything seems perfect:
program:
#include <stdio.h>
#include <math.h>
int main()
{
double a= 14.571205;
double b= 2.837654 ;
printf("%f \t %f \n",ceil(a),ceil(b));
return 0;
}
output: 15.000000 3.000000
I am using GCC in Linux Ubuntu.
==============================================================================
The problem has been solved.
The real fatal problem is here
if ((index>=0) &&(index<=400))
{
grid[index].vcm_x=grid[index].vcm_x+solvent[i].vx;
grid[index].vcm_y=grid[index].vcm_y+solvent[i].vy;
grid[index].number=grid[index].number+1;
}
}
Both solvent[i].vy and solvent[i].vx should be changed i for j.
Just as @Jon and @Blckknght @Eric Lippert have pointed out.
And I obvious get in the wrong trap and mislead @ouah and @Blckknght and In fact, the output sentences do have problem but they do not caused the prorame to crash only the out of boundary will do that.
Thank you ALL!
And I do like to share Eric Lippert 's comment which is powerful and insightful:
More generally, the problem you have is called "select isn't broken" by experienced programmers. That is, you have done something completely wrong, you cannot figure out what you did wrong, and so you conclude that a piece of basic, easy infrastructure written by experts and tested extensively is wrong. I assure you, floor and ceil do exactly what they are supposed to do. They are not wrong. The problem lies somewhere in your code, not in the standard library. – Eric Lippert
