Skip to content

Commit c6b1675

Browse files
llaniewskillaniewski
authored
Merge pull request CFD-GO#518 from CFD-GO/merge/develop
Merging `master` changes to `develop`
2 parents bd81db6 + 9ee8c06 commit c6b1675

25 files changed

+688
-320
lines changed

‎example/particle/3d/cube.xml

+26
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
<?xml version="1.0"?>
2+
<CLBConfig version="2.0" permissive="true" output="output/">
3+
<Units>
4+
<Param name="D" value="1m" gauge="16"/>
5+
<Param name="L" value="1x" gauge="4m"/>
6+
<Param name="DT" value="1" gauge="0.00004s"/>
7+
<Param name="rho" value="1kg/m3" gauge="1"/>
8+
</Units>
9+
<Geometry nx="1x" ny="1x" nz="1x" px="-0.5x" py="-0.5x" pz="-0.5x">
10+
<BGK><Box/></BGK>
11+
</Geometry>
12+
<Model>
13+
<Param name="aX_mean" value="100Pa/m"/>
14+
<Param name="nu" value="1m2/s"/>
15+
<RemoteForceInterface integrator="SIMPLEPART">
16+
<SimplePart>
17+
<Particle x="0" y="0" z="0" r="1" log="y" omegaz="10"/>
18+
<!-- 10m/s -->
19+
<Log Iterations="1" rotation="true"/>
20+
</SimplePart>
21+
</RemoteForceInterface>
22+
</Model>
23+
<VTK Iterations="1000" what="U,Solid"/>
24+
<Log Iterations="100"/>
25+
<Solve Iterations="10000"/>
26+
</CLBConfig>

‎example/particle/3d/roll_lb.xml

+22
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,22 @@
1+
<?xml version="1.0"?>
2+
<CLBConfig version="2.0" permissive="true" output="output/">
3+
<Units>
4+
<Param name="D" value="1m" gauge="64"/>
5+
<Param name="DT" value="1" gauge="0.00004s"/>
6+
<Param name="rho" value="1kg/m3" gauge="1"/>
7+
<!-- 1/s -->
8+
<!-- 40m/s 10m/s-->
9+
</Units>
10+
<Geometry nx="2m" ny="1m+2" nz="2m" px="-1.0m" py="-0.5m-1" pz="-1.0m">
11+
<BGK><Box/></BGK>
12+
<Wall><Box ny="1"/><Box ny="-1"/></Wall>
13+
</Geometry>
14+
<Model>
15+
<Param name="aX_mean" value="100Pa/m"/>
16+
<Param name="nu" value="1m2/s"/>
17+
<RemoteForceInterface integrator="SIMPLEPART"/>
18+
</Model>
19+
<VTK Iterations="100" what="U,Solid"/>
20+
<Log Iterations="100"/>
21+
<Solve Iterations="10000"/>
22+
</CLBConfig>

‎example/particle/3d/roll_sp.xml

+5
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
<SimplePart dt="0.00004">
2+
<Particle x="0" y="-0.4" z="0" r="0.1" log="y" vz="5" omegax="100"/>
3+
<Periodic x="2" z="2"/>
4+
<Log name="output/forces.csv" Iterations="1" rotation="true"/>
5+
</SimplePart>

‎src/CartLatticeAccess.hpp.Rt

+33-9
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,10 @@
1515
#include "CartLatticeContainer.h"
1616
#include "StorageConversions.h"
1717

18+
#ifndef NODE_SYMZ
19+
#define NODE_SYMZ 0
20+
#endif
21+
1822
/// Push all densities
1923

2024
<?R
@@ -398,21 +402,23 @@ public:
398402

399403
si = which(Fields$name == fn)
400404
sf = rows(Fields)[[si]]
401-
sd = d
402-
sd[i] = sd[i]*(-1)
405+
sd_plus = d
406+
sd_plus[i] = autosym_shift-sd_plus[i]
407+
sd_minus = d
408+
sd_minus[i] = -autosym_shift-sd_minus[i]
403409
?>
404410
template < class PARENT >
405411
template <class dx_t, class dy_t, class dz_t>
406412
CudaDeviceFunction real_t SymmetryAccess< PARENT >::<?%s paste0(this_fun, f$nicename) ?> (const dx_t & dx, const dy_t & dy, const dz_t & dz) const
407413
{
408414
<?R if (paste0("SYM",ch[i]) %in% NodeTypes$group) { ?>
409415
if (<?R C(d[i]) ?> > range_int<0>()) {
410-
if ((this->getNodeType() & NODE_SYM<?%s ch[i] ?>) == NODE_Symmetry<?%s ch[i] ?>_plus) {
411-
return <?%s paste0(sig, next_fun, sf$nicename) ?>(<?R C(sd,sep=", ") ?>);
416+
if ((this->getNodeType() & NODE_SYM<?%s ch[i] ?>) == NODE_<?%s autosym_name ?><?%s ch[i] ?>_plus) {
417+
return <?%s paste0(sig, next_fun, sf$nicename) ?>(<?R C(sd_plus,sep=", ",float=FALSE,wrap.const=range_int) ?>);
412418
}
413419
} else if (<?R C(d[i]) ?> < range_int<0>()) {
414-
if ((this->getNodeType() & NODE_SYM<?%s ch[i] ?>) == NODE_Symmetry<?%s ch[i] ?>_minus) {
415-
return <?%s paste0(sig, next_fun, sf$nicename) ?>(<?R C(sd,sep=", ") ?>);
420+
if ((this->getNodeType() & NODE_SYM<?%s ch[i] ?>) == NODE_<?%s autosym_name ?><?%s ch[i] ?>_minus) {
421+
return <?%s paste0(sig, next_fun, sf$nicename) ?>(<?R C(sd_minus,sep=", ",float=FALSE,wrap.const=range_int) ?>);
416422
}
417423
}
418424
<?R } ?>
@@ -425,9 +431,27 @@ CudaDeviceFunction real_t SymmetryAccess< PARENT >::<?%s paste0(this_fun, f$nice
425431
template < class PARENT >
426432
template <class N>
427433
CudaDeviceFunction void SymmetryAccess< PARENT >::pop<?%s s$suffix ?>(N & node) const
428-
{
429-
parent::pop<?%s s$suffix ?>(node);
430-
<?R resolve.symmetries(Density[s$load.densities,,drop=FALSE]) ?>
434+
{ <?R
435+
if (Options$autosym == 0) { ?>
436+
parent::pop<?%s s$suffix ?>(node); <?R
437+
} else if (Options$autosym == 1) { ?>
438+
parent::pop<?%s s$suffix ?>(node); <?R
439+
resolve.symmetries(Density[s$load.densities,,drop=FALSE])
440+
} else if (Options$autosym == 2) { ?>
441+
if (this->getNodeType() & (NODE_SYMX | NODE_SYMY | NODE_SYMZ)) { <?R
442+
dens = Density;
443+
dens$load = s$load.densities;
444+
for (d in rows(dens)) if (d$load) {
445+
f = rows(Fields)[[match(d$field, Fields$name)]]
446+
dp = c(-d$dx, -d$dy, -d$dz) ?>
447+
<?%s paste("node",d$name,sep=".") ?> = load_<?%s f$nicename ?>(range_int< <?%d dp[1] ?> >(),range_int< <?%d dp[2] ?> >(),range_int< <?%d dp[3] ?> >()); <?R
448+
} else if (!is.na(d$default)) { ?>
449+
<?%s paste("node",d$name,sep=".") ?> = <?%f d$default ?>; <?R
450+
} ?>
451+
} else {
452+
parent::pop<?%s s$suffix ?>(node);
453+
} <?R
454+
} else stop("Unknown autosym option") ?>
431455
}
432456
<?R } } ?>
433457

‎src/Handlers/GenericOptimizer.cpp

+1-1
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@ int GenericOptimizer::Init () {
66
int ret;
77
DEBUG_M;
88
ret = OptimizerInit();
9-
MPI_Bcast( &ret, 1, MPI_INT, 0, MPI_COMM_WORLD );
9+
MPI_Bcast( &ret, 1, MPI_INT, 0, solver->mpi_comm );
1010
if (ret) {
1111
ERROR("Failed to initialize Optimizer");
1212
return -1;

‎src/Handlers/OptimalControl.cpp

+2-2
Original file line numberDiff line numberDiff line change
@@ -88,7 +88,7 @@ int OptimalControl::Parameters (int type, double * tab) {
8888
return 0;
8989
case PAR_SET:
9090
output("Setting the params in the zone\n");
91-
MPI_Bcast(tab, Pars, MPI_DOUBLE, 0, MPI_COMM_WORLD);
91+
MPI_Bcast(tab, Pars, MPI_DOUBLE, 0, solver->mpi_comm);
9292
if (f != NULL) {
9393
fprintf(f,"SET");
9494
for (int i=0;i<Pars;i++) fprintf(f,",%lg",(double) tab[i]);
@@ -99,7 +99,7 @@ int OptimalControl::Parameters (int type, double * tab) {
9999
case PAR_GRAD:
100100
output("Getting gradient of a param in zone\n");
101101
solver->lattice->zSet.get_grad(par_index, zone_number, tmptab);
102-
MPI_Reduce(tmptab, tab, Pars, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
102+
MPI_Reduce(tmptab, tab, Pars, MPI_DOUBLE, MPI_SUM, 0, solver->mpi_comm);
103103
if (f != NULL) {
104104
fprintf(f,"GRAD");
105105
for (int i=0;i<Pars;i++) fprintf(f,",%lg",(double) tab[i]);

‎src/Handlers/acAndersen.cpp

+1-1
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@ double acAndersen::skal(real_t * a, real_t * b) {
1818
for (size_t k=0; k<n; k++) sum += a[k]*b[k];
1919
// TODO: MPI scatter gather
2020
double gsum=0;
21-
MPI_Allreduce ( &sum, &gsum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
21+
MPI_Allreduce ( &sum, &gsum, 1, MPI_DOUBLE, MPI_SUM, solver->mpi_comm );
2222
return gsum;
2323
}
2424

‎src/Handlers/acRemoteForceInterface.cpp

+29-1
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,8 @@
22
std::string acRemoteForceInterface::xmlname = "RemoteForceInterface";
33
#include "../HandlerFactory.h"
44

5+
#include <sstream>
6+
57
int acRemoteForceInterface::Init () {
68
Action::Init();
79
pugi::xml_attribute attr = node.attribute("integrator");
@@ -22,6 +24,26 @@ int acRemoteForceInterface::ConnectRemoteForceInterface(std::string integrator_)
2224
solver->lattice->RFI.setUnits(units[0],units[1],units[2]);
2325
solver->lattice->RFI.CanCopeWithUnits(false);
2426

27+
solver->lattice->RFI.setVar("output", solver->outpath);
28+
29+
30+
std::string element_content;
31+
int node_children = 0;
32+
for (pugi::xml_node par = node.first_child(); par; par = par.next_sibling()) {
33+
node_children ++;
34+
if (node_children > 1) {
35+
ERROR("Only a single element/CDATA allowed inside of a RemoteForceInterface xml element\n");
36+
return -1;
37+
}
38+
if ((par.type() == pugi::node_pcdata) || (par.type() == pugi::node_cdata)) {
39+
element_content = par.value();
40+
} else {
41+
std::stringstream ss;
42+
par.print(ss);
43+
element_content = ss.str();
44+
}
45+
}
46+
if (node_children > 0) solver->lattice->RFI.setVar("content", element_content);
2547
bool stats = false;
2648
std::string stats_prefix = solver->outpath;
2749
stats_prefix = stats_prefix + "_RFI";
@@ -56,7 +78,7 @@ int acRemoteForceInterface::ConnectRemoteForceInterface(std::string integrator_)
5678
bool use_box = true;
5779
attr = node.attribute("use_box");
5880
if (attr) use_box = attr.as_bool();
59-
81+
6082
if (use_box) {
6183
lbRegion reg = lattice->getLocalRegion();
6284
double px = lattice->px;
@@ -70,6 +92,12 @@ int acRemoteForceInterface::ConnectRemoteForceInterface(std::string integrator_)
7092
pz + reg.dz - PART_MAR_BOX,
7193
pz + reg.dz + reg.nz + PART_MAR_BOX);
7294
}
95+
96+
attr = node.attribute("omega");
97+
if (attr) solver->lattice->RFI_omega = attr.as_bool();
98+
attr = node.attribute("torque");
99+
if (attr) solver->lattice->RFI_torque = attr.as_bool();
100+
73101
MPI_Barrier(MPMD.local);
74102
lattice->RFI.Connect(MPMD.work,inter.work);
75103

‎src/Handlers/acSolve.cpp

+5-2
Original file line numberDiff line numberDiff line change
@@ -7,15 +7,18 @@ int acSolve::Init () {
77
if (GenericAction::ExecuteInternal()) return -1;
88
int stop=0;
99
do {
10-
int next_it = Next(solver->iter);
10+
int my_next_it = Next(solver->iter);
11+
int next_it = my_next_it;
1112
for (size_t i=0; i<solver->hands.size(); i++) {
1213
int it = solver->hands[i].Next(solver->iter);
1314
if ((it > 0) && (it < next_it)) next_it = it;
1415
}
1516
solver->steps = next_it;
1617
MPI_Bcast(&solver->steps, 1, MPI_INT, 0, MPMD.local);
1718
solver->iter += solver->steps;
18-
solver->lattice->Iterate(solver->steps, solver->iter_type);
19+
int iter_type = solver->iter_type;
20+
if (solver->steps == my_next_it) iter_type |= ITER_LASTGLOB;
21+
solver->lattice->Iterate(solver->steps, iter_type);
1922
CudaDeviceSynchronize();
2023
MPI_Barrier(MPMD.local);
2124
for (size_t i=0; i<solver->hands.size(); i++) {

‎src/Handlers/cbPID.cpp

+1-1
Original file line numberDiff line numberDiff line change
@@ -111,7 +111,7 @@ int cbPID::DoIt () {
111111
old_err = err;
112112

113113
}
114-
MPI_Bcast(&control, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
114+
MPI_Bcast(&control, 1, MPI_DOUBLE, 0, solver->mpi_comm);
115115
double nval;
116116
/* if (zone_number < 0) {
117117
sval = solver->lattice->zSet.get(setting, 0, (size_t) 0);

‎src/Handlers/cbRunR.cpp

+7
Original file line numberDiff line numberDiff line change
@@ -742,6 +742,8 @@ int cbRunR::Init() {
742742
output("%s\n",source.c_str());
743743
output("--------------------------------\n");
744744
}
745+
old_iter_type = solver->iter_type;
746+
solver->iter_type |= ITER_LASTGLOB;
745747
return 0;
746748
}
747749

@@ -788,6 +790,11 @@ int cbRunR::DoIt() {
788790
return 0;
789791
}
790792

793+
int cbRunR::Finish () {
794+
solver->iter_type = old_iter_type;
795+
return Callback::Finish();
796+
}
797+
791798

792799
#endif // WITH_R
793800

‎src/Handlers/cbRunR.h

+2
Original file line numberDiff line numberDiff line change
@@ -29,9 +29,11 @@ class cbRunR : public Callback {
2929
bool python;
3030
static int s_tag;
3131
int tag;
32+
int old_iter_type;
3233
public:
3334
int Init ();
3435
int DoIt ();
36+
int Finish ();
3537
};
3638

3739
#endif // WITH_R

‎src/Lattice.h.Rt

+1
Original file line numberDiff line numberDiff line change
@@ -74,6 +74,7 @@ public:
7474
real_t px, py, pz;
7575
MPIInfo mpi; ///< MPI information
7676
rfi_t RFI;
77+
bool RFI_omega, RFI_torque;
7778
solidcontainer_t SC;
7879
size_t particle_data_size_max;
7980
char snapFileName[STRING_LEN];

‎src/LatticeBase.cpp.Rt

+19-11
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,9 @@ LatticeBase::LatticeBase(int zonesettings, int zones, int num_snaps_, const Unit
1414
data.particle_data_size = 0;
1515
SC.balls = &RFI;
1616
RFI.name = "TCLB";
17+
RFI_omega = true;
18+
RFI_torque = true;
19+
1720
CudaStreamCreate(&kernelStream);
1821
CudaStreamCreate(&inStream);
1922
CudaStreamCreate(&outStream);
@@ -107,29 +110,34 @@ void LatticeBase::CopyInParticles() {
107110
CudaMalloc(&data.particle_data, RFI.mem_size());
108111
}
109112
data.particle_data_size = RFI.size();
110-
for (size_t i = 0; i < RFI.size(); i++) {
111-
RFI.RawData(i, RFI_DATA_FORCE + 0) = 0;
112-
RFI.RawData(i, RFI_DATA_FORCE + 1) = 0;
113-
RFI.RawData(i, RFI_DATA_FORCE + 2) = 0;
113+
for (int j=0; j<6; j++) for (size_t i=0; i<RFI.size(); i++) RFI.RawData(i, RFI_DATA_FORCE+j) = 0;
114+
if (!RFI_omega) for (int j=0; j<3; j++) for (size_t i=0; i<RFI.size(); i++) RFI.RawData(i, RFI_DATA_ANGVEL+j) = 0;
115+
if (RFI.mem_size() > 0) {
116+
CudaMemcpyAsync(data.particle_data, RFI.Particles(), RFI.mem_size(), CudaMemcpyHostToDevice, kernelStream);
114117
}
115-
if (RFI.mem_size() > 0) { CudaMemcpyAsync(data.particle_data, RFI.Particles(), RFI.mem_size(), CudaMemcpyHostToDevice, kernelStream); }
116118
DEBUG_PROF_PUSH("Tree Build");
117119
SC.Build();
118120
DEBUG_PROF_POP();
119121
SC.CopyToGPU(data.solidfinder, kernelStream);
120122
}
121123

122124
void LatticeBase::CopyOutParticles() {
123-
if (RFI.mem_size() > 0) { CudaMemcpyAsync(RFI.Particles(), data.particle_data, RFI.mem_size(), CudaMemcpyDeviceToHost, kernelStream); }
125+
if (RFI.mem_size() > 0) {
126+
CudaMemcpyAsync(RFI.Particles(), data.particle_data, RFI.mem_size(), CudaMemcpyDeviceToHost, kernelStream);
127+
}
124128
CudaStreamSynchronize(kernelStream);
125129
DEBUG_PROF_PUSH("Testing particles for NaNs");
126130
int nans = 0;
127-
for (size_t i = 0; i < RFI.size(); i++) {
128-
for (int j = 0; j < 3; j++) {
129-
if (!isfinite(RFI.RawData(i, RFI_DATA_FORCE + j))) {
130-
nans++;
131-
RFI.RawData(i, RFI_DATA_FORCE + j) = 0.0;
131+
for (int j=0; j<6; j++) {
132+
if (RFI_torque || (j<3)) {
133+
for (size_t i=0; i<RFI.size(); i++){
134+
if (! isfinite(RFI.RawData(i,RFI_DATA_FORCE+j))) {
135+
nans++;
136+
RFI.RawData(i,RFI_DATA_FORCE+j) = 0.0;
137+
}
132138
}
139+
} else {
140+
for (size_t i=0; i<RFI.size(); i++) RFI.RawData(i,RFI_DATA_FORCE+j) = 0.0;
133141
}
134142
}
135143
if (nans > 0) notice("%d NANs in particle forces (overwritten with 0.0)\n", nans);

‎src/LatticeBase.hpp

+1
Original file line numberDiff line numberDiff line change
@@ -64,6 +64,7 @@ class LatticeBase {
6464
solidcontainer_t SC; ///<
6565
size_t particle_data_size_max = 0; ///<
6666
rfi_t RFI; ///<
67+
bool RFI_omega, RFI_torque;
6768
SyntheticTurbulence ST; ///<
6869
std::string snapFileName;
6970

‎src/Particle.hpp

+3-3
Original file line numberDiff line numberDiff line change
@@ -30,9 +30,9 @@ struct ParticleI : Particle {
3030
angvel.x = particle_data[i*RFI_DATA_SIZE+RFI_DATA_ANGVEL+0];
3131
angvel.y = particle_data[i*RFI_DATA_SIZE+RFI_DATA_ANGVEL+1];
3232
angvel.z = particle_data[i*RFI_DATA_SIZE+RFI_DATA_ANGVEL+2];
33-
diff.x = pos.x - node[0];
34-
diff.y = pos.y - node[1];
35-
diff.z = pos.z - node[2];
33+
diff.x = node[0] - pos.x;
34+
diff.y = node[1] - pos.y;
35+
diff.z = node[2] - pos.z;
3636
dist = sqrt(diff.x*diff.x + diff.y*diff.y + diff.z*diff.z);
3737
cvel.x = vel.x + angvel.y*diff.z - angvel.z*diff.y;
3838
cvel.y = vel.y + angvel.z*diff.x - angvel.x*diff.z;

0 commit comments

Comments
 (0)