diff --git a/src/Makefile b/src/Makefile index aa62fc473d..c95e423ba6 100755 --- a/src/Makefile +++ b/src/Makefile @@ -18,7 +18,7 @@ PACKAGE = asphere class2 colloid dipole gpu granular \ shock srd xtc PACKUSER = user-misc user-atc user-awpmd user-cg-cmm \ - user-cuda user-eff user-ewaldn user-omp user-reaxc + user-cuda user-eff user-ewaldn user-omp user-reaxc user-sph PACKALL = $(PACKAGE) $(PACKUSER) diff --git a/src/USER-SPH/Install.sh b/src/USER-SPH/Install.sh new file mode 100644 index 0000000000..56b8b50b5b --- /dev/null +++ b/src/USER-SPH/Install.sh @@ -0,0 +1,59 @@ +# Install/unInstall package files in LAMMPS + +if (test $1 = 1) then + + cp -p atom_vec_meso.cpp .. + cp -p pair_sph_heatconduction.cpp .. + cp -p pair_sph_idealgas.cpp .. + cp -p pair_sph_lj.cpp .. + cp -p pair_sph_rhosum.cpp .. + cp -p pair_sph_taitwater.cpp .. + cp -p pair_sph_taitwater_morris.cpp .. + cp -p compute_meso_e_atom.cpp .. + cp -p compute_meso_rho_atom.cpp .. + cp -p compute_meso_t_atom.cpp .. + cp -p fix_meso.cpp .. + cp -p fix_meso_stationary.cpp .. + + cp -p atom_vec_meso.h .. + cp -p pair_sph_heatconduction.h .. + cp -p pair_sph_idealgas.h .. + cp -p pair_sph_lj.h .. + cp -p pair_sph_rhosum.h .. + cp -p pair_sph_taitwater.h .. + cp -p pair_sph_taitwater_morris.h .. + cp -p compute_meso_e_atom.h .. + cp -p compute_meso_rho_atom.h .. + cp -p compute_meso_t_atom.h .. + cp -p fix_meso.h .. + cp -p fix_meso_stationary.h .. + +elif (test $1 = 0) then + + rm -f ../atom_vec_meso.cpp + rm -f ../pair_sph_heatconduction.cpp + rm -f ../pair_sph_idealgas.cpp + rm -f ../pair_sph_lj.cpp + rm -f ../pair_sph_rhosum.cpp + rm -f ../pair_sph_taitwater.cpp + rm -f ../pair_sph_taitwater_morris.cpp + rm -f ../compute_meso_e_atom.cpp + rm -f ../compute_meso_rho_atom.cpp + rm -f ../compute_meso_t_atom.cpp + rm -f ../fix_meso.cpp + rm -f ../fix_meso_stationary.cpp + + rm -f ../atom_vec_meso.h + rm -f ../pair_sph_heatconduction.h + rm -f ../pair_sph_idealgas.h + rm -f ../pair_sph_lj.h + rm -f ../pair_sph_rhosum.h + rm -f ../pair_sph_taitwater.h + rm -f ../pair_sph_taitwater_morris.h + rm -f ../compute_meso_e_atom.h + rm -f ../compute_meso_rho_atom.h + rm -f ../compute_meso_t_atom.h + rm -f ../fix_meso.h + rm -f ../fix_meso_stationary.h + +fi diff --git a/src/USER-SPH/atom_vec_meso.cpp b/src/USER-SPH/atom_vec_meso.cpp new file mode 100644 index 0000000000..c977a8e324 --- /dev/null +++ b/src/USER-SPH/atom_vec_meso.cpp @@ -0,0 +1,868 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. + ------------------------------------------------------------------------- */ + +#include "stdlib.h" +#include "atom_vec_meso.h" +#include "atom.h" +#include "domain.h" +#include "modify.h" +#include "fix.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +#define DELTA 10000 + +/* ---------------------------------------------------------------------- */ + +AtomVecMeso::AtomVecMeso(LAMMPS *lmp, int narg, char **arg) : + AtomVec(lmp, narg, arg) { + molecular = 0; + mass_type = 1; + + comm_x_only = 0; // we communicate not only x forward but also vest ... + comm_f_only = 0; // we also communicate de and drho in reverse direction + size_forward = 8; // 3 + rho + e + vest[3], that means we may only communicate 5 in hybrid + size_reverse = 5; // 3 + drho + de + size_border = 12; // 6 + rho + e + vest[3] + cv + size_velocity = 3; + size_data_atom = 8; + size_data_vel = 4; + xcol_data = 6; + + atom->e_flag = 1; + atom->rho_flag = 1; + atom->cv_flag = 1; + atom->vest_flag = 1; +} + +/* ---------------------------------------------------------------------- + grow atom arrays + n = 0 grows arrays by DELTA + n > 0 allocates arrays to size n + ------------------------------------------------------------------------- */ + +void AtomVecMeso::grow(int n) { + if (n == 0) + nmax += DELTA; + else + nmax = n; + atom->nmax = nmax; + if (nmax < 0 || nmax > MAXSMALLINT) + error->one("Per-processor system is too big"); + + tag = memory->grow(atom->tag, nmax, "atom:tag"); + type = memory->grow(atom->type, nmax, "atom:type"); + mask = memory->grow(atom->mask, nmax, "atom:mask"); + image = memory->grow(atom->image, nmax, "atom:image"); + x = memory->grow(atom->x, nmax, 3, "atom:x"); + v = memory->grow(atom->v, nmax, 3, "atom:v"); + f = memory->grow(atom->f, nmax, 3, "atom:f"); + + rho = memory->grow(atom->rho, nmax, "atom:rho"); + drho = memory->grow(atom->drho, nmax, "atom:drho"); + e = memory->grow(atom->e, nmax, "atom:e"); + de = memory->grow(atom->de, nmax, "atom:de"); + vest = memory->grow(atom->vest, nmax, 3, "atom:vest"); + cv = memory->grow(atom->cv, nmax, "atom:cv"); + + if (atom->nextra_grow) + for (int iextra = 0; iextra < atom->nextra_grow; iextra++) + modify->fix[atom->extra_grow[iextra]]->grow_arrays(nmax); +} + +/* ---------------------------------------------------------------------- + reset local array ptrs + ------------------------------------------------------------------------- */ + +void AtomVecMeso::grow_reset() { + tag = atom->tag; + type = atom->type; + mask = atom->mask; + image = atom->image; + x = atom->x; + v = atom->v; + f = atom->f; + rho = atom->rho; + drho = atom->drho; + e = atom->e; + de = atom->de; + vest = atom->vest; + cv = atom->cv; +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecMeso::copy(int i, int j, int delflag) { + //printf("in AtomVecMeso::copy\n"); + tag[j] = tag[i]; + type[j] = type[i]; + mask[j] = mask[i]; + image[j] = image[i]; + x[j][0] = x[i][0]; + x[j][1] = x[i][1]; + x[j][2] = x[i][2]; + v[j][0] = v[i][0]; + v[j][1] = v[i][1]; + v[j][2] = v[i][2]; + + rho[j] = rho[i]; + drho[j] = drho[i]; + e[j] = e[i]; + de[j] = de[i]; + cv[j] = cv[i]; + vest[j][0] = vest[i][0]; + vest[j][1] = vest[i][1]; + vest[j][2] = vest[i][2]; + + if (atom->nextra_grow) + for (int iextra = 0; iextra < atom->nextra_grow; iextra++) + modify->fix[atom->extra_grow[iextra]]->copy_arrays(i, j); +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecMeso::pack_comm_hybrid(int n, int *list, double *buf) { + //printf("in AtomVecMeso::pack_comm_hybrid\n"); + int i, j, m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = rho[j]; + buf[m++] = e[j]; + buf[m++] = vest[j][0]; + buf[m++] = vest[j][1]; + buf[m++] = vest[j][2]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecMeso::unpack_comm_hybrid(int n, int first, double *buf) { + //printf("in AtomVecMeso::unpack_comm_hybrid\n"); + int i, m, last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + rho[i] = buf[m++]; + e[i] = buf[m++]; + vest[i][0] = buf[m++]; + vest[i][1] = buf[m++]; + vest[i][2] = buf[m++]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecMeso::pack_border_hybrid(int n, int *list, double *buf) { + //printf("in AtomVecMeso::pack_border_hybrid\n"); + int i, j, m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = rho[j]; + buf[m++] = e[j]; + buf[m++] = vest[j][0]; + buf[m++] = vest[j][1]; + buf[m++] = vest[j][2]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecMeso::unpack_border_hybrid(int n, int first, double *buf) { + //printf("in AtomVecMeso::unpack_border_hybrid\n"); + int i, m, last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + rho[i] = buf[m++]; + e[i] = buf[m++]; + vest[i][0] = buf[m++]; + vest[i][1] = buf[m++]; + vest[i][2] = buf[m++]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecMeso::pack_reverse_hybrid(int n, int first, double *buf) { + //printf("in AtomVecMeso::pack_reverse_hybrid\n"); + int i, m, last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + buf[m++] = drho[i]; + buf[m++] = de[i]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecMeso::unpack_reverse_hybrid(int n, int *list, double *buf) { + //printf("in AtomVecMeso::unpack_reverse_hybrid\n"); + int i, j, m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + drho[j] += buf[m++]; + de[j] += buf[m++]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecMeso::pack_comm(int n, int *list, double *buf, int pbc_flag, + int *pbc) { + //printf("in AtomVecMeso::pack_comm\n"); + int i, j, m; + double dx, dy, dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = rho[j]; + buf[m++] = e[j]; + buf[m++] = vest[j][0]; + buf[m++] = vest[j][1]; + buf[m++] = vest[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0] * domain->xprd; + dy = pbc[1] * domain->yprd; + dz = pbc[2] * domain->zprd; + } else { + dx = pbc[0] * domain->xprd + pbc[5] * domain->xy + pbc[4] * domain->xz; + dy = pbc[1] * domain->yprd + pbc[3] * domain->yz; + dz = pbc[2] * domain->zprd; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = rho[j]; + buf[m++] = e[j]; + buf[m++] = vest[j][0]; + buf[m++] = vest[j][1]; + buf[m++] = vest[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecMeso::pack_comm_vel(int n, int *list, double *buf, int pbc_flag, + int *pbc) { + //printf("in AtomVecMeso::pack_comm_vel\n"); + int i, j, m; + double dx, dy, dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + buf[m++] = rho[j]; + buf[m++] = e[j]; + buf[m++] = vest[j][0]; + buf[m++] = vest[j][1]; + buf[m++] = vest[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0] * domain->xprd; + dy = pbc[1] * domain->yprd; + dz = pbc[2] * domain->zprd; + } else { + dx = pbc[0] * domain->xprd + pbc[5] * domain->xy + pbc[4] * domain->xz; + dy = pbc[1] * domain->yprd + pbc[3] * domain->yz; + dz = pbc[2] * domain->zprd; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + buf[m++] = rho[j]; + buf[m++] = e[j]; + buf[m++] = vest[j][0]; + buf[m++] = vest[j][1]; + buf[m++] = vest[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecMeso::unpack_comm(int n, int first, double *buf) { + //printf("in AtomVecMeso::unpack_comm\n"); + int i, m, last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + rho[i] = buf[m++]; + e[i] = buf[m++]; + vest[i][0] = buf[m++]; + vest[i][1] = buf[m++]; + vest[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecMeso::unpack_comm_vel(int n, int first, double *buf) { + //printf("in AtomVecMeso::unpack_comm_vel\n"); + int i, m, last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + rho[i] = buf[m++]; + e[i] = buf[m++]; + vest[i][0] = buf[m++]; + vest[i][1] = buf[m++]; + vest[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecMeso::pack_reverse(int n, int first, double *buf) { + //printf("in AtomVecMeso::pack_reverse\n"); + int i, m, last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + buf[m++] = f[i][0]; + buf[m++] = f[i][1]; + buf[m++] = f[i][2]; + buf[m++] = drho[i]; + buf[m++] = de[i]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecMeso::unpack_reverse(int n, int *list, double *buf) { + //printf("in AtomVecMeso::unpack_reverse\n"); + int i, j, m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + f[j][0] += buf[m++]; + f[j][1] += buf[m++]; + f[j][2] += buf[m++]; + drho[j] += buf[m++]; + de[j] += buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecMeso::pack_border(int n, int *list, double *buf, int pbc_flag, + int *pbc) { + //printf("in AtomVecMeso::pack_border\n"); + int i, j, m; + double dx, dy, dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = rho[j]; + buf[m++] = e[j]; + buf[m++] = cv[j]; + buf[m++] = vest[j][0]; + buf[m++] = vest[j][1]; + buf[m++] = vest[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0] * domain->xprd; + dy = pbc[1] * domain->yprd; + dz = pbc[2] * domain->zprd; + } else { + dx = pbc[0]; + dy = pbc[1]; + dz = pbc[2]; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = rho[j]; + buf[m++] = e[j]; + buf[m++] = cv[j]; + buf[m++] = vest[j][0]; + buf[m++] = vest[j][1]; + buf[m++] = vest[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecMeso::pack_border_vel(int n, int *list, double *buf, int pbc_flag, + int *pbc) { + //printf("in AtomVecMeso::pack_border_vel\n"); + int i, j, m; + double dx, dy, dz; + + m = 0; + if (pbc_flag == 0) { + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0]; + buf[m++] = x[j][1]; + buf[m++] = x[j][2]; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + buf[m++] = rho[j]; + buf[m++] = e[j]; + buf[m++] = cv[j]; + buf[m++] = vest[j][0]; + buf[m++] = vest[j][1]; + buf[m++] = vest[j][2]; + } + } else { + if (domain->triclinic == 0) { + dx = pbc[0] * domain->xprd; + dy = pbc[1] * domain->yprd; + dz = pbc[2] * domain->zprd; + } else { + dx = pbc[0]; + dy = pbc[1]; + dz = pbc[2]; + } + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = x[j][0] + dx; + buf[m++] = x[j][1] + dy; + buf[m++] = x[j][2] + dz; + buf[m++] = tag[j]; + buf[m++] = type[j]; + buf[m++] = mask[j]; + buf[m++] = v[j][0]; + buf[m++] = v[j][1]; + buf[m++] = v[j][2]; + buf[m++] = rho[j]; + buf[m++] = e[j]; + buf[m++] = cv[j]; + buf[m++] = vest[j][0]; + buf[m++] = vest[j][1]; + buf[m++] = vest[j][2]; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecMeso::unpack_border(int n, int first, double *buf) { + //printf("in AtomVecMeso::unpack_border\n"); + int i, m, last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (i == nmax) + grow(0); + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + tag[i] = static_cast (buf[m++]); + type[i] = static_cast (buf[m++]); + mask[i] = static_cast (buf[m++]); + rho[i] = buf[m++]; + e[i] = buf[m++]; + cv[i] = buf[m++]; + vest[i][0] = buf[m++]; + vest[i][1] = buf[m++]; + vest[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecMeso::unpack_border_vel(int n, int first, double *buf) { + //printf("in AtomVecMeso::unpack_border_vel\n"); + int i, m, last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (i == nmax) + grow(0); + x[i][0] = buf[m++]; + x[i][1] = buf[m++]; + x[i][2] = buf[m++]; + tag[i] = static_cast (buf[m++]); + type[i] = static_cast (buf[m++]); + mask[i] = static_cast (buf[m++]); + v[i][0] = buf[m++]; + v[i][1] = buf[m++]; + v[i][2] = buf[m++]; + rho[i] = buf[m++]; + e[i] = buf[m++]; + cv[i] = buf[m++]; + vest[i][0] = buf[m++]; + vest[i][1] = buf[m++]; + vest[i][2] = buf[m++]; + } +} + +/* ---------------------------------------------------------------------- + pack data for atom I for sending to another proc + xyz must be 1st 3 values, so comm::exchange() can test on them + ------------------------------------------------------------------------- */ + +int AtomVecMeso::pack_exchange(int i, double *buf) { + //printf("in AtomVecMeso::pack_exchange\n"); + int m = 1; + buf[m++] = x[i][0]; + buf[m++] = x[i][1]; + buf[m++] = x[i][2]; + buf[m++] = v[i][0]; + buf[m++] = v[i][1]; + buf[m++] = v[i][2]; + buf[m++] = tag[i]; + buf[m++] = type[i]; + buf[m++] = mask[i]; + buf[m++] = image[i]; + buf[m++] = rho[i]; + buf[m++] = e[i]; + buf[m++] = cv[i]; + buf[m++] = vest[i][0]; + buf[m++] = vest[i][1]; + buf[m++] = vest[i][2]; + + if (atom->nextra_grow) + for (int iextra = 0; iextra < atom->nextra_grow; iextra++) + m += modify->fix[atom->extra_grow[iextra]]->pack_exchange(i, &buf[m]); + + buf[0] = m; + return m; +} + +/* ---------------------------------------------------------------------- */ + +int AtomVecMeso::unpack_exchange(double *buf) { + //printf("in AtomVecMeso::unpack_exchange\n"); + int nlocal = atom->nlocal; + if (nlocal == nmax) + grow(0); + + int m = 1; + x[nlocal][0] = buf[m++]; + x[nlocal][1] = buf[m++]; + x[nlocal][2] = buf[m++]; + v[nlocal][0] = buf[m++]; + v[nlocal][1] = buf[m++]; + v[nlocal][2] = buf[m++]; + tag[nlocal] = static_cast (buf[m++]); + type[nlocal] = static_cast (buf[m++]); + mask[nlocal] = static_cast (buf[m++]); + image[nlocal] = static_cast (buf[m++]); + rho[nlocal] = buf[m++]; + e[nlocal] = buf[m++]; + cv[nlocal] = buf[m++]; + vest[nlocal][0] = buf[m++]; + vest[nlocal][1] = buf[m++]; + vest[nlocal][2] = buf[m++]; + + if (atom->nextra_grow) + for (int iextra = 0; iextra < atom->nextra_grow; iextra++) + m += modify->fix[atom->extra_grow[iextra]]-> unpack_exchange(nlocal, + &buf[m]); + + atom->nlocal++; + return m; +} + +/* ---------------------------------------------------------------------- + size of restart data for all atoms owned by this proc + include extra data stored by fixes + ------------------------------------------------------------------------- */ + +int AtomVecMeso::size_restart() { + int i; + + int nlocal = atom->nlocal; + int n = 17 * nlocal; // 11 + rho + e + cv + vest[3] + + if (atom->nextra_restart) + for (int iextra = 0; iextra < atom->nextra_restart; iextra++) + for (i = 0; i < nlocal; i++) + n += modify->fix[atom->extra_restart[iextra]]->size_restart(i); + + return n; +} + +/* ---------------------------------------------------------------------- + pack atom I's data for restart file including extra quantities + xyz must be 1st 3 values, so that read_restart can test on them + molecular types may be negative, but write as positive + ------------------------------------------------------------------------- */ + +int AtomVecMeso::pack_restart(int i, double *buf) { + int m = 1; + buf[m++] = x[i][0]; + buf[m++] = x[i][1]; + buf[m++] = x[i][2]; + buf[m++] = tag[i]; + buf[m++] = type[i]; + buf[m++] = mask[i]; + buf[m++] = image[i]; + buf[m++] = v[i][0]; + buf[m++] = v[i][1]; + buf[m++] = v[i][2]; + buf[m++] = rho[i]; + buf[m++] = e[i]; + buf[m++] = cv[i]; + buf[m++] = vest[i][0]; + buf[m++] = vest[i][1]; + buf[m++] = vest[i][2]; + + if (atom->nextra_restart) + for (int iextra = 0; iextra < atom->nextra_restart; iextra++) + m += modify->fix[atom->extra_restart[iextra]]->pack_restart(i, &buf[m]); + + buf[0] = m; + return m; +} + +/* ---------------------------------------------------------------------- + unpack data for one atom from restart file including extra quantities + ------------------------------------------------------------------------- */ + +int AtomVecMeso::unpack_restart(double *buf) { + int nlocal = atom->nlocal; + if (nlocal == nmax) { + grow(0); + if (atom->nextra_store) + memory->grow(atom->extra, nmax, atom->nextra_store, "atom:extra"); + } + + int m = 1; + x[nlocal][0] = buf[m++]; + x[nlocal][1] = buf[m++]; + x[nlocal][2] = buf[m++]; + tag[nlocal] = static_cast (buf[m++]); + type[nlocal] = static_cast (buf[m++]); + mask[nlocal] = static_cast (buf[m++]); + image[nlocal] = static_cast (buf[m++]); + v[nlocal][0] = buf[m++]; + v[nlocal][1] = buf[m++]; + v[nlocal][2] = buf[m++]; + rho[nlocal] = buf[m++]; + e[nlocal] = buf[m++]; + cv[nlocal] = buf[m++]; + vest[nlocal][0] = buf[m++]; + vest[nlocal][1] = buf[m++]; + vest[nlocal][2] = buf[m++]; + + double **extra = atom->extra; + if (atom->nextra_store) { + int size = static_cast (buf[0]) - m; + for (int i = 0; i < size; i++) + extra[nlocal][i] = buf[m++]; + } + + atom->nlocal++; + return m; +} + +/* ---------------------------------------------------------------------- + create one atom of itype at coord + set other values to defaults + ------------------------------------------------------------------------- */ + +void AtomVecMeso::create_atom(int itype, double *coord) { + int nlocal = atom->nlocal; + if (nlocal == nmax) + grow(0); + + tag[nlocal] = 0; + type[nlocal] = itype; + x[nlocal][0] = coord[0]; + x[nlocal][1] = coord[1]; + x[nlocal][2] = coord[2]; + mask[nlocal] = 1; + image[nlocal] = (512 << 20) | (512 << 10) | 512; + v[nlocal][0] = 0.0; + v[nlocal][1] = 0.0; + v[nlocal][2] = 0.0; + rho[nlocal] = 0.0; + e[nlocal] = 0.0; + cv[nlocal] = 1.0; + vest[nlocal][0] = 0.0; + vest[nlocal][1] = 0.0; + vest[nlocal][2] = 0.0; + de[nlocal] = 0.0; + drho[nlocal] = 0.0; + + atom->nlocal++; +} + +/* ---------------------------------------------------------------------- + unpack one line from Atoms section of data file + initialize other atom quantities + ------------------------------------------------------------------------- */ + +void AtomVecMeso::data_atom(double *coord, int imagetmp, char **values) { + int nlocal = atom->nlocal; + if (nlocal == nmax) + grow(0); + + tag[nlocal] = atoi(values[0]); + if (tag[nlocal] <= 0) + error->one("Invalid atom ID in Atoms section of data file"); + + type[nlocal] = atoi(values[1]); + if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes) + error->one("Invalid atom type in Atoms section of data file"); + + rho[nlocal] = atof(values[2]); + e[nlocal] = atof(values[3]); + cv[nlocal] = atof(values[4]); + + x[nlocal][0] = coord[0]; + x[nlocal][1] = coord[1]; + x[nlocal][2] = coord[2]; + + //printf("rho=%f, e=%f, cv=%f, x=%f\n", rho[nlocal], e[nlocal], cv[nlocal], x[nlocal][0]); + + image[nlocal] = imagetmp; + + mask[nlocal] = 1; + v[nlocal][0] = 0.0; + v[nlocal][1] = 0.0; + v[nlocal][2] = 0.0; + + vest[nlocal][0] = 0.0; + vest[nlocal][1] = 0.0; + vest[nlocal][2] = 0.0; + + de[nlocal] = 0.0; + drho[nlocal] = 0.0; + + atom->nlocal++; +} + +/* ---------------------------------------------------------------------- + unpack hybrid quantities from one line in Atoms section of data file + initialize other atom quantities for this sub-style + ------------------------------------------------------------------------- */ + +int AtomVecMeso::data_atom_hybrid(int nlocal, char **values) { + + rho[nlocal] = atof(values[0]); + e[nlocal] = atof(values[1]); + cv[nlocal] = atof(values[2]); + + return 3; +} + +/* ---------------------------------------------------------------------- + return # of bytes of allocated memory + ------------------------------------------------------------------------- */ + +bigint AtomVecMeso::memory_usage() { + bigint bytes = 0; + + if (atom->memcheck("tag")) + bytes += memory->usage(tag, nmax); + if (atom->memcheck("type")) + bytes += memory->usage(type, nmax); + if (atom->memcheck("mask")) + bytes += memory->usage(mask, nmax); + if (atom->memcheck("image")) + bytes += memory->usage(image, nmax); + if (atom->memcheck("x")) + bytes += memory->usage(x, nmax); + if (atom->memcheck("v")) + bytes += memory->usage(v, nmax); + if (atom->memcheck("f")) + bytes += memory->usage(f, nmax); + if (atom->memcheck("rho")) + bytes += memory->usage(rho, nmax); + if (atom->memcheck("drho")) + bytes += memory->usage(drho, nmax); + if (atom->memcheck("e")) + bytes += memory->usage(e, nmax); + if (atom->memcheck("de")) + bytes += memory->usage(de, nmax); + if (atom->memcheck("cv")) + bytes += memory->usage(cv, nmax); + if (atom->memcheck("vest")) + bytes += memory->usage(vest, nmax); + + return bytes; +} diff --git a/src/USER-SPH/atom_vec_meso.h b/src/USER-SPH/atom_vec_meso.h new file mode 100644 index 0000000000..b187198c09 --- /dev/null +++ b/src/USER-SPH/atom_vec_meso.h @@ -0,0 +1,74 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef ATOM_CLASS + +AtomStyle(meso,AtomVecMeso) + +#else + +#ifndef LMP_ATOM_VEC_MESO_H +#define LMP_ATOM_VEC_MESO_H + +#include "atom_vec.h" + +namespace LAMMPS_NS { + +class AtomVecMeso : public AtomVec { + public: + AtomVecMeso(class LAMMPS *, int, char **); + ~AtomVecMeso() {} + void grow(int); + void grow_reset(); + void copy(int, int, int); + int pack_comm(int, int *, double *, int, int *); + int pack_comm_vel(int, int *, double *, int, int *); + void unpack_comm(int, int, double *); + void unpack_comm_vel(int, int, double *); + int pack_reverse(int, int, double *); + void unpack_reverse(int, int *, double *); + + int pack_comm_hybrid(int, int *, double *); + int unpack_comm_hybrid(int, int, double *); + + int pack_border_hybrid(int, int *, double *); + int unpack_border_hybrid(int, int, double *); + + int pack_reverse_hybrid(int, int, double *); + int unpack_reverse_hybrid(int, int *, double *); + + int pack_border(int, int *, double *, int, int *); + int pack_border_vel(int, int *, double *, int, int *); + void unpack_border(int, int, double *); + void unpack_border_vel(int, int, double *); + int pack_exchange(int, double *); + int unpack_exchange(double *); + int size_restart(); + int pack_restart(int, double *); + int unpack_restart(double *); + void create_atom(int, double *); + void data_atom(double *, int, char **); + int data_atom_hybrid(int, char **); + bigint memory_usage(); + + private: + int *tag,*type,*mask,*image; + double **x,**v,**f; + double *rho, *drho, *e, *de, *cv; + double **vest; // estimated velocity during force computation +}; + +} + +#endif +#endif diff --git a/src/USER-SPH/compute_meso_e_atom.cpp b/src/USER-SPH/compute_meso_e_atom.cpp new file mode 100644 index 0000000000..934818584b --- /dev/null +++ b/src/USER-SPH/compute_meso_e_atom.cpp @@ -0,0 +1,97 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "string.h" +#include "compute_meso_e_atom.h" +#include "atom.h" +#include "update.h" +#include "modify.h" +#include "comm.h" +#include "force.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeMesoEAtom::ComputeMesoEAtom(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg != 3) error->all("Number of arguments for compute meso_e/atom command != 3"); + if (atom->e_flag != 1) error->all("compute meso_e/atom command requires atom_style with energy (e.g. meso)"); + + peratom_flag = 1; + size_peratom_cols = 0; + + nmax = 0; + evector = NULL; +} + +/* ---------------------------------------------------------------------- */ + +ComputeMesoEAtom::~ComputeMesoEAtom() +{ + memory->sfree(evector); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeMesoEAtom::init() +{ + + int count = 0; + for (int i = 0; i < modify->ncompute; i++) + if (strcmp(modify->compute[i]->style,"evector/atom") == 0) count++; + if (count > 1 && comm->me == 0) + error->warning("More than one compute evector/atom"); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeMesoEAtom::compute_peratom() +{ + invoked_peratom = update->ntimestep; + + // grow evector array if necessary + + if (atom->nlocal > nmax) { + memory->sfree(evector); + nmax = atom->nmax; + evector = (double *) memory->smalloc(nmax*sizeof(double),"evector/atom:evector"); + vector_atom = evector; + } + + double *e = atom->e; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + evector[i] = e[i]; + } + else { + evector[i] = 0.0; + } + } +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +double ComputeMesoEAtom::memory_usage() +{ + double bytes = nmax * sizeof(double); + return bytes; +} diff --git a/src/USER-SPH/compute_meso_e_atom.h b/src/USER-SPH/compute_meso_e_atom.h new file mode 100644 index 0000000000..07fc01c920 --- /dev/null +++ b/src/USER-SPH/compute_meso_e_atom.h @@ -0,0 +1,43 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS + +ComputeStyle(meso_e/atom,ComputeMesoEAtom) + +#else + +#ifndef LMP_COMPUTE_MESO_E_ATOM_H +#define LMP_COMPUTE_MESO_E_ATOM_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeMesoEAtom : public Compute { + public: + ComputeMesoEAtom(class LAMMPS *, int, char **); + ~ComputeMesoEAtom(); + void init(); + void compute_peratom(); + double memory_usage(); + + private: + int nmax; + double *evector; +}; + +} + +#endif +#endif diff --git a/src/USER-SPH/compute_meso_rho_atom.cpp b/src/USER-SPH/compute_meso_rho_atom.cpp new file mode 100644 index 0000000000..98aab3ad72 --- /dev/null +++ b/src/USER-SPH/compute_meso_rho_atom.cpp @@ -0,0 +1,99 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "string.h" +#include "compute_meso_rho_atom.h" +#include "atom.h" +#include "update.h" +#include "modify.h" +#include "comm.h" +#include "force.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeMesoRhoAtom::ComputeMesoRhoAtom(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg != 3) error->all("Illegal compute meso_rho/atom command"); + if (atom->rho_flag != 1) error->all("compute meso_rho/atom command requires atom_style with density (e.g. meso)"); + + peratom_flag = 1; + size_peratom_cols = 0; + + nmax = 0; + rhoVector = NULL; +} + +/* ---------------------------------------------------------------------- */ + +ComputeMesoRhoAtom::~ComputeMesoRhoAtom() +{ + memory->sfree(rhoVector); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeMesoRhoAtom::init() +{ + + int count = 0; + for (int i = 0; i < modify->ncompute; i++) + if (strcmp(modify->compute[i]->style,"rhoVector/atom") == 0) count++; + if (count > 1 && comm->me == 0) + error->warning("More than one compute rhoVector/atom"); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeMesoRhoAtom::compute_peratom() +{ + invoked_peratom = update->ntimestep; + + // grow rhoVector array if necessary + + if (atom->nlocal > nmax) { + memory->sfree(rhoVector); + nmax = atom->nmax; + rhoVector = (double *) memory->smalloc(nmax*sizeof(double),"atom:rhoVector"); + vector_atom = rhoVector; + } + + // compute kinetic energy for each atom in group + + double *rho = atom->rho; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + rhoVector[i] = rho[i]; + } + else { + rhoVector[i] = 0.0; + } + } +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +double ComputeMesoRhoAtom::memory_usage() +{ + double bytes = nmax * sizeof(double); + return bytes; +} diff --git a/src/USER-SPH/compute_meso_rho_atom.h b/src/USER-SPH/compute_meso_rho_atom.h new file mode 100644 index 0000000000..e1100a2f81 --- /dev/null +++ b/src/USER-SPH/compute_meso_rho_atom.h @@ -0,0 +1,43 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS + +ComputeStyle(meso_rho/atom,ComputeMesoRhoAtom) + +#else + +#ifndef LMP_COMPUTE_MESO_RHO_ATOM_H +#define LMP_COMPUTE_MESO_RHO_ATOM_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeMesoRhoAtom : public Compute { + public: + ComputeMesoRhoAtom(class LAMMPS *, int, char **); + ~ComputeMesoRhoAtom(); + void init(); + void compute_peratom(); + double memory_usage(); + + private: + int nmax; + double *rhoVector; +}; + +} + +#endif +#endif diff --git a/src/USER-SPH/compute_meso_t_atom.cpp b/src/USER-SPH/compute_meso_t_atom.cpp new file mode 100644 index 0000000000..21a9a4d660 --- /dev/null +++ b/src/USER-SPH/compute_meso_t_atom.cpp @@ -0,0 +1,101 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "string.h" +#include "compute_meso_t_atom.h" +#include "atom.h" +#include "update.h" +#include "modify.h" +#include "comm.h" +#include "force.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeMesoTAtom::ComputeMesoTAtom(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg != 3) error->all("Number of arguments for compute meso_t/atom command != 3"); + if ((atom->e_flag != 1) || (atom->cv_flag != 1)) + error->all("compute meso_e/atom command requires atom_style with both energy and heat capacity (e.g. meso)"); + + peratom_flag = 1; + size_peratom_cols = 0; + + nmax = 0; + tvector = NULL; +} + +/* ---------------------------------------------------------------------- */ + +ComputeMesoTAtom::~ComputeMesoTAtom() +{ + memory->sfree(tvector); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeMesoTAtom::init() +{ + + int count = 0; + for (int i = 0; i < modify->ncompute; i++) + if (strcmp(modify->compute[i]->style,"meso_t/atom") == 0) count++; + if (count > 1 && comm->me == 0) + error->warning("More than one compute meso_t/atom"); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeMesoTAtom::compute_peratom() +{ + invoked_peratom = update->ntimestep; + + // grow tvector array if necessary + + if (atom->nlocal > nmax) { + memory->sfree(tvector); + nmax = atom->nmax; + tvector = (double *) memory->smalloc(nmax*sizeof(double),"tvector/atom:tvector"); + vector_atom = tvector; + } + + double *e = atom->e; + double *cv = atom->cv; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + if (cv[i] > 0.0) { + tvector[i] = e[i] / cv[i]; + } + } + else { + tvector[i] = 0.0; + } + } +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +double ComputeMesoTAtom::memory_usage() +{ + double bytes = nmax * sizeof(double); + return bytes; +} diff --git a/src/USER-SPH/compute_meso_t_atom.h b/src/USER-SPH/compute_meso_t_atom.h new file mode 100644 index 0000000000..279d749f87 --- /dev/null +++ b/src/USER-SPH/compute_meso_t_atom.h @@ -0,0 +1,43 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS + +ComputeStyle(meso_t/atom,ComputeMesoTAtom) + +#else + +#ifndef LMP_COMPUTE_MESO_T_ATOM_H +#define LMP_COMPUTE_MESO_T_ATOM_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeMesoTAtom : public Compute { + public: + ComputeMesoTAtom(class LAMMPS *, int, char **); + ~ComputeMesoTAtom(); + void init(); + void compute_peratom(); + double memory_usage(); + + private: + int nmax; + double *tvector; +}; + +} + +#endif +#endif diff --git a/src/USER-SPH/fix_meso.cpp b/src/USER-SPH/fix_meso.cpp new file mode 100644 index 0000000000..0b8576f386 --- /dev/null +++ b/src/USER-SPH/fix_meso.cpp @@ -0,0 +1,167 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. + ------------------------------------------------------------------------- */ + +#include "stdio.h" +#include "string.h" +#include "fix_meso.h" +#include "math.h" +#include "stdlib.h" +#include "string.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "update.h" +#include "integrate.h" +#include "respa.h" +#include "memory.h" +#include "error.h" +#include "pair.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +FixMeso::FixMeso(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) { + + if ((atom->e_flag != 1) || (atom->rho_flag != 1)) + error->all( + "fix meso command requires atom_style with both energy and density"); + + if (narg != 3) + error->all("Illegal number of arguments for fix meso command"); + + time_integrate = 1; +} + +/* ---------------------------------------------------------------------- */ + +int FixMeso::setmask() { + int mask = 0; + mask |= INITIAL_INTEGRATE; + mask |= FINAL_INTEGRATE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixMeso::init() { + dtv = update->dt; + dtf = 0.5 * update->dt * force->ftm2v; +} + +/* ---------------------------------------------------------------------- + allow for both per-type and per-atom mass + ------------------------------------------------------------------------- */ + +void FixMeso::initial_integrate(int vflag) { + // update v and x and rho and e of atoms in group + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + double **vest = atom->vest; + double *rho = atom->rho; + double *drho = atom->drho; + double *e = atom->e; + double *de = atom->de; + double *mass = atom->mass; + double *rmass = atom->rmass; + int rmass_flag = atom->rmass_flag; + + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int i; + double dtfm; + + if (igroup == atom->firstgroup) + nlocal = atom->nfirst; + + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + if (rmass_flag) { + dtfm = dtf / rmass[i]; + } else { + dtfm = dtf / mass[type[i]]; + } + + e[i] += dtf * de[i]; // half-step update of particle internal energy + rho[i] += dtf * drho[i]; // ... and density + + // extrapolate velocity for use with velocity-dependent potentials, e.g. SPH + vest[i][0] = v[i][0] + 2.0 * dtfm * f[i][0]; + vest[i][1] = v[i][1] + 2.0 * dtfm * f[i][1]; + vest[i][2] = v[i][2] + 2.0 * dtfm * f[i][2]; + + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixMeso::final_integrate() { + + // update v, rho, and e of atoms in group + + double **v = atom->v; + double **f = atom->f; + double *e = atom->e; + double *de = atom->de; + double *rho = atom->rho; + double *drho = atom->drho; + int *type = atom->type; + int *mask = atom->mask; + double *mass = atom->mass; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) + nlocal = atom->nfirst; + double dtfm; + double *rmass = atom->rmass; + int rmass_flag = atom->rmass_flag; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + + if (rmass_flag) { + dtfm = dtf / rmass[i]; + } else { + dtfm = dtf / mass[type[i]]; + } + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + + e[i] += dtf * de[i]; + rho[i] += dtf * drho[i]; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixMeso::reset_dt() { + dtv = update->dt; + dtf = 0.5 * update->dt * force->ftm2v; +} + diff --git a/src/USER-SPH/fix_meso.h b/src/USER-SPH/fix_meso.h new file mode 100644 index 0000000000..eee5ee1d45 --- /dev/null +++ b/src/USER-SPH/fix_meso.h @@ -0,0 +1,49 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(meso,FixMeso) + +#else + +#ifndef LMP_FIX_MESO_H +#define LMP_FIX_MESO_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixMeso : public Fix { + public: + FixMeso(class LAMMPS *, int, char **); + int setmask(); + virtual void init(); + virtual void initial_integrate(int); + virtual void final_integrate(); + void reset_dt(); + + private: + class NeighList *list; + protected: + double dtv,dtf; + double *step_respa; + int mass_require; + + class Pair *pair; +}; + +} + +#endif +#endif diff --git a/src/USER-SPH/fix_meso_stationary.cpp b/src/USER-SPH/fix_meso_stationary.cpp new file mode 100644 index 0000000000..9a08d28fd9 --- /dev/null +++ b/src/USER-SPH/fix_meso_stationary.cpp @@ -0,0 +1,122 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. + ------------------------------------------------------------------------- */ + +#include "stdio.h" +#include "string.h" +#include "fix_meso_stationary.h" +#include "math.h" +#include "stdlib.h" +#include "string.h" +#include "atom.h" +#include "comm.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "update.h" +#include "integrate.h" +#include "respa.h" +#include "memory.h" +#include "error.h" +#include "pair.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +FixMesoStationary::FixMesoStationary(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) { + + if ((atom->e_flag != 1) || (atom->rho_flag != 1)) + error->all( + "fix meso/stationary command requires atom_style with both energy and density, e.g. meso"); + + if (narg != 3) + error->all("Illegal number of arguments for fix meso/stationary command"); + + time_integrate = 0; +} + +/* ---------------------------------------------------------------------- */ + +int FixMesoStationary::setmask() { + int mask = 0; + mask |= INITIAL_INTEGRATE; + mask |= FINAL_INTEGRATE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixMesoStationary::init() { + dtv = update->dt; + dtf = 0.5 * update->dt * force->ftm2v; +} + +/* ---------------------------------------------------------------------- + allow for both per-type and per-atom mass + ------------------------------------------------------------------------- */ + +void FixMesoStationary::initial_integrate(int vflag) { + + double *rho = atom->rho; + double *drho = atom->drho; + double *e = atom->e; + double *de = atom->de; + + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int i; + + if (igroup == atom->firstgroup) + nlocal = atom->nfirst; + + for (i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + e[i] += dtf * de[i]; // half-step update of particle internal energy + rho[i] += dtf * drho[i]; // ... and density + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixMesoStationary::final_integrate() { + + double *e = atom->e; + double *de = atom->de; + double *rho = atom->rho; + double *drho = atom->drho; + int *type = atom->type; + int *mask = atom->mask; + double *mass = atom->mass; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) + nlocal = atom->nfirst; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + e[i] += dtf * de[i]; + rho[i] += dtf * drho[i]; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixMesoStationary::reset_dt() { + dtv = update->dt; + dtf = 0.5 * update->dt * force->ftm2v; +} + diff --git a/src/USER-SPH/fix_meso_stationary.h b/src/USER-SPH/fix_meso_stationary.h new file mode 100644 index 0000000000..08869d2105 --- /dev/null +++ b/src/USER-SPH/fix_meso_stationary.h @@ -0,0 +1,49 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(meso/stationary,FixMesoStationary) + +#else + +#ifndef LMP_FIX_MESO_STATIONARY_H +#define LMP_FIX_MESO_STATIONARY_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixMesoStationary : public Fix { + public: + FixMesoStationary(class LAMMPS *, int, char **); + int setmask(); + virtual void init(); + virtual void initial_integrate(int); + virtual void final_integrate(); + void reset_dt(); + + private: + class NeighList *list; + protected: + double dtv,dtf; + double *step_respa; + int mass_require; + + class Pair *pair; +}; + +} + +#endif +#endif diff --git a/src/USER-SPH/pair_sph_heatconduction.cpp b/src/USER-SPH/pair_sph_heatconduction.cpp new file mode 100644 index 0000000000..e4af6f7663 --- /dev/null +++ b/src/USER-SPH/pair_sph_heatconduction.cpp @@ -0,0 +1,221 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. + ------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdlib.h" +#include "pair_sph_heatconduction.h" +#include "atom.h" +#include "force.h" +#include "comm.h" +#include "memory.h" +#include "error.h" +#include "neigh_list.h" +#include "domain.h" + +using namespace LAMMPS_NS; + +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#define MAX(a,b) ((a) > (b) ? (a) : (b)) + +/* ---------------------------------------------------------------------- */ + +PairSPHHeatConduction::PairSPHHeatConduction(LAMMPS *lmp) : + Pair(lmp) { +} + +/* ---------------------------------------------------------------------- */ + +PairSPHHeatConduction::~PairSPHHeatConduction() { + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + memory->destroy(cut); + memory->destroy(alpha); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairSPHHeatConduction::compute(int eflag, int vflag) { + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, delx, dely, delz; + + int *ilist, *jlist, *numneigh, **firstneigh; + double imass, jmass, h, ih, ihsq; + double rsq, wfd, D, deltaE; + + if (eflag || vflag) + ev_setup(eflag, vflag); + else + evflag = vflag_fdotr = 0; + + double **x = atom->x; + double *e = atom->e; + double *de = atom->de; + double *mass = atom->mass; + double *rho = atom->rho; + int *type = atom->type; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms and do heat diffusion + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + itype = type[i]; + + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + + jlist = firstneigh[i]; + jnum = numneigh[i]; + + imass = mass[itype]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + jtype = type[j]; + jmass = mass[jtype]; + + if (rsq < cutsq[itype][jtype]) { + h = cut[itype][jtype]; + ih = 1.0 / h; + ihsq = ih * ih; + + // kernel function + wfd = h - sqrt(rsq); + if (domain->dimension == 3) { + // Lucy Kernel, 3d + // Note that wfd, the derivative of the weight function with respect to r, + // is lacking a factor of r. + // The missing factor of r is recovered by + // deltaE, which is missing a factor of 1/r + wfd = -25.066903536973515383e0 * wfd * wfd * ihsq * ihsq * ihsq * ih; + } else { + // Lucy Kernel, 2d + wfd = -19.098593171027440292e0 * wfd * wfd * ihsq * ihsq * ihsq; + } + + jmass = mass[jtype]; + D = alpha[itype][jtype]; // diffusion coefficient + + deltaE = 2.0 * imass * jmass / (imass+jmass); + deltaE *= (rho[i] + rho[j]) / (rho[i] * rho[j]); + deltaE *= D * (e[i] - e[j]) * wfd; + + de[i] += deltaE; + if (newton_pair || j < nlocal) { + de[j] -= deltaE; + } + + } + } + } +} + +/* ---------------------------------------------------------------------- + allocate all arrays + ------------------------------------------------------------------------- */ + +void PairSPHHeatConduction::allocate() { + allocated = 1; + int n = atom->ntypes; + + 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; + + memory->create(cutsq, n + 1, n + 1, "pair:cutsq"); + memory->create(cut, n + 1, n + 1, "pair:cut"); + memory->create(alpha, n + 1, n + 1, "pair:alpha"); +} + +/* ---------------------------------------------------------------------- + global settings + ------------------------------------------------------------------------- */ + +void PairSPHHeatConduction::settings(int narg, char **arg) { + if (narg != 0) + error->all( + "Illegal number of setting arguments for pair_style sph/heatconduction"); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs + ------------------------------------------------------------------------- */ + +void PairSPHHeatConduction::coeff(int narg, char **arg) { + if (narg != 4) + error->all("Incorrect number of args for pair_style sph/heatconduction coefficients"); + if (!allocated) + allocate(); + + int ilo, ihi, jlo, jhi; + force->bounds(arg[0], atom->ntypes, ilo, ihi); + force->bounds(arg[1], atom->ntypes, jlo, jhi); + + double alpha_one = force->numeric(arg[2]); + double cut_one = force->numeric(arg[3]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + //printf("setting cut[%d][%d] = %f\n", i, j, cut_one); + cut[i][j] = cut_one; + alpha[i][j] = alpha_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) + error->all("Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i + ------------------------------------------------------------------------- */ + +double PairSPHHeatConduction::init_one(int i, int j) { + + if (setflag[i][j] == 0) { + error->all("All pair sph/heatconduction coeffs are not set"); + } + + cut[j][i] = cut[i][j]; + alpha[j][i] = alpha[i][j]; + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- */ + +double PairSPHHeatConduction::single(int i, int j, int itype, int jtype, + double rsq, double factor_coul, double factor_lj, double &fforce) { + fforce = 0.0; + + return 0.0; +} diff --git a/src/USER-SPH/pair_sph_heatconduction.h b/src/USER-SPH/pair_sph_heatconduction.h new file mode 100644 index 0000000000..81402fda2e --- /dev/null +++ b/src/USER-SPH/pair_sph_heatconduction.h @@ -0,0 +1,45 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(sph/heatconduction,PairSPHHeatConduction) + +#else + +#ifndef LMP_PAIR_SPH_HEATCONDUCTION_H +#define LMP_PAIR_SPH_HEATCONDUCTION_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairSPHHeatConduction : public Pair { + public: + PairSPHHeatConduction(class LAMMPS *); + virtual ~PairSPHHeatConduction(); + virtual void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + virtual double init_one(int, int); + virtual double single(int, int, int, int, double, double, double, double &); + + protected: + double **cut, **alpha; + void allocate(); +}; + +} + +#endif +#endif diff --git a/src/USER-SPH/pair_sph_idealgas.cpp b/src/USER-SPH/pair_sph_idealgas.cpp new file mode 100644 index 0000000000..35baf218cb --- /dev/null +++ b/src/USER-SPH/pair_sph_idealgas.cpp @@ -0,0 +1,262 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. + ------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdlib.h" +#include "pair_sph_idealgas.h" +#include "atom.h" +#include "force.h" +#include "comm.h" +#include "neigh_list.h" +#include "memory.h" +#include "error.h" +#include "domain.h" + +using namespace LAMMPS_NS; + +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#define MAX(a,b) ((a) > (b) ? (a) : (b)) + +/* ---------------------------------------------------------------------- */ + +PairSPHIdealGas::PairSPHIdealGas(LAMMPS *lmp) : + Pair(lmp) { +} + +/* ---------------------------------------------------------------------- */ + +PairSPHIdealGas::~PairSPHIdealGas() { + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(cut); + memory->destroy(viscosity); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairSPHIdealGas::compute(int eflag, int vflag) { + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, delx, dely, delz, fpair; + + int *ilist, *jlist, *numneigh, **firstneigh; + double vxtmp, vytmp, vztmp, imass, jmass, fi, fj, fvisc, h, ih, ihsq; + double rsq, wfd, delVdotDelR, mu, deltaE, ci, cj; + + if (eflag || vflag) + ev_setup(eflag, vflag); + else + evflag = vflag_fdotr = 0; + + double **v = atom->vest; + double **x = atom->x; + double **f = atom->f; + double *rho = atom->rho; + double *mass = atom->mass; + double *de = atom->de; + double *e = atom->e; + double *drho = atom->drho; + int *type = atom->type; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + vxtmp = v[i][0]; + vytmp = v[i][1]; + vztmp = v[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + imass = mass[itype]; + + fi = 0.4 * e[i] / imass / rho[i]; // ideal gas EOS; this expression is fi = pressure / rho^2 + ci = sqrt(0.4*e[i]/imass); // speed of sound with heat capacity ratio gamma=1.4 + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + jtype = type[j]; + jmass = mass[jtype]; + + if (rsq < cutsq[itype][jtype]) { + h = cut[itype][jtype]; + ih = 1. / h; + ihsq = ih * ih; + + wfd = h - sqrt(rsq); + if (domain->dimension == 3) { + // Lucy Kernel, 3d + // Note that wfd, the derivative of the weight function with respect to r, + // is lacking a factor of r. + // The missing factor of r is recovered by + // (1) using delV . delX instead of delV . (delX/r) and + // (2) using f[i][0] += delx * fpair instead of f[i][0] += (delx/r) * fpair + wfd = -25.066903536973515383e0 * wfd * wfd * ihsq * ihsq * ihsq * ih; + } else { + // Lucy Kernel, 2d + wfd = -19.098593171027440292e0 * wfd * wfd * ihsq * ihsq * ihsq; + } + + fj = 0.4 * e[j] / jmass / rho[j]; + + // dot product of velocity delta and distance vector + delVdotDelR = delx * (vxtmp - v[j][0]) + dely * (vytmp - v[j][1]) + + delz * (vztmp - v[j][2]); + + // artificial viscosity (Monaghan 1992) + if (delVdotDelR < 0.) { + cj = sqrt(0.4*e[j]/jmass); + mu = h * delVdotDelR / (rsq + 0.01 * h * h); + fvisc = -viscosity[itype][jtype] * (ci + cj) * mu / (rho[i] + rho[j]); + } else { + fvisc = 0.; + } + + // total pair force & thermal energy increment + fpair = -imass * jmass * (fi + fj + fvisc) * wfd; + deltaE = -0.5 * fpair * delVdotDelR; + + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; + + // and change in density + drho[i] += jmass * delVdotDelR * wfd; + + // change in thermal energy + de[i] += deltaE; + + if (newton_pair || j < nlocal) { + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; + de[j] += deltaE; + drho[j] += imass * delVdotDelR * wfd; + } + + if (evflag) + ev_tally(i, j, nlocal, newton_pair, 0.0, 0.0, fpair, delx, dely, delz); + + if (evflag) + ev_tally(i, j, nlocal, newton_pair, 0.0, 0.0, fpair, delx, dely, + delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- + allocate all arrays + ------------------------------------------------------------------------- */ + +void PairSPHIdealGas::allocate() { + allocated = 1; + int n = atom->ntypes; + + 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; + + memory->create(cutsq, n + 1, n + 1, "pair:cutsq"); + + memory->create(cut, n + 1, n + 1, "pair:cut"); + memory->create(viscosity, n + 1, n + 1, "pair:viscosity"); +} + +/* ---------------------------------------------------------------------- + global settings + ------------------------------------------------------------------------- */ + +void PairSPHIdealGas::settings(int narg, char **arg) { + if (narg != 0) + error->all( + "Illegal number of setting arguments for pair_style sph/idealgas"); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs + ------------------------------------------------------------------------- */ + +void PairSPHIdealGas::coeff(int narg, char **arg) { + if (narg != 4) + error->all("Incorrect number of args for pair_style sph/idealgas coefficients"); + if (!allocated) + allocate(); + + int ilo, ihi, jlo, jhi; + force->bounds(arg[0], atom->ntypes, ilo, ihi); + force->bounds(arg[1], atom->ntypes, jlo, jhi); + + double viscosity_one = force->numeric(arg[2]); + double cut_one = force->numeric(arg[3]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + viscosity[i][j] = viscosity_one; + //printf("setting cut[%d][%d] = %f\n", i, j, cut_one); + cut[i][j] = cut_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) + error->all("Incorrect args for pair sph/idealgas coefficients"); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i + ------------------------------------------------------------------------- */ + +double PairSPHIdealGas::init_one(int i, int j) { + + if (setflag[i][j] == 0) { + error->all("All pair sph/idealgas coeffs are not set"); + } + + cut[j][i] = cut[i][j]; + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- */ + +double PairSPHIdealGas::single(int i, int j, int itype, int jtype, + double rsq, double factor_coul, double factor_lj, double &fforce) { + fforce = 0.0; + + return 0.0; +} diff --git a/src/USER-SPH/pair_sph_idealgas.h b/src/USER-SPH/pair_sph_idealgas.h new file mode 100644 index 0000000000..540789afeb --- /dev/null +++ b/src/USER-SPH/pair_sph_idealgas.h @@ -0,0 +1,46 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(sph/idealgas,PairSPHIdealGas) + +#else + +#ifndef LMP_PAIR_IDEALGAS_H +#define LMP_PAIR_IDEALGAS_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairSPHIdealGas : public Pair { + public: + PairSPHIdealGas(class LAMMPS *); + virtual ~PairSPHIdealGas(); + virtual void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + virtual double init_one(int, int); + virtual double single(int, int, int, int, double, double, double, double &); + + protected: + double **cut,**viscosity; + + void allocate(); +}; + +} + +#endif +#endif diff --git a/src/USER-SPH/pair_sph_lj.cpp b/src/USER-SPH/pair_sph_lj.cpp new file mode 100644 index 0000000000..fd92cddc8f --- /dev/null +++ b/src/USER-SPH/pair_sph_lj.cpp @@ -0,0 +1,368 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. + ------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdlib.h" +#include "pair_sph_lj.h" +#include "atom.h" +#include "force.h" +#include "comm.h" +#include "neigh_list.h" +#include "memory.h" +#include "error.h" +#include "domain.h" + +using namespace LAMMPS_NS; + +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#define MAX(a,b) ((a) > (b) ? (a) : (b)) + +/* ---------------------------------------------------------------------- */ + +PairSPHLJ::PairSPHLJ(LAMMPS *lmp) : + Pair(lmp) { +} + +/* ---------------------------------------------------------------------- */ + +PairSPHLJ::~PairSPHLJ() { + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(cut); + memory->destroy(viscosity); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairSPHLJ::compute(int eflag, int vflag) { + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, delx, dely, delz, fpair; + + int *ilist, *jlist, *numneigh, **firstneigh; + double vxtmp, vytmp, vztmp, imass, jmass, fi, fj, fvisc, h, ih, ihsq, ihcub; + double rsq, wfd, delVdotDelR, mu, deltaE, ci, cj, lrc; + + if (eflag || vflag) + ev_setup(eflag, vflag); + else + evflag = vflag_fdotr = 0; + + double **v = atom->vest; + double **x = atom->x; + double **f = atom->f; + double *rho = atom->rho; + double *mass = atom->mass; + double *de = atom->de; + double *e = atom->e; + double *cv = atom->cv; + double *drho = atom->drho; + int *type = atom->type; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + vxtmp = v[i][0]; + vytmp = v[i][1]; + vztmp = v[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + imass = mass[itype]; + + // compute pressure of particle i with LJ EOS + LJEOS2(rho[i], e[i], cv[i], &fi, &ci); + fi /= (rho[i] * rho[i]); + //printf("fi = %f\n", fi); + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + jtype = type[j]; + jmass = mass[jtype]; + + if (rsq < cutsq[itype][jtype]) { + h = cut[itype][jtype]; + ih = 1.0 / h; + ihsq = ih * ih; + ihcub = ihsq * ih; + + wfd = h - sqrt(rsq); + if (domain->dimension == 3) { + // Lucy Kernel, 3d + // Note that wfd, the derivative of the weight function with respect to r, + // is lacking a factor of r. + // The missing factor of r is recovered by + // (1) using delV . delX instead of delV . (delX/r) and + // (2) using f[i][0] += delx * fpair instead of f[i][0] += (delx/r) * fpair + wfd = -25.066903536973515383e0 * wfd * wfd * ihsq * ihsq * ihsq * ih; + } else { + // Lucy Kernel, 2d + wfd = -19.098593171027440292e0 * wfd * wfd * ihsq * ihsq * ihsq; + } + + // function call to LJ EOS + LJEOS2(rho[j], e[j], cv[j], &fj, &cj); + fj /= (rho[j] * rho[j]); + + // apply long-range correction to model a LJ fluid with cutoff + // this implies that the modelled LJ fluid has cutoff == SPH cutoff + lrc = - 11.1701 * (ihcub * ihcub * ihcub - 1.5 * ihcub); + fi += lrc; + fj += lrc; + + // dot product of velocity delta and distance vector + delVdotDelR = delx * (vxtmp - v[j][0]) + dely * (vytmp - v[j][1]) + + delz * (vztmp - v[j][2]); + + // artificial viscosity (Monaghan 1992) + if (delVdotDelR < 0.) { + mu = h * delVdotDelR / (rsq + 0.01 * h * h); + fvisc = -viscosity[itype][jtype] * (ci + cj) * mu / (rho[i] + rho[j]); + } else { + fvisc = 0.; + } + + // total pair force & thermal energy increment + fpair = -imass * jmass * (fi + fj + fvisc) * wfd; + deltaE = -0.5 * fpair * delVdotDelR; + + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; + + // and change in density + drho[i] += jmass * delVdotDelR * wfd; + + // change in thermal energy + de[i] += deltaE; + + if (newton_pair || j < nlocal) { + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; + de[j] += deltaE; + drho[j] += imass * delVdotDelR * wfd; + } + + if (evflag) + ev_tally(i, j, nlocal, newton_pair, 0.0, 0.0, fpair, delx, dely, delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- + allocate all arrays + ------------------------------------------------------------------------- */ + +void PairSPHLJ::allocate() { + allocated = 1; + int n = atom->ntypes; + + 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; + + memory->create(cutsq, n + 1, n + 1, "pair:cutsq"); + + memory->create(cut, n + 1, n + 1, "pair:cut"); + memory->create(viscosity, n + 1, n + 1, "pair:viscosity"); +} + +/* ---------------------------------------------------------------------- + global settings + ------------------------------------------------------------------------- */ + +void PairSPHLJ::settings(int narg, char **arg) { + if (narg != 0) + error->all( + "Illegal number of setting arguments for pair_style sph/lj"); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs + ------------------------------------------------------------------------- */ + +void PairSPHLJ::coeff(int narg, char **arg) { + if (narg != 4) + error->all( + "Incorrect args for pair_style sph/lj coefficients"); + if (!allocated) + allocate(); + + int ilo, ihi, jlo, jhi; + force->bounds(arg[0], atom->ntypes, ilo, ihi); + force->bounds(arg[1], atom->ntypes, jlo, jhi); + + double viscosity_one = force->numeric(arg[2]); + double cut_one = force->numeric(arg[3]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + viscosity[i][j] = viscosity_one; + printf("setting cut[%d][%d] = %f\n", i, j, cut_one); + cut[i][j] = cut_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) + error->all("Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i + ------------------------------------------------------------------------- */ + +double PairSPHLJ::init_one(int i, int j) { + + if (setflag[i][j] == 0) { + error->all("All pair sph/lj coeffs are not set"); + } + + cut[j][i] = cut[i][j]; + viscosity[j][i] = viscosity[i][j]; + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- */ + +double PairSPHLJ::single(int i, int j, int itype, int jtype, + double rsq, double factor_coul, double factor_lj, double &fforce) { + fforce = 0.0; + + return 0.0; +} + + +/*double PairSPHLJ::LJEOS2(double rho, double e, double cv) { + + + double T = e / cv; + if (T < 1.e-2) T = 1.e-2; + //printf("%f %f\n", T, rho); + double iT = 0.1e1 / T; + //double itpow1_4 = exp(0.25 * log(iT)); //pow(iT, 0.1e1 / 0.4e1); + double itpow1_4 = pow(iT, 0.1e1 / 0.4e1); + double x = rho * itpow1_4; + double xsq = x * x; + double xpow3 = xsq * x; + double xpow4 = xsq * xsq; + double xpow9 = xpow3 * xpow3 * xpow3; + + + return (0.1e1 + rho * (0.3629e1 + 0.7264e1 * x + 0.104925e2 * xsq + 0.11460e2 + * xpow3 + 0.21760e1 * xpow9 - itpow1_4 * itpow1_4 * (0.5369e1 + 0.13160e2 + * x + 0.18525e2 * xsq - 0.17076e2 * xpow3 + 0.9320e1 * xpow4) + iT + * (-0.3492e1 + 0.18698e2 * x - 0.35505e2 * xsq + 0.31816e2 * xpow3 + - 0.11195e2 * xpow4)) * itpow1_4) * rho * T; +}*/ + + +/* --------------------------------------------------------------------------------------------- */ +/* Lennard-Jones EOS, + Francis H. Ree + "Analytic representation of thermodynamic data for the Lennard‐Jones fluid", + Journal of Chemical Physics 73 pp. 5401-5403 (1980) +*/ + +void PairSPHLJ::LJEOS2(double rho, double e, double cv, double *p, double *c) { + double T = e/cv; + double beta = 1.0 / T; + double beta_sqrt = sqrt(beta); + double x = rho * sqrt(beta_sqrt); + + double xsq = x * x; + double xpow3 = xsq * x; + double xpow4 = xsq * xsq; + + /* differential of Helmholtz free energy w.r.t. x */ + double diff_A_NkT = 3.629 + 7.264*x - beta*(3.492 - 18.698*x + 35.505*xsq - 31.816*xpow3 + 11.195*xpow4) + - beta_sqrt*(5.369 + 13.16*x + 18.525*xsq - 17.076*xpow3 + 9.32*xpow4) + + 10.4925*xsq + 11.46*xpow3 + 2.176*xpow4*xpow4*x; + + /* differential of Helmholtz free energy w.r.t. x^2 */ + double d2A_dx2 = 7.264 + 20.985*x \ + + beta*(18.698 - 71.01*x + 95.448*xsq - 44.78*xpow3)\ + - beta_sqrt*(13.16 + 37.05*x - 51.228*xsq + 37.28*xpow3)\ + + 34.38*xsq + 19.584*xpow4*xpow4; + + // p = rho k T * (1 + rho * d(A/(NkT))/drho) + // dx/drho = rho/x + *p = rho * T * (1.0 + diff_A_NkT * x); // pressure + double csq = T * (1.0 + 2.0 * diff_A_NkT * x + d2A_dx2 * x * x); // soundspeed squared + if (csq > 0.0) { + *c = sqrt(csq); // soundspeed + } else { + *c = 0.0; + } +} + +/* ------------------------------------------------------------------------------ */ + +/* Jirí Kolafa, Ivo Nezbeda + * "The Lennard-Jones fluid: an accurate analytic and theoretically-based equation of state", + * Fluid Phase Equilibria 100 pp. 1-34 (1994) */ +/*double PairSPHLJ::LJEOS2(double rho, double e, double cv) { + double T = e / cv; + + double sT = sqrt(T); + double isT = 1.0 / sT; + double dC = -0.063920968 * log(T) + 0.011117524 / T - 0.076383859 / sT + + 1.080142248 + 0.000693129 * sT; + double eta = 3.141592654 / 6. * rho * (dC * dC * dC); + double zHS = (1 + eta * (1 + eta * (1 - eta / 1.5 * (1 + eta)))) + / ((1. - eta) * (1. - eta) * (1. - eta)); + double BC = (((((-0.58544978 * isT + 0.43102052) * isT + .87361369) * isT + - 4.13749995) * isT + 2.90616279) * isT - 7.02181962) / T + 0.02459877; + double gammaBH = 1.92907278; + + double sum = ((2.01546797 * 2 + rho * ((-28.17881636) * 3 + rho + * (28.28313847 * 4 + rho * (-10.42402873) * 5))) + (-19.58371655 * 2 + + rho * (+75.62340289 * 3 + rho * ((-120.70586598) * 4 + rho + * (+93.92740328 * 5 + rho * (-27.37737354) * 6)))) / sqrt(T) + + ((29.34470520 * 2 + rho * ((-112.35356937) * 3 + rho * (+170.64908980 + * 4 + rho * ((-123.06669187) * 5 + rho * 34.42288969 * 6)))) + + ((-13.37031968) * 2 + rho * (65.38059570 * 3 + rho + * ((-115.09233113) * 4 + rho * (88.91973082 * 5 + rho + * (-25.62099890) * 6)))) / T) / T) * rho * rho; + return ((zHS + BC / exp(gammaBH * rho * rho) * rho * (1 - 2 * gammaBH * rho + * rho)) * T + sum) * rho; + } +*/ diff --git a/src/USER-SPH/pair_sph_lj.h b/src/USER-SPH/pair_sph_lj.h new file mode 100644 index 0000000000..244ee044a5 --- /dev/null +++ b/src/USER-SPH/pair_sph_lj.h @@ -0,0 +1,48 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(sph/lj,PairSPHLJ) + +#else + +#ifndef LMP_PAIR_LJ_H +#define LMP_PAIR_LJ_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairSPHLJ : public Pair { + public: + PairSPHLJ(class LAMMPS *); + virtual ~PairSPHLJ(); + virtual void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + virtual double init_one(int, int); + virtual double single(int, int, int, int, double, double, double, double &); + //double LJEOS(int); + void LJEOS2(double, double, double, double *, double *); + + protected: + double **cut,**viscosity; + + void allocate(); +}; + +} + +#endif +#endif diff --git a/src/USER-SPH/pair_sph_rhosum.cpp b/src/USER-SPH/pair_sph_rhosum.cpp new file mode 100644 index 0000000000..b846b1bb49 --- /dev/null +++ b/src/USER-SPH/pair_sph_rhosum.cpp @@ -0,0 +1,316 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. + ------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdlib.h" +#include "pair_sph_rhosum.h" +#include "atom.h" +#include "force.h" +#include "comm.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "memory.h" +#include "error.h" +#include "neighbor.h" +#include "update.h" +#include "domain.h" + +using namespace LAMMPS_NS; + +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#define MAX(a,b) ((a) > (b) ? (a) : (b)) + +/* ---------------------------------------------------------------------- */ + +PairSPHRhoSum::PairSPHRhoSum(LAMMPS *lmp) : + Pair(lmp) { + + // set comm size needed by this Pair + + comm_forward = 1; + first = 1; +} + +/* ---------------------------------------------------------------------- */ + +PairSPHRhoSum::~PairSPHRhoSum() { + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(cut); + } +} + +/* ---------------------------------------------------------------------- + init specific to this pair style + ------------------------------------------------------------------------- */ + +void PairSPHRhoSum::init_style() { + // need a full neighbor list + int irequest = neighbor->request(this); + neighbor->requests[irequest]->half = 0; + neighbor->requests[irequest]->full = 1; +} + +/* ---------------------------------------------------------------------- */ + +void PairSPHRhoSum::compute(int eflag, int vflag) { + int i, j, ii, jj, jnum, itype, jtype; + double xtmp, ytmp, ztmp, delx, dely, delz; + double r, rsq, imass, h, ih, ihsq; + int *jlist; + double wf; + // neighbor list variables + int inum, *ilist, *numneigh, **firstneigh; + + if (eflag || vflag) + ev_setup(eflag, vflag); + else + evflag = vflag_fdotr = 0; + + double **x = atom->x; + double *rho = atom->rho; + int *type = atom->type; + double *mass = atom->mass; + + // check consistency of pair coefficients + + if (first) { + for (i = 1; i <= atom->ntypes; i++) { + for (j = 1; i <= atom->ntypes; i++) { + if (cutsq[i][j] > 0.0) { + if (!setflag[i][i] || !setflag[j][j]) { + if (comm->me == 0) { + printf( + "SPH particle types %d and %d interact, but not all of their single particle properties are set.\n", + i, j); + } + } + } + } + } + first = 0; + } + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // recompute density + // we use a full neighborlist here + + if (nstep != 0) { + if ((update->ntimestep % nstep) == 0) { + + // initialize density with self-contribution, + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + itype = type[i]; + imass = mass[itype]; + + h = cut[itype][itype]; + if (domain->dimension == 3) { + /* + // Lucy kernel, 3d + wf = 2.0889086280811262819e0 / (h * h * h); + */ + + // quadric kernel, 3d + wf = 2.1541870227086614782 / (h * h * h); + } else { + /* + // Lucy kernel, 2d + wf = 1.5915494309189533576e0 / (h * h); + */ + + // quadric kernel, 2d + wf = 1.5915494309189533576e0 / (h * h); + } + + rho[i] = imass * wf; + } + + // add density at each atom via kernel function overlap + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + jtype = type[j]; + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + + if (rsq < cutsq[itype][jtype]) { + h = cut[itype][jtype]; + ih = 1.0 / h; + ihsq = ih * ih; + + if (domain->dimension == 3) { + /* + // Lucy kernel, 3d + r = sqrt(rsq); + wf = (h - r) * ihsq; + wf = 2.0889086280811262819e0 * (h + 3. * r) * wf * wf * wf * ih; + */ + + // quadric kernel, 3d + wf = 1.0 - rsq * ihsq; + wf = wf * wf; + wf = wf * wf; + wf = 2.1541870227086614782e0 * wf * ihsq * ih; + } else { + // Lucy kernel, 2d + //r = sqrt(rsq); + //wf = (h - r) * ihsq; + //wf = 1.5915494309189533576e0 * (h + 3. * r) * wf * wf * wf; + + // quadric kernel, 2d + wf = 1.0 - rsq * ihsq; + wf = wf * wf; + wf = wf * wf; + wf = 1.5915494309189533576e0 * wf * ihsq; + } + + rho[i] += mass[jtype] * wf; + } + + } + } + } + } + + // communicate densities + comm->forward_comm_pair(this); +} + +/* ---------------------------------------------------------------------- + allocate all arrays + ------------------------------------------------------------------------- */ + +void PairSPHRhoSum::allocate() { + allocated = 1; + int n = atom->ntypes; + + 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; + + memory->create(cutsq, n + 1, n + 1, "pair:cutsq"); + + memory->create(cut, n + 1, n + 1, "pair:cut"); +} + +/* ---------------------------------------------------------------------- + global settings + ------------------------------------------------------------------------- */ + +void PairSPHRhoSum::settings(int narg, char **arg) { + if (narg != 1) + error->all( + "Illegal number of setting arguments for pair_style sph/rhosum"); + nstep = force->inumeric(arg[0]); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs + ------------------------------------------------------------------------- */ + +void PairSPHRhoSum::coeff(int narg, char **arg) { + if (narg != 3) + error->all("Incorrect number of args for sph/rhosum coefficients"); + if (!allocated) + allocate(); + + int ilo, ihi, jlo, jhi; + force->bounds(arg[0], atom->ntypes, ilo, ihi); + force->bounds(arg[1], atom->ntypes, jlo, jhi); + + double cut_one = force->numeric(arg[2]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + //printf("setting cut[%d][%d] = %f\n", i, j, cut_one); + cut[i][j] = cut_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) + error->all("Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i + ------------------------------------------------------------------------- */ + +double PairSPHRhoSum::init_one(int i, int j) { + if (setflag[i][j] == 0) { + error->all("All pair sph/rhosum coeffs are not set"); + } + + cut[j][i] = cut[i][j]; + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- */ + +double PairSPHRhoSum::single(int i, int j, int itype, int jtype, double rsq, + double factor_coul, double factor_lj, double &fforce) { + fforce = 0.0; + + return 0.0; +} + +/* ---------------------------------------------------------------------- */ + +int PairSPHRhoSum::pack_comm(int n, int *list, double *buf, int pbc_flag, + int *pbc) { + int i, j, m; + double *rho = atom->rho; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = rho[j]; + } + return 1; +} + +/* ---------------------------------------------------------------------- */ + +void PairSPHRhoSum::unpack_comm(int n, int first, double *buf) { + int i, m, last; + double *rho = atom->rho; + + m = 0; + last = first + n; + for (i = first; i < last; i++) + rho[i] = buf[m++]; +} + diff --git a/src/USER-SPH/pair_sph_rhosum.h b/src/USER-SPH/pair_sph_rhosum.h new file mode 100644 index 0000000000..227cdef95c --- /dev/null +++ b/src/USER-SPH/pair_sph_rhosum.h @@ -0,0 +1,50 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(sph/rhosum,PairSPHRhoSum) + +#else + +#ifndef LMP_PAIR_SPH_RHOSUM_H +#define LMP_PAIR_SPH_RHOSUM_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairSPHRhoSum : public Pair { + public: + PairSPHRhoSum(class LAMMPS *); + virtual ~PairSPHRhoSum(); + void init_style(); + virtual void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + virtual double init_one(int, int); + virtual double single(int, int, int, int, double, double, double, double &); + int pack_comm(int, int *, double *, int, int *); + void unpack_comm(int, int, double *); + + protected: + double **cut; + int nstep, first; + + void allocate(); +}; + +} + +#endif +#endif diff --git a/src/USER-SPH/pair_sph_taitwater.cpp b/src/USER-SPH/pair_sph_taitwater.cpp new file mode 100644 index 0000000000..0226d8fed4 --- /dev/null +++ b/src/USER-SPH/pair_sph_taitwater.cpp @@ -0,0 +1,303 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. + ------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdlib.h" +#include "pair_sph_taitwater.h" +#include "atom.h" +#include "force.h" +#include "comm.h" +#include "neigh_list.h" +#include "memory.h" +#include "error.h" +#include "domain.h" + +using namespace LAMMPS_NS; + +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#define MAX(a,b) ((a) > (b) ? (a) : (b)) + +/* ---------------------------------------------------------------------- */ + +PairSPHTaitwater::PairSPHTaitwater(LAMMPS *lmp) : + Pair(lmp) { + + first = 1; +} + +/* ---------------------------------------------------------------------- */ + +PairSPHTaitwater::~PairSPHTaitwater() { + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(cut); + memory->destroy(rho0); + memory->destroy(soundspeed); + memory->destroy(B); + memory->destroy(viscosity); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairSPHTaitwater::compute(int eflag, int vflag) { + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, delx, dely, delz, fpair; + + int *ilist, *jlist, *numneigh, **firstneigh; + double vxtmp, vytmp, vztmp, imass, jmass, fi, fj, fvisc, h, ih, ihsq; + double rsq, tmp, wfd, delVdotDelR, mu, deltaE; + + if (eflag || vflag) + ev_setup(eflag, vflag); + else + evflag = vflag_fdotr = 0; + + double **v = atom->vest; + double **x = atom->x; + double **f = atom->f; + double *rho = atom->rho; + double *mass = atom->mass; + double *de = atom->de; + double *drho = atom->drho; + int *type = atom->type; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + + // check consistency of pair coefficients + + if (first) { + for (i = 1; i <= atom->ntypes; i++) { + for (j = 1; i <= atom->ntypes; i++) { + if (cutsq[i][j] > 1.e-32) { + if (!setflag[i][i] || !setflag[j][j]) { + if (comm->me == 0) { + printf( + "SPH particle types %d and %d interact with cutoff=%g, but not all of their single particle properties are set.\n", + i, j, sqrt(cutsq[i][j])); + } + } + } + } + } + first = 0; + } + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + vxtmp = v[i][0]; + vytmp = v[i][1]; + vztmp = v[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + imass = mass[itype]; + + // compute pressure of atom i with Tait EOS + tmp = rho[i] / rho0[itype]; + fi = tmp * tmp * tmp; + fi = B[itype] * (fi * fi * tmp - 1.0) / (rho[i] * rho[i]); + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + jtype = type[j]; + jmass = mass[jtype]; + + if (rsq < cutsq[itype][jtype]) { + + h = cut[itype][jtype]; + ih = 1.0 / h; + ihsq = ih * ih; + + wfd = h - sqrt(rsq); + if (domain->dimension == 3) { + // Lucy Kernel, 3d + // Note that wfd, the derivative of the weight function with respect to r, + // is lacking a factor of r. + // The missing factor of r is recovered by + // (1) using delV . delX instead of delV . (delX/r) and + // (2) using f[i][0] += delx * fpair instead of f[i][0] += (delx/r) * fpair + wfd = -25.066903536973515383e0 * wfd * wfd * ihsq * ihsq * ihsq * ih; + } else { + // Lucy Kernel, 2d + wfd = -19.098593171027440292e0 * wfd * wfd * ihsq * ihsq * ihsq; + } + + // compute pressure of atom j with Tait EOS + tmp = rho[j] / rho0[jtype]; + fj = tmp * tmp * tmp; + fj = B[jtype] * (fj * fj * tmp - 1.0) / (rho[j] * rho[j]); + + // dot product of velocity delta and distance vector + delVdotDelR = delx * (vxtmp - v[j][0]) + dely * (vytmp - v[j][1]) + + delz * (vztmp - v[j][2]); + + // artificial viscosity (Monaghan 1992) + if (delVdotDelR < 0.) { + mu = h * delVdotDelR / (rsq + 0.01 * h * h); + fvisc = -viscosity[itype][jtype] * (soundspeed[itype] + + soundspeed[jtype]) * mu / (rho[i] + rho[j]); + } else { + fvisc = 0.; + } + + // total pair force & thermal energy increment + fpair = -imass * jmass * (fi + fj + fvisc) * wfd; + deltaE = -0.5 * fpair * delVdotDelR; + + f[i][0] += delx * fpair; + f[i][1] += dely * fpair; + f[i][2] += delz * fpair; + + // and change in density + drho[i] += jmass * delVdotDelR * wfd; + + // change in thermal energy + de[i] += deltaE; + + if (newton_pair || j < nlocal) { + f[j][0] -= delx * fpair; + f[j][1] -= dely * fpair; + f[j][2] -= delz * fpair; + de[j] += deltaE; + drho[j] += imass * delVdotDelR * wfd; + } + + if (evflag) + ev_tally(i, j, nlocal, newton_pair, 0.0, 0.0, fpair, delx, dely, delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- + allocate all arrays + ------------------------------------------------------------------------- */ + +void PairSPHTaitwater::allocate() { + allocated = 1; + int n = atom->ntypes; + + 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; + + memory->create(cutsq, n + 1, n + 1, "pair:cutsq"); + + memory->create(rho0, n + 1, "pair:rho0"); + memory->create(soundspeed, n + 1, "pair:soundspeed"); + memory->create(B, n + 1, "pair:B"); + memory->create(cut, n + 1, n + 1, "pair:cut"); + memory->create(viscosity, n + 1, n + 1, "pair:viscosity"); +} + +/* ---------------------------------------------------------------------- + global settings + ------------------------------------------------------------------------- */ + +void PairSPHTaitwater::settings(int narg, char **arg) { + if (narg != 0) + error->all( + "Illegal number of setting arguments for pair_style sph/taitwater"); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs + ------------------------------------------------------------------------- */ + +void PairSPHTaitwater::coeff(int narg, char **arg) { + if (narg != 6) + error->all( + "Incorrect args for pair_style sph/taitwater coefficients"); + if (!allocated) + allocate(); + + int ilo, ihi, jlo, jhi; + force->bounds(arg[0], atom->ntypes, ilo, ihi); + force->bounds(arg[1], atom->ntypes, jlo, jhi); + + double rho0_one = force->numeric(arg[2]); + double soundspeed_one = force->numeric(arg[3]); + double viscosity_one = force->numeric(arg[4]); + double cut_one = force->numeric(arg[5]); + double B_one = soundspeed_one * soundspeed_one * rho0_one / 7.0; + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + rho0[i] = rho0_one; + soundspeed[i] = soundspeed_one; + B[i] = B_one; + for (int j = MAX(jlo,i); j <= jhi; j++) { + viscosity[i][j] = viscosity_one; + //printf("setting cut[%d][%d] = %f\n", i, j, cut_one); + cut[i][j] = cut_one; + + setflag[i][j] = 1; + + //cut[j][i] = cut[i][j]; + //viscosity[j][i] = viscosity[i][j]; + //setflag[j][i] = 1; + count++; + } + } + + if (count == 0) + error->all("Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i + ------------------------------------------------------------------------- */ + +double PairSPHTaitwater::init_one(int i, int j) { + + if (setflag[i][j] == 0) { + error->all("Not all pair sph/taitwater coeffs are set"); + } + + cut[j][i] = cut[i][j]; + viscosity[j][i] = viscosity[i][j]; + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- */ + +double PairSPHTaitwater::single(int i, int j, int itype, int jtype, + double rsq, double factor_coul, double factor_lj, double &fforce) { + fforce = 0.0; + + return 0.0; +} diff --git a/src/USER-SPH/pair_sph_taitwater.h b/src/USER-SPH/pair_sph_taitwater.h new file mode 100644 index 0000000000..de081aa544 --- /dev/null +++ b/src/USER-SPH/pair_sph_taitwater.h @@ -0,0 +1,48 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(sph/taitwater,PairSPHTaitwater) + +#else + +#ifndef LMP_PAIR_TAITWATER_H +#define LMP_PAIR_TAITWATER_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairSPHTaitwater : public Pair { + public: + PairSPHTaitwater(class LAMMPS *); + virtual ~PairSPHTaitwater(); + virtual void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + virtual double init_one(int, int); + virtual double single(int, int, int, int, double, double, double, double &); + + protected: + double *rho0, *soundspeed, *B; + double **cut,**viscosity; + int first; + + void allocate(); +}; + +} + +#endif +#endif diff --git a/src/USER-SPH/pair_sph_taitwater_morris.cpp b/src/USER-SPH/pair_sph_taitwater_morris.cpp new file mode 100644 index 0000000000..5e5eb1d098 --- /dev/null +++ b/src/USER-SPH/pair_sph_taitwater_morris.cpp @@ -0,0 +1,300 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. + ------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdlib.h" +#include "pair_sph_taitwater_morris.h" +#include "atom.h" +#include "force.h" +#include "comm.h" +#include "neigh_list.h" +#include "memory.h" +#include "error.h" +#include "domain.h" + +using namespace LAMMPS_NS; + +#define MIN(a,b) ((a) < (b) ? (a) : (b)) +#define MAX(a,b) ((a) > (b) ? (a) : (b)) + +/* ---------------------------------------------------------------------- */ + +PairSPHTaitwaterMorris::PairSPHTaitwaterMorris(LAMMPS *lmp) : + Pair(lmp) { + + first = 1; +} + +/* ---------------------------------------------------------------------- */ + +PairSPHTaitwaterMorris::~PairSPHTaitwaterMorris() { + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + + memory->destroy(cut); + memory->destroy(rho0); + memory->destroy(soundspeed); + memory->destroy(B); + memory->destroy(viscosity); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairSPHTaitwaterMorris::compute(int eflag, int vflag) { + int i, j, ii, jj, inum, jnum, itype, jtype; + double xtmp, ytmp, ztmp, delx, dely, delz, fpair; + + int *ilist, *jlist, *numneigh, **firstneigh; + double vxtmp, vytmp, vztmp, imass, jmass, fi, fj, fvisc, h, ih, ihsq, velx, vely, velz; + double rsq, tmp, wfd, delVdotDelR, deltaE; + + if (eflag || vflag) + ev_setup(eflag, vflag); + else + evflag = vflag_fdotr = 0; + + double **v = atom->vest; + double **x = atom->x; + double **f = atom->f; + double *rho = atom->rho; + double *mass = atom->mass; + double *de = atom->de; + double *drho = atom->drho; + int *type = atom->type; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + + // check consistency of pair coefficients + + if (first) { + for (i = 1; i <= atom->ntypes; i++) { + for (j = 1; i <= atom->ntypes; i++) { + if (cutsq[i][j] > 1.e-32) { + if (!setflag[i][i] || !setflag[j][j]) { + if (comm->me == 0) { + printf( + "SPH particle types %d and %d interact with cutoff=%g, but not all of their single particle properties are set.\n", + i, j, sqrt(cutsq[i][j])); + } + } + } + } + } + first = 0; + } + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + vxtmp = v[i][0]; + vytmp = v[i][1]; + vztmp = v[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + imass = mass[itype]; + + // compute pressure of atom i with Tait EOS + tmp = rho[i] / rho0[itype]; + fi = tmp * tmp * tmp; + fi = B[itype] * (fi * fi * tmp - 1.0) / (rho[i] * rho[i]); + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + jtype = type[j]; + jmass = mass[jtype]; + + if (rsq < cutsq[itype][jtype]) { + h = cut[itype][jtype]; + ih = 1.0 / h; + ihsq = ih * ih; + + wfd = h - sqrt(rsq); + if (domain->dimension == 3) { + // Lucy Kernel, 3d + // Note that wfd, the derivative of the weight function with respect to r, + // is lacking a factor of r. + // The missing factor of r is recovered by + // (1) using delV . delX instead of delV . (delX/r) and + // (2) using f[i][0] += delx * fpair instead of f[i][0] += (delx/r) * fpair + wfd = -25.066903536973515383e0 * wfd * wfd * ihsq * ihsq * ihsq * ih; + } else { + // Lucy Kernel, 2d + wfd = -19.098593171027440292e0 * wfd * wfd * ihsq * ihsq * ihsq; + } + + // compute pressure of atom j with Tait EOS + tmp = rho[j] / rho0[jtype]; + fj = tmp * tmp * tmp; + fj = B[jtype] * (fj * fj * tmp - 1.0) / (rho[j] * rho[j]); + + velx=vxtmp - v[j][0]; + vely=vytmp - v[j][1]; + velz=vztmp - v[j][2]; + + // dot product of velocity delta and distance vector + delVdotDelR = delx * velx + dely * vely + delz * velz; + + // Morris Viscosity (Morris, 1996) + + fvisc = 2 * viscosity[itype][jtype] / (rho[i] * rho[j]); + + fvisc *= imass * jmass * wfd; + + // total pair force & thermal energy increment + fpair = -imass * jmass * (fi + fj) * wfd; + deltaE = -0.5 *(fpair * delVdotDelR + fvisc * (velx*velx + vely*vely + velz*velz)); + + // printf("testvar= %f, %f \n", delx, dely); + + f[i][0] += delx * fpair + velx * fvisc; + f[i][1] += dely * fpair + vely * fvisc; + f[i][2] += delz * fpair + velz * fvisc; + + // and change in density + drho[i] += jmass * delVdotDelR * wfd; + + // change in thermal energy + de[i] += deltaE; + + if (newton_pair || j < nlocal) { + f[j][0] -= delx * fpair + velx * fvisc; + f[j][1] -= dely * fpair + vely * fvisc; + f[j][2] -= delz * fpair + velz * fvisc; + de[j] += deltaE; + drho[j] += imass * delVdotDelR * wfd; + } + + if (evflag) + ev_tally(i, j, nlocal, newton_pair, 0.0, 0.0, fpair, delx, dely, delz); + } + } + } + + if (vflag_fdotr) virial_fdotr_compute(); +} + +/* ---------------------------------------------------------------------- + allocate all arrays + ------------------------------------------------------------------------- */ + +void PairSPHTaitwaterMorris::allocate() { + allocated = 1; + int n = atom->ntypes; + + 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; + + memory->create(cutsq, n + 1, n + 1, "pair:cutsq"); + + memory->create(rho0, n + 1, "pair:rho0"); + memory->create(soundspeed, n + 1, "pair:soundspeed"); + memory->create(B, n + 1, "pair:B"); + memory->create(cut, n + 1, n + 1, "pair:cut"); + memory->create(viscosity, n + 1, n + 1, "pair:viscosity"); +} + +/* ---------------------------------------------------------------------- + global settings + ------------------------------------------------------------------------- */ + +void PairSPHTaitwaterMorris::settings(int narg, char **arg) { + if (narg != 0) + error->all( + "Illegal number of setting arguments for pair_style sph/taitwater/morris"); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs + ------------------------------------------------------------------------- */ + +void PairSPHTaitwaterMorris::coeff(int narg, char **arg) { + if (narg != 6) + error->all( + "Incorrect args for pair_style sph/taitwater/morris coefficients"); + if (!allocated) + allocate(); + + int ilo, ihi, jlo, jhi; + force->bounds(arg[0], atom->ntypes, ilo, ihi); + force->bounds(arg[1], atom->ntypes, jlo, jhi); + + double rho0_one = force->numeric(arg[2]); + double soundspeed_one = force->numeric(arg[3]); + double viscosity_one = force->numeric(arg[4]); + double cut_one = force->numeric(arg[5]); + double B_one = soundspeed_one * soundspeed_one * rho0_one / 7.0; + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + rho0[i] = rho0_one; + soundspeed[i] = soundspeed_one; + B[i] = B_one; + for (int j = MAX(jlo,i); j <= jhi; j++) { + viscosity[i][j] = viscosity_one; + //printf("setting cut[%d][%d] = %f\n", i, j, cut_one); + cut[i][j] = cut_one; + + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) + error->all("Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i + ------------------------------------------------------------------------- */ + +double PairSPHTaitwaterMorris::init_one(int i, int j) { + + if (setflag[i][j] == 0) { + error->all("Not all pair sph/taitwater/morris coeffs are not set"); + } + + cut[j][i] = cut[i][j]; + viscosity[j][i] = viscosity[i][j]; + + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- */ + +double PairSPHTaitwaterMorris::single(int i, int j, int itype, int jtype, + double rsq, double factor_coul, double factor_lj, double &fforce) { + fforce = 0.0; + + return 0.0; +} diff --git a/src/USER-SPH/pair_sph_taitwater_morris.h b/src/USER-SPH/pair_sph_taitwater_morris.h new file mode 100644 index 0000000000..5390d73f08 --- /dev/null +++ b/src/USER-SPH/pair_sph_taitwater_morris.h @@ -0,0 +1,48 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(sph/taitwater/morris,PairSPHTaitwaterMorris) + +#else + +#ifndef LMP_PAIR_TAITWATER_MORRIS_H +#define LMP_PAIR_TAITWATER_MORRIS_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairSPHTaitwaterMorris : public Pair { + public: + PairSPHTaitwaterMorris(class LAMMPS *); + virtual ~PairSPHTaitwaterMorris(); + virtual void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + virtual double init_one(int, int); + virtual double single(int, int, int, int, double, double, double, double &); + + protected: + double *rho0, *soundspeed, *B; + double **cut,**viscosity; + int first; + + void allocate(); +}; + +} + +#endif +#endif