Merging master changes to develop

example/particle/3d/cube.xml

@@ -0,0 +1,26 @@
+<?xml version="1.0"?>
+<CLBConfig version="2.0" permissive="true" output="output/">
+	<Units>
+		<Param name="D" value="1m" gauge="16"/>
+		<Param name="L" value="1x" gauge="4m"/>
+		<Param name="DT" value="1" gauge="0.00004s"/>
+		<Param name="rho" value="1kg/m3" gauge="1"/>
+	</Units>
+	<Geometry nx="1x" ny="1x" nz="1x" px="-0.5x" py="-0.5x" pz="-0.5x">
+		<BGK><Box/></BGK>
+	</Geometry>
+	<Model>
+		<Param name="aX_mean" value="100Pa/m"/>
+		<Param name="nu" value="1m2/s"/>
+        <RemoteForceInterface integrator="SIMPLEPART">
+			<SimplePart>
+				<Particle x="0" y="0" z="0" r="1" log="y" omegaz="10"/>
+				<!-- 10m/s -->
+				<Log Iterations="1" rotation="true"/>
+			</SimplePart>			
+		</RemoteForceInterface>
+	</Model>
+	<VTK Iterations="1000" what="U,Solid"/>
+	<Log Iterations="100"/>
+	<Solve Iterations="10000"/>
+</CLBConfig>

example/particle/3d/roll_lb.xml

@@ -0,0 +1,22 @@
+<?xml version="1.0"?>
+<CLBConfig version="2.0" permissive="true" output="output/">
+	<Units>
+		<Param name="D" value="1m" gauge="64"/>
+		<Param name="DT" value="1" gauge="0.00004s"/>
+		<Param name="rho" value="1kg/m3" gauge="1"/>
+		<!-- 1/s -->
+		<!-- 40m/s 10m/s-->
+	</Units>
+	<Geometry nx="2m" ny="1m+2" nz="2m" px="-1.0m" py="-0.5m-1" pz="-1.0m">
+		<BGK><Box/></BGK>
+		<Wall><Box ny="1"/><Box ny="-1"/></Wall>
+	</Geometry>
+	<Model>
+		<Param name="aX_mean" value="100Pa/m"/>
+		<Param name="nu" value="1m2/s"/>
+        <RemoteForceInterface integrator="SIMPLEPART"/>
+	</Model>
+	<VTK Iterations="100" what="U,Solid"/>
+	<Log Iterations="100"/>
+	<Solve Iterations="10000"/>
+</CLBConfig>

example/particle/3d/roll_sp.xml

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

src/CartLatticeAccess.hpp.Rt

@@ -15,6 +15,10 @@
 #include "CartLatticeContainer.h"
 #include "StorageConversions.h"
 
+#ifndef NODE_SYMZ
+  #define NODE_SYMZ 0
+#endif
+
 /// Push all densities
 
 <?R
@@ -398,21 +402,23 @@ public:
 
     si = which(Fields$name == fn)
     sf = rows(Fields)[[si]]
-    sd = d
-    sd[i] = sd[i]*(-1)
+    sd_plus = d
+    sd_plus[i] = autosym_shift-sd_plus[i]
+    sd_minus = d
+    sd_minus[i] = -autosym_shift-sd_minus[i]
 ?>
 template < class PARENT >
 template <class dx_t, class dy_t, class dz_t>
 CudaDeviceFunction real_t SymmetryAccess< PARENT >::<?%s paste0(this_fun, f$nicename) ?> (const dx_t & dx, const dy_t & dy, const dz_t & dz) const
 {
   <?R if (paste0("SYM",ch[i]) %in% NodeTypes$group) { ?>
   if (<?R C(d[i]) ?> > range_int<0>()) {
-    if ((this->getNodeType() & NODE_SYM<?%s ch[i] ?>) == NODE_Symmetry<?%s ch[i] ?>_plus) {
-      return <?%s paste0(sig, next_fun, sf$nicename) ?>(<?R C(sd,sep=", ") ?>);
+    if ((this->getNodeType() & NODE_SYM<?%s ch[i] ?>) == NODE_<?%s autosym_name ?><?%s ch[i] ?>_plus) {
+      return <?%s paste0(sig, next_fun, sf$nicename) ?>(<?R C(sd_plus,sep=", ",float=FALSE,wrap.const=range_int) ?>);
     }
   } else if (<?R C(d[i]) ?> < range_int<0>()) {
-    if ((this->getNodeType() & NODE_SYM<?%s ch[i] ?>) == NODE_Symmetry<?%s ch[i] ?>_minus) {
-      return <?%s paste0(sig, next_fun, sf$nicename) ?>(<?R C(sd,sep=", ") ?>);
+    if ((this->getNodeType() & NODE_SYM<?%s ch[i] ?>) == NODE_<?%s autosym_name ?><?%s ch[i] ?>_minus) {
+      return <?%s paste0(sig, next_fun, sf$nicename) ?>(<?R C(sd_minus,sep=", ",float=FALSE,wrap.const=range_int) ?>);
     }
   }
   <?R } ?>
@@ -425,9 +431,27 @@ CudaDeviceFunction real_t SymmetryAccess< PARENT >::<?%s paste0(this_fun, f$nice
 template < class PARENT >
 template <class N>
 CudaDeviceFunction void SymmetryAccess< PARENT >::pop<?%s s$suffix ?>(N & node) const
-{
-  parent::pop<?%s s$suffix ?>(node);
-  <?R resolve.symmetries(Density[s$load.densities,,drop=FALSE]) ?>
+{ <?R
+    if (Options$autosym == 0) { ?>
+      parent::pop<?%s s$suffix ?>(node); <?R
+    } else if (Options$autosym == 1) { ?>
+      parent::pop<?%s s$suffix ?>(node); <?R
+      resolve.symmetries(Density[s$load.densities,,drop=FALSE])
+    } else if (Options$autosym == 2) { ?>
+      if (this->getNodeType() & (NODE_SYMX | NODE_SYMY | NODE_SYMZ)) { <?R
+        dens = Density;
+        dens$load = s$load.densities;
+        for (d in rows(dens)) if (d$load) {
+          f = rows(Fields)[[match(d$field, Fields$name)]]
+          dp = c(-d$dx, -d$dy, -d$dz) ?>
+          <?%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
+        } else if (!is.na(d$default)) { ?>
+          <?%s paste("node",d$name,sep=".") ?> = <?%f d$default ?>; <?R
+        } ?>
+      } else {
+        parent::pop<?%s s$suffix ?>(node);
+      } <?R
+    } else stop("Unknown autosym option") ?>
 }
 <?R } } ?>
 

src/Handlers/GenericOptimizer.cpp

@@ -6,7 +6,7 @@ int GenericOptimizer::Init () {
 		int ret;
 		DEBUG_M;
 		ret = OptimizerInit();
-		MPI_Bcast( &ret, 1, MPI_INT, 0, MPI_COMM_WORLD );
+		MPI_Bcast( &ret, 1, MPI_INT, 0, solver->mpi_comm );
 		if (ret) {
 			ERROR("Failed to initialize Optimizer");
 			return -1;

src/Handlers/OptimalControl.cpp

@@ -88,7 +88,7 @@ int OptimalControl::Parameters (int type, double * tab) {
 			return 0;
 		case PAR_SET:
 			output("Setting the params in the zone\n");
-			MPI_Bcast(tab, Pars, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+			MPI_Bcast(tab, Pars, MPI_DOUBLE, 0, solver->mpi_comm);
 			if (f != NULL) {
 				fprintf(f,"SET");
 				for (int i=0;i<Pars;i++) fprintf(f,",%lg",(double) tab[i]);
@@ -99,7 +99,7 @@ int OptimalControl::Parameters (int type, double * tab) {
 		case PAR_GRAD:
 			output("Getting gradient of a param in zone\n");
 			solver->lattice->zSet.get_grad(par_index, zone_number, tmptab);
-			MPI_Reduce(tmptab, tab, Pars, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+			MPI_Reduce(tmptab, tab, Pars, MPI_DOUBLE, MPI_SUM, 0, solver->mpi_comm);
 			if (f != NULL) {
 				fprintf(f,"GRAD");
 				for (int i=0;i<Pars;i++) fprintf(f,",%lg",(double) tab[i]);

src/Handlers/acAndersen.cpp

@@ -18,7 +18,7 @@ double acAndersen::skal(real_t * a, real_t * b) {
         for (size_t k=0; k<n; k++) sum += a[k]*b[k];
         // TODO: MPI scatter gather
         double gsum=0;
-        MPI_Allreduce ( &sum, &gsum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
+        MPI_Allreduce ( &sum, &gsum, 1, MPI_DOUBLE, MPI_SUM, solver->mpi_comm );
         return gsum;
 }
 

src/Handlers/acRemoteForceInterface.cpp

@@ -2,6 +2,8 @@
 std::string acRemoteForceInterface::xmlname = "RemoteForceInterface";
 #include "../HandlerFactory.h"
 
+#include <sstream>
+
 int acRemoteForceInterface::Init () {
         Action::Init();
         pugi::xml_attribute attr = node.attribute("integrator");
@@ -22,6 +24,26 @@ int acRemoteForceInterface::ConnectRemoteForceInterface(std::string integrator_)
         solver->lattice->RFI.setUnits(units[0],units[1],units[2]);
         solver->lattice->RFI.CanCopeWithUnits(false);
 
+        solver->lattice->RFI.setVar("output", solver->outpath);
+
+
+        std::string element_content;
+        int node_children = 0;
+        for (pugi::xml_node par = node.first_child(); par; par = par.next_sibling()) {
+          node_children ++;
+          if (node_children > 1) {
+              ERROR("Only a single element/CDATA allowed inside of a RemoteForceInterface xml element\n");
+              return -1;
+          }
+		      if ((par.type() == pugi::node_pcdata) || (par.type() == pugi::node_cdata)) {
+            element_content = par.value();
+		      } else {
+            std::stringstream ss;
+            par.print(ss);
+            element_content = ss.str();
+          }
+	      }
+        if (node_children > 0) solver->lattice->RFI.setVar("content", element_content);
         bool stats = false;
         std::string stats_prefix = solver->outpath;
         stats_prefix = stats_prefix + "_RFI";
@@ -56,7 +78,7 @@ int acRemoteForceInterface::ConnectRemoteForceInterface(std::string integrator_)
         bool use_box = true;
         attr = node.attribute("use_box");
         if (attr) use_box = attr.as_bool();
-        
+
         if (use_box) {
           lbRegion reg = lattice->getLocalRegion();
           double px = lattice->px;
@@ -70,6 +92,12 @@ int acRemoteForceInterface::ConnectRemoteForceInterface(std::string integrator_)
             pz + reg.dz - PART_MAR_BOX,
             pz + reg.dz + reg.nz + PART_MAR_BOX);
         }
+
+        attr = node.attribute("omega");
+        if (attr) solver->lattice->RFI_omega = attr.as_bool();
+        attr = node.attribute("torque");
+        if (attr) solver->lattice->RFI_torque = attr.as_bool();
+
         MPI_Barrier(MPMD.local);
         lattice->RFI.Connect(MPMD.work,inter.work);
 

src/Handlers/acSolve.cpp

@@ -7,15 +7,18 @@ int acSolve::Init () {
 		if (GenericAction::ExecuteInternal()) return -1;
 		int stop=0;
 		do {
-			int next_it = Next(solver->iter);
+			int my_next_it = Next(solver->iter);
+			int next_it = my_next_it;
 			for (size_t i=0; i<solver->hands.size(); i++) {
 				int it  = solver->hands[i].Next(solver->iter);
 				if ((it > 0) && (it < next_it)) next_it = it;
 			}
 			solver->steps = next_it;
 			MPI_Bcast(&solver->steps, 1, MPI_INT, 0, MPMD.local);
 			solver->iter += solver->steps;
-			solver->lattice->Iterate(solver->steps, solver->iter_type);
+			int iter_type = solver->iter_type;
+			if (solver->steps == my_next_it) iter_type |= ITER_LASTGLOB;
+			solver->lattice->Iterate(solver->steps, iter_type);
 			CudaDeviceSynchronize();
 			MPI_Barrier(MPMD.local);
 			for (size_t i=0; i<solver->hands.size(); i++) {

src/Handlers/cbPID.cpp

@@ -111,7 +111,7 @@ int cbPID::DoIt () {
 			old_err = err;
 
 		}
-                MPI_Bcast(&control, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
+                MPI_Bcast(&control, 1, MPI_DOUBLE, 0, solver->mpi_comm);
                 double nval;
 /*                if (zone_number < 0) {
                 	sval = solver->lattice->zSet.get(setting, 0, (size_t) 0);

src/Handlers/cbRunR.cpp

@@ -742,6 +742,8 @@ int cbRunR::Init() {
 		output("%s\n",source.c_str());
 		output("--------------------------------\n");
 	}
+	old_iter_type = solver->iter_type;
+	solver->iter_type |= ITER_LASTGLOB;
 	return 0;
 }
 
@@ -788,6 +790,11 @@ int cbRunR::DoIt() {
 	return 0;
 }
 
+int cbRunR::Finish () {
+	solver->iter_type = old_iter_type;
+	return Callback::Finish();
+}
+
 
 #endif // WITH_R
 

src/Handlers/cbRunR.h

@@ -29,9 +29,11 @@ class cbRunR : public  Callback  {
     bool python;
     static int s_tag;
     int tag;
+    int old_iter_type;
 public:
     int Init ();
     int DoIt ();
+    int Finish ();
 };
 
 #endif // WITH_R

src/Lattice.h.Rt

@@ -74,6 +74,7 @@ public:
   real_t px, py, pz; 
   MPIInfo mpi; ///< MPI information
   rfi_t RFI;
+  bool RFI_omega, RFI_torque;
   solidcontainer_t SC;
   size_t particle_data_size_max;
   char snapFileName[STRING_LEN];

src/LatticeBase.cpp.Rt

@@ -14,6 +14,9 @@ LatticeBase::LatticeBase(int zonesettings, int zones, int num_snaps_, const Unit
     data.particle_data_size = 0;
     SC.balls = &RFI;
     RFI.name = "TCLB";
+   	RFI_omega = true;
+	RFI_torque = true;
+
     CudaStreamCreate(&kernelStream);
     CudaStreamCreate(&inStream);
     CudaStreamCreate(&outStream);
@@ -107,29 +110,34 @@ void LatticeBase::CopyInParticles() {
         CudaMalloc(&data.particle_data, RFI.mem_size());
     }
     data.particle_data_size = RFI.size();
-    for (size_t i = 0; i < RFI.size(); i++) {
-        RFI.RawData(i, RFI_DATA_FORCE + 0) = 0;
-        RFI.RawData(i, RFI_DATA_FORCE + 1) = 0;
-        RFI.RawData(i, RFI_DATA_FORCE + 2) = 0;
+    for (int j=0; j<6; j++) for (size_t i=0; i<RFI.size(); i++) RFI.RawData(i, RFI_DATA_FORCE+j) = 0;
+    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;
+    if (RFI.mem_size() > 0) {
+        CudaMemcpyAsync(data.particle_data, RFI.Particles(), RFI.mem_size(), CudaMemcpyHostToDevice, kernelStream);
     }
-    if (RFI.mem_size() > 0) { CudaMemcpyAsync(data.particle_data, RFI.Particles(), RFI.mem_size(), CudaMemcpyHostToDevice, kernelStream); }
     DEBUG_PROF_PUSH("Tree Build");
     SC.Build();
     DEBUG_PROF_POP();
     SC.CopyToGPU(data.solidfinder, kernelStream);
 }
 
 void LatticeBase::CopyOutParticles() {
-    if (RFI.mem_size() > 0) { CudaMemcpyAsync(RFI.Particles(), data.particle_data, RFI.mem_size(), CudaMemcpyDeviceToHost, kernelStream); }
+    if (RFI.mem_size() > 0) {
+        CudaMemcpyAsync(RFI.Particles(), data.particle_data, RFI.mem_size(), CudaMemcpyDeviceToHost, kernelStream);
+    }
     CudaStreamSynchronize(kernelStream);
     DEBUG_PROF_PUSH("Testing particles for NaNs");
     int nans = 0;
-    for (size_t i = 0; i < RFI.size(); i++) {
-        for (int j = 0; j < 3; j++) {
-            if (!isfinite(RFI.RawData(i, RFI_DATA_FORCE + j))) {
-                nans++;
-                RFI.RawData(i, RFI_DATA_FORCE + j) = 0.0;
+    for (int j=0; j<6; j++) {
+        if (RFI_torque || (j<3)) {
+            for (size_t i=0; i<RFI.size(); i++){
+                if (! isfinite(RFI.RawData(i,RFI_DATA_FORCE+j))) {
+                    nans++;
+                    RFI.RawData(i,RFI_DATA_FORCE+j) = 0.0;
+                }
             }
+        } else {
+            for (size_t i=0; i<RFI.size(); i++) RFI.RawData(i,RFI_DATA_FORCE+j) = 0.0;
         }
     }
     if (nans > 0) notice("%d NANs in particle forces (overwritten with 0.0)\n", nans);

src/LatticeBase.hpp

@@ -64,6 +64,7 @@ class LatticeBase {
     solidcontainer_t SC;                      ///<
     size_t particle_data_size_max = 0;        ///<
     rfi_t RFI;                                ///<
+    bool RFI_omega, RFI_torque;
     SyntheticTurbulence ST;                   ///<
     std::string snapFileName;
 

src/Particle.hpp

@@ -30,9 +30,9 @@ struct ParticleI : Particle {
 		angvel.x = particle_data[i*RFI_DATA_SIZE+RFI_DATA_ANGVEL+0];
 		angvel.y = particle_data[i*RFI_DATA_SIZE+RFI_DATA_ANGVEL+1];
 		angvel.z = particle_data[i*RFI_DATA_SIZE+RFI_DATA_ANGVEL+2];
-		diff.x = pos.x - node[0];
-		diff.y = pos.y - node[1];
-		diff.z = pos.z - node[2];
+		diff.x = node[0] - pos.x;
+		diff.y = node[1] - pos.y;
+		diff.z = node[2] - pos.z;
 		dist = sqrt(diff.x*diff.x + diff.y*diff.y + diff.z*diff.z);
 		cvel.x = vel.x + angvel.y*diff.z - angvel.z*diff.y;
 		cvel.y = vel.y + angvel.z*diff.x - angvel.x*diff.z;