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