diff --git a/src/RHEO/atom_vec_rheo.cpp b/src/RHEO/atom_vec_rheo.cpp index 0cbebff008..31f1aa2f7e 100644 --- a/src/RHEO/atom_vec_rheo.cpp +++ b/src/RHEO/atom_vec_rheo.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories diff --git a/src/RHEO/atom_vec_rheo_thermal.cpp b/src/RHEO/atom_vec_rheo_thermal.cpp index 426c059570..09ce6665ff 100644 --- a/src/RHEO/atom_vec_rheo_thermal.cpp +++ b/src/RHEO/atom_vec_rheo_thermal.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -47,8 +46,10 @@ AtomVecRHEOThermal::AtomVecRHEOThermal(LAMMPS *lmp) : AtomVec(lmp) // order of fields in a string does not matter // except: fields_data_atom & fields_data_vel must match data file - fields_grow = {"rheo_status", "rho", "drho", "temperature", "esph", "heatflow", "conductivity", "pressure", "viscosity"}; - fields_copy = {"rheo_status", "rho", "drho", "temperature", "esph", "heatflow", "conductivity", "pressure", "viscosity"}; + fields_grow = {"rheo_status", "rho", "drho", "temperature", "esph", + "heatflow", "conductivity", "pressure", "viscosity"}; + fields_copy = {"rheo_status", "rho", "drho", "temperature", "esph", + "heatflow", "conductivity", "pressure", "viscosity"}; fields_comm = {"rheo_status", "rho", "esph"}; fields_comm_vel = {"rheo_status", "rho", "esph"}; fields_reverse = {"drho", "heatflow"}; @@ -56,7 +57,8 @@ AtomVecRHEOThermal::AtomVecRHEOThermal(LAMMPS *lmp) : AtomVec(lmp) fields_border_vel = {"rheo_status", "rho", "esph"}; fields_exchange = {"rheo_status", "rho", "esph"}; fields_restart = {"rheo_status", "rho", "esph"}; - fields_create = {"rheo_status", "rho", "drho", "temperature", "esph", "heatflow", "conductivity", "pressure", "viscosity"}; + fields_create = {"rheo_status", "rho", "drho", "temperature", "esph", + "heatflow", "conductivity", "pressure", "viscosity"}; fields_data_atom = {"id", "type", "rheo_status", "rho", "esph", "x"}; fields_data_vel = {"id", "v"}; diff --git a/src/RHEO/bond_rheo_shell.cpp b/src/RHEO/bond_rheo_shell.cpp index dbe65e2b2f..7ef05bbc38 100644 --- a/src/RHEO/bond_rheo_shell.cpp +++ b/src/RHEO/bond_rheo_shell.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -37,7 +36,7 @@ #include #include -#define EPSILON 1e-10 +static constexpr double EPSILON = 1e-10; using namespace LAMMPS_NS; using namespace RHEO_NS; @@ -45,8 +44,8 @@ using namespace RHEO_NS; /* ---------------------------------------------------------------------- */ BondRHEOShell::BondRHEOShell(LAMMPS *_lmp) : - BondBPM(_lmp), k(nullptr), ecrit(nullptr), gamma(nullptr), dbond(nullptr), nbond(nullptr), - id_fix(nullptr), compute_surface(nullptr) + BondBPM(_lmp), k(nullptr), ecrit(nullptr), gamma(nullptr), dbond(nullptr), nbond(nullptr), + id_fix(nullptr), compute_surface(nullptr) { partial_flag = 1; comm_reverse = 1; @@ -168,8 +167,7 @@ void BondRHEOShell::compute(int eflag, int vflag) store_data(); } - if (hybrid_flag) - fix_bond_history->compress_history(); + if (hybrid_flag) fix_bond_history->compress_history(); int i1, i2, itmp, n, type; double delx, dely, delz, delvx, delvy, delvz; @@ -191,7 +189,7 @@ void BondRHEOShell::compute(int eflag, int vflag) double **bondstore = fix_bond_history->bondstore; - if (atom->nmax > nmax_store){ + if (atom->nmax > nmax_store) { nmax_store = atom->nmax; memory->destroy(dbond); memory->create(dbond, nmax_store, "rheo/shell:dbond"); @@ -226,8 +224,7 @@ void BondRHEOShell::compute(int eflag, int vflag) r = sqrt(rsq); // If bond hasn't been set - zero data - if (t < EPSILON || std::isnan(t)) - t = store_bond(n, i1, i2); + if (t < EPSILON || std::isnan(t)) t = store_bond(n, i1, i2); delx = x[i1][0] - x[i2][0]; dely = x[i1][1] - x[i2][1]; @@ -297,15 +294,14 @@ void BondRHEOShell::compute(int eflag, int vflag) // Communicate changes in nbond if (newton_bond) comm->reverse_comm(this); - for(int i = 0; i < nlocal; i++) { + for (int i = 0; i < nlocal; i++) { nbond[i] += dbond[i]; // If it has bonds, no shifting if (nbond[i] != 0) status[i] |= STATUS_NO_SHIFT; } - if (hybrid_flag) - fix_bond_history->uncompress_history(); + if (hybrid_flag) fix_bond_history->uncompress_history(); } /* ---------------------------------------------------------------------- */ @@ -368,12 +364,13 @@ void BondRHEOShell::init_style() if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use bond rheo/shell"); class FixRHEO *fix_rheo = dynamic_cast(fixes[0]); - if (!fix_rheo->surface_flag) error->all(FLERR, - "Bond rheo/shell requires surface calculation in fix rheo"); + if (!fix_rheo->surface_flag) + error->all(FLERR, "Bond rheo/shell requires surface calculation in fix rheo"); compute_surface = fix_rheo->compute_surface; fixes = modify->get_fix_by_style("^rheo/oxidation$"); - if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo/oxidation to use bond rheo/shell"); + if (fixes.size() == 0) + error->all(FLERR, "Need to define fix rheo/oxidation to use bond rheo/shell"); class FixRHEOOxidation *fix_rheo_oxidation = dynamic_cast(fixes[0]); rsurf = fix_rheo_oxidation->rsurf; @@ -402,7 +399,6 @@ void BondRHEOShell::settings(int narg, char **arg) error->all(FLERR, "Illegal bond rheo/shell command, must specify positive formation time"); } - /* ---------------------------------------------------------------------- used to check bond communiction cutoff - not perfect, estimates based on local-local only ------------------------------------------------------------------------- */ @@ -464,13 +460,10 @@ void BondRHEOShell::write_restart_settings(FILE *fp) void BondRHEOShell::read_restart_settings(FILE *fp) { - if (comm->me == 0) { - utils::sfread(FLERR, &tform, sizeof(double), 1, fp, nullptr, error); - } + if (comm->me == 0) { utils::sfread(FLERR, &tform, sizeof(double), 1, fp, nullptr, error); } MPI_Bcast(&tform, 1, MPI_DOUBLE, 0, world); } - /* ---------------------------------------------------------------------- */ int BondRHEOShell::pack_reverse_comm(int n, int first, double *buf) @@ -479,9 +472,7 @@ int BondRHEOShell::pack_reverse_comm(int n, int first, double *buf) m = 0; last = first + n; - for (i = first; i < last; i++) { - buf[m++] = dbond[i]; - } + for (i = first; i < last; i++) { buf[m++] = dbond[i]; } return m; } @@ -561,8 +552,8 @@ void BondRHEOShell::process_ineligibility(int i, int j) n = num_bond[i]; bond_type[i][m] = bond_type[i][n - 1]; bond_atom[i][m] = bond_atom[i][n - 1]; - for (auto &ihistory: histories) { - auto fix_bond_history2 = dynamic_cast (ihistory); + for (auto &ihistory : histories) { + auto fix_bond_history2 = dynamic_cast(ihistory); fix_bond_history2->shift_history(i, m, n - 1); fix_bond_history2->delete_history(i, n - 1); } @@ -579,8 +570,8 @@ void BondRHEOShell::process_ineligibility(int i, int j) n = num_bond[j]; bond_type[j][m] = bond_type[j][n - 1]; bond_atom[j][m] = bond_atom[j][n - 1]; - for (auto &ihistory: histories) { - auto fix_bond_history2 = dynamic_cast (ihistory); + for (auto &ihistory : histories) { + auto fix_bond_history2 = dynamic_cast(ihistory); fix_bond_history2->shift_history(j, m, n - 1); fix_bond_history2->delete_history(j, n - 1); } diff --git a/src/RHEO/compute_rheo_grad.cpp b/src/RHEO/compute_rheo_grad.cpp index fdddea18b5..63c77e3f3b 100644 --- a/src/RHEO/compute_rheo_grad.cpp +++ b/src/RHEO/compute_rheo_grad.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -21,15 +20,15 @@ #include "atom.h" #include "comm.h" -#include "compute_rheo_kernel.h" #include "compute_rheo_interface.h" +#include "compute_rheo_kernel.h" #include "domain.h" #include "error.h" #include "fix_rheo.h" #include "force.h" #include "memory.h" -#include "neighbor.h" #include "neigh_list.h" +#include "neighbor.h" #include "update.h" #include @@ -37,24 +36,29 @@ using namespace LAMMPS_NS; using namespace RHEO_NS; -enum{COMMGRAD, COMMFIELD}; +enum { COMMGRAD, COMMFIELD }; /* ---------------------------------------------------------------------- */ ComputeRHEOGrad::ComputeRHEOGrad(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), gradv(nullptr), gradr(nullptr), grade(nullptr), gradn(nullptr), - fix_rheo(nullptr), rho0(nullptr), compute_kernel(nullptr), compute_interface(nullptr), - list(nullptr) + Compute(lmp, narg, arg), gradv(nullptr), gradr(nullptr), grade(nullptr), gradn(nullptr), + fix_rheo(nullptr), rho0(nullptr), compute_kernel(nullptr), compute_interface(nullptr), + list(nullptr) { - if (narg < 4) error->all(FLERR,"Illegal compute rheo/grad command"); + if (narg < 4) error->all(FLERR, "Illegal compute rheo/grad command"); velocity_flag = energy_flag = rho_flag = eta_flag = 0; for (int iarg = 3; iarg < narg; iarg++) { - if (strcmp(arg[iarg], "velocity") == 0) velocity_flag = 1; - else if (strcmp(arg[iarg], "rho") == 0) rho_flag = 1; - else if (strcmp(arg[iarg], "energy") == 0) energy_flag = 1; - else if (strcmp(arg[iarg], "viscosity") == 0) eta_flag = 1; - else error->all(FLERR, "Illegal compute rheo/grad command, {}", arg[iarg]); + if (strcmp(arg[iarg], "velocity") == 0) + velocity_flag = 1; + else if (strcmp(arg[iarg], "rho") == 0) + rho_flag = 1; + else if (strcmp(arg[iarg], "energy") == 0) + energy_flag = 1; + else if (strcmp(arg[iarg], "viscosity") == 0) + eta_flag = 1; + else + error->all(FLERR, "Illegal compute rheo/grad command, {}", arg[iarg]); } ncomm_grad = 0; @@ -89,7 +93,6 @@ ComputeRHEOGrad::ComputeRHEOGrad(LAMMPS *lmp, int narg, char **arg) : nmax_store = 0; grow_arrays(atom->nmax); - } /* ---------------------------------------------------------------------- */ @@ -161,20 +164,16 @@ void ComputeRHEOGrad::compute_peratom() for (i = 0; i < nmax_store; i++) { if (velocity_flag) { - for (k = 0; k < dim * dim; k++) - gradv[i][k] = 0.0; + for (k = 0; k < dim * dim; k++) gradv[i][k] = 0.0; } if (rho_flag) { - for (k = 0; k < dim; k++) - gradr[i][k] = 0.0; + for (k = 0; k < dim; k++) gradr[i][k] = 0.0; } if (energy_flag) { - for (k = 0; k < dim; k++) - grade[i][k] = 0.0; + for (k = 0; k < dim; k++) grade[i][k] = 0.0; } if (eta_flag) { - for (k = 0; k < dim; k++) - gradn[i][k] = 0.0; + for (k = 0; k < dim; k++) gradn[i][k] = 0.0; } } @@ -216,10 +215,10 @@ void ComputeRHEOGrad::compute_peratom() if (interface_flag) { if (fluidi && (!fluidj)) { compute_interface->correct_v(vj, vi, j, i); - rhoj = compute_interface->correct_rho(j, i); + rhoj = compute_interface->correct_rho(j); } else if ((!fluidi) && fluidj) { compute_interface->correct_v(vi, vj, i, j); - rhoi = compute_interface->correct_rho(i, j); + rhoi = compute_interface->correct_rho(i); } else if ((!fluidi) && (!fluidj)) { rhoi = rho0[itype]; rhoj = rho0[jtype]; @@ -248,34 +247,34 @@ void ComputeRHEOGrad::compute_peratom() for (a = 0; a < dim; a++) { for (b = 0; b < dim; b++) { - if (velocity_flag) // uxx uxy uxz uyx uyy uyz uzx uzy uzz + if (velocity_flag) // uxx uxy uxz uyx uyy uyz uzx uzy uzz gradv[i][a * dim + b] -= vij[a] * Volj * dWij[b]; } - if (rho_flag) // P,x P,y P,z + if (rho_flag) // P,x P,y P,z gradr[i][a] -= drho * Volj * dWij[a]; - if (energy_flag) // e,x e,y e,z + if (energy_flag) // e,x e,y e,z grade[i][a] -= de * Volj * dWij[a]; - if (eta_flag) // n,x n,y n,z + if (eta_flag) // n,x n,y n,z gradn[i][a] -= deta * Volj * dWij[a]; } if (newton || j < nlocal) { for (a = 0; a < dim; a++) { for (b = 0; b < dim; b++) { - if (velocity_flag) // uxx uxy uxz uyx uyy uyz uzx uzy uzz + if (velocity_flag) // uxx uxy uxz uyx uyy uyz uzx uzy uzz gradv[j][a * dim + b] += vij[a] * Voli * dWji[b]; } - if (rho_flag) // P,x P,y P,z + if (rho_flag) // P,x P,y P,z gradr[j][a] += drho * Voli * dWji[a]; - if (energy_flag) // e,x e,y e,z + if (energy_flag) // e,x e,y e,z grade[j][a] += de * Voli * dWji[a]; - if (eta_flag) // n,x n,y n,z + if (eta_flag) // n,x n,y n,z gradn[j][a] += deta * Voli * dWji[a]; } } @@ -295,7 +294,6 @@ void ComputeRHEOGrad::forward_gradients() comm->forward_comm(this); } - /* ---------------------------------------------------------------------- */ void ComputeRHEOGrad::forward_fields() @@ -307,10 +305,9 @@ void ComputeRHEOGrad::forward_fields() /* ---------------------------------------------------------------------- */ -int ComputeRHEOGrad::pack_forward_comm(int n, int *list, double *buf, - int pbc_flag, int *pbc) +int ComputeRHEOGrad::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) { - int i,j,k,m; + int i, j, k, m; int *mask = atom->mask; double *rho = atom->rho; double *energy = atom->esph; @@ -333,38 +330,30 @@ int ComputeRHEOGrad::pack_forward_comm(int n, int *list, double *buf, if (comm_stage == COMMGRAD) { if (velocity_flag) - for (k = 0; k < dim * dim; k++) - buf[m++] = gradv[j][k]; + for (k = 0; k < dim * dim; k++) buf[m++] = gradv[j][k]; if (rho_flag) - for (k = 0; k < dim; k++) - buf[m++] = gradr[j][k]; + for (k = 0; k < dim; k++) buf[m++] = gradr[j][k]; if (energy_flag) - for (k = 0; k < dim; k++) - buf[m++] = grade[j][k]; + for (k = 0; k < dim; k++) buf[m++] = grade[j][k]; if (eta_flag) - for (k = 0; k < dim; k++) - buf[m++] = gradn[j][k]; + for (k = 0; k < dim; k++) buf[m++] = gradn[j][k]; } else if (comm_stage == COMMFIELD) { if (velocity_flag) { if (remap_v_flag && pbc_flag && (mask[j] & deform_groupbit)) { - for (k = 0; k < dim; k++) - buf[m++] = v[j][k] + dv[k]; + for (k = 0; k < dim; k++) buf[m++] = v[j][k] + dv[k]; } else { - for (k = 0; k < dim; k++) - buf[m++] = v[j][k]; + for (k = 0; k < dim; k++) buf[m++] = v[j][k]; } } - if (rho_flag) - buf[m++] = rho[j]; + if (rho_flag) buf[m++] = rho[j]; - if (energy_flag) - buf[m++] = energy[j]; + if (energy_flag) buf[m++] = energy[j]; } } return m; @@ -376,7 +365,8 @@ void ComputeRHEOGrad::unpack_forward_comm(int n, int first, double *buf) { int i, k, m, last; double *rho = atom->rho; - double *energy = atom->esph; double **v = atom->v; + double *energy = atom->esph; + double **v = atom->v; int dim = domain->dimension; m = 0; @@ -384,31 +374,24 @@ void ComputeRHEOGrad::unpack_forward_comm(int n, int first, double *buf) for (i = first; i < last; i++) { if (comm_stage == COMMGRAD) { if (velocity_flag) - for (k = 0; k < dim * dim; k++) - gradv[i][k] = buf[m++]; + for (k = 0; k < dim * dim; k++) gradv[i][k] = buf[m++]; if (rho_flag) - for (k = 0; k < dim; k++) - gradr[i][k] = buf[m++]; + for (k = 0; k < dim; k++) gradr[i][k] = buf[m++]; if (energy_flag) - for (k = 0; k < dim; k++) - grade[i][k] = buf[m++]; + for (k = 0; k < dim; k++) grade[i][k] = buf[m++]; if (eta_flag) - for (k = 0; k < dim; k++) - gradn[i][k] = buf[m++]; + for (k = 0; k < dim; k++) gradn[i][k] = buf[m++]; } else if (comm_stage == COMMFIELD) { if (velocity_flag) - for (k = 0; k < dim; k++) - v[i][k] = buf[m++]; + for (k = 0; k < dim; k++) v[i][k] = buf[m++]; - if (rho_flag) - rho[i] = buf[m++]; + if (rho_flag) rho[i] = buf[m++]; - if (energy_flag) - energy[i] = buf[m++]; + if (energy_flag) energy[i] = buf[m++]; } } } @@ -417,27 +400,23 @@ void ComputeRHEOGrad::unpack_forward_comm(int n, int first, double *buf) int ComputeRHEOGrad::pack_reverse_comm(int n, int first, double *buf) { - int i,k,m,last; + int i, k, m, last; int dim = domain->dimension; m = 0; last = first + n; for (i = first; i < last; i++) { if (velocity_flag) - for (k = 0; k < dim * dim; k++) - buf[m++] = gradv[i][k]; + for (k = 0; k < dim * dim; k++) buf[m++] = gradv[i][k]; if (rho_flag) - for (k = 0; k < dim; k++) - buf[m++] = gradr[i][k]; + for (k = 0; k < dim; k++) buf[m++] = gradr[i][k]; if (energy_flag) - for (k = 0; k < dim; k++) - buf[m++] = grade[i][k]; + for (k = 0; k < dim; k++) buf[m++] = grade[i][k]; if (eta_flag) - for (k = 0; k < dim; k++) - buf[m++] = gradn[i][k]; + for (k = 0; k < dim; k++) buf[m++] = gradn[i][k]; } return m; } @@ -446,27 +425,23 @@ int ComputeRHEOGrad::pack_reverse_comm(int n, int first, double *buf) void ComputeRHEOGrad::unpack_reverse_comm(int n, int *list, double *buf) { - int i,k,j,m; + int i, k, j, m; int dim = domain->dimension; m = 0; for (i = 0; i < n; i++) { j = list[i]; if (velocity_flag) - for (k = 0; k < dim * dim; k++) - gradv[j][k] += buf[m++]; + for (k = 0; k < dim * dim; k++) gradv[j][k] += buf[m++]; if (rho_flag) - for (k = 0; k < dim; k++) - gradr[j][k] += buf[m++]; + for (k = 0; k < dim; k++) gradr[j][k] += buf[m++]; if (energy_flag) - for (k = 0; k < dim; k++) - grade[j][k] += buf[m++]; + for (k = 0; k < dim; k++) grade[j][k] += buf[m++]; if (eta_flag) - for (k = 0; k < dim; k++) - gradn[j][k] += buf[m++]; + for (k = 0; k < dim; k++) gradn[j][k] += buf[m++]; } } @@ -475,17 +450,13 @@ void ComputeRHEOGrad::unpack_reverse_comm(int n, int *list, double *buf) void ComputeRHEOGrad::grow_arrays(int nmax) { int dim = domain->dimension; - if (velocity_flag) - memory->grow(gradv, nmax, dim * dim, "rheo:grad_v"); + if (velocity_flag) memory->grow(gradv, nmax, dim * dim, "rheo:grad_v"); - if (rho_flag) - memory->grow(gradr, nmax, dim, "rheo:grad_rho"); + if (rho_flag) memory->grow(gradr, nmax, dim, "rheo:grad_rho"); - if (energy_flag) - memory->grow(grade, nmax, dim, "rheo:grad_energy"); + if (energy_flag) memory->grow(grade, nmax, dim, "rheo:grad_energy"); - if (eta_flag) - memory->grow(gradn, nmax, dim, "rheo:grad_eta"); + if (eta_flag) memory->grow(gradn, nmax, dim, "rheo:grad_eta"); nmax_store = nmax; } @@ -496,17 +467,13 @@ double ComputeRHEOGrad::memory_usage() double bytes = 0.0; int dim = domain->dimension; - if (velocity_flag) - bytes = (size_t) nmax_store * dim * dim * sizeof(double); + if (velocity_flag) bytes = (size_t) nmax_store * dim * dim * sizeof(double); - if (rho_flag) - bytes = (size_t) nmax_store * dim * sizeof(double); + if (rho_flag) bytes = (size_t) nmax_store * dim * sizeof(double); - if (energy_flag) - bytes = (size_t) nmax_store * dim * sizeof(double); + if (energy_flag) bytes = (size_t) nmax_store * dim * sizeof(double); - if (eta_flag) - bytes = (size_t) nmax_store * dim * sizeof(double); + if (eta_flag) bytes = (size_t) nmax_store * dim * sizeof(double); return bytes; } diff --git a/src/RHEO/compute_rheo_interface.cpp b/src/RHEO/compute_rheo_interface.cpp index 2d2a368926..eec5659051 100644 --- a/src/RHEO/compute_rheo_interface.cpp +++ b/src/RHEO/compute_rheo_interface.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -21,18 +20,18 @@ #include "atom.h" #include "comm.h" -#include "domain.h" #include "compute_rheo_kernel.h" +#include "domain.h" #include "error.h" -#include "force.h" #include "fix_rheo.h" #include "fix_rheo_pressure.h" +#include "force.h" #include "math_extra.h" #include "memory.h" #include "modify.h" -#include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" +#include "neighbor.h" #include @@ -45,12 +44,12 @@ static constexpr double EPSILON = 1e-1; /* ---------------------------------------------------------------------- */ ComputeRHEOInterface::ComputeRHEOInterface(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), chi(nullptr), fp_store(nullptr), fix_rheo(nullptr), - rho0(nullptr), norm(nullptr), normwf(nullptr), id_fix_pa(nullptr), list(nullptr), - compute_kernel(nullptr), fix_pressure(nullptr) + Compute(lmp, narg, arg), chi(nullptr), fp_store(nullptr), fix_rheo(nullptr), rho0(nullptr), + norm(nullptr), normwf(nullptr), id_fix_pa(nullptr), list(nullptr), compute_kernel(nullptr), + fix_pressure(nullptr) { - if (narg != 3) error->all(FLERR,"Illegal compute rheo/interface command"); + if (narg != 3) error->all(FLERR, "Illegal compute rheo/interface command"); comm_forward = 3; comm_reverse = 4; @@ -166,20 +165,18 @@ void ComputeRHEOInterface::compute_peratom() if (rsq < cutsq) { jtype = type[j]; fluidj = !(status[j] & PHASECHECK); - w = compute_kernel->calc_w_quintic(i, j, dx[0], dx[1], dx[2], sqrt(rsq)); + w = compute_kernel->calc_w_quintic(sqrt(rsq)); norm[i] += w; status_match = 0; - if ((fluidi && fluidj) || ((!fluidi) && (!fluidj))) - status_match = 1; + if ((fluidi && fluidj) || ((!fluidi) && (!fluidj))) status_match = 1; if (status_match) { chi[i] += w; } else { if (!fluidi) { dot = 0; - for (a = 0; a < 3; a++) - dot += (-fp_store[j][a] + fp_store[i][a]) * dx[a]; + for (a = 0; a < 3; a++) dot += (-fp_store[j][a] + fp_store[i][a]) * dx[a]; rho[i] += w * (fix_pressure->calc_pressure(rho[j], jtype) - rho[j] * dot); normwf[i] += w; @@ -193,8 +190,7 @@ void ComputeRHEOInterface::compute_peratom() } else { if (!fluidj) { dot = 0; - for (a = 0; a < 3; a++) - dot += (-fp_store[i][a] + fp_store[j][a]) * dx[a]; + for (a = 0; a < 3; a++) dot += (-fp_store[i][a] + fp_store[j][a]) * dx[a]; rho[j] += w * (fix_pressure->calc_pressure(rho[i], itype) + rho[i] * dot); normwf[j] += w; @@ -228,7 +224,8 @@ void ComputeRHEOInterface::compute_peratom() /* ---------------------------------------------------------------------- */ -int ComputeRHEOInterface::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/) +int ComputeRHEOInterface::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, + int * /*pbc*/) { int m = 0; double *rho = atom->rho; @@ -293,7 +290,7 @@ void ComputeRHEOInterface::unpack_reverse_comm(int n, int *list, double *buf) int j = list[i]; norm[j] += buf[m++]; chi[j] += buf[m++]; - if (status[j] & PHASECHECK){ + if (status[j] & PHASECHECK) { normwf[j] += buf[m++]; rho[j] += buf[m++]; } else { @@ -321,7 +318,7 @@ void ComputeRHEOInterface::correct_v(double *v_solid, double *v_fluid, int i_sol /* ---------------------------------------------------------------------- */ -double ComputeRHEOInterface::correct_rho(int i_solid, int i_fluid) +double ComputeRHEOInterface::correct_rho(int i_solid) { int itype = atom->type[i_solid]; return MAX(rho0[itype], atom->rho[i_solid]); @@ -347,28 +344,26 @@ void ComputeRHEOInterface::store_forces() if (fixlist.size() != 0) { for (const auto &fix : fixlist) { for (int i = 0; i < atom->nlocal; i++) { - if (rmass) minv = 1.0 / rmass[i]; - else minv = 1.0 / mass[type[i]]; - if (mask[i] & fix->groupbit) - for (int a = 0; a < 3; a++) - fp_store[i][a] = f[i][a] * minv; + if (rmass) + minv = 1.0 / rmass[i]; else - for (int a = 0; a < 3; a++) - fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv; + minv = 1.0 / mass[type[i]]; + if (mask[i] & fix->groupbit) + for (int a = 0; a < 3; a++) fp_store[i][a] = f[i][a] * minv; + else + for (int a = 0; a < 3; a++) fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv; } } } else { if (rmass) { for (int i = 0; i < atom->nlocal; i++) { minv = 1.0 / rmass[i]; - for (int a = 0; a < 3; a++) - fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv; + for (int a = 0; a < 3; a++) fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv; } } else { for (int i = 0; i < atom->nlocal; i++) { minv = 1.0 / mass[type[i]]; - for (int a = 0; a < 3; a++) - fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv; + for (int a = 0; a < 3; a++) fp_store[i][a] = (f[i][a] - fp_store[i][a]) * minv; } } } @@ -388,4 +383,3 @@ double ComputeRHEOInterface::memory_usage() double bytes = 3 * nmax_store * sizeof(double); return bytes; } - diff --git a/src/RHEO/compute_rheo_interface.h b/src/RHEO/compute_rheo_interface.h index a8cd448822..2f88eca50f 100644 --- a/src/RHEO/compute_rheo_interface.h +++ b/src/RHEO/compute_rheo_interface.h @@ -37,7 +37,7 @@ class ComputeRHEOInterface : public Compute { void unpack_reverse_comm(int, int *, double *) override; double memory_usage() override; void correct_v(double *, double *, int, int); - double correct_rho(int, int); + double correct_rho(int); void store_forces(); double *chi, **fp_store; diff --git a/src/RHEO/compute_rheo_kernel.cpp b/src/RHEO/compute_rheo_kernel.cpp index 108ec58b09..dd05901b32 100644 --- a/src/RHEO/compute_rheo_kernel.cpp +++ b/src/RHEO/compute_rheo_kernel.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -30,31 +29,34 @@ #include "math_extra.h" #include "memory.h" #include "modify.h" -#include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" +#include "neighbor.h" #include "pair.h" #include "update.h" #include "utils.h" #include +#include +#include #include #include -#include -#include using namespace LAMMPS_NS; using namespace RHEO_NS; using namespace MathConst; using namespace MathExtra; +// max value of Mdim 1 + dim + dim * (dim + 1) / 2 with dim = 3 +static constexpr int MAX_MDIM = 12; + /* ---------------------------------------------------------------------- */ ComputeRHEOKernel::ComputeRHEOKernel(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), coordination(nullptr), fix_rheo(nullptr), C(nullptr), C0(nullptr), - list(nullptr), compute_interface(nullptr) + Compute(lmp, narg, arg), coordination(nullptr), fix_rheo(nullptr), C(nullptr), C0(nullptr), + list(nullptr), compute_interface(nullptr) { - if (narg != 4) error->all(FLERR,"Illegal compute rheo/kernel command"); + if (narg != 4) error->all(FLERR, "Illegal compute rheo/kernel command"); kernel_style = utils::inumeric(FLERR, arg[3], false, lmp); @@ -114,13 +116,12 @@ void ComputeRHEOKernel::init() cutinv = 1.0 / cut; cutsqinv = cutinv * cutinv; - if (kernel_style != WENDLANDC4) { if (dim == 3) { pre_w = 1.0 / (120.0 * MY_PI) * 27.0 * cutsqinv * cutinv; pre_wp = pre_w * 3.0 * cutinv; } else { - pre_w = 7.0 / (478.0 * MY_PI) * 9 * cutsqinv; + pre_w = 7.0 / (478.0 * MY_PI) * 9 * cutsqinv; pre_wp = pre_w * 3.0 * cutinv; } } else { @@ -157,31 +158,25 @@ int ComputeRHEOKernel::check_corrections(int i) { // Skip if there were gsl errors for this atom if (gsl_error_flag) - if (gsl_error_tags.find(atom->tag[i]) != gsl_error_tags.end()) - return 0; + if (gsl_error_tags.find(atom->tag[i]) != gsl_error_tags.end()) return 0; // Skip if undercoordinated - if (coordination[i] < zmin) - return 0; + if (coordination[i] < zmin) return 0; // Skip if corrections not yet calculated - if (!corrections_calculated) - return 0; + if (!corrections_calculated) return 0; return 1; } /* ---------------------------------------------------------------------- */ -double ComputeRHEOKernel::calc_w_self(int i, int j) +double ComputeRHEOKernel::calc_w_self() { - double w; if (kernel_style == WENDLANDC4) - w = calc_w_wendlandc4(i, j, 0.0, 0.0, 0.0, 0.0); + return calc_w_wendlandc4(0.0); else - w = calc_w_quintic(i, j, 0.0, 0.0, 0.0, 0.0); - - return w; + return calc_w_quintic(0.0); } /* ---------------------------------------------------------------------- */ @@ -191,8 +186,7 @@ double ComputeRHEOKernel::calc_w(int i, int j, double delx, double dely, double double w = 0.0; int corrections_i, corrections_j, corrections; - if (kernel_style == WENDLANDC4) - return calc_w_wendlandc4(i, j, delx, dely, delz, r); + if (kernel_style == WENDLANDC4) return calc_w_wendlandc4(r); if (kernel_style != QUINTIC) { corrections_i = check_corrections(i); @@ -202,10 +196,14 @@ double ComputeRHEOKernel::calc_w(int i, int j, double delx, double dely, double corrections = 0; } - if (!corrections) w = calc_w_quintic(i, j, delx, dely, delz, r); - else if (kernel_style == RK0) w = calc_w_rk0(i, j, delx, dely, delz, r); - else if (kernel_style == RK1) w = calc_w_rk1(i, j, delx, dely, delz, r); - else if (kernel_style == RK2) w = calc_w_rk2(i, j, delx, dely, delz, r); + if (!corrections) + w = calc_w_quintic(r); + else if (kernel_style == RK0) + w = calc_w_rk0(i, j, r); + else if (kernel_style == RK1) + w = calc_w_rk1(i, j, delx, dely, delz, r); + else if (kernel_style == RK2) + w = calc_w_rk2(i, j, delx, dely, delz, r); return w; } @@ -217,8 +215,7 @@ double ComputeRHEOKernel::calc_dw(int i, int j, double delx, double dely, double double wp; int corrections_i, corrections_j; - if (kernel_style == WENDLANDC4) - return calc_dw_wendlandc4(i, j, delx, dely, delz, r, dWij, dWji); + if (kernel_style == WENDLANDC4) return calc_dw_wendlandc4(delx, dely, delz, r, dWij, dWji); if (kernel_style != QUINTIC) { corrections_i = check_corrections(i); @@ -226,15 +223,15 @@ double ComputeRHEOKernel::calc_dw(int i, int j, double delx, double dely, double } // Calc wp and default dW's, a bit inefficient but can redo later - wp = calc_dw_quintic(i, j, delx, dely, delz, r, dWij, dWji); + wp = calc_dw_quintic(delx, dely, delz, r, dWij, dWji); // Overwrite if there are corrections if (kernel_style == RK1) { - if (corrections_i) calc_dw_rk1(i, j, delx, dely, delz, r, dWij); - if (corrections_j) calc_dw_rk1(j, i, -delx, -dely, -delz, r, dWji); + if (corrections_i) calc_dw_rk1(i, delx, dely, delz, r, dWij); + if (corrections_j) calc_dw_rk1(j, -delx, -dely, -delz, r, dWji); } else if (kernel_style == RK2) { - if (corrections_i) calc_dw_rk2(i, j, delx, dely, delz, r, dWij); - if (corrections_j) calc_dw_rk2(j, i, -delx, -dely, -delz, r, dWji); + if (corrections_i) calc_dw_rk2(i, delx, dely, delz, r, dWij); + if (corrections_j) calc_dw_rk2(j, -delx, -dely, -delz, r, dWji); } return wp; @@ -242,30 +239,28 @@ double ComputeRHEOKernel::calc_dw(int i, int j, double delx, double dely, double /* ---------------------------------------------------------------------- */ -double ComputeRHEOKernel::calc_w_quintic(int i, int j, double delx, double dely, double delz, double r) +double ComputeRHEOKernel::calc_w_quintic(double r) { double w, tmp1, tmp2, tmp3, tmp1sq, tmp2sq, tmp3sq, s; s = r * 3.0 * cutinv; - if (s > 3.0) { - w = 0.0; - } + if (s > 3.0) { w = 0.0; } - if (s <= 3.0) { - tmp3 = 3.0 - s; - tmp3sq = tmp3 * tmp3; - w = tmp3sq * tmp3sq * tmp3; - } - if (s <= 2.0) { - tmp2 = 2.0 - s; - tmp2sq = tmp2 * tmp2; - w -= 6.0 * tmp2sq * tmp2sq * tmp2; - } - if (s <= 1.0) { - tmp1 = 1.0 - s; - tmp1sq = tmp1 * tmp1; - w += 15.0 * tmp1sq * tmp1sq * tmp1; - } + if (s <= 3.0) { + tmp3 = 3.0 - s; + tmp3sq = tmp3 * tmp3; + w = tmp3sq * tmp3sq * tmp3; + } + if (s <= 2.0) { + tmp2 = 2.0 - s; + tmp2sq = tmp2 * tmp2; + w -= 6.0 * tmp2sq * tmp2sq * tmp2; + } + if (s <= 1.0) { + tmp1 = 1.0 - s; + tmp1sq = tmp1 * tmp1; + w += 15.0 * tmp1sq * tmp1sq * tmp1; + } w *= pre_w; @@ -277,15 +272,14 @@ double ComputeRHEOKernel::calc_w_quintic(int i, int j, double delx, double dely, /* ---------------------------------------------------------------------- */ -double ComputeRHEOKernel::calc_dw_quintic(int i, int j, double delx, double dely, double delz, double r, double *dW1, double *dW2) +double ComputeRHEOKernel::calc_dw_quintic(double delx, double dely, double delz, double r, + double *dW1, double *dW2) { double wp, tmp1, tmp2, tmp3, tmp1sq, tmp2sq, tmp3sq, s, wprinv; s = r * 3.0 * cutinv; - if (s > 3.0) { - wp = 0.0; - } + if (s > 3.0) { wp = 0.0; } if (s <= 3.0) { tmp3 = 3.0 - s; tmp3sq = tmp3 * tmp3; @@ -317,7 +311,7 @@ double ComputeRHEOKernel::calc_dw_quintic(int i, int j, double delx, double dely /* ---------------------------------------------------------------------- */ -double ComputeRHEOKernel::calc_w_wendlandc4(int i, int j, double delx, double dely, double delz, double r) +double ComputeRHEOKernel::calc_w_wendlandc4(double r) { double w, tmp6, s; s = r * cutinv; @@ -340,7 +334,8 @@ double ComputeRHEOKernel::calc_w_wendlandc4(int i, int j, double delx, double de /* ---------------------------------------------------------------------- */ -double ComputeRHEOKernel::calc_dw_wendlandc4(int i, int j, double delx, double dely, double delz, double r, double *dW1, double *dW2) +double ComputeRHEOKernel::calc_dw_wendlandc4(double delx, double dely, double delz, double r, + double *dW1, double *dW2) { double wp, tmp1, tmp5, tmp6, s, wprinv; @@ -372,11 +367,11 @@ double ComputeRHEOKernel::calc_dw_wendlandc4(int i, int j, double delx, double d /* ---------------------------------------------------------------------- */ -double ComputeRHEOKernel::calc_w_rk0(int i, int j, double delx, double dely, double delz, double r) +double ComputeRHEOKernel::calc_w_rk0(int i, int j, double r) { double w; - w = calc_w_quintic(i, j, delx, dely, delz, r); + w = calc_w_quintic(r); Wij = C0[i] * w; Wji = C0[j] * w; @@ -389,12 +384,12 @@ double ComputeRHEOKernel::calc_w_rk0(int i, int j, double delx, double dely, dou double ComputeRHEOKernel::calc_w_rk1(int i, int j, double delx, double dely, double delz, double r) { int b; - double w, dx[3], H[Mdim]; + double w, dx[3], H[MAX_MDIM]; dx[0] = delx; dx[1] = dely; dx[2] = delz; - w = calc_w_quintic(i, j, delx, dely, delz, r); + w = calc_w_quintic(r); if (dim == 2) { H[0] = 1.0; @@ -408,7 +403,7 @@ double ComputeRHEOKernel::calc_w_rk1(int i, int j, double delx, double dely, dou } Wij = 0; for (b = 0; b < Mdim; b++) { - Wij += C[i][0][b] * H[b]; // C columns: 1 x y (z) xx yy (zz) + Wij += C[i][0][b] * H[b]; // C columns: 1 x y (z) xx yy (zz) } Wij *= w; @@ -419,7 +414,7 @@ double ComputeRHEOKernel::calc_w_rk1(int i, int j, double delx, double dely, dou Wji = 0; for (b = 0; b < Mdim; b++) { - Wji += C[j][0][b] * H[b]; // C columns: 1 x y (z) xx yy (zz) + Wji += C[j][0][b] * H[b]; // C columns: 1 x y (z) xx yy (zz) } Wji *= w; @@ -431,11 +426,11 @@ double ComputeRHEOKernel::calc_w_rk1(int i, int j, double delx, double dely, dou double ComputeRHEOKernel::calc_w_rk2(int i, int j, double delx, double dely, double delz, double r) { int b; - double w, dx[3], H[Mdim]; + double w, dx[3], H[MAX_MDIM]; dx[0] = delx; dx[1] = dely; dx[2] = delz; - w = calc_w_quintic(i, j, delx, dely, delz, r); + w = calc_w_quintic(r); if (dim == 2) { H[0] = 1.0; @@ -458,7 +453,7 @@ double ComputeRHEOKernel::calc_w_rk2(int i, int j, double delx, double dely, dou } Wij = 0; for (b = 0; b < Mdim; b++) { - Wij += C[i][0][b] * H[b]; // C columns: 1 x y (z) xx yy (zz) + Wij += C[i][0][b] * H[b]; // C columns: 1 x y (z) xx yy (zz) } Wij *= w; @@ -469,7 +464,7 @@ double ComputeRHEOKernel::calc_w_rk2(int i, int j, double delx, double dely, dou Wji = 0; for (b = 0; b < Mdim; b++) { - Wji += C[j][0][b] * H[b]; // C columns: 1 x y (z) xx yy (zz) + Wji += C[j][0][b] * H[b]; // C columns: 1 x y (z) xx yy (zz) } Wji *= w; @@ -478,15 +473,16 @@ double ComputeRHEOKernel::calc_w_rk2(int i, int j, double delx, double dely, dou /* ---------------------------------------------------------------------- */ -void ComputeRHEOKernel::calc_dw_rk1(int i, int j, double delx, double dely, double delz, double r, double *dW) +void ComputeRHEOKernel::calc_dw_rk1(int i, double delx, double dely, double delz, double r, + double *dW) { int a, b; - double w, dx[3], H[Mdim]; + double w, dx[3], H[MAX_MDIM]; dx[0] = delx; dx[1] = dely; dx[2] = delz; - w = calc_w_quintic(i, j, delx, dely, delz, r); + w = calc_w_quintic(r); //Populate correction basis if (dim == 2) { @@ -506,24 +502,24 @@ void ComputeRHEOKernel::calc_dw_rk1(int i, int j, double delx, double dely, doub dW[a] = 0.0; for (b = 0; b < Mdim; b++) { //First derivative kernels - dW[a] += C[i][1 + a][b] * H[b]; // C columns: 1 x y (z) + dW[a] += C[i][1 + a][b] * H[b]; // C columns: 1 x y (z) } dW[a] *= w; } } - /* ---------------------------------------------------------------------- */ -void ComputeRHEOKernel::calc_dw_rk2(int i, int j, double delx, double dely, double delz, double r, double *dW) +void ComputeRHEOKernel::calc_dw_rk2(int i, double delx, double dely, double delz, double r, + double *dW) { int a, b; - double w, dx[3], H[Mdim]; + double w, dx[3], H[MAX_MDIM]; dx[0] = delx; dx[1] = dely; dx[2] = delz; - w = calc_w_quintic(i, j, delx, dely, delz, r); + w = calc_w_quintic(r); //Populate correction basis if (dim == 2) { @@ -552,7 +548,7 @@ void ComputeRHEOKernel::calc_dw_rk2(int i, int j, double delx, double dely, doub dW[a] = 0.0; for (b = 0; b < Mdim; b++) { //First derivative kernels - dW[a] += C[i][1 + a][b] * H[b]; // C columns: 1 x y (z) xx yy (zz) + dW[a] += C[i][1 + a][b] * H[b]; // C columns: 1 x y (z) xx yy (zz) } dW[a] *= w; } @@ -619,14 +615,15 @@ void ComputeRHEOKernel::compute_peratom() if (rsq < cutsq) { r = sqrt(rsq); - w = calc_w_quintic(i, j, dx[0], dx[1], dx[2], r); + w = calc_w_quintic(r); rhoj = rho[j]; if (interface_flag) - if (status[j] & PHASECHECK) - rhoj = compute_interface->correct_rho(j, i); + if (status[j] & PHASECHECK) rhoj = compute_interface->correct_rho(j); - if (rmass) vj = rmass[j] / rhoj; - else vj = mass[type[j]] / rhoj; + if (rmass) + vj = rmass[j] / rhoj; + else + vj = mass[type[j]] / rhoj; M += w * vj; } } @@ -637,7 +634,7 @@ void ComputeRHEOKernel::compute_peratom() } else if (correction_order > 0) { // Moment matrix M and polynomial basis vector cut (1d for gsl compatibility) - double H[Mdim], M[Mdim * Mdim]; + double H[MAX_MDIM], M[MAX_MDIM * MAX_MDIM]; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; @@ -650,9 +647,7 @@ void ComputeRHEOKernel::compute_peratom() // Zero upper-triangle M and cut (will be symmetric): for (a = 0; a < Mdim; a++) { - for (b = a; b < Mdim; b++) { - M[a * Mdim + b] = 0; - } + for (b = a; b < Mdim; b++) { M[a * Mdim + b] = 0; } } for (jj = 0; jj < jnum; jj++) { @@ -667,15 +662,16 @@ void ComputeRHEOKernel::compute_peratom() if (rsq < cutsq) { r = sqrt(rsq); - w = calc_w_quintic(i, j, dx[0], dx[1], dx[2], r); + w = calc_w_quintic(r); rhoj = rho[j]; if (interface_flag) - if (status[j] & PHASECHECK) - rhoj = compute_interface->correct_rho(j, i); + if (status[j] & PHASECHECK) rhoj = compute_interface->correct_rho(j); - if (rmass) vj = rmass[j] / rhoj; - else vj = mass[type[j]] / rhoj; + if (rmass) + vj = rmass[j] / rhoj; + else + vj = mass[type[j]] / rhoj; //Populate the H-vector of polynomials (2D) if (dim == 2) { @@ -704,18 +700,14 @@ void ComputeRHEOKernel::compute_peratom() // Populate the upper triangle for (a = 0; a < Mdim; a++) { - for (b = a; b < Mdim; b++) { - M[a * Mdim + b] += H[a] * H[b] * w * vj; - } + for (b = a; b < Mdim; b++) { M[a * Mdim + b] += H[a] * H[b] * w * vj; } } } } // Populate the lower triangle from the symmetric entries of M: for (a = 0; a < Mdim; a++) { - for (b = a; b < Mdim; b++) { - M[b * Mdim + a] = M[a * Mdim + b]; - } + for (b = a; b < Mdim; b++) { M[b * Mdim + a] = M[a * Mdim + b]; } } // Skip if undercoordinated @@ -723,7 +715,7 @@ void ComputeRHEOKernel::compute_peratom() // Use gsl to get Minv, use Cholesky decomposition since the // polynomials are independent, M is symmetrix & positive-definite - gM = gsl_matrix_view_array(M,Mdim,Mdim); + gM = gsl_matrix_view_array(M, Mdim, Mdim); gsl_error = gsl_linalg_cholesky_decomp1(&gM.matrix); if (gsl_error) { @@ -738,7 +730,7 @@ void ComputeRHEOKernel::compute_peratom() continue; } - gsl_linalg_cholesky_invert(&gM.matrix); //M is now M^-1 + gsl_linalg_cholesky_invert(&gM.matrix); //M is now M^-1 // Correction coefficients are columns of M^-1 multiplied by an appropriate coefficient // Solve the linear system several times to get coefficientns @@ -759,16 +751,15 @@ void ComputeRHEOKernel::compute_peratom() // Pack coefficients into C for (a = 0; a < Mdim; a++) { - C[i][0][a] = M[a * Mdim + 0]; // all rows of column 0 + C[i][0][a] = M[a * Mdim + 0]; // all rows of column 0 for (b = 0; b < dim; b++) { //First derivatives C[i][1 + b][a] = -M[a * Mdim + b + 1] * cutinv; // columns 1-2 (2D) or 1-3 (3D) //Second derivatives - if (kernel_style == RK2) - C[i][1 + dim + b][a] = M[a * Mdim + b + 1 + dim] * cutsqinv; - // columns 3-4 (2D) or 4-6 (3D) + if (kernel_style == RK2) C[i][1 + dim + b][a] = M[a * Mdim + b + 1 + dim] * cutsqinv; + // columns 3-4 (2D) or 4-6 (3D) } } } @@ -780,7 +771,6 @@ void ComputeRHEOKernel::compute_peratom() comm->forward_comm(this); } - /* ---------------------------------------------------------------------- */ void ComputeRHEOKernel::compute_coordination() @@ -819,8 +809,7 @@ void ComputeRHEOKernel::compute_coordination() dx[2] = ztmp - x[j][2]; rsq = lensq3(dx); - if (rsq < cutsq) - coordination[i] += 1; + if (rsq < cutsq) coordination[i] += 1; } } @@ -847,8 +836,8 @@ void ComputeRHEOKernel::grow_arrays(int nmax) /* ---------------------------------------------------------------------- */ -int ComputeRHEOKernel::pack_forward_comm(int n, int *list, double *buf, - int /*pbc_flag*/, int * /*pbc*/) +int ComputeRHEOKernel::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, + int * /*pbc*/) { int m = 0; for (int i = 0; i < n; i++) { @@ -860,8 +849,7 @@ int ComputeRHEOKernel::pack_forward_comm(int n, int *list, double *buf, buf[m++] = C0[j]; } else { for (int a = 0; a < ncor; a++) - for (int b = 0; b < Mdim; b++) - buf[m++] = C[j][a][b]; + for (int b = 0; b < Mdim; b++) buf[m++] = C[j][a][b]; } } } @@ -882,8 +870,7 @@ void ComputeRHEOKernel::unpack_forward_comm(int n, int first, double *buf) C0[i] = buf[m++]; } else { for (int a = 0; a < ncor; a++) - for (int b = 0; b < Mdim; b++) - C[i][a][b] = buf[m++]; + for (int b = 0; b < Mdim; b++) C[i][a][b] = buf[m++]; } } } diff --git a/src/RHEO/compute_rheo_kernel.h b/src/RHEO/compute_rheo_kernel.h index d15e8e210a..20516255be 100644 --- a/src/RHEO/compute_rheo_kernel.h +++ b/src/RHEO/compute_rheo_kernel.h @@ -36,13 +36,13 @@ class ComputeRHEOKernel : public Compute { void unpack_forward_comm(int, int, double *) override; double memory_usage() override; void compute_coordination(); - double calc_w_self(int,int); - double calc_w(int,int,double,double,double,double); - double calc_dw(int,int,double,double,double,double); - double calc_w_quintic(int,int,double,double,double,double); - double calc_dw_quintic(int,int,double,double,double,double,double *,double *); - double calc_w_wendlandc4(int,int,double,double,double,double); - double calc_dw_wendlandc4(int,int,double,double,double,double,double *,double *); + double calc_w_self(); + double calc_w(int, int, double, double, double, double); + double calc_dw(int, int, double, double, double, double); + double calc_w_quintic(double); + double calc_dw_quintic(double, double, double, double, double *, double *); + double calc_w_wendlandc4(double); + double calc_dw_wendlandc4(double, double, double, double, double *, double *); void grow_arrays(int); double dWij[3], dWji[3], Wij, Wji; @@ -68,14 +68,12 @@ class ComputeRHEOKernel : public Compute { int check_corrections(int); - double calc_w_rk0(int,int,double,double,double,double); - double calc_w_rk1(int,int,double,double,double,double); - double calc_w_rk2(int,int,double,double,double,double); - void calc_dw_rk1(int,int,double,double,double,double,double *); - void calc_dw_rk2(int,int,double,double,double,double,double *); + double calc_w_rk0(int, int, double); + double calc_w_rk1(int, int, double, double, double, double); + double calc_w_rk2(int, int, double, double, double, double); + void calc_dw_rk1(int, double, double, double, double, double *); + void calc_dw_rk2(int, double, double, double, double, double *); }; - } // namespace LAMMPS_NS - #endif #endif diff --git a/src/RHEO/compute_rheo_property_atom.cpp b/src/RHEO/compute_rheo_property_atom.cpp index cd0ff6692d..fb5848ac75 100644 --- a/src/RHEO/compute_rheo_property_atom.cpp +++ b/src/RHEO/compute_rheo_property_atom.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -21,11 +20,11 @@ #include "atom.h" #include "atom_vec.h" +#include "compute_rheo_grad.h" #include "compute_rheo_interface.h" #include "compute_rheo_kernel.h" #include "compute_rheo_surface.h" #include "compute_rheo_vshift.h" -#include "compute_rheo_grad.h" #include "domain.h" #include "error.h" #include "fix_rheo.h" @@ -46,12 +45,12 @@ using namespace RHEO_NS; /* ---------------------------------------------------------------------- */ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), avec_index(nullptr), col_index(nullptr), col_t_index(nullptr), - buf(nullptr), pack_choice(nullptr), fix_rheo(nullptr), fix_pressure(nullptr), - fix_thermal(nullptr), compute_interface(nullptr), compute_kernel(nullptr), - compute_surface(nullptr), compute_vshift(nullptr), compute_grad(nullptr) + Compute(lmp, narg, arg), avec_index(nullptr), col_index(nullptr), col_t_index(nullptr), + buf(nullptr), pack_choice(nullptr), fix_rheo(nullptr), fix_pressure(nullptr), + fix_thermal(nullptr), compute_interface(nullptr), compute_kernel(nullptr), + compute_surface(nullptr), compute_vshift(nullptr), compute_grad(nullptr) { - if (narg < 4) utils::missing_cmd_args(FLERR, "compute property/atom", error); + if (narg < 4) utils::missing_cmd_args(FLERR, "compute property/atom", error); peratom_flag = 1; int dim = domain->dimension; @@ -74,8 +73,10 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a } } - if (nvalues == 1) size_peratom_cols = 0; - else size_peratom_cols = nvalues; + if (nvalues == 1) + size_peratom_cols = 0; + else + size_peratom_cols = nvalues; pressure_flag = thermal_flag = interface_flag = 0; surface_flag = shift_flag = shell_flag = 0; @@ -173,15 +174,19 @@ ComputeRHEOPropertyAtom::~ComputeRHEOPropertyAtom() void ComputeRHEOPropertyAtom::init() { auto fixes = modify->get_fix_by_style("^rheo$"); - if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use compute rheo/property/atom"); + if (fixes.size() == 0) + error->all(FLERR, "Need to define fix rheo to use compute rheo/property/atom"); fix_rheo = dynamic_cast(fixes[0]); if (interface_flag && !(fix_rheo->interface_flag)) - error->all(FLERR, "Cannot request interfacial property without corresponding option in fix rheo"); + error->all(FLERR, + "Cannot request interfacial property without corresponding option in fix rheo"); if (surface_flag && !(fix_rheo->surface_flag)) error->all(FLERR, "Cannot request surface property without corresponding option in fix rheo"); if (shift_flag && !(fix_rheo->shift_flag)) - error->all(FLERR, "Cannot request velocity shifting property without corresponding option in fix rheo"); + error->all( + FLERR, + "Cannot request velocity shifting property without corresponding option in fix rheo"); if (thermal_flag && !(fix_rheo->thermal_flag)) error->all(FLERR, "Cannot request thermal property without fix rheo/thermal"); @@ -237,10 +242,11 @@ void ComputeRHEOPropertyAtom::compute_peratom() buf = vector_atom; (this->*pack_choice[0])(0); } else { - if (nmax) buf = &array_atom[0][0]; - else buf = nullptr; - for (int n = 0; n < nvalues; n++) - (this->*pack_choice[n])(n); + if (nmax) + buf = &array_atom[0][0]; + else + buf = nullptr; + for (int n = 0; n < nvalues; n++) (this->*pack_choice[n])(n); } } @@ -250,7 +256,7 @@ void ComputeRHEOPropertyAtom::compute_peratom() double ComputeRHEOPropertyAtom::memory_usage() { - double bytes = (double)nmax * nvalues * sizeof(double); + double bytes = (double) nmax * nvalues * sizeof(double); return bytes; } @@ -269,8 +275,10 @@ void ComputeRHEOPropertyAtom::pack_phase(int n) int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) buf[n] = (status[i] & PHASECHECK); - else buf[n] = 0.0; + if (mask[i] & groupbit) + buf[n] = (status[i] & PHASECHECK); + else + buf[n] = 0.0; n += nvalues; } } @@ -284,8 +292,10 @@ void ComputeRHEOPropertyAtom::pack_status(int n) int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) buf[n] = status[i]; - else buf[n] = 0.0; + if (mask[i] & groupbit) + buf[n] = status[i]; + else + buf[n] = 0.0; n += nvalues; } } @@ -299,8 +309,10 @@ void ComputeRHEOPropertyAtom::pack_chi(int n) int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) buf[n] = chi[i]; - else buf[n] = 0.0; + if (mask[i] & groupbit) + buf[n] = chi[i]; + else + buf[n] = 0.0; n += nvalues; } } @@ -334,8 +346,10 @@ void ComputeRHEOPropertyAtom::pack_surface_r(int n) int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) buf[n] = rsurface[i]; - else buf[n] = 0.0; + if (mask[i] & groupbit) + buf[n] = rsurface[i]; + else + buf[n] = 0.0; n += nvalues; } } @@ -349,8 +363,10 @@ void ComputeRHEOPropertyAtom::pack_surface_divr(int n) int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) buf[n] = divr[i]; - else buf[n] = 0.0; + if (mask[i] & groupbit) + buf[n] = divr[i]; + else + buf[n] = 0.0; n += nvalues; } } @@ -365,8 +381,10 @@ void ComputeRHEOPropertyAtom::pack_surface_n(int n) int index = col_index[n]; for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) buf[n] = nsurface[i][index]; - else buf[n] = 0.0; + if (mask[i] & groupbit) + buf[n] = nsurface[i][index]; + else + buf[n] = 0.0; n += nvalues; } } @@ -380,8 +398,10 @@ void ComputeRHEOPropertyAtom::pack_coordination(int n) int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) buf[n] = coordination[i]; - else buf[n] = 0.0; + if (mask[i] & groupbit) + buf[n] = coordination[i]; + else + buf[n] = 0.0; n += nvalues; } } @@ -395,8 +415,10 @@ void ComputeRHEOPropertyAtom::pack_cv(int n) int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) buf[n] = fix_thermal->calc_cv(i, type[i]); - else buf[n] = 0.0; + if (mask[i] & groupbit) + buf[n] = fix_thermal->calc_cv(type[i]); + else + buf[n] = 0.0; n += nvalues; } } @@ -411,8 +433,10 @@ void ComputeRHEOPropertyAtom::pack_pressure(int n) int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) buf[n] = fix_pressure->calc_pressure(rho[i], type[i]); - else buf[n] = 0.0; + if (mask[i] & groupbit) + buf[n] = fix_pressure->calc_pressure(rho[i], type[i]); + else + buf[n] = 0.0; n += nvalues; } } @@ -432,8 +456,10 @@ void ComputeRHEOPropertyAtom::pack_viscous_stress(int n) int index_transpose = b * dim + a; for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) buf[n] = viscosity[i] * (gradv[i][index] + gradv[i][index_transpose]); - else buf[n] = 0.0; + if (mask[i] & groupbit) + buf[n] = viscosity[i] * (gradv[i][index] + gradv[i][index_transpose]); + else + buf[n] = 0.0; n += nvalues; } } @@ -462,7 +488,8 @@ void ComputeRHEOPropertyAtom::pack_total_stress(int n) else p = 0.0; buf[n] = viscosity[i] * (gradv[i][index] + gradv[i][index_transpose]) + p; - } else buf[n] = 0.0; + } else + buf[n] = 0.0; n += nvalues; } } @@ -476,8 +503,10 @@ void ComputeRHEOPropertyAtom::pack_nbond_shell(int n) int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) buf[n] = nbond[i]; - else buf[n] = 0.0; + if (mask[i] & groupbit) + buf[n] = nbond[i]; + else + buf[n] = 0.0; n += nvalues; } } @@ -492,8 +521,10 @@ void ComputeRHEOPropertyAtom::pack_shift_v(int n) int index = col_index[n]; for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) buf[n] = vshift[i][index]; - else buf[n] = 0.0; + if (mask[i] & groupbit) + buf[n] = vshift[i][index]; + else + buf[n] = 0.0; n += nvalues; } } @@ -508,8 +539,10 @@ void ComputeRHEOPropertyAtom::pack_gradv(int n) int index = col_index[n]; for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) buf[n] = gradv[i][index]; - else buf[n] = 0.0; + if (mask[i] & groupbit) + buf[n] = gradv[i][index]; + else + buf[n] = 0.0; n += nvalues; } } @@ -523,7 +556,7 @@ void ComputeRHEOPropertyAtom::pack_atom_style(int n) /* ---------------------------------------------------------------------- */ -int ComputeRHEOPropertyAtom::add_tensor_component(char* option, int i, FnPtrPack pack_function) +int ComputeRHEOPropertyAtom::add_tensor_component(char *option, int i, FnPtrPack pack_function) { int shift; int dim = domain->dimension; @@ -547,11 +580,15 @@ int ComputeRHEOPropertyAtom::add_tensor_component(char* option, int i, FnPtrPack index = 2; if (dim == 2) dim_error = 1; } else if (utils::strmatch(option, "yx$")) { - if (dim == 2) index = 2; - else index = 3; + if (dim == 2) + index = 2; + else + index = 3; } else if (utils::strmatch(option, "yy$")) { - if (dim == 2) index = 3; - else index = 4; + if (dim == 2) + index = 3; + else + index = 4; } else if (utils::strmatch(option, "yz$")) { index = 5; if (dim == 2) dim_error = 1; @@ -581,7 +618,7 @@ int ComputeRHEOPropertyAtom::add_tensor_component(char* option, int i, FnPtrPack /* ---------------------------------------------------------------------- */ -int ComputeRHEOPropertyAtom::add_vector_component(char* option, int i, FnPtrPack pack_function) +int ComputeRHEOPropertyAtom::add_vector_component(char *option, int i, FnPtrPack pack_function) { int shift; int dim = domain->dimension; diff --git a/src/RHEO/compute_rheo_property_atom.h b/src/RHEO/compute_rheo_property_atom.h index 4b1ebf2313..14f28075a3 100644 --- a/src/RHEO/compute_rheo_property_atom.h +++ b/src/RHEO/compute_rheo_property_atom.h @@ -61,8 +61,8 @@ class ComputeRHEOPropertyAtom : public Compute { void pack_nbond_shell(int); void pack_atom_style(int); - int add_vector_component(char*, int, FnPtrPack); - int add_tensor_component(char*, int, FnPtrPack); + int add_vector_component(char *, int, FnPtrPack); + int add_tensor_component(char *, int, FnPtrPack); class FixRHEO *fix_rheo; class FixRHEOPressure *fix_pressure; @@ -73,11 +73,9 @@ class ComputeRHEOPropertyAtom : public Compute { class ComputeRHEOSurface *compute_surface; class ComputeRHEOVShift *compute_vshift; class ComputeRHEOGrad *compute_grad; - }; } // namespace LAMMPS_NS #endif #endif - diff --git a/src/RHEO/compute_rheo_rho_sum.cpp b/src/RHEO/compute_rheo_rho_sum.cpp index 8067bae21b..db85b028fb 100644 --- a/src/RHEO/compute_rheo_rho_sum.cpp +++ b/src/RHEO/compute_rheo_rho_sum.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -25,18 +24,18 @@ #include "error.h" #include "fix_rheo.h" #include "force.h" -#include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" +#include "neighbor.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeRHEORhoSum::ComputeRHEORhoSum(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), fix_rheo(nullptr), compute_kernel(nullptr) + Compute(lmp, narg, arg), fix_rheo(nullptr), compute_kernel(nullptr) { - if (narg != 4) error->all(FLERR,"Illegal compute RHEO/rho command"); + if (narg != 4) error->all(FLERR, "Illegal compute RHEO/rho command"); self_mass_flag = utils::bnumeric(FLERR, arg[3], false, lmp); @@ -69,7 +68,6 @@ void ComputeRHEORhoSum::init_list(int /*id*/, NeighList *ptr) /* ---------------------------------------------------------------------- */ - void ComputeRHEORhoSum::compute_peratom() { int i, j, ii, jj, inum, jnum; @@ -94,9 +92,11 @@ void ComputeRHEORhoSum::compute_peratom() // initialize arrays, local with quintic self-contribution, ghosts are zeroed for (i = 0; i < nlocal; i++) { - w = compute_kernel->calc_w_self(i, i); - if (rmass) rho[i] = w * rmass[i]; - else rho[i] = w * mass[type[i]]; + w = compute_kernel->calc_w_self(); + if (rmass) + rho[i] = w * rmass[i]; + else + rho[i] = w * mass[type[i]]; } for (i = nlocal; i < nall; i++) rho[i] = 0.0; @@ -123,22 +123,18 @@ void ComputeRHEORhoSum::compute_peratom() if (rmass) { if (self_mass_flag) { rho[i] += w * rmass[i]; - if (newton || j < nlocal) - rho[j] += w * rmass[j]; + if (newton || j < nlocal) rho[j] += w * rmass[j]; } else { rho[i] += w * rmass[j]; - if (newton || j < nlocal) - rho[j] += w * rmass[i]; + if (newton || j < nlocal) rho[j] += w * rmass[i]; } } else { if (self_mass_flag) { rho[i] += w * mass[type[i]]; - if (newton || j < nlocal) - rho[j] += w * mass[type[j]]; + if (newton || j < nlocal) rho[j] += w * mass[type[j]]; } else { rho[i] += w * mass[type[j]]; - if (newton || j < nlocal) - rho[j] += w * mass[type[i]]; + if (newton || j < nlocal) rho[j] += w * mass[type[i]]; } } } @@ -151,8 +147,8 @@ void ComputeRHEORhoSum::compute_peratom() /* ---------------------------------------------------------------------- */ -int ComputeRHEORhoSum::pack_forward_comm(int n, int *list, double *buf, - int /*pbc_flag*/, int * /*pbc*/) +int ComputeRHEORhoSum::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, + int * /*pbc*/) { int i, j, m; double *rho = atom->rho; @@ -173,9 +169,7 @@ void ComputeRHEORhoSum::unpack_forward_comm(int n, int first, double *buf) m = 0; last = first + n; - for (i = first; i < last; i++) { - rho[i] = buf[m++]; - } + for (i = first; i < last; i++) { rho[i] = buf[m++]; } } /* ---------------------------------------------------------------------- */ @@ -187,9 +181,7 @@ int ComputeRHEORhoSum::pack_reverse_comm(int n, int first, double *buf) m = 0; last = first + n; - for (i = first; i < last; i++) { - buf[m++] = rho[i]; - } + for (i = first; i < last; i++) { buf[m++] = rho[i]; } return m; } diff --git a/src/RHEO/compute_rheo_surface.cpp b/src/RHEO/compute_rheo_surface.cpp index 314d2d2ef5..4ff9587f6c 100644 --- a/src/RHEO/compute_rheo_surface.cpp +++ b/src/RHEO/compute_rheo_surface.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -29,9 +28,9 @@ #include "force.h" #include "math_extra.h" #include "memory.h" -#include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" +#include "neighbor.h" using namespace LAMMPS_NS; using namespace RHEO_NS; @@ -43,11 +42,11 @@ static constexpr double EPSILON = 1e-10; /* ---------------------------------------------------------------------- */ ComputeRHEOSurface::ComputeRHEOSurface(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), nsurface(nullptr), rsurface(nullptr), divr(nullptr), fix_rheo(nullptr), - rho0(nullptr), B(nullptr), gradC(nullptr), list(nullptr), compute_kernel(nullptr), - compute_interface(nullptr) + Compute(lmp, narg, arg), nsurface(nullptr), rsurface(nullptr), divr(nullptr), fix_rheo(nullptr), + rho0(nullptr), B(nullptr), gradC(nullptr), list(nullptr), compute_kernel(nullptr), + compute_interface(nullptr) { - if (narg != 3) error->all(FLERR,"Illegal compute RHEO/SURFACE command"); + if (narg != 3) error->all(FLERR, "Illegal compute RHEO/SURFACE command"); int dim = domain->dimension; comm_forward = 2; @@ -123,8 +122,7 @@ void ComputeRHEOSurface::compute_peratom() firstneigh = list->firstneigh; // Grow and zero arrays - if (nmax_store < atom->nmax) - grow_arrays(atom->nmax); + if (nmax_store < atom->nmax) grow_arrays(atom->nmax); size_t nbytes = nmax_store * sizeof(double); memset(&divr[0], 0, nbytes); @@ -133,7 +131,6 @@ void ComputeRHEOSurface::compute_peratom() memset(&gradC[0][0], 0, dim * dim * nbytes); memset(&B[0][0], 0, dim * dim * nbytes); - // loop over neighbors to calculate the average orientation of neighbors for (ii = 0; ii < inum; ii++) { i = ilist[ii]; @@ -165,9 +162,9 @@ void ComputeRHEOSurface::compute_peratom() // Add corrections for walls if (interface_flag) { if (fluidi && (!fluidj)) { - rhoj = compute_interface->correct_rho(j, i); + rhoj = compute_interface->correct_rho(j); } else if ((!fluidi) && fluidj) { - rhoi = compute_interface->correct_rho(i, j); + rhoi = compute_interface->correct_rho(i); } else if ((!fluidi) && (!fluidj)) { rhoi = rho0[itype]; rhoj = rho0[jtype]; @@ -181,7 +178,7 @@ void ComputeRHEOSurface::compute_peratom() Voli = mass[itype] / rhoi; Volj = mass[jtype] / rhoj; } - compute_kernel->calc_dw_quintic(i, j, dx[0], dx[1], dx[2], sqrt(rsq), dWij, dWji); + compute_kernel->calc_dw_quintic(dx[0], dx[1], dx[2], sqrt(rsq), dWij, dWji); for (a = 0; a < dim; a++) { divr[i] -= dWij[a] * dx[a] * Volj; @@ -189,7 +186,7 @@ void ComputeRHEOSurface::compute_peratom() } if ((j < nlocal) || newton) { - for (a = 0; a < dim; a++){ + for (a = 0; a < dim; a++) { divr[j] += dWji[a] * dx[a] * Voli; gradC[j][a] += dWji[a] * Voli; } @@ -211,12 +208,10 @@ void ComputeRHEOSurface::compute_peratom() for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { maggC = 0.0; - for (a = 0;a < dim; a++) - maggC += gradC[i][a] * gradC[i][a]; + for (a = 0; a < dim; a++) maggC += gradC[i][a] * gradC[i][a]; maggC = sqrt(maggC) + EPSILON; maggC = 1.0 / maggC; - for (a = 0; a < dim; a++) - nsurface[i][a] = -gradC[i][a] * maggC; + for (a = 0; a < dim; a++) nsurface[i][a] = -gradC[i][a] * maggC; } } @@ -233,8 +228,7 @@ void ComputeRHEOSurface::compute_peratom() test = coordination[i] < threshold_z; // Treat nonfluid particles as bulk - if (status[i] & PHASECHECK) - test = 0; + if (status[i] & PHASECHECK) test = 0; if (test) { if (coordination[i] < threshold_splash) @@ -249,7 +243,6 @@ void ComputeRHEOSurface::compute_peratom() } } - for (ii = 0; ii < inum; ii++) { i = ilist[ii]; xtmp = x[i][0]; @@ -277,18 +270,17 @@ void ComputeRHEOSurface::compute_peratom() status[i] |= STATUS_LAYER; } - if (status[j] & STATUS_SURFACE) - rsurface[i] = MIN(rsurface[i], sqrt(rsq)); + if (status[j] & STATUS_SURFACE) rsurface[i] = MIN(rsurface[i], sqrt(rsq)); } if (fluidj && (j < nlocal || newton)) { - if ((status[j] & STATUS_BULK) && (status[j] & PHASECHECK) && (status[i] & STATUS_SURFACE)) { + if ((status[j] & STATUS_BULK) && (status[j] & PHASECHECK) && + (status[i] & STATUS_SURFACE)) { status[j] &= SURFACEMASK; status[j] |= STATUS_LAYER; } - if (status[i] & STATUS_SURFACE) - rsurface[j] = MIN(rsurface[j], sqrt(rsq)); + if (status[i] & STATUS_SURFACE) rsurface[j] = MIN(rsurface[j], sqrt(rsq)); } } } @@ -299,8 +291,7 @@ void ComputeRHEOSurface::compute_peratom() for (i = 0; i < nall; i++) { if (mask[i] & groupbit) { if (!(status[i] & STATUS_SURFACE)) - for (a = 0; a < dim; a++) - nsurface[i][a] = 0.0; + for (a = 0; a < dim; a++) nsurface[i][a] = 0.0; } } @@ -323,9 +314,8 @@ int ComputeRHEOSurface::pack_reverse_comm(int n, int first, double *buf) for (int i = first; i < last; i++) { if (comm_stage == 0) { buf[m++] = divr[i]; - for (int a = 0; a < dim; a ++ ) - for (int b = 0; b < dim; b ++) - buf[m++] = gradC[i][a * dim + b]; + for (int a = 0; a < dim; a++) + for (int b = 0; b < dim; b++) buf[m++] = gradC[i][a * dim + b]; } else if (comm_stage == 1) { buf[m++] = (double) status[i]; buf[m++] = rsurface[i]; @@ -345,9 +335,8 @@ void ComputeRHEOSurface::unpack_reverse_comm(int n, int *list, double *buf) int j = list[i]; if (comm_stage == 0) { divr[j] += buf[m++]; - for (int a = 0; a < dim; a ++ ) - for (int b = 0; b < dim; b ++) - gradC[j][a * dim + b] += buf[m++]; + for (int a = 0; a < dim; a++) + for (int b = 0; b < dim; b++) gradC[j][a * dim + b] += buf[m++]; } else if (comm_stage == 1) { auto tmp1 = (int) buf[m++]; if ((status[j] & STATUS_BULK) && (tmp1 & STATUS_LAYER)) { @@ -360,11 +349,10 @@ void ComputeRHEOSurface::unpack_reverse_comm(int n, int *list, double *buf) } } - /* ---------------------------------------------------------------------- */ -int ComputeRHEOSurface::pack_forward_comm(int n, int *list, double *buf, - int /*pbc_flag*/, int * /*pbc*/) +int ComputeRHEOSurface::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, + int * /*pbc*/) { int *status = atom->rheo_status; int m = 0; diff --git a/src/RHEO/compute_rheo_vshift.cpp b/src/RHEO/compute_rheo_vshift.cpp index 8d2f04d7f6..8581d9ea8f 100644 --- a/src/RHEO/compute_rheo_vshift.cpp +++ b/src/RHEO/compute_rheo_vshift.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -29,9 +28,9 @@ #include "fix_rheo.h" #include "force.h" #include "memory.h" -#include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" +#include "neighbor.h" #include "update.h" @@ -41,10 +40,10 @@ using namespace RHEO_NS; /* ---------------------------------------------------------------------- */ ComputeRHEOVShift::ComputeRHEOVShift(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), vshift(nullptr), fix_rheo(nullptr), rho0(nullptr), list(nullptr), - compute_interface(nullptr), compute_kernel(nullptr), compute_surface(nullptr) + Compute(lmp, narg, arg), vshift(nullptr), fix_rheo(nullptr), rho0(nullptr), list(nullptr), + compute_interface(nullptr), compute_kernel(nullptr), compute_surface(nullptr) { - if (narg != 3) error->all(FLERR,"Illegal compute RHEO/VShift command"); + if (narg != 3) error->all(FLERR, "Illegal compute RHEO/VShift command"); comm_reverse = 3; surface_flag = 0; @@ -125,8 +124,7 @@ void ComputeRHEOVShift::compute_peratom() } for (i = 0; i < nall; i++) - for (a = 0; a < dim; a++) - vshift[i][a] = 0.0; + for (a = 0; a < dim; a++) vshift[i][a] = 0.0; for (a = 0; a < 3; a++) { vi[a] = 0.0; @@ -141,8 +139,10 @@ void ComputeRHEOVShift::compute_peratom() itype = type[i]; jlist = firstneigh[i]; jnum = numneigh[i]; - if (rmass) imass = rmass[i]; - else imass = mass[itype]; + if (rmass) + imass = rmass[i]; + else + imass = mass[itype]; fluidi = !(status[i] & PHASECHECK); for (jj = 0; jj < jnum; jj++) { @@ -162,8 +162,10 @@ void ComputeRHEOVShift::compute_peratom() if (rsq < cutsq) { jtype = type[j]; - if (rmass) jmass = rmass[j]; - else jmass = mass[jtype]; + if (rmass) + jmass = rmass[j]; + else + jmass = mass[jtype]; r = sqrt(rsq); rinv = 1 / r; @@ -180,10 +182,10 @@ void ComputeRHEOVShift::compute_peratom() if (interface_flag) { if (fluidi && (!fluidj)) { compute_interface->correct_v(vj, vi, j, i); - rhoj = compute_interface->correct_rho(j, i); + rhoj = compute_interface->correct_rho(j); } else if ((!fluidi) && fluidj) { compute_interface->correct_v(vi, vj, i, j); - rhoi = compute_interface->correct_rho(i, j); + rhoi = compute_interface->correct_rho(i); } else if ((!fluidi) && (!fluidj)) { rhoi = rho0[itype]; rhoj = rho0[jtype]; @@ -195,7 +197,7 @@ void ComputeRHEOVShift::compute_peratom() wp = compute_kernel->calc_dw(i, j, dx[0], dx[1], dx[2], r); w = compute_kernel->calc_w(i, j, dx[0], dx[1], dx[2], r); - w0 = compute_kernel->calc_w(i, j, 0, 0, 0, cutthird); // dx, dy, dz irrelevant + w0 = compute_kernel->calc_w(i, j, 0, 0, 0, cutthird); // dx, dy, dz irrelevant w4 = w * w * w * w / (w0 * w0 * w0 * w0); dr = -2 * cutthird * (1 + 0.2 * w4) * wp * rinv; @@ -225,7 +227,6 @@ void ComputeRHEOVShift::compute_peratom() if (newton_pair) comm->reverse_comm(this); } - /* ---------------------------------------------------------------------- */ void ComputeRHEOVShift::correct_surfaces() @@ -300,7 +301,7 @@ int ComputeRHEOVShift::pack_reverse_comm(int n, int first, double *buf) void ComputeRHEOVShift::unpack_reverse_comm(int n, int *list, double *buf) { - int i,j,m; + int i, j, m; m = 0; for (i = 0; i < n; i++) { diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index b093eff681..4a25a46501 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -23,9 +22,9 @@ #include "citeme.h" #include "compute_rheo_grad.h" #include "compute_rheo_interface.h" -#include "compute_rheo_surface.h" #include "compute_rheo_kernel.h" #include "compute_rheo_rho_sum.h" +#include "compute_rheo_surface.h" #include "compute_rheo_vshift.h" #include "domain.h" #include "error.h" @@ -39,20 +38,23 @@ using namespace LAMMPS_NS; using namespace RHEO_NS; using namespace FixConst; +#if 0 +// publication was removed from documentation static const char cite_rheo[] = - "@article{PalermoInPrep,\n" - " journal = {in prep},\n" - " title = {RHEO: A Hybrid Mesh-Free Model Framework for Dynamic Multi-Phase Flows},\n" - " year = {2024},\n" - " author = {Eric T. Palermo and Ki T. Wolf and Joel T. Clemmer and Thomas C. O'Connor},\n" - "}\n\n"; + "@article{PalermoInPrep,\n" + " journal = {in prep},\n" + " title = {RHEO: A Hybrid Mesh-Free Model Framework for Dynamic Multi-Phase Flows},\n" + " year = {2024},\n" + " author = {Eric T. Palermo and Ki T. Wolf and Joel T. Clemmer and Thomas C. O'Connor},\n" + "}\n\n"; +#endif /* ---------------------------------------------------------------------- */ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), rho0(nullptr), csq(nullptr), compute_grad(nullptr), compute_kernel(nullptr), - compute_interface(nullptr), compute_surface(nullptr), compute_rhosum(nullptr), - compute_vshift(nullptr) + Fix(lmp, narg, arg), rho0(nullptr), csq(nullptr), compute_grad(nullptr), + compute_kernel(nullptr), compute_interface(nullptr), compute_surface(nullptr), + compute_rhosum(nullptr), compute_vshift(nullptr) { time_integrate = 1; @@ -78,44 +80,42 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : csq[i] = 1.0; } - if (igroup != 0) - error->all(FLERR, "fix rheo command requires group all"); + if (igroup != 0) error->all(FLERR, "fix rheo command requires group all"); if (atom->pressure_flag != 1) error->all(FLERR, "fix rheo command requires atom_style with pressure"); - if (atom->rho_flag != 1) - error->all(FLERR, "fix rheo command requires atom_style with density"); + if (atom->rho_flag != 1) error->all(FLERR, "fix rheo command requires atom_style with density"); if (atom->viscosity_flag != 1) error->all(FLERR, "fix rheo command requires atom_style with viscosity"); if (atom->rheo_status_flag != 1) error->all(FLERR, "fix rheo command requires atom_style with status"); - if (narg < 5) - utils::missing_cmd_args(FLERR, "fix rheo", error); + if (narg < 5) utils::missing_cmd_args(FLERR, "fix rheo", error); cut = utils::numeric(FLERR, arg[3], false, lmp); if (strcmp(arg[4], "quintic") == 0) { - kernel_style = QUINTIC; + kernel_style = QUINTIC; } else if (strcmp(arg[4], "wendland/c4") == 0) { - kernel_style = WENDLANDC4; + kernel_style = WENDLANDC4; } else if (strcmp(arg[4], "RK0") == 0) { - kernel_style = RK0; + kernel_style = RK0; } else if (strcmp(arg[4], "RK1") == 0) { - kernel_style = RK1; + kernel_style = RK1; } else if (strcmp(arg[4], "RK2") == 0) { - kernel_style = RK2; - } else error->all(FLERR, "Unknown kernel style {} in fix rheo", arg[4]); + kernel_style = RK2; + } else + error->all(FLERR, "Unknown kernel style {} in fix rheo", arg[4]); zmin_kernel = utils::numeric(FLERR, arg[5], false, lmp); int iarg = 6; - while (iarg < narg){ + while (iarg < narg) { if (strcmp(arg[iarg], "shift") == 0) { shift_flag = 1; } else if (strcmp(arg[iarg], "thermal") == 0) { thermal_flag = 1; } else if (strcmp(arg[iarg], "surface/detection") == 0) { surface_flag = 1; - if(iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo surface/detection", error); + if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo surface/detection", error); if (strcmp(arg[iarg + 1], "coordination") == 0) { surface_style = COORDINATION; zmin_surface = utils::inumeric(FLERR, arg[iarg + 2], false, lmp); @@ -137,8 +137,7 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : self_mass_flag = 1; } else if (strcmp(arg[iarg], "density") == 0) { if (iarg + n >= narg) utils::missing_cmd_args(FLERR, "fix rheo density", error); - for (i = 1; i <= n; i++) - rho0[i] = utils::numeric(FLERR, arg[iarg + i], false, lmp); + for (i = 1; i <= n; i++) rho0[i] = utils::numeric(FLERR, arg[iarg + i], false, lmp); iarg += n; } else if (strcmp(arg[iarg], "speed/sound") == 0) { if (iarg + n >= narg) utils::missing_cmd_args(FLERR, "fix rheo speed/sound", error); @@ -156,7 +155,9 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) : if (self_mass_flag && (!rhosum_flag)) error->all(FLERR, "Cannot use self/mass setting without rho/sum"); +#if 0 if (lmp->citeme) lmp->citeme->add(cite_rheo); +#endif } /* ---------------------------------------------------------------------- */ @@ -174,15 +175,14 @@ FixRHEO::~FixRHEO() memory->destroy(rho0); } - /* ---------------------------------------------------------------------- Create necessary internal computes ------------------------------------------------------------------------- */ void FixRHEO::post_constructor() { - compute_kernel = dynamic_cast(modify->add_compute( - fmt::format("rheo_kernel all RHEO/KERNEL {}", kernel_style))); + compute_kernel = dynamic_cast( + modify->add_compute(fmt::format("rheo_kernel all RHEO/KERNEL {}", kernel_style))); compute_kernel->fix_rheo = this; std::string cmd = "rheo_grad all RHEO/GRAD velocity rho viscosity"; @@ -191,26 +191,26 @@ void FixRHEO::post_constructor() compute_grad->fix_rheo = this; if (rhosum_flag) { - compute_rhosum = dynamic_cast(modify->add_compute( - fmt::format("rheo_rhosum all RHEO/RHO/SUM {}", self_mass_flag))); + compute_rhosum = dynamic_cast( + modify->add_compute(fmt::format("rheo_rhosum all RHEO/RHO/SUM {}", self_mass_flag))); compute_rhosum->fix_rheo = this; } if (shift_flag) { - compute_vshift = dynamic_cast(modify->add_compute( - "rheo_vshift all RHEO/VSHIFT")); + compute_vshift = + dynamic_cast(modify->add_compute("rheo_vshift all RHEO/VSHIFT")); compute_vshift->fix_rheo = this; } if (interface_flag) { - compute_interface = dynamic_cast(modify->add_compute( - "rheo_interface all RHEO/INTERFACE")); + compute_interface = dynamic_cast( + modify->add_compute("rheo_interface all RHEO/INTERFACE")); compute_interface->fix_rheo = this; } if (surface_flag) { - compute_surface = dynamic_cast(modify->add_compute( - "rheo_surface all RHEO/SURFACE")); + compute_surface = + dynamic_cast(modify->add_compute("rheo_surface all RHEO/SURFACE")); compute_surface->fix_rheo = this; } } @@ -237,23 +237,23 @@ void FixRHEO::init() error->all(FLERR, "Can only specify one instance of fix rheo"); if (atom->rheo_status_flag != 1) - error->all(FLERR,"fix rheo command requires atom property status"); - if (atom->rho_flag != 1) - error->all(FLERR,"fix rheo command requires atom property rho"); + error->all(FLERR, "fix rheo command requires atom property status"); + if (atom->rho_flag != 1) error->all(FLERR, "fix rheo command requires atom property rho"); if (atom->pressure_flag != 1) - error->all(FLERR,"fix rheo command requires atom property pressure"); + error->all(FLERR, "fix rheo command requires atom property pressure"); if (atom->viscosity_flag != 1) - error->all(FLERR,"fix rheo command requires atom property viscosity"); + error->all(FLERR, "fix rheo command requires atom property viscosity"); if (thermal_flag) { if (atom->esph_flag != 1) - error->all(FLERR,"fix rheo command requires atom property esph with thermal setting"); + error->all(FLERR, "fix rheo command requires atom property esph with thermal setting"); if (atom->temperature_flag != 1) - error->all(FLERR,"fix rheo command requires atom property temperature with thermal setting"); + error->all(FLERR, "fix rheo command requires atom property temperature with thermal setting"); if (atom->heatflow_flag != 1) - error->all(FLERR,"fix rheo command requires atom property heatflow with thermal setting"); + error->all(FLERR, "fix rheo command requires atom property heatflow with thermal setting"); if (atom->conductivity_flag != 1) - error->all(FLERR,"fix rheo command requires atom property conductivity with thermal setting"); + error->all(FLERR, + "fix rheo command requires atom property conductivity with thermal setting"); } } @@ -281,14 +281,11 @@ void FixRHEO::setup(int /*vflag*/) { // Confirm all accessory fixes are defined // Note: fixes set this flag in setup_pre_force() - if (!viscosity_fix_defined) - error->all(FLERR, "Missing fix rheo/viscosity"); + if (!viscosity_fix_defined) error->all(FLERR, "Missing fix rheo/viscosity"); - if (!pressure_fix_defined) - error->all(FLERR, "Missing fix rheo/pressure"); + if (!pressure_fix_defined) error->all(FLERR, "Missing fix rheo/pressure"); - if(thermal_flag && !thermal_fix_defined) - error->all(FLERR, "Missing fix rheo/thermal"); + if (thermal_flag && !thermal_fix_defined) error->all(FLERR, "Missing fix rheo/thermal"); // Reset to zero for future runs thermal_fix_defined = 0; @@ -296,8 +293,7 @@ void FixRHEO::setup(int /*vflag*/) pressure_fix_defined = 0; oxidation_fix_defined = 0; - if (rhosum_flag) - compute_rhosum->compute_peratom(); + if (rhosum_flag) compute_rhosum->compute_peratom(); } /* ---------------------------------------------------------------------- */ @@ -321,15 +317,13 @@ void FixRHEO::initial_integrate(int /*vflag*/) double **gradr = compute_grad->gradr; double **gradv = compute_grad->gradv; double **vshift; - if (shift_flag) - vshift = compute_vshift->vshift; + if (shift_flag) vshift = compute_vshift->vshift; int nlocal = atom->nlocal; int rmass_flag = atom->rmass_flag; int dim = domain->dimension; - if (igroup == atom->firstgroup) - nlocal = atom->nfirst; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; //Density Half-step for (i = 0; i < nlocal; i++) { @@ -349,7 +343,7 @@ void FixRHEO::initial_integrate(int /*vflag*/) } // Update gradients and interpolate solid properties - compute_grad->forward_fields(); // also forwards v and rho for chi + compute_grad->forward_fields(); // also forwards v and rho for chi if (interface_flag) { // Need to save, wiped in exchange compute_interface->store_forces(); @@ -360,9 +354,7 @@ void FixRHEO::initial_integrate(int /*vflag*/) // Position half-step for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - for (a = 0; a < dim; a++) { - x[i][a] += dtv * v[i][a]; - } + for (a = 0; a < dim; a++) { x[i][a] += dtv * v[i][a]; } } } @@ -374,9 +366,7 @@ void FixRHEO::initial_integrate(int /*vflag*/) if (status[i] & PHASECHECK) continue; divu = 0; - for (a = 0; a < dim; a++) { - divu += gradv[i][a * (1 + dim)]; - } + for (a = 0; a < dim; a++) { divu += gradv[i][a * (1 + dim)]; } rho[i] += dtf * (drho[i] - rho[i] * divu); } } @@ -392,16 +382,12 @@ void FixRHEO::initial_integrate(int /*vflag*/) if (mask[i] & groupbit) { for (a = 0; a < dim; a++) { x[i][a] += dtv * vshift[i][a]; - for (b = 0; b < dim; b++) { - v[i][a] += dtv * vshift[i][b] * gradv[i][a * dim + b]; - } + for (b = 0; b < dim; b++) { v[i][a] += dtv * vshift[i][b] * gradv[i][a * dim + b]; } } if (!rhosum_flag) { if (status[i] & PHASECHECK) continue; - for (a = 0; a < dim; a++) { - rho[i] += dtv * vshift[i][a] * gradr[i][a]; - } + for (a = 0; a < dim; a++) { rho[i] += dtv * vshift[i][a] * gradr[i][a]; } } } } @@ -412,10 +398,9 @@ void FixRHEO::initial_integrate(int /*vflag*/) void FixRHEO::pre_force(int /*vflag*/) { - compute_kernel->compute_coordination(); // Needed for rho sum + compute_kernel->compute_coordination(); // Needed for rho sum - if (rhosum_flag) - compute_rhosum->compute_peratom(); + if (rhosum_flag) compute_rhosum->compute_peratom(); compute_kernel->compute_peratom(); @@ -428,22 +413,19 @@ void FixRHEO::pre_force(int /*vflag*/) compute_grad->compute_peratom(); compute_grad->forward_gradients(); - if (shift_flag) - compute_vshift->compute_peratom(); + if (shift_flag) compute_vshift->compute_peratom(); // Remove temporary options int *mask = atom->mask; int *status = atom->rheo_status; int nall = atom->nlocal + atom->nghost; for (int i = 0; i < nall; i++) - if (mask[i] & groupbit) - status[i] &= OPTIONSMASK; + if (mask[i] & groupbit) status[i] &= OPTIONSMASK; // Calculate surfaces, update status if (surface_flag) { compute_surface->compute_peratom(); - if (shift_flag) - compute_vshift->correct_surfaces(); + if (shift_flag) compute_vshift->correct_surfaces(); } } @@ -452,8 +434,7 @@ void FixRHEO::pre_force(int /*vflag*/) void FixRHEO::final_integrate() { int nlocal = atom->nlocal; - if (igroup == atom->firstgroup) - nlocal = atom->nfirst; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; double dtfm, divu; int i, a; @@ -482,9 +463,7 @@ void FixRHEO::final_integrate() dtfm = dtf / mass[type[i]]; } - for (a = 0; a < dim; a++) { - v[i][a] += dtfm * f[i][a]; - } + for (a = 0; a < dim; a++) { v[i][a] += dtfm * f[i][a]; } } } @@ -496,9 +475,7 @@ void FixRHEO::final_integrate() if (status[i] & PHASECHECK) continue; divu = 0; - for (a = 0; a < dim; a++) { - divu += gradv[i][a * (1 + dim)]; - } + for (a = 0; a < dim; a++) { divu += gradv[i][a * (1 + dim)]; } rho[i] += dtf * (drho[i] - rho[i] * divu); } } diff --git a/src/RHEO/fix_rheo.h b/src/RHEO/fix_rheo.h index 8c62197fcd..24e86b21d6 100644 --- a/src/RHEO/fix_rheo.h +++ b/src/RHEO/fix_rheo.h @@ -72,11 +72,11 @@ class FixRHEO : public Fix { namespace RHEO_NS { - enum {QUINTIC, WENDLANDC4, RK0, RK1, RK2}; - enum {COORDINATION, DIVR}; + enum { QUINTIC, WENDLANDC4, RK0, RK1, RK2 }; + enum { COORDINATION, DIVR }; // Status variables - enum Status{ + enum Status { // Phase status STATUS_SOLID = 1 << 0, // Gap for future phase: STATUS_ = 1 << 1, @@ -95,12 +95,13 @@ namespace RHEO_NS { }; // Masks and their inverses - #define PHASEMASK 0xFFFFFFFC // 11111111111111111111111111111100 - #define PHASECHECK 0x00000003 // 00000000000000000000000000000011 - #define SURFACEMASK 0xFFFFFFC3 // 11111111111111111111111111000011 - #define SURFACECHECK 0x0000003C // 00000000000000000000000000111100 - #define OPTIONSMASK 0xFFFFFC3F // 11111111111111111111110000111111 - + enum { + PHASECHECK = 0x00000003, // 00000000000000000000000000000011 + SURFACECHECK = 0x0000003C, // 00000000000000000000000000111100 + OPTIONSMASK = 0xFFFFFC3F, // 11111111111111111111110000111111 + SURFACEMASK = 0xFFFFFFC3, // 11111111111111111111111111000011 + PHASEMASK = 0xFFFFFFFC // 11111111111111111111111111111100 + }; } // namespace RHEO_NS } // namespace LAMMPS_NS diff --git a/src/RHEO/fix_rheo_oxidation.cpp b/src/RHEO/fix_rheo_oxidation.cpp index 699f3428c8..604186ac1b 100644 --- a/src/RHEO/fix_rheo_oxidation.cpp +++ b/src/RHEO/fix_rheo_oxidation.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -28,34 +27,35 @@ #include "fix_rheo.h" #include "force.h" #include "modify.h" -#include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" +#include "neighbor.h" #include "update.h" using namespace LAMMPS_NS; using namespace RHEO_NS; using namespace FixConst; -enum {NONE, CONSTANT}; +enum { NONE, CONSTANT }; static const char cite_rheo_oxide[] = - "@article{ApplMathModel.130.310,\n" - " title = {A hybrid smoothed-particle hydrodynamics model of oxide skins on molten aluminum},\n" - " journal = {Applied Mathematical Modelling},\n" - " volume = {130},\n" - " pages = {310-326},\n" - " year = {2024},\n" - " issn = {0307-904X},\n" - " doi = {https://doi.org/10.1016/j.apm.2024.02.027},\n" - " author = {Joel T. Clemmer and Flint Pierce and Thomas C. O'Connor and Thomas D. Nevins and Elizabeth M.C. Jones and Jeremy B. Lechman and John Tencer},\n" - "}\n\n"; + "@article{ApplMathModel.130.310,\n" + " title = {A hybrid smoothed-particle hydrodynamics model of oxide skins on molten aluminum},\n" + " journal = {Applied Mathematical Modelling},\n" + " volume = {130},\n" + " pages = {310-326},\n" + " year = {2024},\n" + " issn = {0307-904X},\n" + " doi = {https://doi.org/10.1016/j.apm.2024.02.027},\n" + " author = {Joel T. Clemmer and Flint Pierce and Thomas C. O'Connor and Thomas D. Nevins and " + "Elizabeth M.C. Jones and Jeremy B. Lechman and John Tencer},\n" + "}\n\n"; /* ---------------------------------------------------------------------- */ FixRHEOOxidation::FixRHEOOxidation(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), compute_surface(nullptr), fix_rheo(nullptr) + Fix(lmp, narg, arg), compute_surface(nullptr), fix_rheo(nullptr) { - if (narg != 6) error->all(FLERR,"Illegal fix command"); + if (narg != 6) error->all(FLERR, "Illegal fix command"); force_reneighbor = 1; next_reneighbor = -1; @@ -65,7 +65,8 @@ FixRHEOOxidation::FixRHEOOxidation(LAMMPS *lmp, int narg, char **arg) : if (cut <= 0.0) error->all(FLERR, "Illegal bond cutoff {} in fix rheo/oxidation", cut); btype = utils::inumeric(FLERR, arg[4], false, lmp); - if (btype < 1 || btype > atom->nbondtypes) error->all(FLERR, "Illegal value {} for bond type in fix rheo/oxidation", btype); + if (btype < 1 || btype > atom->nbondtypes) + error->all(FLERR, "Illegal value {} for bond type in fix rheo/oxidation", btype); rsurf = utils::numeric(FLERR, arg[5], false, lmp); if (rsurf <= 0.0) error->all(FLERR, "Illegal surface distance {} in fix rheo/oxidation", cut); @@ -77,9 +78,7 @@ FixRHEOOxidation::FixRHEOOxidation(LAMMPS *lmp, int narg, char **arg) : /* ---------------------------------------------------------------------- */ -FixRHEOOxidation::~FixRHEOOxidation() -{ -} +FixRHEOOxidation::~FixRHEOOxidation() {} /* ---------------------------------------------------------------------- */ @@ -100,19 +99,16 @@ void FixRHEOOxidation::init() if (fixes.size() == 0) error->all(FLERR, "Need to define fix rheo to use fix rheo/oxidation"); fix_rheo = dynamic_cast(fixes[0]); - if (cut > fix_rheo->cut) - error->all(FLERR, "Bonding length exceeds kernel cutoff"); + if (cut > fix_rheo->cut) error->all(FLERR, "Bonding length exceeds kernel cutoff"); - if (rsurf >= fix_rheo->cut) - error->all(FLERR, "Surface distance must be less than kernel cutoff"); + if (rsurf >= fix_rheo->cut) error->all(FLERR, "Surface distance must be less than kernel cutoff"); if (!force->bond) error->all(FLERR, "Must define a bond style with fix rheo/oxidation"); if (!atom->avec->bonds_allow) error->all(FLERR, "Fix rheo/oxidation requires atom bonds"); int tmp1, tmp2; index_nb = atom->find_custom("shell_nbond", tmp1, tmp2); - if (index_nb == -1) - error->all(FLERR, "Must use bond style rheo/shell to use fix rheo/oxidation"); + if (index_nb == -1) error->all(FLERR, "Must use bond style rheo/shell to use fix rheo/oxidation"); nbond = atom->ivector[index_nb]; // need a half neighbor list @@ -135,16 +131,14 @@ void FixRHEOOxidation::setup_pre_force(int /*vflag*/) // but enforce to be consistent with other RHEO fixes fix_rheo->oxidation_fix_defined = 1; - if (!fix_rheo->surface_flag) error->all(FLERR, - "fix rheo/oxidation requires surface calculation in fix rheo"); + if (!fix_rheo->surface_flag) + error->all(FLERR, "fix rheo/oxidation requires surface calculation in fix rheo"); compute_surface = fix_rheo->compute_surface; } /* ---------------------------------------------------------------------- */ -void FixRHEOOxidation::pre_force(int /*vflag*/) -{ -} +void FixRHEOOxidation::pre_force(int /*vflag*/) {} /* ---------------------------------------------------------------------- */ @@ -235,7 +229,8 @@ void FixRHEOOxidation::post_integrate() if (!newton_bond || (tagi < tagj)) { if (num_bond[i] == atom->bond_per_atom) - error->one(FLERR,"New bond exceeded bonds per atom in fix rheo/oxidation for atom {}", tagi); + error->one(FLERR, "New bond exceeded bonds per atom in fix rheo/oxidation for atom {}", + tagi); bond_type[i][num_bond[i]] = btype; bond_atom[i][num_bond[i]] = tagj; num_bond[i]++; @@ -246,8 +241,7 @@ void FixRHEOOxidation::post_integrate() int added_bonds_all; MPI_Allreduce(&added_bonds, &added_bonds_all, 1, MPI_INT, MPI_SUM, world); - if (added_bonds_all > 0) - next_reneighbor = update->ntimestep; + if (added_bonds_all > 0) next_reneighbor = update->ntimestep; } /* ---------------------------------------------------------------------- */ @@ -257,14 +251,13 @@ void FixRHEOOxidation::post_force(int /*vflag*/) int *status = atom->rheo_status; int *num_bond = atom->num_bond; for (int i = 0; i < atom->nlocal; i++) - if (num_bond[i] != 0) - status[i] |= STATUS_NO_SHIFT; + if (num_bond[i] != 0) status[i] |= STATUS_NO_SHIFT; } /* ---------------------------------------------------------------------- */ -int FixRHEOOxidation::pack_forward_comm(int n, int *list, double *buf, - int /*pbc_flag*/, int * /*pbc*/) +int FixRHEOOxidation::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, + int * /*pbc*/) { int i, j, m; double **x = atom->x; diff --git a/src/RHEO/fix_rheo_pressure.cpp b/src/RHEO/fix_rheo_pressure.cpp index 84b0825570..10bf6ff2c9 100644 --- a/src/RHEO/fix_rheo_pressure.cpp +++ b/src/RHEO/fix_rheo_pressure.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -32,24 +31,23 @@ using namespace LAMMPS_NS; using namespace FixConst; -enum {NONE, LINEAR, CUBIC, TAITWATER, TAITGENERAL}; +enum { NONE, LINEAR, CUBIC, TAITWATER, TAITGENERAL }; static constexpr double SEVENTH = 1.0 / 7.0; /* ---------------------------------------------------------------------- */ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), c_cubic(nullptr), csq(nullptr), csqinv(nullptr), rho0(nullptr), - rho0inv(nullptr), tpower(nullptr), pbackground(nullptr), pressure_style(nullptr), - fix_rheo(nullptr) + Fix(lmp, narg, arg), c_cubic(nullptr), csq(nullptr), csqinv(nullptr), rho0(nullptr), + rho0inv(nullptr), tpower(nullptr), pbackground(nullptr), pressure_style(nullptr), + fix_rheo(nullptr) { - if (narg < 4) error->all(FLERR,"Illegal fix command"); + if (narg < 4) error->all(FLERR, "Illegal fix command"); comm_forward = 1; // Currently can only have one instance of fix rheo/pressure - if (igroup != 0) - error->all(FLERR,"fix rheo/pressure command requires group all"); + if (igroup != 0) error->all(FLERR, "fix rheo/pressure command requires group all"); int i, nlo, nhi; int n = atom->ntypes; @@ -66,11 +64,9 @@ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) : if (iarg + 1 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/pressure", error); if (strcmp(arg[iarg + 1], "linear") == 0) { - for (i = nlo; i <= nhi; i++) - pressure_style[i] = LINEAR; + for (i = nlo; i <= nhi; i++) pressure_style[i] = LINEAR; } else if (strcmp(arg[iarg + 1], "tait/water") == 0) { - for (i = nlo; i <= nhi; i++) - pressure_style[i] = TAITWATER; + for (i = nlo; i <= nhi; i++) pressure_style[i] = TAITWATER; } else if (strcmp(arg[iarg + 1], "tait/general") == 0) { if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/pressure tait", error); @@ -94,14 +90,14 @@ FixRHEOPressure::FixRHEOPressure(LAMMPS *lmp, int narg, char **arg) : c_cubic[i] = c_cubic_one; } } else { - error->all(FLERR,"Illegal fix command, {}", arg[iarg]); + error->all(FLERR, "Illegal fix command, {}", arg[iarg]); } iarg += 2; } for (i = 1; i <= n; i++) if (pressure_style[i] == NONE) - error->all(FLERR,"Must specify pressure for atom type {} in fix/rheo/pressure", i); + error->all(FLERR, "Must specify pressure for atom type {} in fix/rheo/pressure", i); } /* ---------------------------------------------------------------------- */ @@ -170,16 +166,15 @@ void FixRHEOPressure::pre_force(int /*vflag*/) int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) - pressure[i] = calc_pressure(rho[i], type[i]); + if (mask[i] & groupbit) pressure[i] = calc_pressure(rho[i], type[i]); if (comm_forward) comm->forward_comm(this); } /* ---------------------------------------------------------------------- */ -int FixRHEOPressure::pack_forward_comm(int n, int *list, double *buf, - int /*pbc_flag*/, int * /*pbc*/) +int FixRHEOPressure::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, + int * /*pbc*/) { double *pressure = atom->pressure; int m = 0; @@ -197,9 +192,7 @@ void FixRHEOPressure::unpack_forward_comm(int n, int first, double *buf) double *pressure = atom->pressure; int m = 0; int last = first + n; - for (int i = first; i < last; i++) { - pressure[i] = buf[m++]; - } + for (int i = first; i < last; i++) { pressure[i] = buf[m++]; } } /* ---------------------------------------------------------------------- */ @@ -235,7 +228,8 @@ double FixRHEOPressure::calc_rho(double p, int type) if (pressure_style[type] == LINEAR) { rho = csqinv[type] * p + rho0[type]; } else if (pressure_style[type] == CUBIC) { - error->one(FLERR, "Rho calculation from pressure not yet supported for cubic pressure equation"); + error->one(FLERR, + "Rho calculation from pressure not yet supported for cubic pressure equation"); } else if (pressure_style[type] == TAITWATER) { rho = pow(7.0 * p + csq[type] * rho0[type], SEVENTH); rho *= pow(rho0[type], 6.0 * SEVENTH); diff --git a/src/RHEO/fix_rheo_thermal.cpp b/src/RHEO/fix_rheo_thermal.cpp index a7d9fc8df9..f8449ead70 100644 --- a/src/RHEO/fix_rheo_thermal.cpp +++ b/src/RHEO/fix_rheo_thermal.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -27,44 +26,45 @@ #include "compute_rheo_vshift.h" #include "domain.h" #include "error.h" -#include "fix_rheo.h" #include "fix_bond_history.h" +#include "fix_rheo.h" #include "fix_update_special_bonds.h" #include "force.h" #include "math_extra.h" #include "memory.h" #include "modify.h" -#include "neighbor.h" #include "neigh_list.h" #include "neigh_request.h" +#include "neighbor.h" #include "pair.h" #include "update.h" using namespace LAMMPS_NS; using namespace RHEO_NS; using namespace FixConst; -enum {NONE, CONSTANT}; +enum { NONE, CONSTANT }; static const char cite_rheo_oxide[] = - "@article{ApplMathModel.130.310,\n" - " title = {A hybrid smoothed-particle hydrodynamics model of oxide skins on molten aluminum},\n" - " journal = {Applied Mathematical Modelling},\n" - " volume = {130},\n" - " pages = {310-326},\n" - " year = {2024},\n" - " issn = {0307-904X},\n" - " doi = {https://doi.org/10.1016/j.apm.2024.02.027},\n" - " author = {Joel T. Clemmer and Flint Pierce and Thomas C. O'Connor and Thomas D. Nevins and Elizabeth M.C. Jones and Jeremy B. Lechman and John Tencer},\n" - "}\n\n"; + "@article{ApplMathModel.130.310,\n" + " title = {A hybrid smoothed-particle hydrodynamics model of oxide skins on molten aluminum},\n" + " journal = {Applied Mathematical Modelling},\n" + " volume = {130},\n" + " pages = {310-326},\n" + " year = {2024},\n" + " issn = {0307-904X},\n" + " doi = {https://doi.org/10.1016/j.apm.2024.02.027},\n" + " author = {Joel T. Clemmer and Flint Pierce and Thomas C. O'Connor and Thomas D. Nevins and " + "Elizabeth M.C. Jones and Jeremy B. Lechman and John Tencer},\n" + "}\n\n"; /* ---------------------------------------------------------------------- */ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), cv(nullptr), Tc(nullptr), kappa(nullptr), L(nullptr), cv_style(nullptr), - Tc_style(nullptr), kappa_style(nullptr), L_style(nullptr), list(nullptr), fix_rheo(nullptr), - compute_grad(nullptr), compute_vshift(nullptr), fix_update_special_bonds(nullptr) + Fix(lmp, narg, arg), cv(nullptr), Tc(nullptr), kappa(nullptr), L(nullptr), cv_style(nullptr), + Tc_style(nullptr), kappa_style(nullptr), L_style(nullptr), list(nullptr), fix_rheo(nullptr), + compute_grad(nullptr), compute_vshift(nullptr), fix_update_special_bonds(nullptr) { - if (narg < 4) error->all(FLERR,"Illegal fix command"); + if (narg < 4) error->all(FLERR, "Illegal fix command"); force_reneighbor = 1; next_reneighbor = -1; @@ -72,8 +72,7 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : comm_forward = 0; // Currently can only have one instance of fix rheo/thermal - if (igroup != 0) - error->all(FLERR,"fix rheo/thermal command requires group all"); + if (igroup != 0) error->all(FLERR, "fix rheo/thermal command requires group all"); int i, nlo, nhi; int n = atom->ntypes; @@ -97,13 +96,14 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : int iarg = 3; while (iarg < narg) { - if (strcmp(arg[iarg],"conductivity") == 0) { + if (strcmp(arg[iarg], "conductivity") == 0) { if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal conductivity", error); utils::bounds(FLERR, arg[iarg + 1], 1, n, nlo, nhi, error); // Conductivity arguments if (strcmp(arg[iarg + 2], "constant") == 0) { - if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal conductivity constant", error); + if (iarg + 3 >= narg) + utils::missing_cmd_args(FLERR, "fix rheo/thermal conductivity constant", error); double kappa_one = utils::numeric(FLERR, arg[iarg + 3], false, lmp); if (kappa_one < 0.0) error->all(FLERR, "The conductivity must be positive"); @@ -114,7 +114,7 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : kappa[i] = kappa_one; } } else { - error->all(FLERR,"Illegal fix command, {}", arg[iarg + 2]); + error->all(FLERR, "Illegal fix command, {}", arg[iarg + 2]); } iarg += 2; @@ -124,7 +124,8 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : // Cv arguments if (strcmp(arg[iarg + 2], "constant") == 0) { - if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal specific/heat constant", error); + if (iarg + 3 >= narg) + utils::missing_cmd_args(FLERR, "fix rheo/thermal specific/heat constant", error); double cv_one = utils::numeric(FLERR, arg[iarg + 3], false, lmp); if (cv_one < 0.0) error->all(FLERR, "The specific heat must be positive"); @@ -136,7 +137,7 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : } } else { - error->all(FLERR,"Illegal fix command, {}", arg[iarg + 2]); + error->all(FLERR, "Illegal fix command, {}", arg[iarg + 2]); } iarg += 2; @@ -146,7 +147,8 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : // T freeze arguments if (strcmp(arg[iarg + 2], "constant") == 0) { - if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal Tfreeze constant", error); + if (iarg + 3 >= narg) + utils::missing_cmd_args(FLERR, "fix rheo/thermal Tfreeze constant", error); double Tc_one = utils::numeric(FLERR, arg[iarg + 3], false, lmp); iarg += 2; @@ -167,7 +169,8 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : // Cv arguments if (strcmp(arg[iarg + 2], "constant") == 0) { - if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo/thermal latent/heat constant", error); + if (iarg + 3 >= narg) + utils::missing_cmd_args(FLERR, "fix rheo/thermal latent/heat constant", error); double L_one = utils::numeric(FLERR, arg[iarg + 3], false, lmp); if (L_one < 0.0) error->all(FLERR, "The latent heat must be positive"); @@ -179,7 +182,7 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : } } else { - error->all(FLERR,"Illegal fix command, {}", arg[iarg + 2]); + error->all(FLERR, "Illegal fix command, {}", arg[iarg + 2]); } iarg += 2; @@ -195,18 +198,20 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) : cutsq_bond = cut_bond * cut_bond; iarg += 3; } else { - error->all(FLERR,"Unknown fix rheo/thermal keyword: {}", arg[iarg]); + error->all(FLERR, "Unknown fix rheo/thermal keyword: {}", arg[iarg]); } } - for (i = 1; i <= n; i++) { if (cv_style[i] == NONE) error->all(FLERR, "Must specify specific/heat for atom type {} in fix/rheo/thermal", i); if (kappa_style[i] == NONE) error->all(FLERR, "Must specify conductivity for atom type {} in fix/rheo/thermal", i); if (Tc_style[i] == NONE && L_style[i] != NONE) - error->all(FLERR, "Must specify critical temperature for atom type {} to use latent heat in fix rheo/thermal", i); + error->all(FLERR, + "Must specify critical temperature for atom type {} to use latent heat in fix " + "rheo/thermal", + i); } if (lmp->citeme) lmp->citeme->add(cite_rheo_oxide); @@ -253,11 +258,9 @@ void FixRHEOThermal::init() fix_rheo = dynamic_cast(fixes[0]); cut_kernel = fix_rheo->cut; - if (cut_bond > cut_kernel) - error->all(FLERR, "Bonding length exceeds kernel cutoff"); + if (cut_bond > cut_kernel) error->all(FLERR, "Bonding length exceeds kernel cutoff"); - if (!fix_rheo->thermal_flag) - error->all(FLERR, "Need to define thermal setting in fix rheo"); + if (!fix_rheo->thermal_flag) error->all(FLERR, "Need to define thermal setting in fix rheo"); compute_grad = fix_rheo->compute_grad; compute_vshift = fix_rheo->compute_vshift; @@ -274,19 +277,25 @@ void FixRHEOThermal::init() error->all(FLERR, "fix rheo/thermal command requires atom property conductivity"); if (cut_bond > 0.0) { - if (!force->bond) error->all(FLERR, "Must define a bond style to use reactive bond generation with fix rheo/thermal"); - if (!atom->avec->bonds_allow) error->all(FLERR, "Reactive bond generation in fix rheo/thermal requires atom bonds"); + if (!force->bond) + error->all(FLERR, + "Must define a bond style to use reactive bond generation with fix rheo/thermal"); + if (!atom->avec->bonds_allow) + error->all(FLERR, "Reactive bond generation in fix rheo/thermal requires atom bonds"); // all special weights must be 1.0 (no special neighbors) or there must be an instance of fix update/special/bonds - if (force->special_lj[0] != 1.0 || force->special_lj[1] != 1.0 || force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0) { + if (force->special_lj[0] != 1.0 || force->special_lj[1] != 1.0 || force->special_lj[2] != 1.0 || + force->special_lj[3] != 1.0) { auto fixes = modify->get_fix_by_style("UPDATE_SPECIAL_BONDS"); - if (fixes.size() == 0) error->all(FLERR, "Without fix update/special/bonds, reactive bond generation in fix rheo/thermal requires special weights of 1.0"); + if (fixes.size() == 0) + error->all(FLERR, + "Without fix update/special/bonds, reactive bond generation in fix rheo/thermal " + "requires special weights of 1.0"); fix_update_special_bonds = dynamic_cast(fixes[0]); } // must have newton off so both processors will search nlist to build bonds - if (force->newton_pair) - error->all(FLERR, "Need Newton off for reactive bond generation"); + if (force->newton_pair) error->all(FLERR, "Need Newton off for reactive bond generation"); // need a half neighbor list, built only when particles freeze auto req = neighbor->add_request(this, NeighConst::REQ_OCCASIONAL); @@ -334,13 +343,11 @@ void FixRHEOThermal::initial_integrate(int /*vflag*/) int nlocal = atom->nlocal; int dim = domain->dimension; - if (igroup == atom->firstgroup) - nlocal = atom->nfirst; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; for (i = 0; i < nlocal; i++) { if (status[i] & STATUS_NO_SHIFT) continue; - for (a = 0; a < dim; a++) - energy[i] += dt * vshift[i][a] * grade[i][a]; + for (a = 0; a < dim; a++) energy[i] += dt * vshift[i][a] * grade[i][a]; } } @@ -365,16 +372,16 @@ void FixRHEOThermal::post_integrate() if (status[i] & STATUS_NO_INTEGRATION) continue; itype = type[i]; - cvi = calc_cv(i, itype); + cvi = calc_cv(itype); energy[i] += dth * heatflow[i]; temperature[i] = energy[i] / cvi; if (Tc_style[itype] != NONE) { Ti = temperature[i]; - Tci = calc_Tc(i, itype); + Tci = calc_Tc(itype); if (L_style[itype] != NONE) { - Li = calc_L(i, itype); + Li = calc_L(itype); if (Ti > Tci) Ti = MAX(Tci, (energy[i] - Li) / cvi); temperature[i] = Ti; } @@ -448,8 +455,7 @@ void FixRHEOThermal::post_neighbor() for (i = 0; i < nall; i++) { itype = type[i]; - if (kappa_style[itype] == CONSTANT) - conductivity[i] = kappa[itype]; + if (kappa_style[itype] == CONSTANT) conductivity[i] = kappa[itype]; } } @@ -470,15 +476,15 @@ void FixRHEOThermal::pre_force(int /*vflag*/) // Calculate temperature for (int i = 0; i < nall; i++) { int itype = type[i]; - cvi = calc_cv(i, itype); + cvi = calc_cv(itype); temperature[i] = energy[i] / cvi; if (Tc_style[itype] != NONE) { Ti = temperature[i]; - Tci = calc_Tc(i, itype); + Tci = calc_Tc(itype); if (L_style[itype] != NONE) { - Li = calc_L(i, itype); + Li = calc_L(itype); if (Ti > Tci) Ti = MAX(Tci, (energy[i] - Li) / cvi); temperature[i] = Ti; } @@ -539,14 +545,14 @@ void FixRHEOThermal::break_bonds() nmax = num_bond[i] - 1; if (m == nmax) { if (n_histories > 0) - for (auto &ihistory: histories) + for (auto &ihistory : histories) dynamic_cast(ihistory)->delete_history(i, m); } else { bond_type[i][m] = bond_type[i][nmax]; bond_atom[i][m] = bond_atom[i][nmax]; if (n_histories > 0) { - for (auto &ihistory: histories) { - auto fix_bond_history = dynamic_cast (ihistory); + for (auto &ihistory : histories) { + auto fix_bond_history = dynamic_cast(ihistory); fix_bond_history->shift_history(i, m, nmax); fix_bond_history->delete_history(i, nmax); } @@ -559,9 +565,7 @@ void FixRHEOThermal::break_bonds() // only update for atom with lower tag if (fix_update_special_bonds) { if ((i < nlocal) && (j < nlocal) && melti && meltj) { - if (tag[i] < tag[j]) { - fix_update_special_bonds->add_broken_bond(i, j); - } + if (tag[i] < tag[j]) { fix_update_special_bonds->add_broken_bond(i, j); } } else { fix_update_special_bonds->add_broken_bond(i, j); } @@ -592,8 +596,8 @@ void FixRHEOThermal::break_bonds() bond_type[i][m] = bond_type[i][nmax]; bond_atom[i][m] = bond_atom[i][nmax]; if (n_histories > 0) - for (auto &ihistory: histories) { - auto fix_bond_history = dynamic_cast (ihistory); + for (auto &ihistory : histories) { + auto fix_bond_history = dynamic_cast(ihistory); fix_bond_history->shift_history(i, m, nmax); fix_bond_history->delete_history(i, nmax); } @@ -611,8 +615,8 @@ void FixRHEOThermal::break_bonds() bond_type[j][m] = bond_type[j][nmax]; bond_atom[j][m] = bond_atom[j][nmax]; if (n_histories > 0) - for (auto &ihistory: histories) { - auto fix_bond_history = dynamic_cast (ihistory); + for (auto &ihistory : histories) { + auto fix_bond_history = dynamic_cast(ihistory); fix_bond_history->shift_history(j, m, nmax); fix_bond_history->delete_history(j, nmax); } @@ -691,7 +695,7 @@ void FixRHEOThermal::create_bonds() // If newton bond off, add to both, otherwise add to whichever has a smaller tag if ((i < nlocal) && (!newton_bond || (tag[i] < tag[j]))) { if (num_bond[i] == atom->bond_per_atom) - error->one(FLERR,"New bond exceeded bonds per atom in fix rheo/thermal"); + error->one(FLERR, "New bond exceeded bonds per atom in fix rheo/thermal"); bond_type[i][num_bond[i]] = btype; bond_atom[i][num_bond[i]] = tag[j]; num_bond[i]++; @@ -699,7 +703,7 @@ void FixRHEOThermal::create_bonds() if ((j < nlocal) && (!newton_bond || (tag[j] < tag[i]))) { if (num_bond[j] == atom->bond_per_atom) - error->one(FLERR,"New bond exceeded bonds per atom in fix rheo/thermal"); + error->one(FLERR, "New bond exceeded bonds per atom in fix rheo/thermal"); bond_type[j][num_bond[j]] = btype; bond_atom[j][num_bond[j]] = tag[i]; num_bond[j]++; @@ -712,41 +716,35 @@ void FixRHEOThermal::create_bonds() /* ---------------------------------------------------------------------- */ -double FixRHEOThermal::calc_cv(int i, int itype) +double FixRHEOThermal::calc_cv(int itype) { - if (cv_style[itype] == CONSTANT) { - return cv[itype]; - } + if (cv_style[itype] == CONSTANT) { return cv[itype]; } return 0.0; } /* ---------------------------------------------------------------------- */ -double FixRHEOThermal::calc_Tc(int i, int itype) +double FixRHEOThermal::calc_Tc(int itype) { - if (Tc_style[itype] == CONSTANT) { - return Tc[itype]; - } + if (Tc_style[itype] == CONSTANT) { return Tc[itype]; } return 0.0; } /* ---------------------------------------------------------------------- */ -double FixRHEOThermal::calc_L(int i, int itype) +double FixRHEOThermal::calc_L(int itype) { - if (L_style[itype] == CONSTANT) { - return L[itype]; - } + if (L_style[itype] == CONSTANT) { return L[itype]; } return 0.0; } /* ---------------------------------------------------------------------- */ -int FixRHEOThermal::pack_forward_comm(int n, int *list, double *buf, - int /*pbc_flag*/, int * /*pbc*/) +int FixRHEOThermal::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, + int * /*pbc*/) { int *status = atom->rheo_status; double **x = atom->x; diff --git a/src/RHEO/fix_rheo_thermal.h b/src/RHEO/fix_rheo_thermal.h index 264bbdbce3..a957d59fc8 100644 --- a/src/RHEO/fix_rheo_thermal.h +++ b/src/RHEO/fix_rheo_thermal.h @@ -42,9 +42,9 @@ class FixRHEOThermal : public Fix { int pack_forward_comm(int, int *, double *, int, int *) override; void unpack_forward_comm(int, int, double *) override; void reset_dt() override; - double calc_cv(int, int); - double calc_Tc(int, int); - double calc_L(int, int); + double calc_cv(int); + double calc_Tc(int); + double calc_L(int); private: double *cv, *Tc, *kappa, *L; diff --git a/src/RHEO/fix_rheo_viscosity.cpp b/src/RHEO/fix_rheo_viscosity.cpp index 5f7a1208c5..88a4929ed8 100644 --- a/src/RHEO/fix_rheo_viscosity.cpp +++ b/src/RHEO/fix_rheo_viscosity.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -33,13 +32,13 @@ using namespace LAMMPS_NS; using namespace FixConst; -enum {NONE, CONSTANT, POWER}; +enum { NONE, CONSTANT, POWER }; /* ---------------------------------------------------------------------- */ FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), eta(nullptr), npow(nullptr), K(nullptr), gd0(nullptr), tau0(nullptr), - viscosity_style(nullptr), fix_rheo(nullptr), compute_grad(nullptr) + Fix(lmp, narg, arg), eta(nullptr), npow(nullptr), K(nullptr), gd0(nullptr), tau0(nullptr), + viscosity_style(nullptr), fix_rheo(nullptr), compute_grad(nullptr) { if (narg < 4) error->all(FLERR, "Illegal fix command"); @@ -48,8 +47,7 @@ FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) : evolve_flag = 0; // Currently can only have one instance of fix rheo/viscosity - if (igroup != 0) - error->all(FLERR,"fix rheo/viscosity command requires group all"); + if (igroup != 0) error->all(FLERR, "fix rheo/viscosity command requires group all"); int i, nlo, nhi; int n = atom->ntypes; @@ -107,7 +105,7 @@ FixRHEOViscosity::FixRHEOViscosity(LAMMPS *lmp, int narg, char **arg) : for (i = 1; i <= n; i++) if (viscosity_style[i] == NONE) - error->all(FLERR,"Must specify viscosity for atom type {} in fix/rheo/viscosity", i); + error->all(FLERR, "Must specify viscosity for atom type {} in fix/rheo/viscosity", i); } /* ---------------------------------------------------------------------- */ @@ -172,8 +170,7 @@ void FixRHEOViscosity::post_neighbor() for (i = 0; i < nall; i++) { itype = type[i]; - if (viscosity_style[itype]) - viscosity[i] = eta[itype]; + if (viscosity_style[itype]) viscosity[i] = eta[itype]; } } @@ -195,7 +192,6 @@ void FixRHEOViscosity::pre_force(int /*vflag*/) int nlocal = atom->nlocal; int dim = domain->dimension; - for (i = 0; i < nlocal; i++) { itype = type[i]; if (viscosity_style[itype] == POWER) { @@ -222,8 +218,8 @@ void FixRHEOViscosity::pre_force(int /*vflag*/) /* ---------------------------------------------------------------------- */ -int FixRHEOViscosity::pack_forward_comm(int n, int *list, double *buf, - int /*pbc_flag*/, int * /*pbc*/) +int FixRHEOViscosity::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, + int * /*pbc*/) { double *viscosity = atom->viscosity; int m = 0; @@ -241,7 +237,5 @@ void FixRHEOViscosity::unpack_forward_comm(int n, int first, double *buf) double *viscosity = atom->viscosity; int m = 0; int last = first + n; - for (int i = first; i < last; i++) { - viscosity[i] = buf[m++]; - } + for (int i = first; i < last; i++) { viscosity[i] = buf[m++]; } } diff --git a/src/RHEO/pair_rheo.cpp b/src/RHEO/pair_rheo.cpp index 8feac7ddcd..aa0c685aaa 100644 --- a/src/RHEO/pair_rheo.cpp +++ b/src/RHEO/pair_rheo.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -21,9 +20,9 @@ #include "atom.h" #include "comm.h" -#include "compute_rheo_kernel.h" #include "compute_rheo_grad.h" #include "compute_rheo_interface.h" +#include "compute_rheo_kernel.h" #include "domain.h" #include "error.h" #include "fix_rheo.h" @@ -32,8 +31,8 @@ #include "math_extra.h" #include "memory.h" #include "modify.h" -#include "neighbor.h" #include "neigh_list.h" +#include "neighbor.h" #include "utils.h" #include @@ -47,9 +46,8 @@ static constexpr double EPSILON = 1e-2; /* ---------------------------------------------------------------------- */ PairRHEO::PairRHEO(LAMMPS *lmp) : - Pair(lmp), csq(nullptr), rho0(nullptr), cs(nullptr), compute_kernel(nullptr), - compute_grad(nullptr), compute_interface(nullptr), fix_rheo(nullptr), - fix_pressure(nullptr) + Pair(lmp), csq(nullptr), rho0(nullptr), cs(nullptr), compute_kernel(nullptr), + compute_grad(nullptr), compute_interface(nullptr), fix_rheo(nullptr), fix_pressure(nullptr) { restartinfo = 0; single_enable = 0; @@ -82,7 +80,8 @@ void PairRHEO::compute(int eflag, int vflag) int pair_force_flag, pair_rho_flag, pair_avisc_flag; int fluidi, fluidj; double xtmp, ytmp, ztmp, wp, Ti, Tj, dT, csq_ave, cs_ave; - double rhoi, rhoj, rho0i, rho0j, voli, volj, Pi, Pj, etai, etaj, kappai, kappaj, eta_ave, kappa_ave, dT_prefactor; + double rhoi, rhoj, rho0i, rho0j, voli, volj, Pi, Pj, etai, etaj, kappai, kappaj, eta_ave, + kappa_ave, dT_prefactor; double mu, q, fp_prefactor, drho_damp, fmag, psi_ij, Fij; double *dWij, *dWji; double dx[3], du[3], dv[3], fv[3], dfp[3], fsolid[3], ft[3], vi[3], vj[3]; @@ -150,8 +149,10 @@ void PairRHEO::compute(int eflag, int vflag) itype = type[i]; jlist = firstneigh[i]; jnum = numneigh[i]; - if (rmass) imass = rmass[i]; - else imass = mass[itype]; + if (rmass) + imass = rmass[i]; + else + imass = mass[itype]; etai = viscosity[i]; fluidi = !(status[i] & PHASECHECK); if (thermal_flag) { @@ -173,8 +174,10 @@ void PairRHEO::compute(int eflag, int vflag) r = sqrt(rsq); rinv = 1 / r; - if (rmass) jmass = rmass[i]; - else jmass = mass[jtype]; + if (rmass) + jmass = rmass[i]; + else + jmass = mass[jtype]; etaj = viscosity[j]; fluidj = !(status[j] & PHASECHECK); if (thermal_flag) { @@ -217,7 +220,7 @@ void PairRHEO::compute(int eflag, int vflag) if (interface_flag) { if (fluidi && (!fluidj)) { compute_interface->correct_v(vj, vi, j, i); - rhoj = compute_interface->correct_rho(j, i); + rhoj = compute_interface->correct_rho(j); Pj = fix_pressure->calc_pressure(rhoj, jtype); if ((chi[j] > 0.9) && (r < (cutk * 0.5))) @@ -225,7 +228,7 @@ void PairRHEO::compute(int eflag, int vflag) } else if ((!fluidi) && fluidj) { compute_interface->correct_v(vi, vj, i, j); - rhoi = compute_interface->correct_rho(i, j); + rhoi = compute_interface->correct_rho(i); Pi = fix_pressure->calc_pressure(rhoi, itype); if (chi[i] > 0.9 && r < (cutk * 0.5)) @@ -251,7 +254,8 @@ void PairRHEO::compute(int eflag, int vflag) } else { kappa_ave = 0.5 * (kappai + kappaj); } - dT_prefactor = 2.0 * kappa_ave * (Ti - Tj) * rinv * rinv * voli * volj * 2.0 / (rhoi + rhoj); + dT_prefactor = + 2.0 * kappa_ave * (Ti - Tj) * rinv * rinv * voli * volj * 2.0 / (rhoi + rhoj); dT = dot3(dx, dWij); heatflow[i] += dT * dT_prefactor; @@ -295,8 +299,7 @@ void PairRHEO::compute(int eflag, int vflag) // Now compute viscous eta*Lap[v] terms for (a = 0; a < dim; a++) { fv[a] = 0.0; - for (b = 0; b < dim; b++) - fv[a] += dv[a] * dx[b] * dWij[b]; + for (b = 0; b < dim; b++) fv[a] += dv[a] * dx[b] * dWij[b]; fv[a] *= 2.0 * eta_ave * voli * volj * rinv * rinv; } @@ -309,13 +312,13 @@ void PairRHEO::compute(int eflag, int vflag) // Note the virial's definition is hazy, e.g. viscous contributions will depend on rotation if (evflag) - ev_tally_xyz(i, j, nlocal, newton_pair, 0.0, 0.0, ft[0], ft[1], ft[2], dx[0], dx[1], dx[2]); + ev_tally_xyz(i, j, nlocal, newton_pair, 0.0, 0.0, ft[0], ft[1], ft[2], dx[0], dx[1], + dx[2]); if (newton_pair || j < nlocal) { - for (a = 0; a < dim; a ++) { + for (a = 0; a < dim; a++) { fv[a] = 0.0; - for (b = 0; b < dim; b++) - fv[a] += (vi[a] - vj[a]) * dx[b] * dWji[b]; + for (b = 0; b < dim; b++) fv[a] += (vi[a] - vj[a]) * dx[b] * dWji[b]; fv[a] *= -2.0 * eta_ave * voli * volj * rinv * rinv; // flip sign here b/c -= at accummulator } @@ -329,7 +332,8 @@ void PairRHEO::compute(int eflag, int vflag) f[j][2] -= ft[2]; if (evflag) - ev_tally_xyz(i, j, nlocal, newton_pair, 0.0, 0.0, ft[0], ft[1], ft[2], -dx[0], -dx[1], -dx[2]); + ev_tally_xyz(i, j, nlocal, newton_pair, 0.0, 0.0, ft[0], ft[1], ft[2], -dx[0], -dx[1], + -dx[2]); } if (compute_interface) { @@ -352,8 +356,7 @@ void PairRHEO::compute(int eflag, int vflag) if (laplacian_order >= 1) { psi_ij = rhoj - rhoi; Fij = -rinv * rinv * dot3(dx, dWij); - for (a = 0; a < dim; a++) - psi_ij += 0.5 * (gradr[i][a] + gradr[j][a]) * dx[a]; + for (a = 0; a < dim; a++) psi_ij += 0.5 * (gradr[i][a] + gradr[j][a]) * dx[a]; drho[i] += 2 * rho_damp * psi_ij * Fij * volj; } else { drho_damp = 2 * rho_damp * ((rhoj - rho0[jtype]) - (rhoi - rho0[itype])) * rinv * wp; @@ -393,8 +396,7 @@ void PairRHEO::allocate() memory->create(setflag, n + 1, n + 1, "pair:setflag"); for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - setflag[i][j] = 0; + for (int j = i; j <= n; j++) setflag[i][j] = 0; memory->create(cutsq, n + 1, n + 1, "pair:cutsq"); } @@ -423,7 +425,8 @@ void PairRHEO::settings(int narg, char **arg) iarg++; } else if (strcmp(arg[iarg], "harmonic/means") == 0) { harmonic_means_flag = 1; - } else error->all(FLERR, "Illegal pair_style command, {}", arg[iarg]); + } else + error->all(FLERR, "Illegal pair_style command, {}", arg[iarg]); iarg++; } } @@ -434,10 +437,8 @@ void PairRHEO::settings(int narg, char **arg) void PairRHEO::coeff(int narg, char **arg) { - if (narg != 2) - error->all(FLERR, "Incorrect number of args for pair_style rheo coefficients"); - if (!allocated) - allocate(); + if (narg != 2) error->all(FLERR, "Incorrect number of args for pair_style rheo coefficients"); + if (!allocated) allocate(); int ilo, ihi, jlo, jhi; utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error); @@ -451,8 +452,7 @@ void PairRHEO::coeff(int narg, char **arg) } } - if (count == 0) - error->all(FLERR, "Incorrect args for pair rheo coefficients"); + if (count == 0) error->all(FLERR, "Incorrect args for pair rheo coefficients"); } /* ---------------------------------------------------------------------- @@ -481,7 +481,8 @@ void PairRHEO::setup() rho0 = fix_rheo->rho0; if (cutk != fix_rheo->cut) - error->all(FLERR, "Pair rheo cutoff {} does not agree with fix rheo cutoff {}", cutk, fix_rheo->cut); + error->all(FLERR, "Pair rheo cutoff {} does not agree with fix rheo cutoff {}", cutk, + fix_rheo->cut); cutksq = cutk * cutk; cutkinv = 1.0 / cutk; @@ -490,11 +491,9 @@ void PairRHEO::setup() int n = atom->ntypes; memory->create(cs, n + 1, "rheo:cs"); - for (int i = 1; i <= n; i++) - cs[i] = sqrt(csq[i]); + for (int i = 1; i <= n; i++) cs[i] = sqrt(csq[i]); - if (comm->ghost_velocity == 0) - error->all(FLERR, "Pair RHEO requires ghost atoms store velocity"); + if (comm->ghost_velocity == 0) error->all(FLERR, "Pair RHEO requires ghost atoms store velocity"); if (laplacian_order == -1) { if (fix_rheo->kernel_style == RK2) @@ -512,8 +511,7 @@ void PairRHEO::setup() double PairRHEO::init_one(int i, int j) { - if (setflag[i][j] == 0) - error->all(FLERR, "All pair rheo coeffs are not set"); + if (setflag[i][j] == 0) error->all(FLERR, "All pair rheo coeffs are not set"); return cutk; } diff --git a/src/RHEO/pair_rheo.h b/src/RHEO/pair_rheo.h index eba3a70eea..444fcc2cb4 100644 --- a/src/RHEO/pair_rheo.h +++ b/src/RHEO/pair_rheo.h @@ -37,7 +37,7 @@ class PairRHEO : public Pair { void unpack_reverse_comm(int, int *, double *) override; protected: - double cutk, *csq, *rho0; // From fix RHEO + double cutk, *csq, *rho0; // From fix RHEO double *cs, cutksq, cutkinv, cutkinv3, av, rho_damp; int laplacian_order; diff --git a/src/RHEO/pair_rheo_solid.cpp b/src/RHEO/pair_rheo_solid.cpp index 070cabaf86..6f6f7949ae 100644 --- a/src/RHEO/pair_rheo_solid.cpp +++ b/src/RHEO/pair_rheo_solid.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -221,15 +220,13 @@ void PairRHEOSolid::coeff(int narg, char **arg) void PairRHEOSolid::init_style() { if (comm->ghost_velocity == 0) - error->all(FLERR,"Pair rheo/solid requires ghost atoms store velocity"); + error->all(FLERR, "Pair rheo/solid requires ghost atoms store velocity"); - if (!atom->rheo_status_flag) - error->all(FLERR,"Pair rheo/solid requires atom_style rheo"); + if (!atom->rheo_status_flag) error->all(FLERR, "Pair rheo/solid requires atom_style rheo"); neighbor->add_request(this); } - /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */