Commit b9230376 authored by Jibril B. Coulibaly's avatar Jibril B. Coulibaly
Browse files

change and restore mass using scaling for floating point accuracy

parent 4a2ab6d2
Loading
Loading
Loading
Loading
+10 −13
Original line number Diff line number Diff line
@@ -602,13 +602,11 @@ void FixAdapt::change_settings()
	
	for (i = 0; i < nall; i++) {
	  if (mask[i] & groupbit) {
	    if (discflag) density = rmass[i] / (MY_PI * radius[i]*radius[i]);
	    else density = rmass[i] / (4.0*MY_PI/3.0 *
				       radius[i]*radius[i]*radius[i]);
      if(!scaleflag) scale = 0.5*value / radius[i];
	    if (scaleflag) radius[i] *= scale;
	    else radius[i] = 0.5*value;
	    if (discflag) rmass[i] = MY_PI * radius[i]*radius[i] * density;
	    else rmass[i] = 4.0*MY_PI/3.0 * radius[i]*radius[i]*radius[i] * density;
	    if (discflag) rmass[i] *= scale*scale;
	    else rmass[i] *= scale*scale*scale;
	  }
	}
	  
@@ -694,7 +692,7 @@ void FixAdapt::restore_settings()

    } else if (ad->which == ATOM) {
      if (diamflag) {
        double density;
        double density,scale;

        double *vec = fix_diam->vstore;
        double *radius = atom->radius;
@@ -702,15 +700,14 @@ void FixAdapt::restore_settings()
        int *mask = atom->mask;
        int nlocal = atom->nlocal;
        
        if (scaleflag) scale = previous_diam_scale;

        for (int i = 0; i < nlocal; i++)
          if (mask[i] & groupbit) {
            if (discflag) density = rmass[i] / (MY_PI * radius[i]*radius[i]);
            else density = rmass[i] / (4.0*MY_PI/3.0 *
                                       radius[i]*radius[i]*radius[i]);
            if(!scaleflag) scale = vec[i] / radius[i];
            radius[i] = vec[i];
            if (discflag) rmass[i] = MY_PI * radius[i]*radius[i] * density;
            else rmass[i] = 4.0*MY_PI/3.0 *
                            radius[i]*radius[i]*radius[i] * density;
            if (discflag) rmass[i] *= scale*scale;
            else rmass[i] *= scale*scale*scale;
          }
      }
      if (chgflag) {