1 // clang-format off
2 /* ----------------------------------------------------------------------
3    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4    https://www.lammps.org/, Sandia National Laboratories
5    Steve Plimpton, sjplimp@sandia.gov
6 
7    Copyright (2003) Sandia Corporation.  Under the terms of Contract
8    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9    certain rights in this software.  This software is distributed under
10    the GNU General Public License.
11 
12    See the README file in the top-level LAMMPS directory.
13 ------------------------------------------------------------------------- */
14 
15 #include "compute_reduce_region.h"
16 
17 #include "arg_info.h"
18 #include "atom.h"
19 #include "domain.h"
20 #include "error.h"
21 #include "fix.h"
22 #include "group.h"
23 #include "input.h"
24 #include "memory.h"
25 #include "modify.h"
26 #include "region.h"
27 #include "update.h"
28 #include "variable.h"
29 
30 using namespace LAMMPS_NS;
31 
32 #define BIG 1.0e20
33 
34 /* ---------------------------------------------------------------------- */
35 
ComputeReduceRegion(LAMMPS * lmp,int narg,char ** arg)36 ComputeReduceRegion::ComputeReduceRegion(LAMMPS *lmp, int narg, char **arg) :
37   ComputeReduce(lmp, narg, arg) {}
38 
39 /* ----------------------------------------------------------------------
40    calculate reduced value for one input M and return it
41    if flag = -1:
42      sum/min/max/ave all values in vector
43      for per-atom quantities, limit to atoms in group and region
44      if mode = MIN or MAX, also set index to which vector value wins
45    if flag >= 0: simply return vector[flag]
46 ------------------------------------------------------------------------- */
47 
compute_one(int m,int flag)48 double ComputeReduceRegion::compute_one(int m, int flag)
49 {
50   int i;
51 
52   Region *region = domain->regions[iregion];
53   region->prematch();
54 
55   // invoke the appropriate attribute,compute,fix,variable
56   // compute scalar quantity by summing over atom scalars
57   // only include atoms in group
58 
59   index = -1;
60   double **x = atom->x;
61   int *mask = atom->mask;
62   int nlocal = atom->nlocal;
63 
64   int n = value2index[m];
65 
66   // initialization in case it has not yet been run,
67   // e.g. when invoked
68   if (n == ArgInfo::UNKNOWN) {
69     init();
70     n = value2index[m];
71   }
72 
73   int j = argindex[m];
74 
75   double one = 0.0;
76   if (mode == MINN) one = BIG;
77   if (mode == MAXX) one = -BIG;
78 
79   if (which[m] == ArgInfo::X) {
80     if (flag < 0) {
81       for (i = 0; i < nlocal; i++)
82         if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
83           combine(one,x[i][j],i);
84     } else one = x[flag][j];
85   } else if (which[m] == ArgInfo::V) {
86     double **v = atom->v;
87     if (flag < 0) {
88       for (i = 0; i < nlocal; i++)
89         if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
90           combine(one,v[i][j],i);
91     } else one = v[flag][j];
92   } else if (which[m] == ArgInfo::F) {
93     double **f = atom->f;
94     if (flag < 0) {
95       for (i = 0; i < nlocal; i++)
96         if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
97           combine(one,f[i][j],i);
98     } else one = f[flag][j];
99 
100   // invoke compute if not previously invoked
101 
102   } else if (which[m] == ArgInfo::COMPUTE) {
103     Compute *compute = modify->compute[n];
104 
105     if (flavor[m] == PERATOM) {
106       if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) {
107         compute->compute_peratom();
108         compute->invoked_flag |= Compute::INVOKED_PERATOM;
109       }
110 
111       if (j == 0) {
112         double *compute_vector = compute->vector_atom;
113         if (flag < 0) {
114           for (i = 0; i < nlocal; i++)
115             if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
116               combine(one,compute_vector[i],i);
117         } else one = compute_vector[flag];
118       } else {
119         double **compute_array = compute->array_atom;
120         int jm1 = j - 1;
121         if (flag < 0) {
122           for (i = 0; i < nlocal; i++)
123             if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
124               combine(one,compute_array[i][jm1],i);
125         } else one = compute_array[flag][jm1];
126       }
127 
128     } else if (flavor[m] == LOCAL) {
129       if (!(compute->invoked_flag & Compute::INVOKED_LOCAL)) {
130         compute->compute_local();
131         compute->invoked_flag |= Compute::INVOKED_LOCAL;
132       }
133 
134       if (j == 0) {
135         double *compute_vector = compute->vector_local;
136         if (flag < 0)
137           for (i = 0; i < compute->size_local_rows; i++)
138             combine(one,compute_vector[i],i);
139         else one = compute_vector[flag];
140       } else {
141         double **compute_array = compute->array_local;
142         int jm1 = j - 1;
143         if (flag < 0)
144           for (i = 0; i < compute->size_local_rows; i++)
145             combine(one,compute_array[i][jm1],i);
146         else one = compute_array[flag][jm1];
147       }
148     }
149 
150   // check if fix frequency is a match
151 
152   } else if (which[m] == ArgInfo::FIX) {
153     if (update->ntimestep % modify->fix[n]->peratom_freq)
154       error->all(FLERR,"Fix used in compute reduce not computed at "
155                  "compatible time");
156     Fix *fix = modify->fix[n];
157 
158     if (flavor[m] == PERATOM) {
159       if (j == 0) {
160         double *fix_vector = fix->vector_atom;
161         if (flag < 0) {
162           for (i = 0; i < nlocal; i++)
163             if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
164               combine(one,fix_vector[i],i);
165         } else one = fix_vector[flag];
166       } else {
167         double **fix_array = fix->array_atom;
168         int jm1 = j - 1;
169         if (flag < 0) {
170           for (i = 0; i < nlocal; i++)
171             if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
172               combine(one,fix_array[i][jm1],i);
173         } else one = fix_array[flag][jm1];
174       }
175 
176     } else if (flavor[m] == LOCAL) {
177       if (j == 0) {
178         double *fix_vector = fix->vector_local;
179         if (flag < 0)
180           for (i = 0; i < fix->size_local_rows; i++)
181             combine(one,fix_vector[i],i);
182         else one = fix_vector[flag];
183       } else {
184         double **fix_array = fix->array_local;
185         int jm1 = j - 1;
186         if (flag < 0)
187           for (i = 0; i < fix->size_local_rows; i++)
188             combine(one,fix_array[i][jm1],i);
189         else one = fix_array[flag][jm1];
190       }
191     }
192 
193   // evaluate atom-style variable
194 
195   } else if (which[m] == ArgInfo::VARIABLE) {
196     if (atom->nmax > maxatom) {
197       maxatom = atom->nmax;
198       memory->destroy(varatom);
199       memory->create(varatom,maxatom,"reduce/region:varatom");
200     }
201 
202     input->variable->compute_atom(n,igroup,varatom,1,0);
203     if (flag < 0) {
204       for (i = 0; i < nlocal; i++)
205         if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2]))
206           combine(one,varatom[i],i);
207     } else one = varatom[flag];
208   }
209 
210   return one;
211 }
212 
213 /* ---------------------------------------------------------------------- */
214 
count(int m)215 bigint ComputeReduceRegion::count(int m)
216 {
217   int n = value2index[m];
218 
219   if (which[m] == ArgInfo::X || which[m] == ArgInfo::V || which[m] == ArgInfo::F)
220     return group->count(igroup,iregion);
221   else if (which[m] == ArgInfo::COMPUTE) {
222     Compute *compute = modify->compute[n];
223     if (flavor[m] == PERATOM) {
224       return group->count(igroup,iregion);
225     } else if (flavor[m] == LOCAL) {
226       bigint ncount = compute->size_local_rows;
227       bigint ncountall;
228       MPI_Allreduce(&ncount,&ncountall,1,MPI_DOUBLE,MPI_SUM,world);
229       return ncountall;
230     }
231   } else if (which[m] == ArgInfo::FIX) {
232     Fix *fix = modify->fix[n];
233     if (flavor[m] == PERATOM) {
234       return group->count(igroup,iregion);
235     } else if (flavor[m] == LOCAL) {
236       bigint ncount = fix->size_local_rows;
237       bigint ncountall;
238       MPI_Allreduce(&ncount,&ncountall,1,MPI_DOUBLE,MPI_SUM,world);
239       return ncountall;
240     }
241   } else if (which[m] == ArgInfo::VARIABLE)
242     return group->count(igroup,iregion);
243 
244   bigint dummy = 0;
245   return dummy;
246 }
247