Merge pull request #512 from llaniewski/feature/torque

Adding torque exchange to lammps/liggghts integration

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/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->info.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->info.outpath;
         stats_prefix = stats_prefix + "_RFI";
@@ -55,7 +77,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 = solver->lattice->region;
           double px = solver->lattice->px;
@@ -69,6 +91,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);
         solver->lattice->RFI.Connect(MPMD.work,inter.work);
 

src/Lattice.cu.Rt

@@ -146,6 +146,8 @@ Lattice::Lattice(lbRegion _region, MPIInfo mpi_, int ns): zSet(ZONESETTINGS, ZON
 	container->particle_data_size = 0;
 	SC.balls = &RFI;
 	RFI.name = "TCLB";
+	RFI_omega = true;
+	RFI_torque = true;
 }
 
 /// Initialization of MPI buffors
@@ -398,11 +400,9 @@ void Lattice::CopyInParticles() {
 		CudaMalloc(&container->particle_data, RFI.mem_size());
 	}
 	container->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(container->particle_data, RFI.Particles(), RFI.mem_size(), CudaMemcpyHostToDevice, kernelStream);
 	}
@@ -419,12 +419,16 @@ void Lattice::CopyOutParticles() {
     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/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/Particle.hpp

@@ -28,9 +28,9 @@ struct ParticleI : Particle {
 		angvel.x = constContainer.particle_data[i*RFI_DATA_SIZE+RFI_DATA_ANGVEL+0];
 		angvel.y = constContainer.particle_data[i*RFI_DATA_SIZE+RFI_DATA_ANGVEL+1];
 		angvel.z = constContainer.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;

src/RemoteForceInterface.h

@@ -6,6 +6,7 @@
 #include <stdio.h>
 #include <vector>
 #include <string>
+#include <map>
 
 namespace rfi {
 
@@ -94,6 +95,10 @@ class RemoteForceInterface {
   std::vector< rfi_real_t > unit;
   bool non_trivial_units;
   bool can_cope_with_units;
+  typedef std::string vars_name_t;
+  typedef std::string vars_value_t;
+  typedef std::map< vars_name_t, vars_value_t > vars_t;
+  vars_t vars;
   void ISendSizes();
   void WSendSizes();
   void ISendParticles();
@@ -114,6 +119,7 @@ class RemoteForceInterface {
 public:
   int particle_size;
   std::string name;
+  rfi_real_t auto_timestep;
   RemoteForceInterface();
   ~RemoteForceInterface();
   void MakeTypes(bool,bool);  
@@ -143,6 +149,9 @@ class RemoteForceInterface {
   template <class T> inline std::vector<T> Exchange(std::vector<T> out);
   template <class T> inline std::basic_string<T> Exchange(std::basic_string<T> out);
   void setUnits(rfi_real_t meter, rfi_real_t second, rfi_real_t kilogram);
+  void setVar(const vars_name_t& name, const vars_value_t& value);
+  bool hasVar(const vars_name_t& name) { return vars.find(name) != vars.end(); };
+  const vars_value_t& getVar(const vars_name_t& name) { return vars[name]; };
   inline rfi_real_t& RawData(size_t i, int j) {
     if (STORAGE == ArrayOfStructures) {
       return tab[i*particle_size + j];

src/RemoteForceInterface.hpp

@@ -18,7 +18,7 @@
 
 namespace rfi {
 
-const int version = 0x000104;
+const int version = 0x000105;
 
 #define safe_MPI_Type_free(datatype) { if ((*datatype) != NULL) MPI_Type_free(datatype); }
 
@@ -117,6 +117,7 @@ RemoteForceInterface < TYPE, ROT, STORAGE, rfi_real_t, tab_allocator >::RemoteFo
    base_units[0] = 1.0;
    base_units[1] = 1.0;
    base_units[2] = 1.0;
+   auto_timestep = 1.0;
    non_trivial_units = false;
    can_cope_with_units = true;
 }
@@ -211,6 +212,15 @@ inline void RemoteForceInterface < TYPE, ROT, STORAGE, rfi_real_t, tab_allocator
  non_trivial_units = true;
 }
 
+template < rfi_type_t TYPE, rfi_rot_t ROT, rfi_storage_t STORAGE, typename rfi_real_t, typename tab_allocator >
+inline void RemoteForceInterface < TYPE, ROT, STORAGE, rfi_real_t, tab_allocator >::setVar(const vars_name_t& name, const vars_value_t& value) {
+ if (Connected()) {
+   ERROR("Vars can be set only before connection is established\n");
+   exit(-1);
+ }
+ vars[name] = value;
+}
+
 template < rfi_type_t TYPE, rfi_rot_t ROT, rfi_storage_t STORAGE, typename rfi_real_t, typename tab_allocator >
 inline void RemoteForceInterface < TYPE, ROT, STORAGE, rfi_real_t, tab_allocator >::CanCopeWithUnits(bool ccwu_) {
    if (Connected()) {
@@ -362,6 +372,43 @@ int RemoteForceInterface < TYPE, ROT, STORAGE, rfi_real_t, tab_allocator >::Nego
         unit[RFI_DATA_FORCE+i] = kilogram*meter/(second*second);
         unit[RFI_DATA_MOMENT+i] = kilogram*meter*meter/(second*second);
       }
+      auto_timestep = 1.0/second;
+    }
+
+    typedef std::vector<std::string::value_type> vars_pack_t;
+    vars_pack_t my_vars, other_vars;
+    for (vars_t::const_iterator it = vars.begin(); it != vars.end(); it++) {
+      std::string v;
+      v = it->first;
+      for (std::string::const_iterator it2 = v.begin(); it2 != v.end(); it2++) my_vars.push_back(*it2);
+      my_vars.push_back(0);
+      v = it->second;
+      for (std::string::const_iterator it2 = v.begin(); it2 != v.end(); it2++) my_vars.push_back(*it2);
+      my_vars.push_back(0);
+    }
+    other_vars = Exchange(my_vars);
+    bool is_name = true;
+    std::string buf;
+    std::string name_buf;
+    std::string value_buf;
+    for (vars_pack_t::const_iterator it = other_vars.begin(); it != other_vars.end(); it++) {
+      if (*it == 0) {
+        if (is_name) {
+          name_buf = buf;
+          is_name = false;
+        } else {
+          value_buf = buf;
+          is_name = true;
+          debug1("RFI: %s: Received: %s = %s\n",name.c_str(), name_buf.c_str(), value_buf.c_str());
+          if (vars.find(name_buf) != vars.end()) {
+            debug1("RFI: %s: Variable overwritten.", name.c_str());
+          }
+          vars[name_buf] = value_buf;
+        }
+        buf.clear();
+      } else {
+        buf.push_back(*it);
+      }
     }
   }
 

src/SolidGrid.h

@@ -167,7 +167,7 @@ class SolidGrid {
 public:
     BALLS* balls;
     inline SolidGrid() {
-        depth = 0;
+        depth = 4;
     }
     void Build();
     void InitFinder(finder_t&);

src/SolidGrid.hpp

@@ -58,7 +58,6 @@ void SolidGrid<BALLS>::Build () {
         }
         size_t grid_size = 1;
         for (int k=0; k<3; k++) grid_size = grid_size * (maxs[k]-mins[k]+1);
-        depth = 4;
         while (TryDepth(grid_size)) {
             depth = depth*2;
             output("Too many particles per bin in SolidGrid. Increasing bin size to %d\n", depth);