Skip to content
This repository

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse code

Add more configuration files

  • Loading branch information...
commit 88acc932271ebb438226876808e8b86e6dca61a3 1 parent 9151670
Xie Qingguang authored

Showing 20 changed files with 673 additions and 106 deletions. Show diff stats Hide diff stats

  1. +109 0 examples/USER/sph/sdpd-polymer-3D-test/addpolymer.awk
  2. +6 0 examples/USER/sph/sdpd-polymer-3D-test/post.sh
  3. +17 0 examples/USER/sph/sdpd-polymer-3D-test/run3d.sh
  4. +29 0 examples/USER/sph/sdpd-polymer-3D-test/sdpd-polymer3D-inti.lmp
  5. +63 0 examples/USER/sph/sdpd-polymer-3D-test/sdpd-polymer3D-run.lmp
  6. +119 0 examples/USER/sph/sdpd-polymer-auto/addpolymer.awk
  7. +6 0 examples/USER/sph/sdpd-polymer-auto/fitcorr.gp
  8. +8 0 examples/USER/sph/sdpd-polymer-auto/post.sh
  9. +59 0 examples/USER/sph/sdpd-polymer-auto/relaxtime.m
  10. +17 0 examples/USER/sph/sdpd-polymer-auto/run.sh
  11. +65 0 examples/USER/sph/sdpd-polymer-auto/sdpd-polymer-inti.lmp
  12. +70 0 examples/USER/sph/sdpd-polymer-auto/sdpd-polymer-run.lmp
  13. +18 0 examples/USER/sph/sdpd-polymer-auto/vautocor.m
  14. +24 0 examples/USER/sph/sdpd-polymer-auto/vautocov.m
  15. +7 5 examples/USER/sph/sdpd-polymer-wall/run.sh
  16. +10 45 examples/USER/sph/sdpd-polymer-wall/sdpd-polymer-inti.lmp
  17. +28 10 examples/USER/sph/sdpd-polymer-wall/sdpd-polymer-run.lmp
  18. +7 7 examples/USER/sph/sdpd-polymer/sdpd-polymer-run.lmp
  19. +4 4 examples/USER/sph/sdpd-test/sdpd_test_3d.lmp
  20. +7 35 src/USER-SPH/pair_sdpd.cpp
109 examples/USER/sph/sdpd-polymer-3D-test/addpolymer.awk
... ... @@ -0,0 +1,109 @@
  1 +function fabs(var) {
  2 + return var>0?var:-var
  3 +}
  4 +
  5 +function isbound(atom_number, period, rem, current_npoly) {
  6 + period = Nbeads + Nsolvent
  7 + rem = (atom_number-1)%(period) # from 0 to period-1
  8 + current_npoly = int(atom_number/period) + 1
  9 + return (rem<Nbeads-1) && (atom_number<iatom) && (current_npoly<=Npoly)
  10 +}
  11 +
  12 +BEGIN {
  13 + inatoms=0
  14 + lo=1; hi=2
  15 + x=1; y=2; z=3
  16 + iatom=0
  17 + if (Npoly=="full") {
  18 + Npoly = 1e22
  19 + }
  20 +}
  21 +
  22 +/LAMMPS/{
  23 + print
  24 + next
  25 +}
  26 +
  27 +/xlo xhi/{
  28 + box[x,lo]=$1
  29 + box[x,hi]=$1
  30 +}
  31 +
  32 +/ylo yhi/{
  33 + box[y,lo]=$1
  34 + box[y,hi]=$2
  35 +}
  36 +
  37 +/zlo zhi/{
  38 + box[z,lo]=$1
  39 + box[z,hi]=$2
  40 +}
  41 +
  42 +/atom types/{
  43 + print
  44 + print "1 bond types"
  45 + next
  46 +}
  47 +
  48 +/atoms/{
  49 + natoms=$1
  50 + print
  51 + printf("%s bonds\n", "_NUMBER_OF_BOUNDS_")
  52 + next
  53 +}
  54 +
  55 +(NF>0)&&($1=="Atoms"){
  56 + inatoms=1
  57 + print
  58 + # skip empty line
  59 + getline
  60 + printf "\n"
  61 + next
  62 +}
  63 +
  64 +inatoms && (NF==0) {
  65 + inatoms = 0
  66 + print
  67 + next
  68 +}
  69 +
  70 +inatoms{
  71 + # here I get one atom
  72 + iatom++
  73 + R[x]=$4; R[y]=$5; R[z]=$6
  74 + if (iatom>1) {
  75 + for (idim=1; idim<=3; idim++) {
  76 + if (fabs(R[idim]- prevR[idim])>cutoff) {
  77 + if (R[idim]<prevR[idim]) image[idim]++; else image[idim]--
  78 + }
  79 + }
  80 + } else {
  81 + image[x]=0; image[y]=0; image[z]=0
  82 + }
  83 + prevR[x]=R[x]; prevR[y]=R[y]; prevR[z]=R[z]
  84 + # change image field
  85 + $(NF-2)=image[x]; $(NF-1)=image[y]; $(NF)=image[z];
  86 + # add molecule ID
  87 + # $6=$6 " 0"
  88 + print $0
  89 + next
  90 +}
  91 +
  92 +!inatoms {
  93 + print
  94 + next
  95 +}
  96 +
  97 +END {
  98 + printf("\nBonds\n\n")
  99 + ibond = 0
  100 + for (q=1; q<iatom; q++) {
  101 + if (isbound(q)) {
  102 + ibond++
  103 + ip = q
  104 + jp = q+1
  105 + bondtype=1
  106 + print ibond, bondtype, ip, jp
  107 + }
  108 + }
  109 +}
6 examples/USER/sph/sdpd-polymer-3D-test/post.sh
... ... @@ -0,0 +1,6 @@
  1 +#! /bin/bash
  2 +rm -rf punto.dat
  3 +awk 'fl{print $3, $4,$5, $6, $7, $8} /ITEM: ATOMS/{fl=1}' *.dat > punto.dat
  4 +#for filename in dump*.dat; do echo $filename; sed '1,8d' $filename > ${filename/.dat/.del};done
  5 +#awk 'fl{print $1, $2, $3, $4} /ITEM: ATOMS/{fl=1}' *.del > punto.dat
  6 +
17 examples/USER/sph/sdpd-polymer-3D-test/run3d.sh
... ... @@ -0,0 +1,17 @@
  1 +#! /bin/bash
  2 +
  3 +
  4 +
  5 +
  6 +../../../../src/lmp_linux -in sdpd-polymer3D-inti.lmp
  7 +../../../../tools/restart2data poly3d.restart poly3d.txt
  8 +
  9 +
  10 + awk -v cutoff=3.0 -v Nbeads=1 -v Nsolvent=1 -v Npoly=full \
  11 + -f addpolymer.awk poly3d.txt > poly3.txt
  12 + nbound=$(tail -n 1 poly3.txt | awk '{print $1}')
  13 + sed "s/_NUMBER_OF_BOUNDS_/$nbound/1" poly3.txt > poly3d.txt
  14 +
  15 +time /scratch/qingguang/prefix-nana/bin/mpirun -np 1 ../../../../src/lmp_linux -in sdpd-polymer3D-run.lmp
  16 +#mpirun -np 2 ../../../../src/lmp_linux -in sdpd-polymer3D-run.lmp
  17 +#../../../../src/lmp_linux -in sdpd-polymer-run.lmp
29 examples/USER/sph/sdpd-polymer-3D-test/sdpd-polymer3D-inti.lmp
... ... @@ -0,0 +1,29 @@
  1 +
  2 +echo both
  3 +dimension 3
  4 +units si
  5 +atom_style hybrid bond meso
  6 +boundary p p p
  7 +
  8 +
  9 +# create simulation box
  10 +#3D box
  11 +variable L equal 2.0e-4
  12 +region box block 0.0 ${L} 0.0 ${L} 0.0 ${L} units box
  13 +create_box 1 box
  14 +
  15 +# create fluid particles
  16 +lattice fcc 0.025e-3
  17 +#lattice sc 0.025e-3
  18 +create_atoms 1 region box
  19 +
  20 +variable sdpd_mass equal 3.90625e-12
  21 +mass 1 ${sdpd_mass}
  22 +
  23 +variable sdpd_temp index 0.0
  24 +
  25 +
  26 +timestep 0.0e-5
  27 +
  28 +write_restart poly3d.restart
  29 +
63 examples/USER/sph/sdpd-polymer-3D-test/sdpd-polymer3D-run.lmp
... ... @@ -0,0 +1,63 @@
  1 +echo both
  2 +dimension 3
  3 +units si
  4 +atom_style hybrid bond meso
  5 +boundary p p p
  6 +# create simulation box
  7 +#3D box
  8 +variable L equal 2.0e-4
  9 +region box block 0.0 ${L} 0.0 ${L} 0.0 ${L} units box
  10 +read_data poly3d.txt
  11 +
  12 +# create fluid particles
  13 +
  14 +#lattice fcc 0.025e-3
  15 +lattice sc 0.025e-3
  16 +create_atoms 1 region box
  17 +
  18 +#special_bonds lj 1 1 1
  19 +#bond_style fene
  20 +#bond_coeff * 3e-2 5e-5 0 0
  21 +
  22 +
  23 +mass 1 3.90625e-12
  24 +set group all meso_rho 1000.0
  25 +
  26 +# smothing length
  27 +variable h equal 6.5e-5
  28 +
  29 +
  30 +# use Tait's EOS in combination with Morris' laminar viscosity.
  31 +# We set rho_0 = 1000 kg/m^3, c = 0.1 m/s, h = 6.5e-5 m.
  32 +# The dynamic viscosity is set to 1.0e-3 Pa s, corresponding to a kinematic viscosity of 1.0e-6 m^2/s
  33 +
  34 +variable sdpd_temp index 1e1
  35 +pair_style hybrid sdpd
  36 +pair_coeff * * sdpd 1000 0.1 5.0e-3 6.5e-5 ${sdpd_temp}
  37 +
  38 +compute rho_peratom all meso_rho/atom
  39 +compute e_peratom all meso_e/atom
  40 +compute ke_peratom all ke/atom
  41 +compute esph all reduce sum c_e_peratom
  42 +compute ke all ke
  43 +compute sdpd_kin all temp
  44 +variable etot equal c_ke+c_esph
  45 +
  46 +#compute sdpd_kin all temp
  47 +#compute rho_peratom all meso_rho/atom
  48 +
  49 +# do full time integration for shear driver and fluid, but keep walls stationary
  50 +fix integrate_fix_full all meso
  51 +
  52 +dump dump_id all custom 10 dump*.dat id type x y z vx vy vy c_rho_peratom
  53 +dump_modify dump_id first yes
  54 +dump_modify dump_id sort id
  55 +dump_modify dump_id pad 8
  56 +thermo_style custom step c_sdpd_kin
  57 +thermo_modify norm no
  58 +thermo 1
  59 +
  60 +
  61 +neighbor 3.0e-6 bin
  62 +timestep 1e-5
  63 +run 40000
119 examples/USER/sph/sdpd-polymer-auto/addpolymer.awk
... ... @@ -0,0 +1,119 @@
  1 +function fabs(var) {
  2 + return var>0?var:-var
  3 +}
  4 +
  5 +function isbound(atom_number, period, rem, current_npoly) {
  6 + period = Nbeads + Nsolvent
  7 + rem = (atom_number-1)%(period) # from 0 to period-1
  8 + current_npoly = int(atom_number/period) + 1
  9 + return (rem<Nbeads-1) && (atom_number<iatom) && (current_npoly<=Npoly)
  10 +}
  11 +
  12 +BEGIN {
  13 + inatoms=0
  14 + lo=1; hi=2
  15 + x=1; y=2; z=3
  16 + iatom=0
  17 + if (Npoly=="full") {
  18 + Npoly = 1e22
  19 + }
  20 +}
  21 +
  22 +/LAMMPS/{
  23 + print
  24 + next
  25 +}
  26 +
  27 +/xlo xhi/{
  28 + box[x,lo]=$1
  29 + box[x,hi]=$1
  30 +}
  31 +
  32 +/ylo yhi/{
  33 + box[y,lo]=$1
  34 + box[y,hi]=$2
  35 +}
  36 +
  37 +/zlo zhi/{
  38 + box[z,lo]=$1
  39 + box[z,hi]=$2
  40 +}
  41 +
  42 +/atom types/{
  43 + print
  44 + print "1 bond types"
  45 + next
  46 +}
  47 +
  48 +/atoms/{
  49 + natoms=$1
  50 + print
  51 + printf("%s bonds\n", "_NUMBER_OF_BOUNDS_")
  52 + next
  53 +}
  54 +
  55 +(NF>0)&&($1=="Atoms"){
  56 + inatoms=1
  57 + print
  58 + # skip empty line
  59 + getline
  60 + printf "\n"
  61 + next
  62 +}
  63 +
  64 +inatoms && (NF==0) {
  65 + inatoms = 0
  66 + print
  67 + next
  68 +}
  69 +
  70 +inatoms{
  71 + # here I get one atom
  72 + iatom++
  73 + R[x]=$4; R[y]=$5; R[z]=$6
  74 + if (iatom>1) {
  75 + for (idim=1; idim<=3; idim++) {
  76 + if (fabs(R[idim]- prevR[idim])>cutoff) {
  77 + if (R[idim]<prevR[idim]) image[idim]++; else image[idim]--
  78 + }
  79 + }
  80 + } else {
  81 + image[x]=0; image[y]=0; image[z]=0
  82 + }
  83 + prevR[x]=R[x]; prevR[y]=R[y]; prevR[z]=R[z]
  84 + # change image field
  85 + $(NF-2)=image[x]; $(NF-1)=image[y]; $(NF)=image[z];
  86 + # add molecule ID
  87 + # $6=$6 " 0"
  88 + print $0
  89 + next
  90 +}
  91 +
  92 +!inatoms {
  93 + print
  94 + next
  95 +}
  96 +
  97 +END {
  98 + printf("\nBonds\n\n")
  99 + ibond = 0
  100 +ipoly=0
  101 +printf("") > "poly.id"
  102 + for (q=1; q<iatom; q++) {
  103 + if (isbound(q)) {
  104 + ibond++
  105 + ip = q
  106 + jp = q+1
  107 + bondtype=1
  108 + print ibond, bondtype, ip, jp
  109 +if (ip != prev) {
  110 +if (prev>0) {
  111 +print prev, ipoly >> "poly.id"
  112 +}
  113 + ipoly++
  114 + }
  115 + print ip, ipoly >> "poly.id"
  116 + prev=jp
  117 +}
  118 + }
  119 +}
6 examples/USER/sph/sdpd-polymer-auto/fitcorr.gp
... ... @@ -0,0 +1,6 @@
  1 +# gnuplot script to fit autocorrelation data
  2 +f(x) = exp(-x/tau)
  3 +
  4 +tlim=30
  5 +fit [0:tlim] f(x) "corr.dat" via tau
  6 +plot "corr.dat", f(x)
8 examples/USER/sph/sdpd-polymer-auto/post.sh
... ... @@ -0,0 +1,8 @@
  1 +#! /bin/bash
  2 +rm -rf punto.dat
  3 +rm -rf pdata
  4 +awk 'fl{print $3, $4, $5, $6} /ITEM: ATOMS/{fl=1}' *.dat > punto.dat
  5 +awk -f extpolymer.awk dump*.dat
  6 +# write auto-correlation to corr.dat
  7 +octave -q --eval relaxtime
  8 +gnuplot fitcorr.gp
59 examples/USER/sph/sdpd-polymer-auto/relaxtime.m
... ... @@ -0,0 +1,59 @@
  1 +function relaxtime()
  2 +% get polymer configurations
  3 + flist=dir("pdata/poly.*");
  4 +
  5 + nfile = size(flist,1);
  6 + for kk=1:nfile-1
  7 + fname=fullfile("pdata", flist(kk).name);
  8 +
  9 + if (kk==1)
  10 + cr = getonecorr(fname);
  11 + else
  12 + printf("kk=%d\n", kk);
  13 + cr = cr + getonecorr(fname);
  14 + endif
  15 + end
  16 + dtime = 0:size(cr, 1)-1;
  17 + plot (dtime, cr/nfile);
  18 +
  19 + dlmwrite( "corr.dat", [dtime', cr/nfile], ' ');
  20 +
  21 +endfunction
  22 +
  23 +function [corrfun] = getonecorr(fname)
  24 +% get number of beads
  25 + nb=getnbead(fname);
  26 +
  27 + % read data
  28 + data = dlmread(fname);
  29 +
  30 + % keep only x and y
  31 + data = data(:, 1:2);
  32 + data = reshape(data'(:), nb, 2, []);
  33 + R = end2end(data);
  34 +
  35 + corrfun = vautocor(R);
  36 +endfunction
  37 +
  38 +% get the number of beads
  39 +function [nbead] = getnbead(fname)
  40 +fid = fopen (fname);
  41 + nbead=-1;
  42 + do
  43 + line = strtrim(fgetl(fid));
  44 + nbead=nbead+1;
  45 + until strcmp(line, "")
  46 + fclose (fid);
  47 +endfunction
  48 +
  49 +function [R] = end2end(data)
  50 +np=size(data, 3);
  51 + Rhead = zeros(2, np);
  52 + Rtail = zeros(2, np);
  53 + Rhead(:, :) = data(1, :, :);
  54 + Rtail(:, :) = data(end, :, :);
  55 +
  56 + R = Rtail' - Rhead';
  57 +endfunction
  58 +
  59 +
17 examples/USER/sph/sdpd-polymer-auto/run.sh
... ... @@ -0,0 +1,17 @@
  1 +#! /bin/bash
  2 +
  3 +
  4 +rm -rf dump*
  5 +rm -rf image*
  6 +
  7 +../../../../src/lmp_linux -in sdpd-polymer-inti.lmp
  8 +../../../../tools/restart2data poly.restart poly.txt
  9 +
  10 +
  11 + awk -v cutoff=3.0 -v Nbeads=18 -v Nsolvent=18 -v Npoly=full \
  12 + -f addpolymer.awk poly.txt > poly2.txt
  13 + nbound=$(tail -n 1 poly2.txt | awk '{print $1}')
  14 + sed "s/_NUMBER_OF_BOUNDS_/$nbound/1" poly2.txt > poly.txt
  15 +
  16 +time /scratch/qingguang/prefix-nana/bin/mpirun -np 6 ../../../../src/lmp_linux -in sdpd-polymer-run.lmp
  17 +#../../../../src/lmp_linux -in sdpd-polymer-run.lmp
65 examples/USER/sph/sdpd-polymer-auto/sdpd-polymer-inti.lmp
... ... @@ -0,0 +1,65 @@
  1 +echo both
  2 +#atom_modify sort 3.0e-6 3.0e-6
  3 +dimension 2
  4 +units si
  5 +atom_style hybrid bond meso
  6 +boundary p p p
  7 +
  8 +# create simulation box
  9 +#2D box
  10 +variable Lx equal 2.000e-3
  11 +region box block 0.0000e-3 ${Lx} 0.0000e-3 ${Lx} -1.0e-6 1.0e-6 units box
  12 +create_box 1 box
  13 +
  14 +# create fluid particles
  15 +lattice sq 0.025e-3
  16 +create_atoms 1 region box
  17 +
  18 +#set bead number of polymers and free particles
  19 +#attention the relation between sum(Sum) of beadnum with free and Lx/latticelength(N) satisfy N%Sum=0
  20 +#beadnum 5
  21 +#freenum 5
  22 +#string coefficient H
  23 +#Maxiaml distance R0
  24 +
  25 +variable sdpd_mass equal 2.776E-6
  26 +mass 1 ${sdpd_mass}
  27 +#set group all meso_rho 1000.0
  28 +
  29 +# use Tait's EOS in combination with Morris' laminar viscosity.
  30 +# We set rho_0 = 1000 kg/m^3, c = 0.1 m/s, h = 6.5e-5 m.
  31 +# The dynamic viscosity is set to 1.0e-3 Pa s, corresponding to a kinematic viscosity of 1.0e-6 m^2/s
  32 +#The temperature is set to 1e5C.
  33 +variable sdpd_temp index 0.0
  34 +#pair_style hybrid sdpd
  35 +#pair_coeff * * sdpd 1000 0.1 5.0e-3 6.5e-5 ${sdpd_temp}
  36 +
  37 +# compute rho_peratom all meso_rho/atom
  38 +# compute e_peratom all meso_e/atom
  39 +# compute ke_peratom all ke/atom
  40 +# compute esph all reduce sum c_e_peratom
  41 +# compute ke all ke
  42 +# compute sdpd_kin all temp
  43 +# variable etot equal c_ke+c_esph
  44 +
  45 +
  46 +# do full time integration for shear driver and fluid, but keep walls stationary
  47 +#fix integrate_fix_full all meso
  48 +#variable A_kol equal 1.0
  49 +#variable var_kol atom ${A_kol}*${sdpd_mass}*sin(2*PI*y/${Lx})
  50 +#fix kol_force all addforce v_var_kol 0.0 0.0
  51 +
  52 +
  53 +# dump dump_id all custom 100 dump*.dat id type x y vx vy c_rho_peratom c_e_peratom
  54 +# dump_modify dump_id first yes
  55 +# dump_modify dump_id sort id
  56 +# dump_modify dump_id pad 8
  57 +# thermo 10
  58 +# thermo_style custom step c_sdpd_kin
  59 +# thermo_modify norm no
  60 +
  61 +#neighbor 1.0e-6 bin
  62 +timestep 0.0e-5
  63 +write_restart poly.restart
  64 +
  65 +#run 0
70 examples/USER/sph/sdpd-polymer-auto/sdpd-polymer-run.lmp
... ... @@ -0,0 +1,70 @@
  1 +echo both
  2 +#atom_modify sort 3.0e-6 3.0e-6
  3 +dimension 2
  4 +units si
  5 +atom_style hybrid bond meso
  6 +boundary p p p
  7 +
  8 +# create simulation box
  9 +#2D box
  10 +variable Lx equal 2.00e-3
  11 +
  12 +region box block 0.0000e-3 ${Lx} 0.0000e-3 ${Lx} -1.0e-6 1.0e-6 units box
  13 +read_data poly.txt
  14 +
  15 +# create fluid particles
  16 +lattice sq 0.025e-3
  17 +bond_style fene
  18 +special_bonds lj 1 1 1
  19 +#H*ro^2/(kT)=100
  20 +variable H_fene equal 1
  21 +#variable H_fene equal 0.5e-1
  22 +variable r0_fene equal 5e-5
  23 +bond_coeff 1 ${H_fene} ${r0_fene} 0.0 0.0
  24 +
  25 +
  26 +variable sdpd_mass equal 2.778E-6
  27 +mass 1 ${sdpd_mass}
  28 +set group all meso_rho 1000.0
  29 +
  30 +# use Tait's EOS in combination with Morris' laminar viscosity.
  31 +# We set rho_0 = 1000 kg/m^3, c = 0.1 m/s, h = 6.5e-5 m.
  32 +# The dynamic viscosity is set to 1.0e-3 Pa s, corresponding to a kinematic viscosity of 1.0e-6 m^2/s
  33 +#The temperature is set to 1e14.
  34 +variable sdpd_temp index 1e13
  35 +pair_style hybrid sdpd
  36 +pair_coeff * * sdpd 1000 0.1 1.0e-3 6.5e-5 ${sdpd_temp}
  37 +
  38 +compute rho_peratom all meso_rho/atom
  39 +compute e_peratom all meso_e/atom
  40 +compute ke_peratom all ke/atom
  41 +compute esph all reduce sum c_e_peratom
  42 +compute ke all ke
  43 +compute sdpd_kin all temp
  44 +variable etot equal c_ke+c_esph
  45 +
  46 +
  47 +# do full time integration for shear driver and fluid, but keep walls stationary
  48 +fix integrate_fix_full all meso
  49 +#variable A_kol equal 0
  50 +#variable var_kol atom ${A_kol}*${sdpd_mass}*cos(2*PI*y/${Lx})
  51 +#fix kol_force all addforce v_var_kol 0.0 0.0
  52 +
  53 +
  54 +dump dump_id all custom 10000 dump*.dat id type xu yu vx vy c_rho_peratom c_e_peratom
  55 +dump_modify dump_id first yes
  56 +dump_modify dump_id sort id
  57 +dump_modify dump_id pad 8
  58 +thermo 10
  59 +thermo_style custom step c_sdpd_kin
  60 +thermo_modify norm no
  61 +
  62 +dump imgDump all image 10000 image.*.jpg type type bond atom 1e-5 adiam 1e-5 view 0 0 zoom 1.75
  63 +dump_modify imgDump pad 9
  64 +
  65 +neighbor 1.0e-6 bin
  66 +#timestep <0.5*delta/Cs
  67 +#include detime
  68 +timestep 1.0e-6
  69 +
  70 +run 10000000
18 examples/USER/sph/sdpd-polymer-auto/vautocor.m
... ... @@ -0,0 +1,18 @@
  1 +function retval = vautocor (X, h)
  2 +if (nargin == 1)
  3 + retval = vautocov (X);
  4 + elseif (nargin == 2)
  5 + retval = vautocov (X, h);
  6 + else
  7 + print_usage ();
  8 + endif
  9 +
  10 + if (min (retval (1,:)) != 0)
  11 + retval = retval ./ (ones (rows (retval), 1) * retval(1,:));
  12 + endif
  13 +
  14 +endfunction
  15 +
  16 +
  17 +
  18 +
24 examples/USER/sph/sdpd-polymer-auto/vautocov.m
... ... @@ -0,0 +1,24 @@
  1 +function retval = vautocov (X, h)
  2 +[n, c] = size (X);
  3 +
  4 + if (isvector (X))
  5 + n = length (X);
  6 + c = 1;
  7 + X = reshape (X, n, 1);
  8 + endif
  9 +
  10 + X = center (X);
  11 +
  12 + if (nargin == 1)
  13 + h = n - 1;
  14 + endif
  15 +
  16 + retval = zeros (h + 1, 1);
  17 +
  18 + for i = 0 : h
  19 + retval(i+1) = sum(( X(i+1:n, :) .* X(1:n-i, :) )(:)) / n;
  20 + endfor
  21 +
  22 +endfunction
  23 +
  24 +
12 examples/USER/sph/sdpd-polymer-wall/run.sh
@@ -4,14 +4,16 @@
4 4 rm -rf dump*
5 5 rm -rf image*
6 6
7   -../../../../src/lmp_linux -in sdpd-polymer-inti.lmp
  7 +dx=0.025e-3
  8 +lmp=../../../../src/lmp_linux
  9 +${lmp} -var dx ${dx} -in sdpd-polymer-inti.lmp
8 10 ../../../../tools/restart2data poly.restart poly.txt
9 11
10   -
11   - awk -v cutoff=3.0 -v Nbeads=18 -v Nsolvent=18 -v Npoly=full \
  12 + awk -v cutoff=3.0 -v Nbeads=9 -v Nsolvent=9 -v Npoly=full \
12 13 -f addpolymer.awk poly.txt > poly2.txt
13 14 nbound=$(tail -n 1 poly2.txt | awk '{print $1}')
14 15 sed "s/_NUMBER_OF_BOUNDS_/$nbound/1" poly2.txt > poly.txt
15 16
16   -time /scratch/qingguang/prefix-nana/bin/mpirun -np 6 ../../../../src/lmp_linux -in sdpd-polymer-run.lmp
17   -#../../../../src/lmp_linux -in sdpd-polymer-run.lmp
  17 +time /scratch/qingguang/prefix-nana/bin/mpirun -np 1 \
  18 + ${lmp} -var dx ${dx} -in sdpd-polymer-run.lmp
  19 +
55 examples/USER/sph/sdpd-polymer-wall/sdpd-polymer-inti.lmp
@@ -3,63 +3,28 @@ echo both
3 3 dimension 2
4 4 units si
5 5 atom_style hybrid bond meso
6   -boundary p p p
  6 +boundary p f p
7 7
8 8 # create simulation box
9 9 #2D box
10 10 variable Lx equal 2.000e-3
11   -region box block 0.0000e-3 ${Lx} 0.0000e-3 ${Lx} -1.0e-6 1.0e-6 units box
12   -create_box 1 box
  11 +variable ub equal ${Lx}-3.5*${dx}
  12 +variable lb equal 2.5*${dx}
  13 +region box block 0.0000e-3 ${Lx} 0.0000e-3 ${Lx} -1.0e-6 1.0e-6 units box
  14 +region fluid block 0.0000e-3 ${Lx} ${lb} ${ub} -1.0e-6 1.0e-6 units box
  15 +create_box 2 box
13 16
14   -# create fluid particles
15   -lattice sq 0.025e-3
16   -create_atoms 1 region box
  17 +#create fluid particles
  18 +lattice sq ${dx}
  19 +create_atoms 1 region fluid
17 20
18   -#set bead number of polymers and free particles
19   -#attention the relation between sum(Sum) of beadnum with free and Lx/latticelength(N) satisfy N%Sum=0
20   -#beadnum 5
21   -#freenum 5
22   -#string coefficient H
23   -#Maxiaml distance R0
24 21
25 22 variable sdpd_mass equal 2.776E-6
26 23 mass 1 ${sdpd_mass}
  24 +mass 2 ${sdpd_mass}
27 25 #set group all meso_rho 1000.0
28 26
29   -# use Tait's EOS in combination with Morris' laminar viscosity.
30   -# We set rho_0 = 1000 kg/m^3, c = 0.1 m/s, h = 6.5e-5 m.
31   -# The dynamic viscosity is set to 1.0e-3 Pa s, corresponding to a kinematic viscosity of 1.0e-6 m^2/s
32   -#The temperature is set to 1e5C.
33   -variable sdpd_temp index 0.0
34   -#pair_style hybrid sdpd
35   -#pair_coeff * * sdpd 1000 0.1 5.0e-3 6.5e-5 ${sdpd_temp}
36 27
37   -# compute rho_peratom all meso_rho/atom
38   -# compute e_peratom all meso_e/atom
39   -# compute ke_peratom all ke/atom
40   -# compute esph all reduce sum c_e_peratom
41   -# compute ke all ke
42   -# compute sdpd_kin all temp
43   -# variable etot equal c_ke+c_esph
44   -
45   -
46   -# do full time integration for shear driver and fluid, but keep walls stationary
47   -#fix integrate_fix_full all meso
48   -#variable A_kol equal 1.0
49   -#variable var_kol atom ${A_kol}*${sdpd_mass}*sin(2*PI*y/${Lx})
50   -#fix kol_force all addforce v_var_kol 0.0 0.0
51   -
52   -
53   -# dump dump_id all custom 100 dump*.dat id type x y vx vy c_rho_peratom c_e_peratom
54   -# dump_modify dump_id first yes
55   -# dump_modify dump_id sort id
56   -# dump_modify dump_id pad 8
57   -# thermo 10
58   -# thermo_style custom step c_sdpd_kin
59   -# thermo_modify norm no
60   -
61   -#neighbor 1.0e-6 bin
62 28 timestep 0.0e-5
63 29 write_restart poly.restart
64 30
65   -#run 0
38 examples/USER/sph/sdpd-polymer-wall/sdpd-polymer-run.lmp
@@ -3,31 +3,45 @@ echo both
3 3 dimension 2
4 4 units si
5 5 atom_style hybrid bond meso
6   -boundary p p p
  6 +boundary p f p
7 7
8 8 # create simulation box
9 9 #2D box
10 10 variable Lx equal 2.00e-3
11 11 read_data poly.txt
12 12
  13 +# create wall particles
  14 +variable ub equal ${Lx}-3.5*${dx}
  15 +variable lb equal 2.5*${dx}
  16 +region uwall block 0.0000e-3 ${Lx} ${ub} ${Lx} -1.0e-6 1.0e-6 units box
  17 +region lwall block 0.0000e-3 ${Lx} 0.0000e-3 ${lb} -1.0e-6 1.0e-6 units box
  18 +lattice sq ${dx}
  19 +create_atoms 2 region uwall
  20 +create_atoms 2 region lwall
  21 +
13 22 # create fluid particles
14 23 bond_style fene
15   -special_bonds lj 1 1 1i
  24 +special_bonds lj 1 1 1
16 25 #H*ro^2/(kT)=100
17 26 variable H_fene equal 5
18   -#variable H_fene equal 0.5e-1
19 27 variable r0_fene equal 5e-5
20 28 bond_coeff 1 ${H_fene} ${r0_fene} 0.0 0.0
21 29
22 30
  31 +group uwall region uwall
  32 +group lwall region lwall
  33 +#group fluid region fluid
  34 +#group integrate_fulll union fluid lwall
  35 +
23 36 variable sdpd_mass equal 2.778E-6
24 37 mass 1 ${sdpd_mass}
  38 +mass 2 ${sdpd_mass}
25 39 set group all meso_rho 1000.0
26 40
27 41 # use Tait's EOS in combination with Morris' laminar viscosity.
28 42 # We set rho_0 = 1000 kg/m^3, c = 0.1 m/s, h = 6.5e-5 m.
29 43 # The dynamic viscosity is set to 1.0e-3 Pa s, corresponding to a kinematic viscosity of 1.0e-6 m^2/s
30   -#The temperature is set to 1e14.
  44 +#The temperature is set to .
31 45 variable sdpd_temp index 1e13
32 46 pair_style hybrid sdpd
33 47 pair_coeff * * sdpd 1000 0.1 1.0e-3 6.5e-5 ${sdpd_temp}
@@ -43,13 +57,17 @@ variable etot equal c_ke+c_esph
43 57
44 58 # do full time integration for shear driver and fluid, but keep walls stationary
45 59 fix integrate_fix_full all meso
46   -variable A_kol equal 0.1
47   -variable var_kol atom ${A_kol}*${sdpd_mass}*cos(2*PI*y/${Lx})
48   -fix kol_force all addforce v_var_kol 0.0 0.0
  60 +#variable A_kol equal 0.1
  61 +#variable var_kol atom ${A_kol}*${sdpd_mass}*cos(2*PI*y/${Lx})
  62 +#fix kol_force all addforce v_var_kol 0.0 0.0
  63 +velocity uwall set 0.01 0.0 0.0 units box
  64 +velocity lwall set -0.01 0.0 0.0 units box
  65 +fix freeze_fixu uwall setforce 0.0 0.0 0.0
  66 +fix freeze_fixl lwall setforce 0.0 0.0 0.0
49 67
50 68
51 69 #dump dump_id all custom 10000 dump*.dat id type x y vx vy c_rho_peratom c_e_peratom
52   -dump dump_id all custom 1000 dump*.dat x y vx vy id type c_rho_peratom
  70 +dump dump_id all custom 100 dump*.dat x y vx vy id type c_rho_peratom
53 71 dump_modify dump_id first yes
54 72 dump_modify dump_id sort id
55 73 dump_modify dump_id pad 8
@@ -57,7 +75,7 @@ thermo 10
57 75 thermo_style custom step c_sdpd_kin
58 76 thermo_modify norm no
59 77
60   -dump imgDump all image 1000 image.*.jpg type type bond atom 1e-5 adiam 1e-5 view 0 0 zoom 1.75
  78 +dump imgDump all image 100 image.*.jpg type type bond atom 1e-5 adiam 1e-5 view 0 0 zoom 1.75 box no 0.01
61 79 dump_modify imgDump pad 9
62 80
63 81 neighbor 1.0e-6 bin
@@ -65,4 +83,4 @@ neighbor 1.0e-6 bin
65 83 #include detime
66 84 timestep 1.0e-5
67 85
68   -run 1000000000
  86 +run 10000000
14 examples/USER/sph/sdpd-polymer/sdpd-polymer-run.lmp
@@ -15,7 +15,7 @@ read_data poly.txt
15 15 # create fluid particles
16 16 lattice sq 0.025e-3
17 17 bond_style fene
18   -special_bonds lj 1 1 1i
  18 +special_bonds lj 1 1 1
19 19 #H*ro^2/(kT)=100
20 20 variable H_fene equal 5
21 21 #variable H_fene equal 0.5e-1
@@ -31,7 +31,7 @@ set group all meso_rho 1000.0
31 31 # We set rho_0 = 1000 kg/m^3, c = 0.1 m/s, h = 6.5e-5 m.
32 32 # The dynamic viscosity is set to 1.0e-3 Pa s, corresponding to a kinematic viscosity of 1.0e-6 m^2/s
33 33 #The temperature is set to 1e14.
34   -variable sdpd_temp index 1e13
  34 +variable sdpd_temp index 1e14
35 35 pair_style hybrid sdpd
36 36 pair_coeff * * sdpd 1000 0.1 1.0e-3 6.5e-5 ${sdpd_temp}
37 37
@@ -46,13 +46,13 @@ variable etot equal c_ke+c_esph
46 46
47 47 # do full time integration for shear driver and fluid, but keep walls stationary
48 48 fix integrate_fix_full all meso
49   -variable A_kol equal 0.1
  49 +variable A_kol equal 5
50 50 variable var_kol atom ${A_kol}*${sdpd_mass}*cos(2*PI*y/${Lx})
51 51 fix kol_force all addforce v_var_kol 0.0 0.0
52 52
53 53
54 54 #dump dump_id all custom 10000 dump*.dat id type x y vx vy c_rho_peratom c_e_peratom
55   -dump dump_id all custom 1000 dump*.dat x y vx vy id type c_rho_peratom
  55 +dump dump_id all custom 100000 dump*.dat x y vx vy id type c_rho_peratom
56 56 dump_modify dump_id first yes
57 57 dump_modify dump_id sort id
58 58 dump_modify dump_id pad 8
@@ -60,12 +60,12 @@ thermo 10
60 60 thermo_style custom step c_sdpd_kin
61 61 thermo_modify norm no
62 62
63   -dump imgDump all image 1000 image.*.jpg type type bond atom 1e-5 adiam 1e-5 view 0 0 zoom 1.75
  63 +dump imgDump all image 100000 image.*.jpg type type bond atom 1e-5 adiam 1e-5 view 0 0 zoom 1.75
64 64 dump_modify imgDump pad 9
65 65
66 66 neighbor 1.0e-6 bin
67 67 #timestep <0.5*delta/Cs
68 68 #include detime
69   -timestep 1.0e-5
  69 +timestep 1.0e-8
70 70
71   -run 1000000000
  71 +run 100000000
8 examples/USER/sph/sdpd-test/sdpd_test_3d.lmp
@@ -15,8 +15,8 @@ lattice fcc 0.025e-3
15 15 create_atoms 1 region box
16 16
17 17
18   -#mass 1 3.90625e-12
19   -mass 1 3.90625e-17
  18 +mass 1 3.90625e-12
  19 +#mass 1 3.90625e-17
20 20 set group all meso_rho 1000.0
21 21
22 22 # smothing length
@@ -27,7 +27,7 @@ variable h equal 6.5e-5
27 27 # We set rho_0 = 1000 kg/m^3, c = 0.1 m/s, h = 6.5e-5 m.
28 28 # The dynamic viscosity is set to 1.0e-3 Pa s, corresponding to a kinematic viscosity of 1.0e-6 m^2/s
29 29 pair_style hybrid/overlay sph/rhosum 1 sdpd
30   -variable sdpd_temp index 1.0
  30 +variable sdpd_temp index 1e10
31 31 pair_coeff * * sph/rhosum ${h}
32 32 pair_coeff * * sdpd 1000 1e-4 1e-4 ${h} ${sdpd_temp}
33 33
@@ -46,5 +46,5 @@ thermo 1
46 46
47 47
48 48 neighbor 3.0e-6 bin
49   -timestep 5e-4
  49 +timestep 1e-6
50 50 run 40000
42 src/USER-SPH/pair_sdpd.cpp
@@ -123,6 +123,8 @@ if (first) {
123 123
124 124 // compute pressure of atom i with Tait EOS
125 125 tmp = rho[i] / rho0[itype];
  126 +//std::cerr<<"rho[i]"<<rho[i]<<'\n';
  127 +
126 128 fi = tmp * tmp * tmp;
127 129 fi = B[itype] * (fi * fi * tmp - 1.0) / (rho[i] * rho[i]);
128 130
@@ -166,39 +168,8 @@ if (first) {
166 168
167 169 // dot product of velocity delta and distance vector
168 170 delVdotDelR = delx * velx + dely * vely + delz * velz;
169   -/*---------------
170   -//create polymer ID in Lammps
171   -constant int beadnum=sdpd_bead[itype][jtype];
172   -constant int freenum=sdpd_free[itype][jtype];
173   -
174   -if(list->inum%(beadnum+freenum)<(beadnum+1)&&list->inum%(beadnum+freenum)>0)
175   -polyID->i=atom->ID;
176   -else
177   -{polyID->i=0;}
178   -
179   -//FENE force between polymer beads
180   -double fene[ndim];
181   -if (PolyID->i>0&& PolyID->j>0)
182   -{
183   -if (PolyID->i-polyID->j==1)
184   -{
185   -if (ndim==2)
186   -{
187   -fene[0]=polymet_H/(1-rsq)*(delx);
188   -fene[1]=polymet_H/(1-rsq)*(dely);
189   -}
190   -else
191   -{
192   -fene[0]=polymet_H/(1-(rsq/R)*(rsq/R))*delx;
193   -fene[1]=polymet_H/(1-(rsq/R)*(rsq/R))*dely;
194   -fene[2]=polymet_H/(1-(rsq/R)*(rsq/R))*delz;
195   -}
196   -}
197   -}
198   -*/
199   -
200   -
201   - // Morris Viscosity (Morris, 1996)
  171 +
  172 + // Morris Viscosity (Morris, 1996)
202 173
203 174 if (ndim==2)
204 175 {
@@ -262,8 +233,9 @@ fene[2]=polymet_H/(1-(rsq/R)*(rsq/R))*delz;
262 233 f[j][0] -= delx*fpair + velx*fvisc + _dUi[0];
263 234 f[j][1] -= dely*fpair + vely*fvisc + _dUi[1];
264 235 if (domain->dimension ==3 ) {
265   - f[j][2] -= delz*fpair + velz*fvisc + _dUi[2];
266   -
  236 +
  237 + f[j][2] -= delz*fpair + velz*fvisc + _dUi[2];
  238 +
267 239 }
268 240 de[j] += deltaE;
269 241 drho[j] += imass * delVdotDelR * wfd;

0 comments on commit 88acc93

Please sign in to comment.
Something went wrong with that request. Please try again.