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

Fixing order of correct_v

Browse Source
This commit is contained in:
jtclemm
2023-06-23 14:12:59 -06:00
parent 6f59b7c5e0
commit 4ae41edee7
8 changed files with 177 additions and 12 deletions
Download Patch File Download Diff File Expand all files Collapse all files

87
examples/PACKAGES/rheo/dam-break/in.rheo.dam.break Normal file
View File

@ -0,0 +1,87 @@
# ------ 2D dam break ------ #
dimension 2
units lj
atom_style rheo
boundary f s p
comm_modify vel yes
newton off
comm_style tiled
# ------ Create simulation box ------ #
variable sf equal 0.05
variable n equal 1.0/(${sf}^2)
variable cut equal 3.0*${sf}
region box block -1 55 -1 20 -0.01 0.01 units box
create_box 2 box
lattice sq ${n}
region left_wall block -1 0 EDGE EDGE -0.01 0.01 units box
region right_wall block 53.66 EDGE EDGE EDGE -0.01 0.01 units box
region bot_wall block 0.01 53.65 -1 0 -0.01 0.01 units box
region interior block 0.01 10 0.01 20 -0.01 0.01 units box
create_atoms 1 region interior
create_atoms 2 region left_wall
create_atoms 2 region right_wall
create_atoms 2 region bot_wall
group fluid type 1
group static_wall type 2
group dyn_wall type 3
group rig union static_wall dyn_wall
# ------ Model parameters ------ #
variable rho0 equal 1.0
variable mp equal ${rho0}/${n}
variable cs equal 10
variable zeta equal 1
variable Dr equal 0.1*${cut}*${cs}
variable dt_max equal 0.1*${cut}/${cs}/3
variable g equal 0.0245
variable fext equal ${g}/${n}
variable eta equal 0.05
mass * ${mp}
variable d0 atom (${rho0}*${g}*(20-y)/${cs}/${cs})+${rho0}
set group all rheo/rho ${rho0}
set group fluid rheo/rho v_d0
set group all rheo/status 0
set group rig rheo/status 2
timestep ${dt_max}
pair_style rheo ${cut} artificial/visc ${zeta} rho/damp ${Dr}
pair_coeff * *
# ------ Fixes & computes ------ #
fix 1 all rheo ${cut} RK1 32 shift surface/detection coordination 26
fix 2 all rheo/viscosity constant ${eta}
fix 3 all rheo/pressure linear
fix 4 rig setforce 0.0 0.0 0.0
fix 5 fluid addforce 0.0 -${fext} 0.0
fix 6 all balance 1000 1.1 rcb
compute 1 all rheo/property/atom rho phase surface
# ------ Output & Run ------ #
thermo 10
thermo_style custom step time ke press
variable skin equal 0.2*${cut}
neighbor ${skin} bin
neigh_modify one 5000
dump 1 all custom 200 atomDump id type x y vx vy fx fy c_1[*]
run 400000

16
examples/PACKAGES/rheo/poiseuille/in.rheo.poiseuille
View File

@ -1,3 +1,5 @@
# ------ 2D Poiseuille flow ------ #
dimension 2
units lj
atom_style rheo
@ -5,16 +7,12 @@ boundary p p p
comm_modify vel yes
# ------ Particle Lattice/Resolution Parameters ------ #
# ------ Create simulation box ------ #
variable L equal 10
variable sf equal 0.2
variable n equal 1.0/(${sf}^2)
variable cut equal 3.5*${sf}
# ------ Create simulation box ------ #
region box block 0 20 -10 10 -0.01 0.01 units box
create_box 2 box
lattice sq ${n}
@ -46,11 +44,11 @@ variable eta equal 0.1
variable dt_max equal 0.1*${cut}/${cs}/3
variable Dr equal 0.05*${cut}*${cs}
mass 1 ${mp}
mass 2 ${mp}
mass * ${mp}
set group all rheo/rho ${rho0}
set group all rheo/status 0
set group rig rheo/status 2
timestep ${dt_max}
pair_style rheo ${cut} artificial/visc ${zeta} rho/damp ${Dr}
@ -71,13 +69,13 @@ compute 1 all rheo/property/atom rho phase
# ------ Output & Run ------ #
thermo 200
thermo_style custom step time temp press
thermo_style custom step time ke press
variable skin equal 0.2*${cut}
neighbor ${skin} bin
neigh_modify one 5000
dump 1 all custom 200 atomDump id type x y vx vy fx fy c_1[*]
#dump 1 all custom 200 atomDump id type x y vx vy fx fy c_1[*]
run 20000

74
examples/PACKAGES/rheo/taylor-green/in.rheo.taylor.green Normal file
View File

@ -0,0 +1,74 @@
# ------ 2D Taylor Green vortex ------ #
dimension 2
units lj
atom_style rheo
boundary p p p
comm_modify vel yes
newton off
# ------ Create simulation box ------ #
variable sf equal 0.1
variable n equal 1.0/(${sf}^2)
variable cut equal 3.5*${sf}
region box block 0 10 0 10 -0.01 0.01 units box
create_box 1 box
lattice sq ${n}
create_atoms 1 region box
variable dr equal 0.1*${sf}
displace_atoms all random ${dr} ${dr} 0 135414 units box
# ------ Model parameters ------ #
variable rho0 equal 1.0
variable mp equal ${rho0}/${n}
variable cs equal 1.0
variable eta equal 0.05
variable zeta equal 1
variable dt_max equal 0.1*${cut}/${cs}/3
variable Dr equal 0.1*${cut}*${cs}
mass * ${mp}
set group all rheo/rho ${rho0}
set group all rheo/status 0
variable u0 equal 0.05
variable uy atom ${u0}*sin(2*PI*x/lx)*cos(2*PI*y/ly)
variable ux atom -${u0}*sin(2*PI*y/ly)*cos(2*PI*x/ly)
variable d0 atom ${rho0}-${u0}*${u0}*${rho0}*0.25*(cos(4*PI*x/lx)+cos(4*PI*y/ly))/${cs}/${cs}
velocity all set v_ux v_uy 0.0 units box
timestep ${dt_max}
pair_style rheo ${cut} artificial/visc ${zeta} rho/damp ${Dr}
pair_coeff * *
# ------ Fixes & computes ------ #
fix 1 all rheo ${cut} quintic 0 shift
fix 2 all rheo/viscosity constant ${eta}
fix 3 all rheo/pressure linear
compute 1 all rheo/property/atom rho phase
# ------ Output & Run ------ #
thermo 200
thermo_style custom step time ke press
variable skin equal 0.2*${cut}
neighbor ${skin} bin
neigh_modify one 10000 #increase number of allowed neighbors
#dump 1 all custom 200 atomDump id type x y vx vy fx fy c_1[*]
run 10000

2
src/RHEO/compute_rheo_grad.cpp
View File

@ -212,9 +212,11 @@ void ComputeRHEOGrad::compute_peratom()
if (interface_flag) {
if (fluidi && (!fluidj)) {
compute_interface->correct_v(vi, vj, i, j);
//compute_interface->correct_v(vj, vi, j, i);
rhoj = compute_interface->correct_rho(j, i);
} else if ((!fluidi) && fluidj) {
compute_interface->correct_v(vj, vi, j, i);
//compute_interface->correct_v(vi, vj, i, j);
rhoi = compute_interface->correct_rho(i, j);
} else if ((!fluidi) && (!fluidj)) {
rhoi = rho0;

2
src/RHEO/compute_rheo_surface.cpp
View File

@ -126,7 +126,7 @@ void ComputeRHEOSurface::compute_peratom()
{
int i, j, ii, jj, inum, jnum, a, b, itype, jtype, fluidi, fluidj;
double xtmp, ytmp, ztmp, rsq, Voli, Volj, rhoi, rhoj, wp;
double *dWij, *dWji, dx[3];
double dWij[3], dWji[3], dx[3];
int *ilist, *jlist, *numneigh, **firstneigh;
int nlocal = atom->nlocal;

2
src/RHEO/compute_rheo_vshift.cpp
View File

@ -171,9 +171,11 @@ void ComputeRHEOVShift::compute_peratom()
if (interface_flag) {
if (fluidi && (!fluidj)) {
compute_interface->correct_v(vi, vj, i, j);
//compute_interface->correct_v(vj, vi, j, i);
rhoj = compute_interface->correct_rho(j,i);
} else if ((!fluidi) && fluidj) {
compute_interface->correct_v(vj, vi, j, i);
//compute_interface->correct_v(vi, vj, i, j);
rhoi = compute_interface->correct_rho(i,j);
} else if ((!fluidi) && (!fluidj)) {
rhoi = 1.0;

2
src/RHEO/fix_rheo.cpp
View File

@ -154,7 +154,7 @@ void FixRHEO::post_constructor()
if (rhosum_flag) {
compute_rhosum = dynamic_cast<ComputeRHEORhoSum *>(modify->add_compute(
"rheo_rho_sum all RHEO/RHO/SUM"));
"rheo_rhosum all RHEO/RHO/SUM"));
compute_rhosum->fix_rheo = this;
}

2
src/RHEO/pair_rheo.cpp
View File

@ -207,6 +207,7 @@ void PairRHEO::compute(int eflag, int vflag)
if (interface_flag) {
if (fluidi && (!fluidj)) {
compute_interface->correct_v(vi, vj, i, j);
//compute_interface->correct_v(vj, vi, j, i);
rhoj = compute_interface->correct_rho(j, i);
Pj = fix_pressure->calc_pressure(rhoj);
@ -215,6 +216,7 @@ void PairRHEO::compute(int eflag, int vflag)
} else if ((!fluidi) && fluidj) {
compute_interface->correct_v(vj, vi, j, i);
//compute_interface->correct_v(vi, vj, i, j);
rhoi = compute_interface->correct_rho(i, j);
Pi = fix_pressure->calc_pressure(rhoi);