Skip to content

Commit

Permalink
Merge pull request #519 from llaniewski/feature/torque
Browse files Browse the repository at this point in the history 
Variable passing to LAMMPS/LIGGGHTS integration
  • Loading branch information
@llaniewski
llaniewski authored Jun 27, 2024
2 parents 68e4df6 + 3c1f54c commit 10e3d2c
Show file tree
Hide file tree
Showing 8 changed files with 279 additions and 41 deletions.
Binary file added example/data/pipe2x1.stl
View file
Binary file not shown.
Binary file added example/data/plane1x1.stl
View file
Binary file not shown.
85 changes: 85 additions & 0 deletions example/particle/3d/cylinder.xml
View file
Original file line number Diff line number Diff line change
@@ -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>
94 changes: 94 additions & 0 deletions example/particle/3d/shear1.xml
View file
Original file line number Diff line number Diff line change
@@ -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>
55 changes: 30 additions & 25 deletions src/Handlers/acRemoteForceInterface.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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");
Expand All @@ -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;
Expand All @@ -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);
Expand All @@ -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;
Expand All @@ -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;
}


Expand Down
5 changes: 5 additions & 0 deletions src/RemoteForceInterface.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down
40 changes: 29 additions & 11 deletions src/lammps.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
Loading

0 comments on commit 10e3d2c

Please sign in to comment.