CFD-GO · llaniewski · Jun 27, 2024 · Jun 5, 2024 · Jun 5, 2024 · Jun 5, 2024
diff --git a/example/data/pipe2x1.stl b/example/data/pipe2x1.stl
diff --git a/example/data/plane1x1.stl b/example/data/plane1x1.stl
diff --git a/example/particle/3d/cylinder.xml b/example/particle/3d/cylinder.xml
@@ -0,0 +1,85 @@
+<?xml version="1.0"?>
+<CLBConfig version="2.0" permissive="true" output="output/">
+	<Units>
+        <!-- particle diameter 0.0150m -->
+        <Param name="X1" value="0.180m" gauge="1x"/>
+        <Param name="X2" value="1x" gauge="64"/>
+        <Param name="Y" value="0.090m" gauge="1y"/>
+        <Param name="nu" value="0.01m2/s" gauge="0.16666"/>
+        <!-- time step 1/s-->
+		<Param name="rho" value="1kg/m3" gauge="1"/>
+	</Units>
+	<Geometry nx="1x" ny="1y+2" nz="1y+2" py="-0.5y-1" pz="-0.5y-1">
+		<BGK><Box/></BGK>
+        <Wall mask="ALL">
+            <STL file="example/data/pipe2x1.stl" scale="1y" y="0.5y+1" z="0.5y+1" side="out" ray_axis="y"/>
+        </Wall>
+	</Geometry>
+	<Model>
+        <Param name="nu" value="0.01m2/s"/>
+		<Param name="aX_mean" value="1Pa/m"/>
+        <RemoteForceInterface integrator="LAMMPS">
+units 		cgs
+boundary 	p f f
+newton		off # required off for tangential history
+atom_style 	sphere
+atom_modify	map array
+atom_modify	sort 1 0.4
+communicate single vel yes
+processors	* 1 1
+
+neighbor      	0.006 bin # ensure skin distance + rp_lrg + rp_sml > dp_lrg
+neigh_modify  	delay 0
+
+# Declare domain
+region		domain block 0 0.18 -0.045 0.045 -0.045 0.045
+create_box	1 domain
+
+# Specify particle groups
+group  particle_group type 1
+
+# Define region for particle insertion
+region pack cylinder x 0 0 0.0375 0 0.12
+
+# Insert particles
+fix     part_1 particle_group particletemplate/sphere 17891 atom_type 1 density constant 1.0 radius constant 0.0075000000
+fix     dist particle_group particledistribution/discrete 18143 1 part_1 1
+fix 	ins particle_group insert/pack seed 100003 distributiontemplate dist maxattempt 500 insert_every once overlapcheck yes all_in yes  region pack volumefraction_region 0.050000 check_dist_from_subdomain_border no
+run	1
+
+# Specify particle groups
+group  particle_group type 1
+
+# Define material properties (from which kn kt etc. are calculated for hertz interactions)
+soft_particles yes
+fix	m1 all property/global youngsModulus peratomtype 3000.000000 # defines kn, kt, gamma_n, gamma_t
+fix     m2 all property/global poissonsRatio peratomtype 0.5 # defines kn, kt, gamma_n, gamma_t
+fix	m3 all property/global coefficientRestitution peratomtypepair 1 0.8 # defines damping, must be >0.05
+fix	m4 all property/global coefficientFriction peratomtypepair 1 0.5 # defines friction
+
+fix	frac_wall all mesh/surface file example/data/pipe2x1.stl type 1 scale 0.09 move 0 0 0
+
+# Define physics for particle interactions
+pair_style 	gran model hertz tangential history # 'tangential off' sets Ft=0; 'tangential no_history' incorporates damping to Ft, sets kt=0; 'tangential history' incorporate kt and damping into Ft
+pair_coeff	* *
+
+fix granwalls all wall/gran model hertz tangential history mesh n_meshes 1 meshes frac_wall
+
+# Apply integration
+fix	integr particle_group nve/sphere
+
+# Couple to TCLB
+fix 	tclb all external pf/callback 1 1
+
+dump vtk_dump all atom/vtk 1000 ${output}_part_*.vtu
+
+timestep ${timestep}
+
+run	50000
+        </RemoteForceInterface>
+	</Model>
+	<!-- <VTK Iterations="1000" what="U,Solid"/> -->
+	<Log Iterations="100"/>
+    <VTK Iterations="1000"/>
+	<Solve Iterations="50000"/>
+</CLBConfig>
diff --git a/example/particle/3d/shear1.xml b/example/particle/3d/shear1.xml
@@ -0,0 +1,94 @@
+<?xml version="1.0"?>
+<CLBConfig version="2.0" permissive="true" output="output/">
+	<Units>
+		<Param name="H" value="1m" gauge="32"/>
+		<Param name="L" value="1x" gauge="1m"/>
+		<Param name="nu" value="1m2/s" gauge="0.16"/><!-- 1s -->
+		<Param name="rho" value="1kg/m3" gauge="1"/>
+	</Units>
+	<Geometry nx="1x" ny="1m+2" nz="1x" px="-0.5x" py="-1" pz="-0.5x">
+		<BGK><Box/></BGK>
+		<MovingWall_N mask="ALL" name="topwall"><Box dy="-1"/></MovingWall_N>
+		<Wall mask="ALL"><Box ny="1"/></Wall>
+	</Geometry>
+	<Model>
+		<Param name="VelocityX" value="1m/s" zone="topwall"/>
+		<RemoteForceInterface integrator="LAMMPS" radius="3/m" height="1m/m" length="1x/m" iterations="5s" vtk_it="1s" log_it="100">
+			units 		cgs
+			boundary 	p f f
+			newton		off # required off for tangential history
+			atom_style 	sphere
+			atom_modify	map array
+			atom_modify	sort 1 0.4
+			communicate single vel yes
+			processors	* 1 1
+
+			neighbor      	0.006 bin # ensure skin distance + rp_lrg + rp_sml > dp_lrg
+			neigh_modify  	delay 0
+
+			# Declare domain
+			variable len2 equal ${length}*0.5
+			region		domain block -${len2} ${len2} 0 ${height} -${len2} ${len2}
+			create_box	1 domain
+
+			# Specify particle groups
+			group  particle_group type 1
+
+			# Define region for particle insertion
+			region pack block -${len2} ${len2} 0 ${height} -${len2} ${len2}
+
+			# Insert particles
+			fix     part_1 particle_group particletemplate/sphere 17891 atom_type 1 density constant 1.0 radius constant ${radius}
+			fix     dist particle_group particledistribution/discrete 18143 1 part_1 1
+			fix 	ins particle_group insert/pack seed 100003 distributiontemplate dist maxattempt 500 insert_every once overlapcheck yes all_in yes  region pack volumefraction_region 0.30000 check_dist_from_subdomain_border no
+			run	1
+
+			# Specify particle groups
+			group  particle_group type 1
+
+			# Define material properties (from which kn kt etc. are calculated for hertz interactions)
+			soft_particles yes
+			fix	m1 all property/global youngsModulus peratomtype 30000.000000 # defines kn, kt, gamma_n, gamma_t
+			fix     m2 all property/global poissonsRatio peratomtype 0.5 # defines kn, kt, gamma_n, gamma_t
+			fix	m3 all property/global coefficientRestitution peratomtypepair 1 0.8 # defines damping, must be >0.05
+			fix	m4 all property/global coefficientFriction peratomtypepair 1 0.5 # defines friction
+
+			fix	topwall all mesh/surface/stress file example/data/plane1x1.stl type 1 scale ${length} rotate axis 1 0 0 angle 90 move -${len2} ${height} -${len2} surface_vel 1 0 0
+			fix	bottomwall all mesh/surface/stress file example/data/plane1x1.stl type 1 scale ${length} rotate axis 1 0 0 angle 90 move -${len2} 0 -${len2}
+
+			# Define physics for particle interactions
+			pair_style 	gran model hertz tangential history # 'tangential off' sets Ft=0; 'tangential no_history' incorporates damping to Ft, sets kt=0; 'tangential history' incorporate kt and damping into Ft
+			pair_coeff	* *
+
+			fix granwalls all wall/gran model hertz tangential history mesh n_meshes 2 meshes topwall bottomwall
+
+			# Apply integration
+			fix	integr particle_group nve/sphere
+
+			# Couple to TCLB
+			fix 	tclb all external pf/callback 1 1
+
+			variable time equal step*dt
+			variable tfx equal f_topwall[1]
+			variable tfy equal f_topwall[2]
+			variable tfz equal f_topwall[3]
+			variable bfx equal f_bottomwall[1]
+			variable bfy equal f_bottomwall[2]
+			variable bfz equal f_bottomwall[3]
+			dump forces all mesh/vtk ${vtk_it} ${output}_wall_*.vtk output interpolate id stress stresscomponents
+			fix forceslog all print ${log_it} "${time},${tfx},${tfy},${tfz},${bfx},${bfy},${bfz}" file ${output}_forces.csv title "t,tFx,tFy,tFz,bFx,bFy,bFz" screen no
+
+			dump vtk_dump all atom/vtk ${vtk_it} ${output}_part_*.vtu
+
+
+			timestep ${timestep}
+
+			run	${iterations}
+		</RemoteForceInterface>
+	</Model>
+	<!-- <VTK Iterations="1000" what="U,Solid"/> -->
+	<!-- 10s 1.8m/s 0.01Pa 1/Pa-->
+	<Log Iterations="100"/>
+	<VTK Iterations="5s"/>
+	<Solve Iterations="5s"/>
+</CLBConfig>
diff --git a/src/Handlers/acRemoteForceInterface.cpp b/src/Handlers/acRemoteForceInterface.cpp
@@ -15,7 +15,6 @@ int acRemoteForceInterface::Init () {
 
 int acRemoteForceInterface::ConnectRemoteForceInterface(std::string integrator_) {
         output("Connecting RFI to %s\n",integrator_.c_str());
-        pugi::xml_attribute attr;
         double units[3];
         units[0] = solver->units.alt("1m");
         units[1] = solver->units.alt("1s");
@@ -25,7 +24,6 @@ int acRemoteForceInterface::ConnectRemoteForceInterface(std::string integrator_)
         solver->lattice->RFI.CanCopeWithUnits(false);
 
         solver->lattice->RFI.setVar("output", solver->info.outpath);
-
 
         std::string element_content;
         int node_children = 0;
@@ -44,24 +42,40 @@ int acRemoteForceInterface::ConnectRemoteForceInterface(std::string integrator_)
           }
 	      }
         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";
         int stats_iter = 200;
-
-        attr = node.attribute("stats");
-        if (attr) stats = attr.as_bool();
-        attr = node.attribute("stats_iter");
-        if (attr) {
-          stats_iter = solver->units.alt(attr.value());
-          stats = true;
-        }
-        attr = node.attribute("stats_prefix");
-        if (attr) {
-          stats_prefix = attr.value();
-          stats = true;
+        bool use_box = true;
+
+
+        for (pugi::xml_attribute attr = node.first_attribute(); attr; attr = attr.next_attribute()) {
+          std::string attr_name = attr.name();
+          if (attr_name == "integrator") {
+            // ignore
+          } else if (attr_name == "stats") {
+            stats = attr.as_bool();
+          } else if (attr_name == "stats_iter") {
+            stats_iter = solver->units.alt(attr.value());
+            stats = true;
+          } else if (attr_name == "stats_prefix") {
+            stats_prefix = attr.value();
+            stats = true;
+          } else if (attr_name == "use_box") {
+            use_box = attr.as_bool();
+          } else if (attr_name == "omega") {
+            solver->lattice->RFI_omega = attr.as_bool();
+          } else if (attr_name == "torque") {
+            solver->lattice->RFI_torque = attr.as_bool();
+          } else {
+            double val = solver->units.alt(attr.value());
+            char str[STRING_LEN];
+            sprintf(str, "%.15lg", val);
+            solver->lattice->RFI.setVar(attr.name(), str);
+          }
         }
-        
+
         if (stats) {
           output("Asking for stats on RFI ( %s every %d it)\n", stats_prefix.c_str(), stats_iter);
           solver->lattice->RFI.enableStats(stats_prefix.c_str(), stats_iter);
@@ -74,10 +88,6 @@ int acRemoteForceInterface::ConnectRemoteForceInterface(std::string integrator_)
         }
         integrator = 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;
@@ -92,15 +102,10 @@ int acRemoteForceInterface::ConnectRemoteForceInterface(std::string integrator_)
             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);
 
-	return 0;
+        return 0;
 }
 
 

diff --git a/src/RemoteForceInterface.h b/src/RemoteForceInterface.h
@@ -151,6 +151,11 @@ class RemoteForceInterface {
   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(); };
+  std::vector<vars_name_t> listVars() {
+    std::vector<vars_name_t> names;
+    for (auto it : vars) names.push_back(it.first);
+    return names;
+  }
   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) {

diff --git a/src/lammps.cpp b/src/lammps.cpp
@@ -113,18 +113,36 @@ int main(int argc, char* argv[]) {
                 MPI_Abort(MPI_COMM_WORLD, 1);
                 exit(1);
             }
-            FILE* fp = NULL;
-            fp = fopen(infile.c_str(), "w");
-            if (fp == NULL) {
-                printf("ERROR: failed to write to %s\n",infile.c_str());
-                MPI_Abort(MPI_COMM_WORLD, 1);
-                exit(1);
-            }
-            if (RFI.hasVar("output")) {
-                fprintf(fp, "variable output string %s\n", RFI.getVar("output").c_str());
+            if (MPMD.local_rank == 0) {
+                FILE* fp = NULL;
+                fp = fopen(infile.c_str(), "w");
+                if (fp == NULL) {
+                    printf("ERROR: failed to write to %s\n",infile.c_str());
+                    MPI_Abort(MPI_COMM_WORLD, 1);
+                    exit(1);
+                }
+                const std::vector<std::string> var_names = RFI.listVars();
+                for (const auto& v : var_names) {
+                    if (v == "content") continue;
+                    auto& value = RFI.getVar(v);
+                    bool is_numeric = false;
+                    if (v != "output") {
+                        double val;
+                        int ret, len;
+                        ret = sscanf(value.c_str(),"%lf%n", &val, &len);
+                        if ((ret > 0) && (len == value.size())) {
+                            is_numeric = true;
+                            fprintf(fp, "variable %s equal %.13lg\n", v.c_str(), val);
+                        }
+                    }
+                    if (!is_numeric) fprintf(fp, "variable %s string %s\n", v.c_str(), value.c_str());
+                }
+                fprintf(fp, "variable timestep equal %.13lg\n", RFI.auto_timestep);
+                fprintf(fp, "\n");
+                fprintf(fp, "%s\n", RFI.getVar("content").c_str());
+                fprintf(fp, "\n");
+                fclose(fp);
             }
-            fprintf(fp, "%s\n", RFI.getVar("content").c_str());
-            fclose(fp);
         } else {
             printf("ERROR: No configuration provided (either xml file or content from force calculator\n");
             MPI_Abort(MPI_COMM_WORLD,1);