Skip to content

Commit 10e3d2c

Browse files
llaniewskillaniewski
authored
Merge pull request CFD-GO#519 from llaniewski/feature/torque
Variable passing to LAMMPS/LIGGGHTS integration
2 parents 68e4df6 + 3c1f54c commit 10e3d2c

File tree

8 files changed

+279
-41
lines changed

8 files changed

+279
-41
lines changed

‎example/data/pipe2x1.stl

3.01 KB
Binary file not shown.

‎example/data/plane1x1.stl

184 Bytes
Binary file not shown.

‎example/particle/3d/cylinder.xml

+85
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,85 @@
1+
<?xml version="1.0"?>
2+
<CLBConfig version="2.0" permissive="true" output="output/">
3+
<Units>
4+
<!-- particle diameter 0.0150m -->
5+
<Param name="X1" value="0.180m" gauge="1x"/>
6+
<Param name="X2" value="1x" gauge="64"/>
7+
<Param name="Y" value="0.090m" gauge="1y"/>
8+
<Param name="nu" value="0.01m2/s" gauge="0.16666"/>
9+
<!-- time step 1/s-->
10+
<Param name="rho" value="1kg/m3" gauge="1"/>
11+
</Units>
12+
<Geometry nx="1x" ny="1y+2" nz="1y+2" py="-0.5y-1" pz="-0.5y-1">
13+
<BGK><Box/></BGK>
14+
<Wall mask="ALL">
15+
<STL file="example/data/pipe2x1.stl" scale="1y" y="0.5y+1" z="0.5y+1" side="out" ray_axis="y"/>
16+
</Wall>
17+
</Geometry>
18+
<Model>
19+
<Param name="nu" value="0.01m2/s"/>
20+
<Param name="aX_mean" value="1Pa/m"/>
21+
<RemoteForceInterface integrator="LAMMPS">
22+
units cgs
23+
boundary p f f
24+
newton off # required off for tangential history
25+
atom_style sphere
26+
atom_modify map array
27+
atom_modify sort 1 0.4
28+
communicate single vel yes
29+
processors * 1 1
30+
31+
neighbor 0.006 bin # ensure skin distance + rp_lrg + rp_sml > dp_lrg
32+
neigh_modify delay 0
33+
34+
# Declare domain
35+
region domain block 0 0.18 -0.045 0.045 -0.045 0.045
36+
create_box 1 domain
37+
38+
# Specify particle groups
39+
group particle_group type 1
40+
41+
# Define region for particle insertion
42+
region pack cylinder x 0 0 0.0375 0 0.12
43+
44+
# Insert particles
45+
fix part_1 particle_group particletemplate/sphere 17891 atom_type 1 density constant 1.0 radius constant 0.0075000000
46+
fix dist particle_group particledistribution/discrete 18143 1 part_1 1
47+
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
48+
run 1
49+
50+
# Specify particle groups
51+
group particle_group type 1
52+
53+
# Define material properties (from which kn kt etc. are calculated for hertz interactions)
54+
soft_particles yes
55+
fix m1 all property/global youngsModulus peratomtype 3000.000000 # defines kn, kt, gamma_n, gamma_t
56+
fix m2 all property/global poissonsRatio peratomtype 0.5 # defines kn, kt, gamma_n, gamma_t
57+
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.8 # defines damping, must be >0.05
58+
fix m4 all property/global coefficientFriction peratomtypepair 1 0.5 # defines friction
59+
60+
fix frac_wall all mesh/surface file example/data/pipe2x1.stl type 1 scale 0.09 move 0 0 0
61+
62+
# Define physics for particle interactions
63+
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
64+
pair_coeff * *
65+
66+
fix granwalls all wall/gran model hertz tangential history mesh n_meshes 1 meshes frac_wall
67+
68+
# Apply integration
69+
fix integr particle_group nve/sphere
70+
71+
# Couple to TCLB
72+
fix tclb all external pf/callback 1 1
73+
74+
dump vtk_dump all atom/vtk 1000 ${output}_part_*.vtu
75+
76+
timestep ${timestep}
77+
78+
run 50000
79+
</RemoteForceInterface>
80+
</Model>
81+
<!-- <VTK Iterations="1000" what="U,Solid"/> -->
82+
<Log Iterations="100"/>
83+
<VTK Iterations="1000"/>
84+
<Solve Iterations="50000"/>
85+
</CLBConfig>

‎example/particle/3d/shear1.xml

+94
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,94 @@
1+
<?xml version="1.0"?>
2+
<CLBConfig version="2.0" permissive="true" output="output/">
3+
<Units>
4+
<Param name="H" value="1m" gauge="32"/>
5+
<Param name="L" value="1x" gauge="1m"/>
6+
<Param name="nu" value="1m2/s" gauge="0.16"/><!-- 1s -->
7+
<Param name="rho" value="1kg/m3" gauge="1"/>
8+
</Units>
9+
<Geometry nx="1x" ny="1m+2" nz="1x" px="-0.5x" py="-1" pz="-0.5x">
10+
<BGK><Box/></BGK>
11+
<MovingWall_N mask="ALL" name="topwall"><Box dy="-1"/></MovingWall_N>
12+
<Wall mask="ALL"><Box ny="1"/></Wall>
13+
</Geometry>
14+
<Model>
15+
<Param name="VelocityX" value="1m/s" zone="topwall"/>
16+
<RemoteForceInterface integrator="LAMMPS" radius="3/m" height="1m/m" length="1x/m" iterations="5s" vtk_it="1s" log_it="100">
17+
units cgs
18+
boundary p f f
19+
newton off # required off for tangential history
20+
atom_style sphere
21+
atom_modify map array
22+
atom_modify sort 1 0.4
23+
communicate single vel yes
24+
processors * 1 1
25+
26+
neighbor 0.006 bin # ensure skin distance + rp_lrg + rp_sml > dp_lrg
27+
neigh_modify delay 0
28+
29+
# Declare domain
30+
variable len2 equal ${length}*0.5
31+
region domain block -${len2} ${len2} 0 ${height} -${len2} ${len2}
32+
create_box 1 domain
33+
34+
# Specify particle groups
35+
group particle_group type 1
36+
37+
# Define region for particle insertion
38+
region pack block -${len2} ${len2} 0 ${height} -${len2} ${len2}
39+
40+
# Insert particles
41+
fix part_1 particle_group particletemplate/sphere 17891 atom_type 1 density constant 1.0 radius constant ${radius}
42+
fix dist particle_group particledistribution/discrete 18143 1 part_1 1
43+
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
44+
run 1
45+
46+
# Specify particle groups
47+
group particle_group type 1
48+
49+
# Define material properties (from which kn kt etc. are calculated for hertz interactions)
50+
soft_particles yes
51+
fix m1 all property/global youngsModulus peratomtype 30000.000000 # defines kn, kt, gamma_n, gamma_t
52+
fix m2 all property/global poissonsRatio peratomtype 0.5 # defines kn, kt, gamma_n, gamma_t
53+
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.8 # defines damping, must be >0.05
54+
fix m4 all property/global coefficientFriction peratomtypepair 1 0.5 # defines friction
55+
56+
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
57+
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}
58+
59+
# Define physics for particle interactions
60+
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
61+
pair_coeff * *
62+
63+
fix granwalls all wall/gran model hertz tangential history mesh n_meshes 2 meshes topwall bottomwall
64+
65+
# Apply integration
66+
fix integr particle_group nve/sphere
67+
68+
# Couple to TCLB
69+
fix tclb all external pf/callback 1 1
70+
71+
variable time equal step*dt
72+
variable tfx equal f_topwall[1]
73+
variable tfy equal f_topwall[2]
74+
variable tfz equal f_topwall[3]
75+
variable bfx equal f_bottomwall[1]
76+
variable bfy equal f_bottomwall[2]
77+
variable bfz equal f_bottomwall[3]
78+
dump forces all mesh/vtk ${vtk_it} ${output}_wall_*.vtk output interpolate id stress stresscomponents
79+
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
80+
81+
dump vtk_dump all atom/vtk ${vtk_it} ${output}_part_*.vtu
82+
83+
84+
timestep ${timestep}
85+
86+
run ${iterations}
87+
</RemoteForceInterface>
88+
</Model>
89+
<!-- <VTK Iterations="1000" what="U,Solid"/> -->
90+
<!-- 10s 1.8m/s 0.01Pa 1/Pa-->
91+
<Log Iterations="100"/>
92+
<VTK Iterations="5s"/>
93+
<Solve Iterations="5s"/>
94+
</CLBConfig>

‎src/Handlers/acRemoteForceInterface.cpp

+30-25
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,6 @@ int acRemoteForceInterface::Init () {
1515

1616
int acRemoteForceInterface::ConnectRemoteForceInterface(std::string integrator_) {
1717
output("Connecting RFI to %s\n",integrator_.c_str());
18-
pugi::xml_attribute attr;
1918
double units[3];
2019
units[0] = solver->units.alt("1m");
2120
units[1] = solver->units.alt("1s");
@@ -25,7 +24,6 @@ int acRemoteForceInterface::ConnectRemoteForceInterface(std::string integrator_)
2524
solver->lattice->RFI.CanCopeWithUnits(false);
2625

2726
solver->lattice->RFI.setVar("output", solver->info.outpath);
28-
2927

3028
std::string element_content;
3129
int node_children = 0;
@@ -44,24 +42,40 @@ int acRemoteForceInterface::ConnectRemoteForceInterface(std::string integrator_)
4442
}
4543
}
4644
if (node_children > 0) solver->lattice->RFI.setVar("content", element_content);
45+
4746
bool stats = false;
4847
std::string stats_prefix = solver->info.outpath;
4948
stats_prefix = stats_prefix + "_RFI";
5049
int stats_iter = 200;
51-
52-
attr = node.attribute("stats");
53-
if (attr) stats = attr.as_bool();
54-
attr = node.attribute("stats_iter");
55-
if (attr) {
56-
stats_iter = solver->units.alt(attr.value());
57-
stats = true;
58-
}
59-
attr = node.attribute("stats_prefix");
60-
if (attr) {
61-
stats_prefix = attr.value();
62-
stats = true;
50+
bool use_box = true;
51+
52+
53+
for (pugi::xml_attribute attr = node.first_attribute(); attr; attr = attr.next_attribute()) {
54+
std::string attr_name = attr.name();
55+
if (attr_name == "integrator") {
56+
// ignore
57+
} else if (attr_name == "stats") {
58+
stats = attr.as_bool();
59+
} else if (attr_name == "stats_iter") {
60+
stats_iter = solver->units.alt(attr.value());
61+
stats = true;
62+
} else if (attr_name == "stats_prefix") {
63+
stats_prefix = attr.value();
64+
stats = true;
65+
} else if (attr_name == "use_box") {
66+
use_box = attr.as_bool();
67+
} else if (attr_name == "omega") {
68+
solver->lattice->RFI_omega = attr.as_bool();
69+
} else if (attr_name == "torque") {
70+
solver->lattice->RFI_torque = attr.as_bool();
71+
} else {
72+
double val = solver->units.alt(attr.value());
73+
char str[STRING_LEN];
74+
sprintf(str, "%.15lg", val);
75+
solver->lattice->RFI.setVar(attr.name(), str);
76+
}
6377
}
64-
78+
6579
if (stats) {
6680
output("Asking for stats on RFI ( %s every %d it)\n", stats_prefix.c_str(), stats_iter);
6781
solver->lattice->RFI.enableStats(stats_prefix.c_str(), stats_iter);
@@ -74,10 +88,6 @@ int acRemoteForceInterface::ConnectRemoteForceInterface(std::string integrator_)
7488
}
7589
integrator = integrator_;
7690

77-
bool use_box = true;
78-
attr = node.attribute("use_box");
79-
if (attr) use_box = attr.as_bool();
80-
8191
if (use_box) {
8292
lbRegion reg = solver->lattice->region;
8393
double px = solver->lattice->px;
@@ -92,15 +102,10 @@ int acRemoteForceInterface::ConnectRemoteForceInterface(std::string integrator_)
92102
pz + reg.dz + reg.nz + PART_MAR_BOX);
93103
}
94104

95-
attr = node.attribute("omega");
96-
if (attr) solver->lattice->RFI_omega = attr.as_bool();
97-
attr = node.attribute("torque");
98-
if (attr) solver->lattice->RFI_torque = attr.as_bool();
99-
100105
MPI_Barrier(MPMD.local);
101106
solver->lattice->RFI.Connect(MPMD.work,inter.work);
102107

103-
return 0;
108+
return 0;
104109
}
105110

106111

‎src/RemoteForceInterface.h

+5
Original file line numberDiff line numberDiff line change
@@ -151,6 +151,11 @@ class RemoteForceInterface {
151151
void setUnits(rfi_real_t meter, rfi_real_t second, rfi_real_t kilogram);
152152
void setVar(const vars_name_t& name, const vars_value_t& value);
153153
bool hasVar(const vars_name_t& name) { return vars.find(name) != vars.end(); };
154+
std::vector<vars_name_t> listVars() {
155+
std::vector<vars_name_t> names;
156+
for (auto it : vars) names.push_back(it.first);
157+
return names;
158+
}
154159
const vars_value_t& getVar(const vars_name_t& name) { return vars[name]; };
155160
inline rfi_real_t& RawData(size_t i, int j) {
156161
if (STORAGE == ArrayOfStructures) {

‎src/lammps.cpp

+29-11
Original file line numberDiff line numberDiff line change
@@ -113,18 +113,36 @@ int main(int argc, char* argv[]) {
113113
MPI_Abort(MPI_COMM_WORLD, 1);
114114
exit(1);
115115
}
116-
FILE* fp = NULL;
117-
fp = fopen(infile.c_str(), "w");
118-
if (fp == NULL) {
119-
printf("ERROR: failed to write to %s\n",infile.c_str());
120-
MPI_Abort(MPI_COMM_WORLD, 1);
121-
exit(1);
122-
}
123-
if (RFI.hasVar("output")) {
124-
fprintf(fp, "variable output string %s\n", RFI.getVar("output").c_str());
116+
if (MPMD.local_rank == 0) {
117+
FILE* fp = NULL;
118+
fp = fopen(infile.c_str(), "w");
119+
if (fp == NULL) {
120+
printf("ERROR: failed to write to %s\n",infile.c_str());
121+
MPI_Abort(MPI_COMM_WORLD, 1);
122+
exit(1);
123+
}
124+
const std::vector<std::string> var_names = RFI.listVars();
125+
for (const auto& v : var_names) {
126+
if (v == "content") continue;
127+
auto& value = RFI.getVar(v);
128+
bool is_numeric = false;
129+
if (v != "output") {
130+
double val;
131+
int ret, len;
132+
ret = sscanf(value.c_str(),"%lf%n", &val, &len);
133+
if ((ret > 0) && (len == value.size())) {
134+
is_numeric = true;
135+
fprintf(fp, "variable %s equal %.13lg\n", v.c_str(), val);
136+
}
137+
}
138+
if (!is_numeric) fprintf(fp, "variable %s string %s\n", v.c_str(), value.c_str());
139+
}
140+
fprintf(fp, "variable timestep equal %.13lg\n", RFI.auto_timestep);
141+
fprintf(fp, "\n");
142+
fprintf(fp, "%s\n", RFI.getVar("content").c_str());
143+
fprintf(fp, "\n");
144+
fclose(fp);
125145
}
126-
fprintf(fp, "%s\n", RFI.getVar("content").c_str());
127-
fclose(fp);
128146
} else {
129147
printf("ERROR: No configuration provided (either xml file or content from force calculator\n");
130148
MPI_Abort(MPI_COMM_WORLD,1);

0 commit comments

Comments
 (0)