ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

change to clang-format for LATBOLTZ routines

Browse Source
This commit is contained in:
Colin Denniston
2022-03-09 20:41:28 -05:00
parent 888db4b0ef
commit 9b26726f51
5 changed files with 3251 additions and 3058 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5611
src/LATBOLTZ/fix_lb_fluid.cpp
View File

File diff suppressed because it is too large Load Diff

482
src/LATBOLTZ/fix_lb_fluid.h
View File

@ -28,36 +28,24 @@ FixStyle(lb/fluid,FixLbFluid)
namespace LAMMPS_NS {
static const double kappa_lb=0.0;
static const double kappa_lb = 0.0;
// 15-velocity lattice propogation vectors
static const int e15[15][3] = {
{ 0, 0, 0},
{ 1, 0, 0},
{ 0, 1, 0},
{-1, 0, 0},
{ 0,-1, 0},
{ 0, 0, 1},
{ 0, 0,-1},
{ 1, 1, 1},
{-1, 1, 1},
{-1,-1, 1},
{ 1,-1, 1},
{ 1, 1,-1},
{-1, 1,-1},
{-1,-1,-1},
{ 1,-1,-1}
};
// 15-velocity lattice propogation vectors
static const int e15[15][3] = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {-1, 0, 0}, {0, -1, 0},
{0, 0, 1}, {0, 0, -1}, {1, 1, 1}, {-1, 1, 1}, {-1, -1, 1},
{1, -1, 1}, {1, 1, -1}, {-1, 1, -1}, {-1, -1, -1}, {1, -1, -1}};
// 15-velocity weights
static const double w_lb15[15] =
{ 2./9., 1./9., 1./9., 1./9., 1./9., 1./9., 1./9., 1./72., 1./72., 1./72., 1./72., 1./72., 1./72., 1./72., 1./72.};
// 15-velocity weights
static const double w_lb15[15] = {2. / 9., 1. / 9., 1. / 9., 1. / 9., 1. / 9.,
1. / 9., 1. / 9., 1. / 72., 1. / 72., 1. / 72.,
1. / 72., 1. / 72., 1. / 72., 1. / 72., 1. / 72.};
// 15-velocity normalizations
static const double Ng_lb15[15] =
{ 1., 3., 3., 3., 9./2., 9./2., 9./2., 9., 9., 9., 27./2., 27./2., 27./2., 9., 1.};
// 15-velocity normalizations
static const double Ng_lb15[15] = {1., 3., 3., 3., 9. / 2., 9. / 2., 9. / 2., 9.,
9., 9., 27. / 2., 27. / 2., 27. / 2., 9., 1.};
// 15-velcity transformation matrix for f_i to moments
// 15-velcity transformation matrix for f_i to moments
// clang-format off
static const double mg_lb15[15][15] = {
{ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.},
{ 0., 1., 0., -1., 0., 0., 0., 1., -1., -1., 1., 1., -1., -1., 1.},
@ -75,54 +63,55 @@ static const double mg_lb15[15][15] = {
{ 0., 0., 0., 0., 0., 0., 0., 1., -1., 1., -1., -1., 1., -1., 1.},
{M_SQRT2,-M_SQRT1_2,-M_SQRT1_2,-M_SQRT1_2,-M_SQRT1_2,-M_SQRT1_2,-M_SQRT1_2,M_SQRT2,M_SQRT2,M_SQRT2,M_SQRT2,M_SQRT2,M_SQRT2,M_SQRT2,M_SQRT2}
};
// clang-format on
// 15-velocity opposite lattice directions for bounce-back, i.e. od[i] = j such that e15[j]=-e15[i]
static const int od[15] =
{ 0, 3, 4, 1, 2, 6, 5, 13, 14, 11, 12, 9, 10, 7, 8};
// 15-velocity opposite lattice directions for bounce-back, i.e. od[i] = j such that e15[j]=-e15[i]
static const int od[15] = {0, 3, 4, 1, 2, 6, 5, 13, 14, 11, 12, 9, 10, 7, 8};
// 15-velocity bounce-back list
// bbl[i][0] = number of bounce-back directions for orientation i
// bbl[i][j]...bbl[i][bbl[i][0]] directions that would be coming from inside the wall so need to come from bounce-back
// bbl[i][[bbl[i][0]+1]...bbl[i][16] directions where standard propogation can proceed (pointing into or along wall)
// inside edge has 1/4 inside domain, 3/4 outside domain
// outside edge has 3/4 outside domain, 1/4 inside domain
// Note: 1. This list is not exhaustive (eg. there should be 12 inside and 12 outside edges possible, it just covers cases
// accessible in the pit routines. Could be generalized to include other geometries
// 2. Need better labelling for corners (particularly in-out) that distinguishes different cases (e.g. 10 and 29 are NOT same, also 11,31)
// ori wall normals (point into domain)
// 0 not relevent, ori==0 only for lattice type 0 (standard bulk fluid) and 2 (outside domain)
// 1 wall +x
// 2 wall +y
// 3 wall +z
// 4 outside edge +x,+z
// 5 inside edge +x,+z
// 6 inside edge +y,+z
// 7 outside edge +y,+z
// 8 inside edge -x,-y
// 9 inside edge -x,+y
// 10 in-out corner +x,+y,+z
// 11 in-out corner +x,-y,+z
// 12 inside corner -x,+y,+z
// 13 inside corner -x,-y,+z
// 14 wall -x
// 15 wall -y
// 16 wall -z
// 17 outside edge -x,+z
// 18 inside edge -x,+z
// 19 inside edge -y,+z
// 20 outside edge -y,+z
// 21 inside edge +x,-y
// 22 inside edge +x,+y
// 23 in-out corner -x,+y,+z
// 24 in-out corner -x,-y,+z
// 25 inside corner +x,+y,+z
// 26 inside corner +x,-y,+z
// 27 inside edge +y,-z
// 28 inside edge -y,-z
// 29 in-out corner +x,+y,+z
// 30 in-out corner +x,-y,+z
// 31 in-out corner -x,+y,+z
// 32 in-out corner -x,-y,+z
// 15-velocity bounce-back list
// bbl[i][0] = number of bounce-back directions for orientation i
// bbl[i][j]...bbl[i][bbl[i][0]] directions that would be coming from inside the wall so need to come from bounce-back
// bbl[i][[bbl[i][0]+1]...bbl[i][16] directions where standard propogation can proceed (pointing into or along wall)
// inside edge has 1/4 inside domain, 3/4 outside domain
// outside edge has 3/4 outside domain, 1/4 inside domain
// Note: 1. This list is not exhaustive (eg. there should be 12 inside and 12 outside edges possible, it just covers cases
// accessible in the pit routines. Could be generalized to include other geometries
// 2. Need better labelling for corners (particularly in-out) that distinguishes different cases (e.g. 10 and 29 are NOT same, also 11,31)
// ori wall normals (point into domain)
// 0 not relevent, ori==0 only for lattice type 0 (standard bulk fluid) and 2 (outside domain)
// 1 wall +x
// 2 wall +y
// 3 wall +z
// 4 outside edge +x,+z
// 5 inside edge +x,+z
// 6 inside edge +y,+z
// 7 outside edge +y,+z
// 8 inside edge -x,-y
// 9 inside edge -x,+y
// 10 in-out corner +x,+y,+z
// 11 in-out corner +x,-y,+z
// 12 inside corner -x,+y,+z
// 13 inside corner -x,-y,+z
// 14 wall -x
// 15 wall -y
// 16 wall -z
// 17 outside edge -x,+z
// 18 inside edge -x,+z
// 19 inside edge -y,+z
// 20 outside edge -y,+z
// 21 inside edge +x,-y
// 22 inside edge +x,+y
// 23 in-out corner -x,+y,+z
// 24 in-out corner -x,-y,+z
// 25 inside corner +x,+y,+z
// 26 inside corner +x,-y,+z
// 27 inside edge +y,-z
// 28 inside edge -y,-z
// 29 in-out corner +x,+y,+z
// 30 in-out corner +x,-y,+z
// 31 in-out corner -x,+y,+z
// 32 in-out corner -x,-y,+z
// clang-format off
static const int bbl[33][16] = {
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{ 5, 1, 7, 10, 11, 14, 0, 2, 3, 4, 5, 6, 8, 9, 12, 13},
@ -158,38 +147,24 @@ static const int bbl[33][16] = {
{ 6, 2, 7, 8, 12, 9, 11, 0, 1, 3, 4, 5, 6, 10, 13, 14},
{ 6, 4, 9, 10, 13, 8, 14, 0, 1, 2, 3, 5, 6, 7, 11, 12}
};
//clang-format on
// 19-velocity lattice propogation vectors
static const int e19[19][3] = {
{ 0, 0, 0},
{ 1, 0, 0},
{ 0, 1, 0},
{-1, 0, 0},
{ 0,-1, 0},
{ 0, 0, 1},
{ 0, 0,-1},
{ 1, 1, 0},
{1, -1, 0},
{-1, 1, 0},
{-1,-1, 0},
{ 1, 0, 1},
{ 1, 0,-1},
{-1, 0, 1},
{-1, 0,-1},
{ 0, 1, 1},
{ 0, 1,-1},
{ 0,-1, 1},
{ 0,-1,-1}
};
// 19-velocity lattice propogation vectors
static const int e19[19][3] = {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {-1, 0, 0}, {0, -1, 0},
{0, 0, 1}, {0, 0, -1}, {1, 1, 0}, {1, -1, 0}, {-1, 1, 0},
{-1, -1, 0}, {1, 0, 1}, {1, 0, -1}, {-1, 0, 1}, {-1, 0, -1},
{0, 1, 1}, {0, 1, -1}, {0, -1, 1}, {0, -1, -1}};
static const double w_lb19[19] =
{1./3.,1./18.,1./18.,1./18.,1./18.,1./18.,1./18.,1./36.,1./36.,1./36.,1./36.,1./36.,1./36.,1./36.,1./36.,1./36., 1./36., 1./36.,1./36.};
static const double w_lb19[19] = {1. / 3., 1. / 18., 1. / 18., 1. / 18., 1. / 18.,
1. / 18., 1. / 18., 1. / 36., 1. / 36., 1. / 36.,
1. / 36., 1. / 36., 1. / 36., 1. / 36., 1. / 36.,
1. / 36., 1. / 36., 1. / 36., 1. / 36.};
static const double Ng_lb19[19] =
{ 1., 3., 3., 3., 9./2., 9./2., 9./2., 9., 9., 9.,27./2.,27./2.,27./2., 18., 18., 18.,162./7.,126./5., 30.};
static const double Ng_lb19[19] = {1., 3., 3., 3., 9. / 2., 9. / 2., 9. / 2.,
9., 9., 9., 27. / 2., 27. / 2., 27. / 2., 18.,
18., 18., 162. / 7., 126. / 5., 30.};
// clang-format off
static const double mg_lb19[19][19] = {
{ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.},
{ 0., 1., 0., -1., 0., 0., 0., 1., 1., -1., -1., 1., 1., -1., -1., 0., 0., 0., 0.},
@ -211,194 +186,191 @@ static const double mg_lb19[19][19] = {
{1./14.,-5./14., 1./7.,-5./14., 1./7.,-3./14.,-3./14., 0., 0., 0., 0.,5./14.,5./14.,5./14.,5./14.,-1./7.,-1./7.,-1./7.,-1./7.},
{1./10., 0.,-3./10., 0.,-3./10.,-3./10.,-3./10., 0., 0., 0., 0., 0., 0., 0., 0.,3./10.,3./10.,3./10.,3./10.}
};
// clang-format on
class Site {
public:
int type;
int orientation;
};
public:
int type;
int orientation;
};
class FixLbFluid : public Fix {
friend class FixLbMomentum;
friend class FixLbViscous;
friend class FixLbMomentum;
friend class FixLbViscous;
public:
FixLbFluid(class LAMMPS *, int, char **);
~FixLbFluid() override;
int setmask() override;
void init() override;
void initial_integrate(int) override;
void setup(int) override;
void pre_force(int) override;
void post_force(int) override;
void final_integrate() override;
void end_of_step() override;
public:
FixLbFluid(class LAMMPS *, int, char **);
~FixLbFluid() override;
int setmask() override;
void init() override;
void initial_integrate(int) override;
void setup(int) override;
void pre_force(int) override;
void post_force(int) override;
void final_integrate() override;
void end_of_step() override;
void grow_arrays(int) override;
void copy_arrays(int, int, int) override;
int pack_exchange(int, double *) override;
int unpack_exchange(int, double *) override;
void grow_arrays(int) override;
void copy_arrays(int, int, int) override;
int pack_exchange(int, double *) override;
int unpack_exchange(int, double *) override;
double compute_scalar() override;
double compute_vector(int) override;
double compute_scalar() override;
double compute_vector(int) override;
void dump(int);
private:
double viscosity,densityinit_real,a_0_real,T;
int setdx,seta0,setdof;
int numvel;
void dump(int);
double dm_lb,dx_lb,dt_lb; // Lattice units for mass, distance, time.
private:
double viscosity, densityinit_real, a_0_real, T;
int setdx, seta0, setdof;
int numvel;
int Nbx,Nby,Nbz; // Total # of x,y,z grid points.
int subNbx,subNby,subNbz; // # of x,y,z, grid points (including buffer)
// on local processor.
int fluid_local_n0[3]; // Size of local including both lower and upper end points
int fluid_global_o0[3]; // Offset of local in global from lower end point
int fluid_global_n0[3]; // Size of global including both lower and upper end points
int me, nprocs; // MPI variables: processor ID, # of processors
MPI_Datatype oneslice; // MPI datatypes to pass arrays.
MPI_Datatype passxu,passyu,passzu;
MPI_Datatype passxf,passyf,passzf;
MPI_Datatype passxrho,passyrho,passzrho;
MPI_Datatype passxtemp,passytemp,passztemp;
MPI_Datatype passxWtemp,passyWtemp,passzWtemp;
MPI_Datatype fluid_density_2_mpitype;
MPI_Datatype fluid_velocity_2_mpitype;
double dm_lb, dx_lb, dt_lb; // Lattice units for mass, distance, time.
MPI_Datatype realType3_mpitype; // MPI type for a 3-vector
int Nbx, Nby, Nbz; // Total # of x,y,z grid points.
int subNbx, subNby, subNbz; // # of x,y,z, grid points (including buffer)
// on local processor.
int fluid_local_n0[3]; // Size of local including both lower and upper end points
int fluid_global_o0[3]; // Offset of local in global from lower end point
int fluid_global_n0[3]; // Size of global including both lower and upper end points
double kB,densityinit,a_0; // Boltzmann constant, initial density,
// and a_0 all in lattice units.
double *Gamma;
double **hydroF;
double *massp;
int me, nprocs; // MPI variables: processor ID, # of processors
MPI_Datatype oneslice; // MPI datatypes to pass arrays.
MPI_Datatype passxu, passyu, passzu;
MPI_Datatype passxf, passyf, passzf;
MPI_Datatype passxrho, passyrho, passzrho;
MPI_Datatype passxtemp, passytemp, passztemp;
MPI_Datatype passxWtemp, passyWtemp, passzWtemp;
MPI_Datatype fluid_density_2_mpitype;
MPI_Datatype fluid_velocity_2_mpitype;
MPI_Datatype realType3_mpitype; // MPI type for a 3-vector
double kB, densityinit, a_0; // Boltzmann constant, initial density,
// and a_0 all in lattice units.
double *Gamma;
double **hydroF;
double *massp;
int dump_interval, dump_time_index;
std::string dump_file_name_xdmf;
std::string dump_file_name_raw;
FILE* dump_file_handle_xdmf;
MPI_File dump_file_handle_raw;
MPI_Datatype dump_file_mpitype;
int groupbit_viscouslb;
std::string dump_file_name_xdmf;
std::string dump_file_name_raw;
FILE *dump_file_handle_xdmf;
MPI_File dump_file_handle_raw;
MPI_Datatype dump_file_mpitype;
double ***density_lb; // fluid density
double ****u_lb; // fluid velocity
double ****f_lb; // distributions
double ****fnew; // used in the calculation of the new
// distributions.
double ****feq; // equilibrium distributions
int groupbit_viscouslb;
double ****Ff,***Wf; // Force, total weight, from the MD particles on the fluid.
double ****Fftempx,***Wftempx;
double ****Fftempy,***Wftempy;
double ****Fftempz,***Wftempz;
double ***density_lb; // fluid density
double ****u_lb; // fluid velocity
double ****f_lb; // distributions
double ****fnew; // used in the calculation of the new
// distributions.
double ****feq; // equilibrium distributions
double tau; // Lattice Boltzmann variables.
double K_0;
double ****Ff, ***Wf; // Force, total weight, from the MD particles on the fluid.
double ****Fftempx, ***Wftempx;
double ****Fftempy, ***Wftempy;
double ****Fftempz, ***Wftempz;
int step;
double tau; // Lattice Boltzmann variables.
double K_0;
int n_stencil=2; // Number of points for spread/interpolate stencil
int step;
double bodyforcex,bodyforcey,bodyforcez; // Body Forces acting on the fluid (default=0)
double vwtp,vwbt; // Velocities of the z walls in the y
// direction. (must have fixed boundary
// conditions in z)
int n_stencil = 2; // Number of points for spread/interpolate stencil
int pressure; // pressure boundary conditions on/off
double rhofactor; // factor for density/pressure jump at boundary
double rhoH,rhoL; // density target for high/low side of jump
double meanrho1,meanrhoLx; // current densities at boundary on high/low side of jump
int noisestress; // 1 to include noise in the system,
// 0 otherwise.
int lin_init; // 1 to initialize with linear interpolation
// between boundaries.
// 0 initialize to uniform density, 0.0 velocities.
double namp,noisefactor;
int seed;
class RanMars *random;
double bodyforcex, bodyforcey, bodyforcez; // Body Forces acting on the fluid (default=0)
double vwtp, vwbt; // Velocities of the z walls in the y
// direction. (must have fixed boundary
// conditions in z)
int readrestart; // 1 to read in data from a restart file.
MPI_File pFileRead;
int pressure; // pressure boundary conditions on/off
double rhofactor; // factor for density/pressure jump at boundary
double rhoH, rhoL; // density target for high/low side of jump
double meanrho1, meanrhoLx; // current densities at boundary on high/low side of jump
int printrestart; // 1 to write data to a restart file.
MPI_File pFileWrite;
int noisestress; // 1 to include noise in the system,
// 0 otherwise.
int lin_init; // 1 to initialize with linear interpolation
// between boundaries.
// 0 initialize to uniform density, 0.0 velocities.
double timeEqb,timeUpdate,timePCalc,timefluidForce,timeCorrectU;
double namp, noisefactor;
int seed;
class RanMars *random;
double dof_lb=0;
int readrestart; // 1 to read in data from a restart file.
MPI_File pFileRead;
int fixviscouslb;
int printrestart; // 1 to write data to a restart file.
MPI_File pFileWrite;
void rescale(void);
void SetupBuffers(void);
void InitializeFirstRun(void);
void initializeLB(void);
void initialize_feq(void);
void (FixLbFluid::*equilibriumdist)(int,int,int,int,int,int);
void equilibriumdist15(int,int,int,int,int,int);
void equilibriumdist19(int,int,int,int,int,int);
void parametercalc_part(int,int,int,int,int,int);
void correctu_part(int,int,int,int,int,int);
void parametercalc_full(void);
void update_periodic(int,int,int,int,int,int);
void correctu_full(void);
double timeEqb, timeUpdate, timePCalc, timefluidForce, timeCorrectU;
void calc_mass_momentum(double& totalmass, double totalmomentum[3]);
void calc_MPT(double& totalmass, double totalmomentum[3], double& Tave);
void (FixLbFluid::*update_full)(void);
void update_full15(void);
void update_full19(void);
double dof_lb = 0;
void read_restartfile(void);
void write_restartfile(void);
int fixviscouslb;
void (FixLbFluid::*interpolate)(int, int);
void (FixLbFluid::*interpolationweight)(int);
void keys_interpolation(int, int);
void keys_interpolationweight(int);
void trilinear_interpolation(int, int);
void trilinear_interpolationweight(int);
void IBM3_interpolation(int, int);
void IBM3_interpolationweight(int);
void calc_fluidforceI(void);
void calc_fluidforceII(void);
void calc_fluidforceweight(void);
void rescale(void);
void SetupBuffers(void);
void InitializeFirstRun(void);
void initializeLB(void);
void initialize_feq(void);
void (FixLbFluid::*equilibriumdist)(int, int, int, int, int, int);
void equilibriumdist15(int, int, int, int, int, int);
void equilibriumdist19(int, int, int, int, int, int);
void parametercalc_part(int, int, int, int, int, int);
void correctu_part(int, int, int, int, int, int);
void parametercalc_full(void);
void update_periodic(int, int, int, int, int, int);
void correctu_full(void);
int adjust_dof_fix();
double dof_compute();
void calc_mass_momentum(double &totalmass, double totalmomentum[3]);
void calc_MPT(double &totalmass, double totalmomentum[3], double &Tave);
/* nanopit parameters */
int npits; // number of nanopits
int h_s; // slit height
int h_p; // pit height
int w_p; // pit width
int l_p; // pit length
int l_pp; // distance between consecutive pits
int l_e; // length of end segments
int sw; // side walls on/off
int openingsites; // Number of active fluid sites at x=0
Site ***sublattice, ***wholelattice; // lattice geometry
void (FixLbFluid::*update_full)(void);
void update_full15(void);
void update_full19(void);
/* nanopit routines */
void addslit(int &,int,int,int,int);
void addpit(int &,int,int,int,int,int);
void initializeGeometry(void);
void initializeGlobalGeometry(void);
};
void read_restartfile(void);
void write_restartfile(void);
void (FixLbFluid::*interpolate)(int, int);
void (FixLbFluid::*interpolationweight)(int);
void keys_interpolation(int, int);
void keys_interpolationweight(int);
void trilinear_interpolation(int, int);
void trilinear_interpolationweight(int);
void IBM3_interpolation(int, int);
void IBM3_interpolationweight(int);
void calc_fluidforceI(void);
void calc_fluidforceII(void);
void calc_fluidforceweight(void);
int adjust_dof_fix();
double dof_compute();
/* nanopit parameters */
int npits; // number of nanopits
int h_s; // slit height
int h_p; // pit height
int w_p; // pit width
int l_p; // pit length
int l_pp; // distance between consecutive pits
int l_e; // length of end segments
int sw; // side walls on/off
int openingsites; // Number of active fluid sites at x=0
Site ***sublattice, ***wholelattice; // lattice geometry
/* nanopit routines */
void addslit(int &, int, int, int, int);
void addpit(int &, int, int, int, int, int);
void initializeGeometry(void);
void initializeGlobalGeometry(void);
};
} // namespace LAMMPS_NS
#endif
#endif

143
src/LATBOLTZ/fix_lb_momentum.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -21,25 +20,24 @@
#include "fix_lb_momentum.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "group.h"
#include "error.h"
#include "fix_lb_fluid.h"
#include "group.h"
#include "modify.h"
#include "comm.h"
using namespace LAMMPS_NS;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixLbMomentum::FixLbMomentum(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
FixLbMomentum::FixLbMomentum(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
{
if (narg < 4) error->all(FLERR,"Illegal fix lb/momentum command");
nevery = utils::inumeric(FLERR,arg[3],false,lmp);;
if (nevery <= 0) error->all(FLERR,"Illegal fix lb/momentum command");
if (narg < 4) error->all(FLERR, "Illegal fix lb/momentum command");
nevery = utils::inumeric(FLERR, arg[3], false, lmp);
;
if (nevery <= 0) error->all(FLERR, "Illegal fix lb/momentum command");
linear = 1;
xflag = 1;
@ -48,23 +46,22 @@ FixLbMomentum::FixLbMomentum(LAMMPS *lmp, int narg, char **arg) :
int iarg = 4;
while (iarg < narg) {
if (strcmp(arg[iarg],"linear") == 0) {
if (iarg+4 > narg) error->all(FLERR,"Illegal fix lb/momentum command");
if (strcmp(arg[iarg], "linear") == 0) {
if (iarg + 4 > narg) error->all(FLERR, "Illegal fix lb/momentum command");
linear = 1;
xflag = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
yflag = utils::inumeric(FLERR,arg[iarg+2],false,lmp);
zflag = utils::inumeric(FLERR,arg[iarg+3],false,lmp);
xflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
yflag = utils::inumeric(FLERR, arg[iarg + 2], false, lmp);
zflag = utils::inumeric(FLERR, arg[iarg + 3], false, lmp);
iarg += 4;
} else error->all(FLERR,"Illegal fix lb/momentum command");
} else
error->all(FLERR, "Illegal fix lb/momentum command");
}
if (linear == 0)
error->all(FLERR,"Illegal fix lb/momentum command");
if (linear == 0) error->all(FLERR, "Illegal fix lb/momentum command");
if (linear)
if (xflag < 0 || xflag > 1 || yflag < 0 || yflag > 1 ||
zflag < 0 || zflag > 1)
error->all(FLERR,"Illegal fix lb/momentum command");
if (xflag < 0 || xflag > 1 || yflag < 0 || yflag > 1 || zflag < 0 || zflag > 1)
error->all(FLERR, "Illegal fix lb/momentum command");
}
/* ---------------------------------------------------------------------- */
@ -81,13 +78,13 @@ int FixLbMomentum::setmask()
void FixLbMomentum::init()
{
// warn if 0 atoms in group
if ((count = group->count(igroup)) == 0)
error->warning(FLERR,"Fix lb/momentum group has no atoms: Only fluid momentum affected");
if ((count = group->count(igroup)) == 0)
error->warning(FLERR, "Fix lb/momentum group has no atoms: Only fluid momentum affected");
for (int ifix = 0; ifix < modify->nfix; ifix++)
if (strcmp(modify->fix[ifix]->style, "lb/fluid") == 0)
fix_lb_fluid = (FixLbFluid *) modify->fix[ifix];
for(int ifix=0; ifix<modify->nfix; ifix++)
if(strcmp(modify->fix[ifix]->style,"lb/fluid")==0)
fix_lb_fluid = (FixLbFluid *)modify->fix[ifix];
count ? masstotal = group->mass(igroup) : 0;
}
@ -97,17 +94,17 @@ void FixLbMomentum::end_of_step()
{
// particle com velcity
double vcmp[3] = {0, 0, 0};
if (count) group->vcm(igroup,masstotal,vcmp);
if (count) group->vcm(igroup, masstotal, vcmp);
// total fluid mass and momentum.
double masslb,momentumlb[3];
double masslb, momentumlb[3];
fix_lb_fluid->calc_mass_momentum(masslb, momentumlb);
//Calculate the center of mass velocity of the particles + fluid.
double vcmall[3];
vcmall[0] = (masstotal*vcmp[0] + momentumlb[0])/(masslb + masstotal);
vcmall[1] = (masstotal*vcmp[1] + momentumlb[1])/(masslb + masstotal);
vcmall[2] = (masstotal*vcmp[2] + momentumlb[2])/(masslb + masstotal);
vcmall[0] = (masstotal * vcmp[0] + momentumlb[0]) / (masslb + masstotal);
vcmall[1] = (masstotal * vcmp[1] + momentumlb[1]) / (masslb + masstotal);
vcmall[2] = (masstotal * vcmp[2] + momentumlb[2]) / (masslb + masstotal);
// adjust velocities to zero net linear momentum
// only adjust a component if flag is set
@ -119,9 +116,9 @@ void FixLbMomentum::end_of_step()
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (xflag) v[i][0] -= vcmall[0];
if (yflag) v[i][1] -= vcmall[1];
if (zflag) v[i][2] -= vcmall[2];
if (xflag) v[i][0] -= vcmall[0];
if (yflag) v[i][1] -= vcmall[1];
if (zflag) v[i][2] -= vcmall[2];
}
}
@ -130,9 +127,9 @@ void FixLbMomentum::end_of_step()
// convert vcmall to lattice units used in lb/fluid
double dx_lb = fix_lb_fluid->dx_lb;
double dt_lb = fix_lb_fluid->dt_lb;
vcmall[0] *= dt_lb/dx_lb;
vcmall[1] *= dt_lb/dx_lb;
vcmall[2] *= dt_lb/dx_lb;
vcmall[0] *= dt_lb / dx_lb;
vcmall[1] *= dt_lb / dx_lb;
vcmall[2] *= dt_lb / dx_lb;
double ucmx, ucmy, ucmz;
ucmx = xflag ? vcmall[0] : 0.0;
@ -140,7 +137,7 @@ void FixLbMomentum::end_of_step()
ucmz = zflag ? vcmall[2] : 0.0;
int numvel = fix_lb_fluid->numvel;
double etacov[14]; // This is for changes and only moments up to 13 are affected
double etacov[14]; // This is for changes and only moments up to 13 are affected
double rho;
int subNbx = fix_lb_fluid->subNbx;
int subNby = fix_lb_fluid->subNby;
@ -150,41 +147,41 @@ void FixLbMomentum::end_of_step()
double ****u_lb = fix_lb_fluid->u_lb;
etacov[0] = etacov[10] = etacov[11] = etacov[12] = etacov[13] = 0.0;
for(int i=0; i<subNbx; i++)
for(int j=0; j<subNby; j++)
for(int k=0; k<subNbz; k++){
rho = density_lb[i][j][k];
etacov[1]=rho*ucmx;
etacov[2]=rho*ucmy;
etacov[3]=rho*ucmz;
etacov[4]=rho*(2.*u_lb[i][j][k][0]*ucmx-ucmx*ucmx);
etacov[5]=rho*(2.*u_lb[i][j][k][1]*ucmy-ucmy*ucmy);
etacov[6]=rho*(2.*u_lb[i][j][k][2]*ucmz-ucmz*ucmz);
etacov[7]=rho*(u_lb[i][j][k][0]*ucmy+u_lb[i][j][k][1]*ucmx-ucmx*ucmy);
etacov[8]=rho*(u_lb[i][j][k][1]*ucmz+u_lb[i][j][k][2]*ucmy-ucmy*ucmz);
etacov[9]=rho*(u_lb[i][j][k][0]*ucmz+u_lb[i][j][k][2]*ucmx-ucmx*ucmz);
for (int i = 0; i < subNbx; i++)
for (int j = 0; j < subNby; j++)
for (int k = 0; k < subNbz; k++) {
rho = density_lb[i][j][k];
if (numvel==15) {
etacov[13]=rho*(u_lb[i][j][k][0]*u_lb[i][j][k][1]*ucmz+u_lb[i][j][k][0]*ucmy*u_lb[i][j][k][2]-
u_lb[i][j][k][0]*ucmy*ucmz+ucmx*u_lb[i][j][k][1]*u_lb[i][j][k][2]-
ucmx*u_lb[i][j][k][1]*ucmz-ucmx*ucmy*u_lb[i][j][k][2]+
ucmx*ucmy*ucmz);
for(int l=0; l<15; l++)
for(int ii=1; ii<14; ii++){
f_lb[i][j][k][l] -= w_lb15[l]*mg_lb15[ii][l]*etacov[ii]*Ng_lb15[ii];
}
}
else // 19-velocity model
for(int l=0; l<19; l++)
for(int ii=1; ii<14; ii++){
f_lb[i][j][k][l] -= w_lb19[l]*mg_lb19[ii][l]*etacov[ii]*Ng_lb19[ii];
}
if(xflag) u_lb[i][j][k][0] -= vcmall[0];
if(yflag) u_lb[i][j][k][1] -= vcmall[1];
if(zflag) u_lb[i][j][k][2] -= vcmall[2];
etacov[1] = rho * ucmx;
etacov[2] = rho * ucmy;
etacov[3] = rho * ucmz;
etacov[4] = rho * (2. * u_lb[i][j][k][0] * ucmx - ucmx * ucmx);
etacov[5] = rho * (2. * u_lb[i][j][k][1] * ucmy - ucmy * ucmy);
etacov[6] = rho * (2. * u_lb[i][j][k][2] * ucmz - ucmz * ucmz);
etacov[7] = rho * (u_lb[i][j][k][0] * ucmy + u_lb[i][j][k][1] * ucmx - ucmx * ucmy);
etacov[8] = rho * (u_lb[i][j][k][1] * ucmz + u_lb[i][j][k][2] * ucmy - ucmy * ucmz);
etacov[9] = rho * (u_lb[i][j][k][0] * ucmz + u_lb[i][j][k][2] * ucmx - ucmx * ucmz);
if (numvel == 15) {
etacov[13] = rho *
(u_lb[i][j][k][0] * u_lb[i][j][k][1] * ucmz +
u_lb[i][j][k][0] * ucmy * u_lb[i][j][k][2] - u_lb[i][j][k][0] * ucmy * ucmz +
ucmx * u_lb[i][j][k][1] * u_lb[i][j][k][2] - ucmx * u_lb[i][j][k][1] * ucmz -
ucmx * ucmy * u_lb[i][j][k][2] + ucmx * ucmy * ucmz);
for (int l = 0; l < 15; l++)
for (int ii = 1; ii < 14; ii++) {
f_lb[i][j][k][l] -= w_lb15[l] * mg_lb15[ii][l] * etacov[ii] * Ng_lb15[ii];
}
} else // 19-velocity model
for (int l = 0; l < 19; l++)
for (int ii = 1; ii < 14; ii++) {
f_lb[i][j][k][l] -= w_lb19[l] * mg_lb19[ii][l] * etacov[ii] * Ng_lb19[ii];
}
if (xflag) u_lb[i][j][k][0] -= vcmall[0];
if (yflag) u_lb[i][j][k][1] -= vcmall[1];
if (zflag) u_lb[i][j][k][2] -= vcmall[2];
}
}

6
src/LATBOLTZ/fix_lb_momentum.h
View File

@ -23,7 +23,7 @@ FixStyle(lb/momentum,FixLbMomentum)
#include "fix.h"
namespace LAMMPS_NS {
class FixLbMomentum : public Fix {
public:
FixLbMomentum(class LAMMPS *, int, char **);
@ -33,14 +33,14 @@ class FixLbMomentum : public Fix {
private:
int linear;
int xflag,yflag,zflag;
int xflag, yflag, zflag;
double masstotal;
int count;
class FixLbFluid *fix_lb_fluid;
};
} // namespace LAMMPS_NS
} // namespace LAMMPS_NS
#endif
#endif

67
src/LATBOLTZ/fix_lb_viscous.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -18,48 +17,46 @@
#include "fix_lb_viscous.h"
#include "atom.h"
#include "update.h"
#include "respa.h"
#include "error.h"
#include "fix_lb_fluid.h"
#include "modify.h"
#include "group.h"
#include "modify.h"
#include "respa.h"
#include "update.h"
using namespace LAMMPS_NS;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixLbViscous::FixLbViscous(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
FixLbViscous::FixLbViscous(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
{
if (narg < 3) error->all(FLERR,"Illegal fix lb/viscous command");
if (narg < 3) error->all(FLERR, "Illegal fix lb/viscous command");
int groupbit_lb_fluid = 0;
for(int ifix=0; ifix<modify->nfix; ifix++)
if(strcmp(modify->fix[ifix]->style,"lb/fluid")==0){
fix_lb_fluid = (FixLbFluid *)modify->fix[ifix];
for (int ifix = 0; ifix < modify->nfix; ifix++)
if (strcmp(modify->fix[ifix]->style, "lb/fluid") == 0) {
fix_lb_fluid = (FixLbFluid *) modify->fix[ifix];
groupbit_lb_fluid = group->bitmask[modify->fix[ifix]->igroup];
}
if(groupbit_lb_fluid == 0)
error->all(FLERR,"the lb/fluid fix must also be used if using the lb/viscous fix");
if (groupbit_lb_fluid == 0)
error->all(FLERR, "the lb/fluid fix must also be used if using the lb/viscous fix");
int *mask = atom->mask;
int nlocal = atom->nlocal;
for(int j=0; j<nlocal; j++){
if((mask[j] & groupbit) && !(mask[j] & groupbit_lb_fluid))
error->one(FLERR,"to apply a fluid force onto an atom, the lb/fluid fix must be used for that atom.");
for (int j = 0; j < nlocal; j++) {
if ((mask[j] & groupbit) && !(mask[j] & groupbit_lb_fluid))
error->one(
FLERR,
"to apply a fluid force onto an atom, the lb/fluid fix must be used for that atom.");
}
}
/* ---------------------------------------------------------------------- */
FixLbViscous::~FixLbViscous()
{
}
FixLbViscous::~FixLbViscous() {}
/* ---------------------------------------------------------------------- */
@ -77,21 +74,20 @@ int FixLbViscous::setmask()
void FixLbViscous::init()
{
if (strcmp(update->integrate_style,"respa") == 0)
if (strcmp(update->integrate_style, "respa") == 0)
nlevels_respa = ((Respa *) update->integrate)->nlevels;
}
/* ---------------------------------------------------------------------- */
void FixLbViscous::setup(int vflag)
{
if (utils::strmatch(update->integrate_style,"^verlet"))
if (utils::strmatch(update->integrate_style, "^verlet"))
post_force(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa - 1);
post_force_respa(vflag, nlevels_respa - 1, 0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa - 1);
}
}
@ -124,31 +120,30 @@ void FixLbViscous::post_force(int /*vflag*/)
if (rmass) {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
massfactor = rmass[i]/(rmass[i]+massp[i]);
massfactor = rmass[i] / (rmass[i] + massp[i]);
f[i][0] = hydroF[i][0] + massfactor*f[i][0];
f[i][1] = hydroF[i][1] + massfactor*f[i][1];
f[i][2] = hydroF[i][2] + massfactor*f[i][2];
f[i][0] = hydroF[i][0] + massfactor * f[i][0];
f[i][1] = hydroF[i][1] + massfactor * f[i][1];
f[i][2] = hydroF[i][2] + massfactor * f[i][2];
}
} else {
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
massfactor = mass[type[i]]/(mass[type[i]]+massp[i]);
massfactor = mass[type[i]] / (mass[type[i]] + massp[i]);
f[i][0] = hydroF[i][0] + massfactor*f[i][0];
f[i][1] = hydroF[i][1] + massfactor*f[i][1];
f[i][2] = hydroF[i][2] + massfactor*f[i][2];
f[i][0] = hydroF[i][0] + massfactor * f[i][0];
f[i][1] = hydroF[i][1] + massfactor * f[i][1];
f[i][2] = hydroF[i][2] + massfactor * f[i][2];
}
}
}
/* ---------------------------------------------------------------------- */
void FixLbViscous::post_force_respa(int vflag, int ilevel, int /*iloop*/)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
if (ilevel == nlevels_respa - 1) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
@ -157,5 +152,3 @@ void FixLbViscous::min_post_force(int vflag)
{
post_force(vflag);
}