more normalization

This commit is contained in:
Steve Plimpton
2022-08-19 14:05:41 -06:00
parent 77c0ad4d26
commit 06e6a168f6

View File

@ -646,7 +646,7 @@ void FixAveGrid::end_of_step()
// for ATOM mode, perform ghost to owned grid communication
// done every sample for norm = SAMPLE, else once per Nfreq epoch
// nvalues + 1 for atom count
// nvalues+1 includes atom count
if (modeatom) {
if (dimension == 2)
@ -658,15 +658,11 @@ void FixAveGrid::end_of_step()
}
// for norm = SAMPLE, sum per-sample values and count to per-epoch
// for ATOM mode, also need to normalize per-sample values
// for ATOM mode, also need to normalize for single sample
// then return if irepeat < nrepeat
if (normflag == SAMPLE) {
if (modeatom) normalize_atom(1);
sum_sample_to_epoch();
if (irepeat < nrepeat && (normflag == ALL || normflag == NONORM)) {
@ -686,149 +682,24 @@ void FixAveGrid::end_of_step()
if (ngridout == 0) return;
// average the final results across Nrepeat samples
// for ATOM mode, result = total_value / total_count
// exception is DENSITY_NUMBER:
// result = value / (current binvol * Nrepeat)
// exception is DENSITY_MASS:
// result = (value * mv2d) / (current binvol * Nrepeat)
// exception is TEMPERATURE:
// result = (value * mvv2e) / (Nrepeat*cdof + adof*count) * boltz)
// for GRID mode, final is total value / Nrepeat
// average the final results for the entire epoch
// for norm = ALL, normalize values_sample by Nrepeat and counts
// for norm = SAMPLE, normalize values_epoch by Nrepeat
if (modeatom) {
double mvv2e = force->mvv2e;
double mv2d = force->mv2d;
double boltz = force->boltz;
double count,norm;
double repeat = nrepeat;
double invrepeat = 1.0/nrepeat;
double *prd = domain->prd;
double dx = prd[0]/nxgrid;
double dy = prd[1]/nygrid;
double dz = prd[2]/nzgrid;
double binvol;
if (dimension == 2) binvol = dx*dy;
else binvol = dx*dy*dz;
if (dimension == 2) {
if (nvalues == 1) {
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++) {
count = count2d[iy][ix];
if (count) {
if (which[0] == ArgInfo::DENSITY_NUMBER)
norm = 1.0 / (binvol * repeat);
else if (which[0] == ArgInfo::DENSITY_MASS)
norm = mv2d / (binvol * nrepeat);
else if (which[0] == ArgInfo::TEMPERATURE)
norm = mvv2e /((repeat*cdof + adof*count) * boltz);
else
norm = 1.0/count;
vec2d[iy][ix] *= norm;
count2d[iy][ix] *= invrepeat;
}
}
} else {
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++) {
count = count2d[iy][ix];
if (count) {
for (m = 0; m <= nvalues; m++) {
if (which[m] == ArgInfo::DENSITY_NUMBER)
norm = 1.0 / (binvol * repeat);
else if (which[m] == ArgInfo::DENSITY_MASS)
norm = mv2d / (binvol * nrepeat);
else if (which[m] == ArgInfo::TEMPERATURE)
norm = mvv2e /((repeat*cdof + adof*count) * boltz);
else
norm = 1.0/count;
array2d[iy][ix][m] *= norm;
}
count2d[iy][ix] *= invrepeat;
}
}
}
} else {
if (nvalues == 1) {
for (iz = nzlo_in; iz <= nzhi_in; iz++)
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++) {
count = count3d[iz][iy][ix];
if (count) {
if (which[0] == ArgInfo::DENSITY_NUMBER)
norm = 1.0 / (binvol * repeat);
else if (which[0] == ArgInfo::DENSITY_MASS)
norm = mv2d / (binvol * nrepeat);
else if (which[0] == ArgInfo::TEMPERATURE)
norm = mvv2e /((repeat*cdof + adof*count) * boltz);
else
norm = 1.0/count;
vec3d[iz][iy][ix] *= norm;
count3d[iz][iy][ix] *= invrepeat;
}
}
} else {
for (iz = nzlo_in; iz <= nzhi_in; iz++)
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++) {
count = count3d[iz][iy][ix];
if (count) {
for (m = 0; m <= nvalues; m++) {
if (which[m] == ArgInfo::DENSITY_NUMBER)
norm = 1.0 / (binvol * repeat);
else if (which[m] == ArgInfo::DENSITY_MASS)
norm = mv2d / (binvol * nrepeat);
else if (which[m] == ArgInfo::TEMPERATURE)
norm = mvv2e /((repeat*cdof + adof*count) * boltz);
else
norm = 1.0/count;
array3d[iz][iy][ix][m] *= norm;
}
count3d[iz][iy][ix] *= invrepeat;
}
}
}
}
if (normflag == ALL) normalize_atom(nrepeat);
else if (normflag == SAMPLE)
normalize_grid(vec2d_epoch,array2d_epoch,vec3d_epoch,array3d_epoch);
else if (normflag == NONORM) normalize_atom();
}
if (modegrid) {
double invrepeat = 1.0/nrepeat;
if (modegrid)
normalize_grid(vec2d_sample,array2d_sample,vec3d_sample,array3d_sample);
// create Nfreq output
if (dimension == 2) {
if (nvalues == 1) {
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++)
vec2d[iy][ix] *= invrepeat;
} else {
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++)
for (m = 0; m <= nvalues; m++)
array2d[iy][ix][m] *= invrepeat;
}
} else {
if (nvalues == 1) {
for (iz = nzlo_in; iz <= nzhi_in; iz++)
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++)
vec3d[iz][iy][ix] *= invrepeat;
} else {
for (iz = nzlo_in; iz <= nzhi_in; iz++)
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++)
for (m = 0; m <= nvalues; m++)
array3d[iz][iy][ix][m] *= invrepeat;
}
}
}
}
/* ----------------------------------------------------------------------
@ -1430,6 +1301,167 @@ void FixAveGrid::sum_sample_to_epoch()
}
}
/* ----------------------------------------------------------------------
normalize grid values and per-cell counts for ATOM mode for owned cells
flag = SAMPLE or EPOCH
set repeat = 1 for SAMPLE, Nrepeat for EPOCH
sample_value & count are either for one sample or summed over epoch
result = sample_value / count
exception is DENSITY_NUMBER:
result = value / (current binvol * Nrepeat)
exception is DENSITY_MASS:
result = (value * mv2d) / (current binvol * Nrepeat)
exception is TEMPERATURE:
result = (value * mvv2e) / (Nrepeat*cdof + adof*count) * boltz)
------------------------------------------------------------------------- */
void FixAveGrid::normalize_atom(int flag)
{
int ix,iy,iz,m;
double count,norm;
double repeat,invrepeat;
double mvv2e = force->mvv2e;
double mv2d = force->mv2d;
double boltz = force->boltz;
if (flag == EPOCH) {
repeat = nrepeat;
invrepeat = 1.0/nrepeat;
} else if (flag == SAMPLE) {
repeat = invrepeat = 1.0;
}
double *prd = domain->prd;
double dx = prd[0]/nxgrid;
double dy = prd[1]/nygrid;
double dz = prd[2]/nzgrid;
double binvol;
if (dimension == 2) binvol = dx*dy;
else binvol = dx*dy*dz;
if (dimension == 2) {
if (nvalues == 1) {
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++) {
count = count2d[iy][ix];
if (count) {
if (which[0] == ArgInfo::DENSITY_NUMBER)
norm = 1.0 / (binvol * repeat);
else if (which[0] == ArgInfo::DENSITY_MASS)
norm = mv2d / (binvol * nrepeat);
else if (which[0] == ArgInfo::TEMPERATURE)
norm = mvv2e /((repeat*cdof + adof*count) * boltz);
else
norm = 1.0/count;
vec2d_sample[iy][ix] *= norm;
count2d_sample[iy][ix] *= invrepeat;
}
}
} else {
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++) {
count = count2d[iy][ix];
if (count) {
for (m = 0; m <= nvalues; m++) {
if (which[m] == ArgInfo::DENSITY_NUMBER)
norm = 1.0 / (binvol * repeat);
else if (which[m] == ArgInfo::DENSITY_MASS)
norm = mv2d / (binvol * nrepeat);
else if (which[m] == ArgInfo::TEMPERATURE)
norm = mvv2e /((repeat*cdof + adof*count) * boltz);
else
norm = 1.0/count;
array2d_sample[iy][ix][m] *= norm;
}
count2d_sample[iy][ix] *= invrepeat;
}
}
}
} else {
if (nvalues == 1) {
for (iz = nzlo_in; iz <= nzhi_in; iz++)
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++) {
count = count3d[iz][iy][ix];
if (count) {
if (which[0] == ArgInfo::DENSITY_NUMBER)
norm = 1.0 / (binvol * repeat);
else if (which[0] == ArgInfo::DENSITY_MASS)
norm = mv2d / (binvol * nrepeat);
else if (which[0] == ArgInfo::TEMPERATURE)
norm = mvv2e /((repeat*cdof + adof*count) * boltz);
else
norm = 1.0/count;
vec3d_sample[iz][iy][ix] *= norm;
count3d_value[iz][iy][ix] *= invrepeat;
}
}
} else {
for (iz = nzlo_in; iz <= nzhi_in; iz++)
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++) {
count = count3d[iz][iy][ix];
if (count) {
for (m = 0; m <= nvalues; m++) {
if (which[m] == ArgInfo::DENSITY_NUMBER)
norm = 1.0 / (binvol * repeat);
else if (which[m] == ArgInfo::DENSITY_MASS)
norm = mv2d / (binvol * nrepeat);
else if (which[m] == ArgInfo::TEMPERATURE)
norm = mvv2e /((repeat*cdof + adof*count) * boltz);
else
norm = 1.0/count;
array3d_sample[iz][iy][ix][m] *= norm;
}
count3d_sample[iz][iy][ix] *= invrepeat;
}
}
}
}
}
/* ----------------------------------------------------------------------
normalize grid values by Nrepeat
used for GRID mode
used for ATOM MODE with norm = SAMPLE
------------------------------------------------------------------------- */
void FixAveGrid::normalize_grid(double **v2d, double ***a2d,
double ***v3d, double ****a3d)
{
int ix,iy,iz,m;
double invrepeat = 1.0/nrepeat;
if (dimension == 2) {
if (nvalues == 1) {
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++)
v2d[iy][ix] *= invrepeat;
} else {
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++)
for (m = 0; m <= nvalues; m++)
a2d[iy][ix][m] *= invrepeat;
}
} else {
if (nvalues == 1) {
for (iz = nzlo_in; iz <= nzhi_in; iz++)
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++)
v3d[iz][iy][ix] *= invrepeat;
} else {
for (iz = nzlo_in; iz <= nzhi_in; iz++)
for (iy = nylo_in; iy <= nyhi_in; iy++)
for (ix = nxlo_in; ix <= nxhi_in; ix++)
for (m = 0; m <= nvalues; m++)
a3d[iz][iy][ix][m] *= invrepeat;
}
}
}
/* ----------------------------------------------------------------------
pack ghost values into buf to send to another proc
nvalues per grid point + count