diff --git a/example/particle/3d/cube.xml b/example/particle/3d/cube.xml new file mode 100644 index 000000000..d0739c874 --- /dev/null +++ b/example/particle/3d/cube.xml @@ -0,0 +1,26 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/example/particle/3d/roll_lb.xml b/example/particle/3d/roll_lb.xml new file mode 100644 index 000000000..748b4c06b --- /dev/null +++ b/example/particle/3d/roll_lb.xml @@ -0,0 +1,22 @@ + + + + + + + + + + + + + + + + + + + + + + diff --git a/example/particle/3d/roll_sp.xml b/example/particle/3d/roll_sp.xml new file mode 100644 index 000000000..1e4035306 --- /dev/null +++ b/example/particle/3d/roll_sp.xml @@ -0,0 +1,5 @@ + + + + + diff --git a/src/CartLatticeAccess.hpp.Rt b/src/CartLatticeAccess.hpp.Rt index 7b8911922..af0d74d93 100644 --- a/src/CartLatticeAccess.hpp.Rt +++ b/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 template < class PARENT > template @@ -407,12 +413,12 @@ CudaDeviceFunction real_t SymmetryAccess< PARENT >:: if ( > range_int<0>()) { - if ((this->getNodeType() & NODE_SYM) == NODE_Symmetry_plus) { - return (); + if ((this->getNodeType() & NODE_SYM) == NODE__plus) { + return (); } } else if ( < range_int<0>()) { - if ((this->getNodeType() & NODE_SYM) == NODE_Symmetry_minus) { - return (); + if ((this->getNodeType() & NODE_SYM) == NODE__minus) { + return (); } } @@ -425,9 +431,27 @@ CudaDeviceFunction real_t SymmetryAccess< PARENT >:: template CudaDeviceFunction void SymmetryAccess< PARENT >::pop(N & node) const -{ - parent::pop(node); - +{ + parent::pop(node); + parent::pop(node); + if (this->getNodeType() & (NODE_SYMX | NODE_SYMY | NODE_SYMZ)) { + = load_(range_int< >(),range_int< >(),range_int< >()); + = ; + } else { + parent::pop(node); + } } diff --git a/src/Handlers/GenericOptimizer.cpp b/src/Handlers/GenericOptimizer.cpp index ad78d632a..ccb97ac49 100644 --- a/src/Handlers/GenericOptimizer.cpp +++ b/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; diff --git a/src/Handlers/OptimalControl.cpp b/src/Handlers/OptimalControl.cpp index 92ad09185..039b020bd 100644 --- a/src/Handlers/OptimalControl.cpp +++ b/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;ilattice->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;impi_comm ); return gsum; } diff --git a/src/Handlers/acRemoteForceInterface.cpp b/src/Handlers/acRemoteForceInterface.cpp index 993fa127c..a3b26d9dd 100644 --- a/src/Handlers/acRemoteForceInterface.cpp +++ b/src/Handlers/acRemoteForceInterface.cpp @@ -2,6 +2,8 @@ std::string acRemoteForceInterface::xmlname = "RemoteForceInterface"; #include "../HandlerFactory.h" +#include + 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->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); diff --git a/src/Handlers/acSolve.cpp b/src/Handlers/acSolve.cpp index f66f841c6..030a316c4 100644 --- a/src/Handlers/acSolve.cpp +++ b/src/Handlers/acSolve.cpp @@ -7,7 +7,8 @@ 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; ihands.size(); i++) { int it = solver->hands[i].Next(solver->iter); if ((it > 0) && (it < next_it)) next_it = it; @@ -15,7 +16,9 @@ int acSolve::Init () { 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; ihands.size(); i++) { diff --git a/src/Handlers/cbPID.cpp b/src/Handlers/cbPID.cpp index 6b6fc05bb..013dc1482 100644 --- a/src/Handlers/cbPID.cpp +++ b/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); diff --git a/src/Handlers/cbRunR.cpp b/src/Handlers/cbRunR.cpp index b73eb39e9..42ba5d6b0 100644 --- a/src/Handlers/cbRunR.cpp +++ b/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 diff --git a/src/Handlers/cbRunR.h b/src/Handlers/cbRunR.h index a62fa8b3a..d9330a653 100644 --- a/src/Handlers/cbRunR.h +++ b/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 diff --git a/src/Lattice.h.Rt b/src/Lattice.h.Rt index 6cdde4882..220b997a5 100644 --- a/src/Lattice.h.Rt +++ b/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]; diff --git a/src/LatticeBase.cpp.Rt b/src/LatticeBase.cpp.Rt index ef3e075c3..0460f2904 100644 --- a/src/LatticeBase.cpp.Rt +++ b/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,12 +110,11 @@ 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 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(); @@ -120,16 +122,22 @@ void LatticeBase::CopyInParticles() { } 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 0) notice("%d NANs in particle forces (overwritten with 0.0)\n", nans); diff --git a/src/LatticeBase.hpp b/src/LatticeBase.hpp index 8604790f5..275d3ec70 100644 --- a/src/LatticeBase.hpp +++ b/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; diff --git a/src/Particle.hpp b/src/Particle.hpp index 6285aebe0..91b70841c 100644 --- a/src/Particle.hpp +++ b/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; diff --git a/src/RemoteForceInterface.h b/src/RemoteForceInterface.h index af70e81ce..1dc109122 100644 --- a/src/RemoteForceInterface.h +++ b/src/RemoteForceInterface.h @@ -6,6 +6,7 @@ #include #include #include +#include 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 inline std::vector Exchange(std::vector out); template inline std::basic_string Exchange(std::basic_string 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]; diff --git a/src/RemoteForceInterface.hpp b/src/RemoteForceInterface.hpp index 6fc034766..7eadf302f 100644 --- a/src/RemoteForceInterface.hpp +++ b/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 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); + } } } diff --git a/src/SolidGrid.h b/src/SolidGrid.h index 7b5520ed2..ebaf89a7b 100644 --- a/src/SolidGrid.h +++ b/src/SolidGrid.h @@ -169,6 +169,9 @@ class SolidGrid { inline int TryDepth(size_t grid_size); public: BALLS* balls; + inline SolidGrid() { + depth = 4; + } inline void Build(); inline void InitFinder(finder_t&); inline void CleanFinder(finder_t&); diff --git a/src/SolidGrid.hpp b/src/SolidGrid.hpp new file mode 100644 index 000000000..e69de29bb diff --git a/src/conf.R b/src/conf.R index e71557187..ddc509ce2 100644 --- a/src/conf.R +++ b/src/conf.R @@ -32,7 +32,7 @@ if (! SYMALGEBRA) { library(symAlgebra,quietly=TRUE,warn.conflicts=FALSE) } -if (is.null(Options$autosym)) Options$autosym = FALSE +if (is.null(Options$autosym)) Options$autosym = 0 #source("linemark.R") @@ -369,7 +369,7 @@ if (is.null(Description)) { } -if (Options$autosym) { ## Automatic symmetries +if (Options$autosym > 0) { ## Automatic symmetries symmetries = data.frame(symX=c(-1,1,1),symY=c(1,-1,1),symZ=c(1,1,-1)) for (g in unique(DensityAll$group)) { @@ -398,15 +398,60 @@ if (Options$autosym) { ## Automatic symmetries Fields[sel,s] = Fields$name[sel] } - AddNodeType("SymmetryX_plus", group="SYMX") - AddNodeType("SymmetryX_minus", group="SYMX") - AddNodeType("SymmetryY_plus", group="SYMY") - AddNodeType("SymmetryY_minus", group="SYMY") + rownames(Fields) = Fields$name + + if (Options$autosym == 1) { + autosym_shift = 0 + autosym_name = "Symmetry" + } else if (Options$autosym == 2) { + autosym_shift = 1 + autosym_name = "SymmetryEdge" + } else stop("unknown autosym value") + + directions = lapply(rows(Fields), function(f) expand.grid(dx=f$minx:f$maxx,dy=f$miny:f$maxy,dz=f$minz:f$maxz)) + names(directions) = Fields$name + dir.sort = function(d) { + d = unique(d) + d[order(d[,3],d[,2],d[,1]),] + } + directions = lapply(directions,dir.sort) + tmp = NULL + while (!identical(directions, tmp)) { + tmp = directions + for (f in rows(Fields)) { + d = directions[[f$name]] + for (i in 1:3) { + nfn = f[[names(symmetries)[i]]] + od = directions[[nfn]] + cr = od + nd = d[d[,i] < 0,, drop=FALSE] + nd[,i] = -nd[,i] - autosym_shift + cr = rbind(cr,nd) + nd = d[d[,i] > 0,, drop=FALSE] + nd[,i] = -nd[,i] + autosym_shift + cr = rbind(cr,nd) + directions[[nfn]] = cr + } + } + directions = lapply(directions,dir.sort) + } + Fields$minx = sapply(directions, function(x) min(x$dx)) + Fields$maxx = sapply(directions, function(x) max(x$dx)) + Fields$miny = sapply(directions, function(x) min(x$dy)) + Fields$maxy = sapply(directions, function(x) max(x$dy)) + Fields$minz = sapply(directions, function(x) min(x$dz)) + Fields$maxz = sapply(directions, function(x) max(x$dz)) + + + AddNodeType(paste0(autosym_name, "X_plus"), group="SYMX") + AddNodeType(paste0(autosym_name, "X_minus"), group="SYMX") + AddNodeType(paste0(autosym_name, "Y_plus"), group="SYMY") + AddNodeType(paste0(autosym_name, "Y_minus"), group="SYMY") if (all(range(Fields$minz,Fields$maxz) == c(0,0))) { # we're in 2D } else { - AddNodeType("SymmetryZ_plus", group="SYMZ") - AddNodeType("SymmetryZ_minus", group="SYMZ") + AddNodeType(paste0(autosym_name, "Z_plus"), group="SYMZ") + AddNodeType(paste0(autosym_name, "Z_minus"), group="SYMZ") } } diff --git a/src/lammps.cpp b/src/lammps.cpp index 0b2cdd79d..987ec7966 100644 --- a/src/lammps.cpp +++ b/src/lammps.cpp @@ -1,268 +1,325 @@ -#include "MPMD.hpp" -#include -#include -#include -#include -#include #include +#include #include #include +#include +#include +#include +#include + +#include + #include "Global.h" +#include "MPMD.hpp" #include "RemoteForceInterface.hpp" -#include using namespace LAMMPS_NS; - -typedef rfi::RemoteForceInterface< rfi::ForceIntegrator, rfi::RotParticle, rfi::ArrayOfStructures, real_t > RFI_t; +typedef rfi::RemoteForceInterface RFI_t; struct Info { - MPMDHelper *MPMD; - Memory *memory; - LAMMPS *lmp; - RFI_t * RFI; - std::vector wsize; - std::vector windex; + MPMDHelper* MPMD; + LAMMPS* lmp; + RFI_t* RFI; + std::vector wsize; + std::vector windex; + bool atom_force; + bool first_print; }; - -void tclb_callback(void *, bigint, int, int *, double **, double **); +void tclb_callback(void*, bigint, int, int*, double**, double**); int match_pattern(char* str, char* pattern) { - bool fit = true; - while (fit) { - if (pattern[0] == '\0') break; - if (pattern[0] == ' ') { - if (str[0] == ' ') str++; - else if (str[0] == '\t') str++; - else pattern++; - } else if (pattern[0] == '*') { - if (str[0] == ' ') pattern++; - else if (str[0] == '\t') pattern++; - else str++; - } else if (pattern[0] == str[0]) { - pattern++; str++; - } else fit = false; - } - return fit; + bool fit = true; + while (fit) { + if (pattern[0] == '\0') break; + if (pattern[0] == ' ') { + if (str[0] == ' ') str++; + else if (str[0] == '\t') + str++; + else + pattern++; + } else if (pattern[0] == '*') { + if (str[0] == ' ') pattern++; + else if (str[0] == '\t') + pattern++; + else + str++; + } else if (pattern[0] == str[0]) { + pattern++; + str++; + } else + fit = false; + } + return fit; } +int main(int argc, char* argv[]) { + int ret; + Info info; + MPMDHelper MPMD; + MPI_Init(&argc, &argv); + MPMD.Init(MPI_COMM_WORLD, "LAMMPS"); + DEBUG_SETRANK(MPMD.local_rank); + InitPrint(DEBUG_LEVEL, 6, 8); + MPMD.Identify(); + RFI_t RFI; + RFI.name = "LAMMPS"; -int main(int argc, char *argv[]) -{ - int ret; - Info info; - MPMDHelper MPMD; - MPI_Init(&argc, &argv); - MPMD.Init(MPI_COMM_WORLD, "LAMMPS"); - DEBUG_SETRANK(MPMD.local_rank); - InitPrint(DEBUG_LEVEL, 6, 8); - MPMD.Identify(); - RFI_t RFI; - RFI.name = "LAMMPS"; - - if (argc < 2) { - printf("Syntax: lammps in.lammps [args]\n"); - MPI_Abort(MPI_COMM_WORLD,1); - exit(1); - } - std::vector lammps_args; - char * infile = NULL; - bool logset = false; - for (int i=0; i lammps_args; + std::string infile = ""; + bool logset = false; + for (int i=0; i %s",line); fflush(stdout); - } + FILE* fp = NULL; + if (MPMD.local_rank == 0) { + printf("LAMMPS: running input file: %s\n", infile.c_str()); + printf("LAMMPS: arguments:"); + for (int i = 1; i < lammps_args.size(); i++) printf(" %s", lammps_args[i]); + printf("\n"); + fp = fopen(infile.c_str(), "r"); + if (fp == NULL) { + printf("ERROR: Could not open LAMMPS input script %s\n", infile.c_str()); + MPI_Abort(MPI_COMM_WORLD, 1); + exit(1); + } + } - lammps_command(lmp,line); + LAMMPS* lmp = NULL; + lmp = new LAMMPS(lammps_args.size(), &lammps_args[0], MPMD.local); - if (match_pattern(line, " fix tclb * external")) { - printf("LAMMPS: Added fix tclb external! (%s) Adding callback.\n", line); - MPMDIntercomm inter = MPMD["TCLB"]; - if (!inter) { - fprintf(stderr,"Didn't find TCLB in MPMD\n"); - return -1; - } - ret = RFI.Connect(MPMD.work,inter.work); - if (ret) return ret; - assert(RFI.Connected()); + int n; + char line[1024]; + while (1) { + if (MPMD.local_rank == 0) { + if (fgets(line, 1024, fp) == NULL) n = 0; + else + n = strlen(line) + 1; + if (n == 0) fclose(fp); + } + MPI_Bcast(&n, 1, MPI_INT, 0, MPMD.local); + if (n == 0) break; + MPI_Bcast(line, n, MPI_CHAR, 0, MPMD.local); + if (MPMD.local_rank == 0) { + fprintf(stdout, "LAMMPS> %s", line); + fflush(stdout); + } - int ifix = lmp->modify->find_fix("tclb"); - printf("ifix: %d\n",ifix); - FixExternal *fix = (FixExternal *) lmp->modify->fix[ifix]; - printf("fix: %p\n",fix); + lammps_command(lmp, line); - info.MPMD = &MPMD; - info.memory = NULL; - info.lmp = lmp; - info.RFI = &RFI; - info.wsize.resize(RFI.Workers()); - info.windex.resize(RFI.Workers()); + if (match_pattern(line, " fix tclb * external")) { + printf("LAMMPS: Added fix tclb external! (%s) Adding callback.\n", line); + assert(RFI.Connected()); - fix->set_callback(tclb_callback,&info); - } - } + int ifix = lmp->modify->find_fix("tclb"); + FixExternal* fix = (FixExternal*)lmp->modify->fix[ifix]; - lammps_close(lmp); - if (RFI.Connected()) { - RFI.Close(); - } - MPI_Finalize(); - return 0; -} + info.MPMD = &MPMD; + info.lmp = lmp; + info.RFI = &RFI; + info.atom_force = false; + info.first_print = true; + info.wsize.resize(RFI.Workers()); + info.windex.resize(RFI.Workers()); + fix->set_callback(tclb_callback, &info); + } + } -void tclb_callback(void *ptr, bigint ntimestep, int nlocal, int *id, double **x_, double **f) -{ - Info *info = (Info *) ptr; - if (! info->RFI->Active()) return; + lammps_close(lmp); + if (RFI.Connected()) { RFI.Close(); } + MPI_Finalize(); + return 0; +} - double ** v = info->lmp->atom->v; - double ** x = info->lmp->atom->x; - double * r = info->lmp->atom->radius; +void tclb_callback(void* ptr, bigint ntimestep, int nlocal, int* id, double** x_, double** f_) { + Info* info = (Info*)ptr; + if (!info->RFI->Active()) return; - for (int phase = 0; phase < 3; phase++) { - if (phase == 0) { - for (int i = 0; i < info->RFI->Workers(); i++) info->wsize[i] = 0; + double** v = info->lmp->atom->v; + double** x = info->lmp->atom->x; + double** omega = info->lmp->atom->omega; + double** f = NULL; + double** torque = NULL; + for (size_t k=0; katom_force) { + f = info->lmp->atom->f; } else { - for (int i = 0; i < info->RFI->Workers(); i++) info->windex[i] = 0; + f = f_; } + torque = info->lmp->atom->torque; + double* r = info->lmp->atom->radius; - for (size_t k = 0; k < nlocal; k++) { - if (phase == 2) { - f[k][0] = 0; - f[k][1] = 0; - f[k][2] = 0; - } - int minper[3], maxper[3], d[3]; - size_t offset = 0; - for (int worker = 0; worker < info->RFI->Workers(); worker++) { - for (int j=0; j<3; j++) { - double prd = info->lmp->domain->prd[j]; - double lower = 0; - double upper = info->lmp->domain->prd[j]; - if (info->RFI->WorkerBox(worker).declared) { - lower = info->RFI->WorkerBox(worker).lower[j]; - upper = info->RFI->WorkerBox(worker).upper[j]; - } - if (info->lmp->domain->periodicity[j]) { - maxper[j] = floor((upper - x[k][j] + r[k]) / prd); - minper[j] = ceil((lower - x[k][j] - r[k]) / prd); - } else { - if ((x[k][j] + r[k] >= lower) && (x[k][j] - r[k] <= upper)) { - minper[j] = 0; maxper[j] = 0; - } else { - minper[j] = 0; maxper[j] = -1; // no balls - } - } -// printf("particle %ld dimenstion %d in %d worker interval [%lg %lg] and periodicity %lg: %lg copied %d:%d\n", k, j, worker, lower, upper, prd, x[k][j], minper[j], maxper[j]); - } - - int copies = (maxper[0]-minper[0]+1)*(maxper[1]-minper[1]+1)*(maxper[2]-minper[2]+1); - // if (copies > 1) printf("particle %ld is copied %d times (%d %d)x(%d %d)x(%d %d)\n", k, copies, minper[0], maxper[0], minper[1], maxper[1], minper[2], maxper[2]); - for (d[0]=minper[0]; d[0]<=maxper[0]; d[0]++) { - for (d[1]=minper[1]; d[1]<=maxper[1]; d[1]++) { - for (d[2]=minper[2]; d[2]<=maxper[2]; d[2]++) { - double px[3]; - for (int j=0; j<3; j++) px[j] = x[k][j] + d[j] * info->lmp->domain->prd[j]; - if (phase == 0) { - info->wsize[worker]++; - } else { - size_t i = offset + info->windex[worker]; - if (phase == 1) { - //printf("particle %ld sent %d at index %ld\n", k, worker, i); - info->RFI->setData(i, RFI_DATA_R, r[k]); - info->RFI->setData(i, RFI_DATA_POS+0, px[0]); - info->RFI->setData(i, RFI_DATA_POS+1, px[1]); - info->RFI->setData(i, RFI_DATA_POS+2, px[2]); - info->RFI->setData(i, RFI_DATA_VEL+0, v[k][0]); - info->RFI->setData(i, RFI_DATA_VEL+1, v[k][1]); - info->RFI->setData(i, RFI_DATA_VEL+2, v[k][2]); - if (info->RFI->Rot()) { - info->RFI->setData(i, RFI_DATA_ANGVEL+0, 0.0); - info->RFI->setData(i, RFI_DATA_ANGVEL+1, 0.0); - info->RFI->setData(i, RFI_DATA_ANGVEL+2, 0.0); - } - } else { - f[k][0] += info->RFI->getData(i, RFI_DATA_FORCE+0); - f[k][1] += info->RFI->getData(i, RFI_DATA_FORCE+1); - f[k][2] += info->RFI->getData(i, RFI_DATA_FORCE+2); - } - info->windex[worker]++; - } - } - } - } - offset += info->wsize[worker]; - } + if (info->first_print) { + printf("LAMMPS: RFI.Rot:%s atom_force:%s f:%s torque:%s\n", + info->RFI->Rot() ? "true" : "false", + info->atom_force ? "true" : "false", + (f != NULL) ? "true" : "false", + (torque != NULL) ? "true" : "false" + ); + info->first_print = false; } - if (phase == 0) { - for (int worker = 0; worker < info->RFI->Workers(); worker++) info->RFI->Size(worker) = info->wsize[worker]; - //printf("sizes:"); for (int worker = 0; worker < info->RFI->Workers(); worker++) printf(" %ld", (size_t) info->wsize[worker]); printf("\n"); - info->RFI->SendSizes(); - info->RFI->Alloc(); - } else if (phase == 1) { - //printf("indexes:"); for (int worker = 0; worker < info->RFI->Workers(); worker++) printf(" %ld", (size_t) info->windex[worker]); printf("\n"); - info->RFI->SendParticles(); - info->RFI->SendForces(); - } else { + + for (int phase=0; phase<3; phase++) { + if (phase == 0) { + for (int i = 0; i < info->RFI->Workers(); i++) info->wsize[i] = 0; + } else { + for (int i = 0; i < info->RFI->Workers(); i++) info->windex[i] = 0; + } + + for (size_t k=0; kRFI->Workers(); worker++) { + for (int j=0; j<3; j++) { + double prd = info->lmp->domain->prd[j]; + double lower = 0; + double upper = info->lmp->domain->prd[j]; + if (info->RFI->WorkerBox(worker).declared) { + lower = info->RFI->WorkerBox(worker).lower[j]; + upper = info->RFI->WorkerBox(worker).upper[j]; + } + if (info->lmp->domain->periodicity[j]) { + maxper[j] = floor((upper - x[k][j] + r[k])/prd); + minper[j] = ceil((lower - x[k][j] - r[k])/prd); + } else { + if ((x[k][j] + r[k] >= lower) && (x[k][j] - r[k] <= upper)) { + minper[j] = 0; + maxper[j] = 0; + } else { + minper[j] = 0; + maxper[j] = -1; // no balls + } + } + // printf("particle %ld dimenstion %d in %d worker interval [%lg %lg] and periodicity %lg: %lg copied %d:%d\n", k, j, worker, lower, upper, prd, x[k][j], minper[j], maxper[j]); + } + + int copies = (maxper[0] - minper[0] + 1)*(maxper[1] - minper[1] + 1)*(maxper[2] - minper[2] + 1); + // if (copies > 1) printf("particle %ld is copied %d times (%d %d)x(%d %d)x(%d %d)\n", k, copies, minper[0], maxper[0], minper[1], maxper[1], minper[2], maxper[2]); + for (d[0] = minper[0]; d[0] <= maxper[0]; d[0]++) { + for (d[1] = minper[1]; d[1] <= maxper[1]; d[1]++) { + for (d[2] = minper[2]; d[2] <= maxper[2]; d[2]++) { + double px[3]; + for (int j=0; j<3; j++) px[j] = x[k][j] + d[j]*info->lmp->domain->prd[j]; + if (phase == 0) { + info->wsize[worker]++; + } else { + size_t i = offset + info->windex[worker]; + if (phase == 1) { + //printf("particle %ld sent %d at index %ld\n", k, worker, i); + info->RFI->setData(i, RFI_DATA_R, r[k]); + info->RFI->setData(i, RFI_DATA_POS + 0, px[0]); + info->RFI->setData(i, RFI_DATA_POS + 1, px[1]); + info->RFI->setData(i, RFI_DATA_POS + 2, px[2]); + info->RFI->setData(i, RFI_DATA_VEL + 0, v[k][0]); + info->RFI->setData(i, RFI_DATA_VEL + 1, v[k][1]); + info->RFI->setData(i, RFI_DATA_VEL + 2, v[k][2]); + if (info->RFI->Rot()) { + if (omega) { + info->RFI->setData(i, RFI_DATA_ANGVEL + 0, omega[k][0]); + info->RFI->setData(i, RFI_DATA_ANGVEL + 1, omega[k][1]); + info->RFI->setData(i, RFI_DATA_ANGVEL + 2, omega[k][2]); + } else { + info->RFI->setData(i, RFI_DATA_ANGVEL + 0, 0.0); + info->RFI->setData(i, RFI_DATA_ANGVEL + 1, 0.0); + info->RFI->setData(i, RFI_DATA_ANGVEL + 2, 0.0); + } + } + } else { + if (f) { + f[k][0] += info->RFI->getData(i, RFI_DATA_FORCE + 0); + f[k][1] += info->RFI->getData(i, RFI_DATA_FORCE + 1); + f[k][2] += info->RFI->getData(i, RFI_DATA_FORCE + 2); + } + if (torque) { + torque[k][0] += info->RFI->getData(i, RFI_DATA_MOMENT + 0); + torque[k][1] += info->RFI->getData(i, RFI_DATA_MOMENT + 1); + torque[k][2] += info->RFI->getData(i, RFI_DATA_MOMENT + 2); + } + } + info->windex[worker]++; + } + } + } + } + offset += info->wsize[worker]; + } + } + if (phase == 0) { + for (int worker = 0; worker < info->RFI->Workers(); worker++) info->RFI->Size(worker) = info->wsize[worker]; + //printf("sizes:"); for (int worker = 0; worker < info->RFI->Workers(); worker++) printf(" %ld", (size_t) info->wsize[worker]); printf("\n"); + info->RFI->SendSizes(); + info->RFI->Alloc(); + } else if (phase == 1) { + //printf("indexes:"); for (int worker = 0; worker < info->RFI->Workers(); worker++) printf(" %ld", (size_t) info->windex[worker]); printf("\n"); + info->RFI->SendParticles(); + info->RFI->SendForces(); + } else { + } } - } } - diff --git a/src/makefile.main.Rt b/src/makefile.main.Rt index c40b42a80..877c6e997 100644 --- a/src/makefile.main.Rt +++ b/src/makefile.main.Rt @@ -148,7 +148,7 @@ for (d in destinations) { paste0(" ", names(opts), " = ", - ifelse(opts==0,"FALSE","TRUE"), + opts, ifelse(seq_along(opts) != length(opts),",","") ) ) diff --git a/src/models.R b/src/models.R index b4db6571f..98cf69650 100644 --- a/src/models.R +++ b/src/models.R @@ -47,11 +47,23 @@ get.models = function() { opts_terms = terms(opts) opts = attr(opts_terms,"factors") opts = data.frame(t(opts)) - rownames(opts) = paste(name,gsub(":","_",rownames(opts)),sep="_") + opts[] = opts > 0 if (attr(opts_terms, "intercept") == 1) { - opts[name,]=0 + opts = rbind(opts, FALSE) } -# opts = apply(opts, 2, function(x) x > 0) + if ("autosym" %in% names(opts)) { + opts$autosym = ifelse(opts$autosym, 1, 0) + x = opts[opts$autosym > 0,,drop=FALSE] + x$autosym = 2 + opts = rbind(opts, x) + } + rownames(opts) = sapply(seq_len(nrow(opts)), function(i) { + x = opts[i,] + w = names(opts) + w = paste0(w, ifelse(x>1,x,"")) + w = c(name, w[x>0]) + paste(w,collapse="_") + }) } else { opts = data.frame(row.names=name) } diff --git a/src/simplepart.cpp b/src/simplepart.cpp index c9a579f51..e5f5b5820 100644 --- a/src/simplepart.cpp +++ b/src/simplepart.cpp @@ -14,6 +14,8 @@ struct Particle { double v[3]; double f[3]; double favg[3]; + double omega[3]; + double torque[3]; size_t n; bool logging; Particle() { @@ -23,6 +25,8 @@ struct Particle { v[i] = 0; f[i] = 0; favg[i] = 0; + omega[i] = 0; + torque[i] = 0; } m = 0; r = 0; @@ -30,6 +34,33 @@ struct Particle { } }; +struct attr_name_t { + std::string vector; + std::string non_vector; + int d; + attr_name_t(const std::string& name) { + bool vec = true; + d = -1; + auto w = name.back(); + if (w == 'x') { d = 0; } + else if (w == 'y') { d = 1; } + else if (w == 'z') { d = 2; } + else { vec = false; } + if (vec) { + vector = name; + vector.pop_back(); + non_vector = "="; + } else { + non_vector = name; + vector = "="; + } + } + attr_name_t(const char* name) : attr_name_t(std::string(name)) {}; + bool operator==(const std::string& name) { + return non_vector == name; + } +}; + typedef std::vector Particles; int main(int argc, char *argv[]) { @@ -44,7 +75,7 @@ int main(int argc, char *argv[]) { return -1; } MPMD.Identify(); - rfi::RemoteForceInterface RFI; + rfi::RemoteForceInterface RFI; RFI.name = "SIMPLEPART"; MPMDIntercomm inter = MPMD["TCLB"]; @@ -57,19 +88,27 @@ int main(int argc, char *argv[]) { wsize.resize(RFI.Workers()); windex.resize(RFI.Workers()); Particles particles; - double dt = 0; + double dt = RFI.auto_timestep; bool logging = false; - std::string logging_filename; + std::string logging_filename = ""; + if (RFI.hasVar("output")) logging_filename = RFI.getVar("output") + "_SP_Log.csv"; int logging_iter = 1; FILE* logging_f = NULL; bool avg = false; - double periodicity[3]; + bool log_position = true; + bool log_velocity = true; + bool log_force = true; + bool log_omega = false; + bool log_torque = false; + + double periodicity[3], periodic_origin[3]; bool periodic[3]; for (int i = 0; i < 3; i++) { periodic[i] = false; periodicity[i] = 0.0; + periodic_origin[i] = 0.0; } double acc_vec[3]; @@ -77,14 +116,33 @@ int main(int argc, char *argv[]) { for (int i = 0; i < 3; i++) acc_vec[i] = 0.0; acc_freq = 0.0; - if (argc != 2) { + if (argc < 1 || argc > 2) { printf("Syntax: simplepart config.xml\n"); + printf(" You can omit config.xml if configuration is provided by the force calculator (eg. TCLB xml)\n"); MPI_Abort(MPI_COMM_WORLD,1); exit(1); } - char * filename = argv[1]; + + char * filename = NULL; + if (argc > 1) { + filename = argv[1]; + } pugi::xml_document config; - pugi::xml_parse_result result = config.load_file(filename, pugi::parse_default | pugi::parse_comments); + pugi::xml_parse_result result; + if (filename != NULL) { + if (RFI.hasVar("content")) { + WARNING("Ignoring content (configuration) sent by calculator"); + } + result = config.load_file(filename, pugi::parse_default | pugi::parse_comments); + } else { + if (RFI.hasVar("content")) { + result = config.load_string(RFI.getVar("content").c_str(), pugi::parse_default | pugi::parse_comments); + } else { + printf("No configuration provided (either xml file or content from force calculator\n"); + MPI_Abort(MPI_COMM_WORLD,1); + exit(1); + } + } if (!result) { ERROR("Error while parsing %s: %s\n", filename, result.description()); return -1; @@ -118,19 +176,13 @@ int main(int argc, char *argv[]) { if (node_name == "Particle") { Particle p; for (pugi::xml_attribute attr = node.first_attribute(); attr; attr = attr.next_attribute()) { - std::string attr_name = attr.name(); - if (attr_name == "x") { - p.x[0] = attr.as_double(); - } else if (attr_name == "y") { - p.x[1] = attr.as_double(); - } else if (attr_name == "z") { - p.x[2] = attr.as_double(); - } else if (attr_name == "vx") { - p.v[0] = attr.as_double(); - } else if (attr_name == "vy") { - p.v[1] = attr.as_double(); - } else if (attr_name == "vz") { - p.v[2] = attr.as_double(); + attr_name_t attr_name = attr.name(); + if (attr_name.vector == "") { + p.x[attr_name.d] = attr.as_double(); + } else if (attr_name.vector == "v") { + p.v[attr_name.d] = attr.as_double(); + } else if (attr_name.vector == "omega") { + p.omega[attr_name.d] = attr.as_double(); } else if (attr_name == "r") { p.r = attr.as_double(); } else if (attr_name == "m") { @@ -150,16 +202,12 @@ int main(int argc, char *argv[]) { particles.push_back(p); } else if (node_name == "Periodic") { for (pugi::xml_attribute attr = node.first_attribute(); attr; attr = attr.next_attribute()) { - std::string attr_name = attr.name(); - if (attr_name == "x") { - periodic[0] = true; - periodicity[0] = attr.as_double(); - } else if (attr_name == "y") { - periodic[1] = true; - periodicity[1] = attr.as_double(); - } else if (attr_name == "z") { - periodic[2] = true; - periodicity[2] = attr.as_double(); + attr_name_t attr_name = attr.name(); + if (attr_name.vector == "") { + periodic[attr_name.d] = true; + periodicity[attr_name.d] = attr.as_double(); + } else if (attr_name.vector == "p") { + periodic_origin[attr_name.d] = attr.as_double(); } else { ERROR("Unknown atribute '%s' in '%s'", attr.name(), node.name()); return -1; @@ -172,8 +220,8 @@ int main(int argc, char *argv[]) { } for (pugi::xml_attribute attr = node.first_attribute(); attr; attr = attr.next_attribute()) { std::string attr_name = attr.name(); + logging = true; if (attr_name == "name") { - logging = true; logging_filename = attr.value(); } else if (attr_name == "Iterations") { logging_iter = attr.as_int(); @@ -184,13 +232,17 @@ int main(int argc, char *argv[]) { } else if (attr_name == "average") { avg = attr.as_bool(); if (avg) notice("SIMPLEPART: Particle force averaging is ON"); + } else if (attr_name == "rotation") { + log_omega = attr.as_bool(); + log_torque = log_omega; + if (log_omega) notice("SIMPLEPART: Particle omega and torque is logged"); } else { ERROR("Unknown atribute '%s' in '%s'", attr.name(), node.name()); return -1; } } - if (!logging) { - ERROR("Name not set in '%s' element", node.name()); + if (logging && logging_filename == "") { + ERROR("Loggin file name not set in '%s' element", node.name()); return -1; } } else { @@ -207,7 +259,11 @@ int main(int argc, char *argv[]) { fprintf(logging_f, "Iteration,Time"); for (Particles::iterator p = particles.begin(); p != particles.end(); p++) if (p->logging) { size_t n = p->n; - fprintf(logging_f, ",p%ld_x,p%ld_y,p%ld_z,p%ld_vx,p%ld_vy,p%ld_vz,p%ld_fx,p%ld_fy,p%ld_fz",n,n,n,n,n,n,n,n,n); + if (log_position) fprintf(logging_f, ",p%2$ld_%1$sx,p%2$ld_%1$sy,p%2$ld_%1$sz","",n); + if (log_velocity) fprintf(logging_f, ",p%2$ld_%1$sx,p%2$ld_%1$sy,p%2$ld_%1$sz","v",n); + if (log_force) fprintf(logging_f, ",p%2$ld_%1$sx,p%2$ld_%1$sy,p%2$ld_%1$sz","f",n); + if (log_omega) fprintf(logging_f, ",p%2$ld_%1$sx,p%2$ld_%1$sy,p%2$ld_%1$sz","o",n); + if (log_torque) fprintf(logging_f, ",p%2$ld_%1$sx,p%2$ld_%1$sy,p%2$ld_%1$sz","t",n); } fprintf(logging_f, "\n"); } @@ -230,6 +286,9 @@ int main(int argc, char *argv[]) { p->f[0] = 0; p->f[1] = 0; p->f[2] = 0; + p->torque[0] = 0; + p->torque[1] = 0; + p->torque[2] = 0; } int minper[3], maxper[3], d[3]; size_t offset = 0; @@ -276,14 +335,19 @@ int main(int argc, char *argv[]) { RFI.setData(i, RFI_DATA_VEL + 1, p->v[1]); RFI.setData(i, RFI_DATA_VEL + 2, p->v[2]); if (RFI.Rot()) { - RFI.setData(i, RFI_DATA_ANGVEL + 0, 0.0); - RFI.setData(i, RFI_DATA_ANGVEL + 1, 0.0); - RFI.setData(i, RFI_DATA_ANGVEL + 2, 0.0); + RFI.setData(i, RFI_DATA_ANGVEL + 0, p->omega[0]); + RFI.setData(i, RFI_DATA_ANGVEL + 1, p->omega[1]); + RFI.setData(i, RFI_DATA_ANGVEL + 2, p->omega[2]); } } else { p->f[0] += RFI.getData(i, RFI_DATA_FORCE + 0); p->f[1] += RFI.getData(i, RFI_DATA_FORCE + 1); p->f[2] += RFI.getData(i, RFI_DATA_FORCE + 2); + if (RFI.Rot()) { + p->torque[0] += RFI.getData(i, RFI_DATA_MOMENT + 0); + p->torque[1] += RFI.getData(i, RFI_DATA_MOMENT + 1); + p->torque[2] += RFI.getData(i, RFI_DATA_MOMENT + 2); + } } windex[worker]++; } @@ -307,14 +371,18 @@ int main(int argc, char *argv[]) { if (logging && (iter % logging_iter == 0)) { fprintf(logging_f, "%d,%.15lg", iter, dt*iter); for (Particles::iterator p = particles.begin(); p != particles.end(); p++) if (p->logging) { - for (int i=0; i<3; i++) fprintf(logging_f, ",%.15lg", p->x[i]); - for (int i=0; i<3; i++) fprintf(logging_f, ",%.15lg", p->v[i]); - if (avg) { - for (int i=0; i<3; i++) fprintf(logging_f, ",%.15lg", p->favg[i]/logging_iter); - } else { - for (int i=0; i<3; i++) fprintf(logging_f, ",%.15lg", p->f[i]); + if (log_position) for (int i=0; i<3; i++) fprintf(logging_f, ",%.15lg", p->x[i]); + if (log_velocity) for (int i=0; i<3; i++) fprintf(logging_f, ",%.15lg", p->v[i]); + if (log_force) { + if (avg) { + for (int i=0; i<3; i++) fprintf(logging_f, ",%.15lg", p->favg[i]/logging_iter); + for (int i=0; i<3; i++) p->favg[i] = 0; + } else { + for (int i=0; i<3; i++) fprintf(logging_f, ",%.15lg", p->f[i]); + } } - for (int i=0; i<3; i++) p->favg[i] = 0; + if (log_omega) for (int i=0; i<3; i++) fprintf(logging_f, ",%.15lg", p->omega[i]); + if (log_torque) for (int i=0; i<3; i++) fprintf(logging_f, ",%.15lg", p->torque[i]); } fprintf(logging_f, "\n"); }