forked from lijiext/lammps
clean up whitespace, indentation, isolate all cauchystat-related code, remove extraneous examples
This commit is contained in:
Ronald E. Miller
Axel Kohlmeyer
parent
43f83abe5e
commit
5c686be33d
|
|
@ -1,69 +0,0 @@
|
|||
units metal
|
||||
atom_style atomic
|
||||
atom_modify map array
|
||||
|
||||
processors 1 1 1
|
||||
|
||||
# Box and atom positions:
|
||||
boundary p p p
|
||||
|
||||
# Defining lattice and creating simulation
|
||||
# box with atoms inside
|
||||
lattice fcc 4.05
|
||||
region simbox prism 0 5 0 5 0 5 0 0 0 units lattice
|
||||
create_box 2 simbox
|
||||
create_atoms 2 box
|
||||
|
||||
# Atomic mass:
|
||||
mass 1 58.69
|
||||
mass 2 26.98154
|
||||
|
||||
# Potential, Al fcc crystal
|
||||
pair_style eam/alloy
|
||||
pair_coeff * * ../Mishin-Ni-Al-2009.eam.alloy Ni Al
|
||||
|
||||
|
||||
thermo 100
|
||||
thermo_style custom step temp pxx pyy pzz pxy pxz pyz
|
||||
compute cna all cna/atom 2.8
|
||||
|
||||
fix 1 all npt temp 600.0 600.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 0.0 0.1 couple none
|
||||
|
||||
dump 1 all cfg 1000 test*.cfg mass type xs ys zs type c_cna
|
||||
|
||||
timestep 0.002
|
||||
|
||||
variable px equal pxx
|
||||
variable py equal pyy
|
||||
variable pz equal pzz
|
||||
variable sxy equal pxy
|
||||
variable sxz equal pxz
|
||||
variable syz equal pyz
|
||||
variable t equal temp
|
||||
|
||||
fix avg all ave/time 1 100 100 v_t v_px v_py v_pz v_sxy v_sxz v_syz file avg.txt
|
||||
|
||||
variable lx equal lx
|
||||
variable ly equal ly
|
||||
variable lz equal ly
|
||||
variable xy equal xy
|
||||
variable xz equal xz
|
||||
variable yz equal yz
|
||||
|
||||
fix box all ave/time 1 100 100 v_lx v_ly v_lz v_xy v_xz v_yz file box.txt
|
||||
|
||||
velocity all create 1200 4928459 rot yes dist gaussian
|
||||
|
||||
run 10000
|
||||
|
||||
unfix 1
|
||||
|
||||
fix 1 all npt temp 600.0 600.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 0.0 0.1 xy -10000.0 -10000.0 0.1 couple none
|
||||
|
||||
run 10000
|
||||
|
||||
unfix 1
|
||||
|
||||
fix 1 all npt temp 600.0 600.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 0.0 0.1 xy -10000.0 -10000.0 0.1 couple none
|
||||
|
||||
run 10000
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -1,45 +0,0 @@
|
|||
units metal
|
||||
atom_style atomic
|
||||
atom_modify map array
|
||||
|
||||
processors 1 1 1
|
||||
|
||||
# Box and atom positions:
|
||||
boundary p p p
|
||||
read_data input.dat
|
||||
|
||||
# Atomic mass:
|
||||
mass 1 58.69
|
||||
mass 2 26.98154
|
||||
|
||||
# Potential, Al fcc crystal
|
||||
pair_style eam/alloy
|
||||
pair_coeff * * ../Mishin-Ni-Al-2009.eam.alloy Ni Al
|
||||
|
||||
thermo 100
|
||||
thermo_style custom step temp pxx pyy pzz lx ly lz
|
||||
compute cna all cna/atom 2.8
|
||||
|
||||
velocity all create 2400.0 4928459 rot yes dist gaussian
|
||||
|
||||
fix 1 all npt temp 1200.0 1200.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z -2800.0 -2800.0 0.1 couple none cauchystat 0.001 no
|
||||
|
||||
dump 1 all cfg 1000 test*.cfg mass type xs ys zs type c_cna
|
||||
|
||||
timestep 0.002
|
||||
|
||||
variable l equal lz
|
||||
variable p equal pzz
|
||||
variable t equal temp
|
||||
|
||||
fix avg all ave/time 1 1000 1000 v_t v_l v_p file avg.txt
|
||||
|
||||
velocity all create 2400 4928459 rot yes dist gaussian
|
||||
|
||||
run 10000
|
||||
|
||||
unfix 1
|
||||
|
||||
fix 1 all npt temp 1200.0 300.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z -2800.0 -2800.0 0.1 couple none cauchystat 0.001 yes
|
||||
|
||||
run 100000
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -1,45 +0,0 @@
|
|||
units metal
|
||||
atom_style atomic
|
||||
atom_modify map array
|
||||
|
||||
processors 1 1 1
|
||||
|
||||
# Box and atom positions:
|
||||
boundary p p p
|
||||
read_data input.dat
|
||||
|
||||
# Atomic mass:
|
||||
mass 1 58.69
|
||||
mass 2 26.98154
|
||||
|
||||
# Potential, Al fcc crystal
|
||||
pair_style eam/alloy
|
||||
pair_coeff * * ../Mishin-Ni-Al-2009.eam.alloy Ni Al
|
||||
|
||||
thermo 100
|
||||
thermo_style custom step temp pxx pyy pzz lx ly lz
|
||||
compute cna all cna/atom 2.8
|
||||
|
||||
velocity all create 1000.0 4928459 rot yes dist gaussian
|
||||
|
||||
fix 1 all npt temp 500.0 500.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 0.0 0.1 couple none cauchystat 0.001 no
|
||||
|
||||
dump 1 all cfg 5000 test*.cfg mass type xs ys zs type c_cna
|
||||
|
||||
timestep 0.002
|
||||
|
||||
variable l equal lz
|
||||
variable p equal pzz
|
||||
variable t equal temp
|
||||
|
||||
fix avg all ave/time 1 1000 1000 v_t v_l v_p file avg.txt
|
||||
|
||||
velocity all create 2400 4928459 rot yes dist gaussian
|
||||
|
||||
run 10000
|
||||
|
||||
unfix 1
|
||||
|
||||
fix 1 all npt temp 500.0 500.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 -14000.0 0.1 couple none cauchystat 0.001 yes
|
||||
|
||||
run 350000
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -1,45 +0,0 @@
|
|||
units metal
|
||||
atom_style atomic
|
||||
atom_modify map array
|
||||
|
||||
processors 1 1 1
|
||||
|
||||
# Box and atom positions:
|
||||
boundary p p p
|
||||
read_data input.dat
|
||||
|
||||
# Atomic mass:
|
||||
mass 1 58.69
|
||||
mass 2 26.98154
|
||||
|
||||
# Potential, Al fcc crystal
|
||||
pair_style eam/alloy
|
||||
pair_coeff * * ../Mishin-Ni-Al-2009.eam.alloy Ni Al
|
||||
|
||||
thermo 100
|
||||
thermo_style custom step temp pxx pyy pzz lx ly lz
|
||||
compute cna all cna/atom 2.8
|
||||
|
||||
velocity all create 1000.0 4928459 rot yes dist gaussian
|
||||
|
||||
fix 1 all npt temp 500.0 500.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 0.0 0.1 couple none
|
||||
|
||||
dump 1 all cfg 5000 test*.cfg mass type xs ys zs type c_cna
|
||||
|
||||
timestep 0.002
|
||||
|
||||
variable l equal lz
|
||||
variable p equal pzz
|
||||
variable t equal temp
|
||||
|
||||
fix avg all ave/time 1 1000 1000 v_t v_l v_p file avg.txt
|
||||
|
||||
velocity all create 2400 4928459 rot yes dist gaussian
|
||||
|
||||
run 10000
|
||||
|
||||
unfix 1
|
||||
|
||||
fix 1 all npt temp 500.0 500.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 -14000.0 0.1 couple none
|
||||
|
||||
run 350000
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -1,48 +0,0 @@
|
|||
units metal
|
||||
atom_style atomic
|
||||
atom_modify map array
|
||||
|
||||
processors 1 1 1
|
||||
|
||||
# Box and atom positions:
|
||||
boundary p p p
|
||||
read_data input.dat
|
||||
|
||||
# Atomic mass:
|
||||
mass 1 58.69
|
||||
mass 2 26.98154
|
||||
|
||||
# Potential, Al fcc crystal
|
||||
pair_style eam/alloy
|
||||
pair_coeff * * ../Mishin-Ni-Al-2009.eam.alloy Ni Al
|
||||
|
||||
|
||||
thermo 100
|
||||
thermo_style custom step temp pxx pyy pzz lx ly lz
|
||||
compute cna all cna/atom 2.8
|
||||
|
||||
velocity all create 2400.0 4928459 rot yes dist gaussian
|
||||
|
||||
fix 1 all npt temp 600.0 600.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 0.0 0.1 couple none cauchystat 0.001 no
|
||||
|
||||
dump 1 all cfg 1000 test*.cfg mass type xs ys zs type c_cna
|
||||
|
||||
timestep 0.002
|
||||
|
||||
variable l equal lz
|
||||
variable p equal pzz
|
||||
variable t equal temp
|
||||
|
||||
fix avg all ave/time 1 1000 1000 v_t v_l v_p file avg.txt
|
||||
|
||||
#restart 1000 restart.*.test
|
||||
|
||||
velocity all create 2400 4928459 rot yes dist gaussian
|
||||
|
||||
run 10000
|
||||
|
||||
unfix 1
|
||||
|
||||
fix 1 all npt temp 600.0 600.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 -6000.0 0.1 couple none cauchystat 0.001 yes
|
||||
|
||||
run 100000
|
||||
|
|
@ -1,42 +1,10 @@
|
|||
There are 6 examples in this directory. All are run by executing:
|
||||
Run this example by executing:
|
||||
|
||||
% lmp_serial < lammps.in
|
||||
|
||||
from the respective home directory.
|
||||
|
||||
Note that these examples use an EAM potential, and therefore must be
|
||||
Note that this example use an EAM potential, and therefore must be
|
||||
run with a LAMMPS executable built with the MANYBODY package.
|
||||
|
||||
Al_shear_PR: Application of a simple shear stress, illustrating that
|
||||
there is a non-zero direct stress in the equilibrated system due to
|
||||
the fact that Parinello-Rahman controls the second Piola-Kirchhoff
|
||||
stress instead of the Cauchy stress. First fix equilibrates the
|
||||
temperature at zero stress, second fix applies the shear. The third
|
||||
fix corrects the problem because each new run resets the reference
|
||||
simulation cell (H0) to the current cell shape.
|
||||
|
||||
Al_shear_CS: Cauchystat example. Same as above, using the Cauchystat
|
||||
instead of PR control. In this case, the equilibrated stresses are as
|
||||
prescribed (ie., the set stresses control the Cauchy stress, not Piola
|
||||
Kirchhoff)
|
||||
|
||||
NiAl_cool: Cauchystat example. First fix equilibrates the system to
|
||||
uniaxial tension of 2800 bars, 1200K. Second fix cools to 300 K,
|
||||
causing a phase transformation. Cauchy stress remains at the set
|
||||
level after the transformation
|
||||
|
||||
NiAl_stretch: Cauchystat example. After a brief equilibration
|
||||
at 600 K and zero stress, quickly ramp up the applied tensile stress
|
||||
from 0 to 6000 bars over 100000 time steps. After the phase
|
||||
transformation, Cauchy stress continues to follow the set value.
|
||||
|
||||
NiAl_slow_stretch_CS: Cauchystat example. After a brief equilibration
|
||||
at 500 K and zero stress, gradually ramp up the applied tensile stress
|
||||
from 0 to 14000 bars over 350000 time steps. After the phase
|
||||
transformation, Cauchy stress continues to follow the set value.
|
||||
|
||||
NiAl_slow_stretch_PR: After a brief equilibration at 500 K and zero
|
||||
stress, gradually ramp up the applied tensile stress from 0 to 14000
|
||||
bars over 350000 time steps. After the phase transformation, Cauchy
|
||||
stress departs substantially from the set value, because the PR
|
||||
control is used instead of the Cauchystat.
|
||||
The first cauchystat fix equilibrates the temperature at zero stress,
|
||||
the second fix applies a shear stress. Output in avg.txt shows
|
||||
convergence to correct values.
|
||||
|
|
|
|||
|
|
@ -20,7 +20,7 @@ mass 2 26.98154
|
|||
|
||||
# Potential, Al fcc crystal
|
||||
pair_style eam/alloy
|
||||
pair_coeff * * ../Mishin-Ni-Al-2009.eam.alloy Ni Al
|
||||
pair_coeff * * Mishin-Ni-Al-2009.eam.alloy Ni Al
|
||||
|
||||
|
||||
thermo 100
|
||||
|
|
@ -61,9 +61,3 @@ unfix 1
|
|||
fix 1 all npt temp 600.0 600.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 0.0 0.1 xy -10000.0 -10000.0 0.1 couple none cauchystat 0.001 yes
|
||||
|
||||
run 10000
|
||||
|
||||
unfix 1
|
||||
|
||||
fix 1 all npt temp 600.0 600.0 1.0 x 0.0 0.0 0.1 y 0.0 0.0 0.1 z 0.0 0.0 0.1 xy -10000.0 -10000.0 0.1 couple none
|
||||
|
||||
run 10000
|
||||
622
src/fix_nh.cpp
622
src/fix_nh.cpp
|
|
@ -347,77 +347,15 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) :
|
|||
} else error->all(FLERR,"Illegal fix nvt/npt/nph command");
|
||||
iarg += 2;
|
||||
} else if (strcmp(arg[iarg],"cauchystat") == 0) {
|
||||
if (iarg+3 > narg) error->all(FLERR,"Illegal fix npt cauchystat command: wrong number of arguments");
|
||||
if (strcmp(arg[iarg+2],"yes") != 0 && strcmp(arg[iarg+2],"no") != 0) error->all(FLERR,"Illegal cauchystat continue value. Must be 'yes' or 'no'");
|
||||
if (iarg+3 > narg) error->all(FLERR,
|
||||
"Illegal fix npt cauchystat command:"
|
||||
" wrong number of arguments");
|
||||
if (strcmp(arg[iarg+2],"yes") != 0 && strcmp(arg[iarg+2],"no") != 0)
|
||||
error->all(FLERR,"Illegal cauchystat continue value. Must be 'yes' or 'no'");
|
||||
alpha = force->numeric(FLERR,arg[iarg+1]);
|
||||
restartPK = !strcmp(arg[iarg+2],"yes");
|
||||
|
||||
|
||||
if (comm->me == 0) {
|
||||
if (screen) {
|
||||
fprintf(screen,"Using the Cauchystat fix with alpha=%f\n",alpha);
|
||||
if(restartPK==1) {
|
||||
fprintf(screen," (this is a continuation fix)\n");
|
||||
}
|
||||
else {
|
||||
fprintf(screen," (this is NOT a continuation fix)\n");
|
||||
}
|
||||
}
|
||||
if (logfile) {
|
||||
fprintf(logfile,"Using the Cauchystat with alpha=%f\n",alpha);
|
||||
if(restartPK==1) {
|
||||
fprintf(logfile," this is a continuation run\n");
|
||||
}
|
||||
else {
|
||||
fprintf(logfile," this is NOT a continuation run\n");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
if(restartPK==1 && restart_stored==0) {
|
||||
error->all(FLERR,"Illegal cauchystat command. Continuation run must follow a previously equilibrated Cauchystat run");
|
||||
}
|
||||
|
||||
if(alpha<=0.0) {
|
||||
error->all(FLERR,"Illegal cauchystat command. Alpha cannot be zero or negative");
|
||||
}
|
||||
|
||||
initRUN = 0;
|
||||
initPK = 1;
|
||||
usePK = 0;
|
||||
|
||||
#define H0(row,col) (H0[(row-1)][(col-1)])
|
||||
#define invH0(row,col) (invH0[(row-1)][(col-1)])
|
||||
|
||||
|
||||
// initialize H0 to current box shape
|
||||
|
||||
int triclinic = domain->triclinic;
|
||||
double *h = domain->h;
|
||||
double *h_inv = domain->h_inv;
|
||||
double *boxhi = domain->boxhi;
|
||||
double *boxlo = domain->boxlo;
|
||||
double yz = domain->yz;
|
||||
double xz = domain->xz;
|
||||
double xy = domain->xy;
|
||||
|
||||
|
||||
H0(1,1)=h[0]; H0(1,2)=h[5]; H0(1,3)=h[4];
|
||||
H0(2,1)=0.0; H0(2,2)=h[1]; H0(2,3)=h[3];
|
||||
H0(3,1)=0.0; H0(3,2)=0.0; H0(3,3)=h[2];
|
||||
|
||||
invH0(1,1)=h_inv[0]; invH0(1,2)=h_inv[5]; invH0(1,3)=h_inv[4];
|
||||
invH0(2,1)=0.0; invH0(2,2)=h_inv[1]; invH0(2,3)=h_inv[3];
|
||||
invH0(3,1)=0.0; invH0(3,2)=0.0; invH0(3,3)=h_inv[2];
|
||||
|
||||
myvol0 = abs(MathExtra::det3(H0)); //find reference volume
|
||||
|
||||
#undef H0
|
||||
#undef invH0
|
||||
|
||||
CauchyStat_init();
|
||||
iarg += 3;
|
||||
|
||||
} else if (strcmp(arg[iarg],"fixedpoint") == 0) {
|
||||
if (iarg+4 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
|
||||
fixedpoint[0] = force->numeric(FLERR,arg[iarg+1]);
|
||||
|
|
@ -2307,170 +2245,19 @@ void FixNH::compute_press_target()
|
|||
for (int i = 3; i < 6; i++)
|
||||
p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]);
|
||||
|
||||
// CauchyStat: call CauchyStat to modify p_target[i] and p_hydro, they are used in compute_sigma()
|
||||
// CauchyStat: call CauchyStat to modify p_target[i] and p_hydro,
|
||||
// if CauchyStat enabled and pressure->vector computation has been initiated
|
||||
if( (usePK==0) && (initRUN==1) ) {
|
||||
|
||||
double* h = domain->h; // shape matrix in Voigt notation
|
||||
double* h_rate = domain->h_rate; // rate of box size/shape change in Voigt notation
|
||||
double H[3][3];
|
||||
double Hdot[3][3];
|
||||
|
||||
#define H(row,col) (H[(row-1)][(col-1)])
|
||||
#define Hdot(row,col) (Hdot[(row-1)][(col-1)])
|
||||
#define F(row,col) (F[(row-1)][(col-1)])
|
||||
#define Fi(row,col) (Fi[(row-1)][(col-1)])
|
||||
#define Fdot(row,col) (Fdot[(row-1)][(col-1)])
|
||||
#define invH0(row,col) (invH0[(row-1)][(col-1)])
|
||||
#define cauchy(row,col) (cauchy[(row-1)][(col-1)])
|
||||
#define setcauchy(row,col) (setcauchy[(row-1)][(col-1)])
|
||||
#define setPK(row,col) (setPK[(row-1)][(col-1)])
|
||||
#define sigmahat(row,col) (sigmahat[(row-1)][(col-1)])
|
||||
|
||||
H(1,1)=h[0]; H(1,2)=0.0; H(1,3)=0.0;
|
||||
H(2,1)=0.0; H(2,2)=h[1]; H(2,3)=0.0;
|
||||
H(3,1)=0.0; H(3,2)=0.0; H(3,3)=h[2];
|
||||
|
||||
for (int i=0;i<6;i++) {
|
||||
h_rate[i]=(h[i]-h_old[i])/update->dt;
|
||||
h_old[i]=h[i];
|
||||
}
|
||||
|
||||
Hdot(1,1)=h_rate[0]; Hdot(1,2)=0.0; Hdot(1,3)=0.0;
|
||||
Hdot(2,1)=0.0; Hdot(2,2)=h_rate[1]; Hdot(2,3)=0.0;
|
||||
Hdot(3,1)=0.0; Hdot(3,2)=0.0; Hdot(3,3)=h_rate[2];
|
||||
|
||||
if (domain->triclinic) {
|
||||
H(1,2)=h[5]; H(1,3)=h[4]; H(2,3)=h[3];
|
||||
Hdot(1,2)=h_rate[5]; Hdot(1,3)=h_rate[4]; Hdot(2,3)=h_rate[3];
|
||||
}
|
||||
|
||||
double F[3][3]={0.0};
|
||||
double FT[3][3]={0.0};
|
||||
double Fi[3][3]={0.0};
|
||||
double Fdot[3][3]={0.0};
|
||||
double J,vol;
|
||||
|
||||
MathExtra::times3(H,invH0,F); //find F
|
||||
MathExtra::times3(Hdot,invH0,Fdot); //find Fdot
|
||||
|
||||
MathExtra::invert3(F,Fi);
|
||||
J = MathExtra::det3(F); //find J
|
||||
vol=myvol0*J; // actual volume
|
||||
|
||||
double deltat;
|
||||
deltat = update->dt; //increment of time step
|
||||
|
||||
// Current pressure on the cell boundaries:
|
||||
//tensor: 0 1 2 3 4 5
|
||||
// x y z xy xz yz
|
||||
|
||||
double *tensor = pressure->vector;
|
||||
|
||||
double cauchy[3][3]; //stress
|
||||
cauchy(1,1)=-tensor[0]; cauchy(1,2)=0.0; cauchy(1,3)=0.0;
|
||||
cauchy(2,1)=0.0; cauchy(2,2)=-tensor[1]; cauchy(2,3)=0.0;
|
||||
cauchy(3,1)=0.0; cauchy(3,2)=0.0; cauchy(3,3)=-tensor[2];
|
||||
|
||||
if (domain->triclinic) {
|
||||
cauchy(1,2)=-tensor[3]; cauchy(1,3)=-tensor[4];
|
||||
cauchy(2,1)=-tensor[3]; cauchy(2,3)=-tensor[5];
|
||||
cauchy(3,1)=-tensor[4]; cauchy(3,2)=-tensor[5];
|
||||
}
|
||||
|
||||
// target pressure on the cell boundaries:
|
||||
//p_target: 0 1 2 3 4 5
|
||||
// x y z yz xz xy
|
||||
|
||||
double setcauchy[3][3]; //stress
|
||||
setcauchy(1,1)=-p_target[0]; setcauchy(1,2)=0.0; setcauchy(1,3)=0.0;
|
||||
setcauchy(2,1)=0.0; setcauchy(2,2)=-p_target[1]; setcauchy(2,3)=0.0;
|
||||
setcauchy(3,1)=0.0; setcauchy(3,2)=0.0; setcauchy(3,3)=-p_target[2];
|
||||
|
||||
if (domain->triclinic) {
|
||||
setcauchy(1,2)=-p_target[5]; setcauchy(1,3)=-p_target[4];
|
||||
setcauchy(2,1)=-p_target[5]; setcauchy(2,3)=-p_target[3];
|
||||
setcauchy(3,1)=-p_target[4]; setcauchy(3,2)=-p_target[3];
|
||||
}
|
||||
|
||||
//initialize:
|
||||
if(initPK==1) {
|
||||
if(restartPK==1) {
|
||||
setPK(1,1)=setPKinit[0]; setPK(1,2)=setPKinit[1]; setPK(1,3)=setPKinit[2];
|
||||
setPK(2,1)=setPKinit[1]; setPK(2,2)=setPKinit[3]; setPK(2,3)=setPKinit[4];
|
||||
setPK(3,1)=setPKinit[2]; setPK(3,2)=setPKinit[4]; setPK(3,3)=setPKinit[5];
|
||||
}else {
|
||||
setPK(1,1)=cauchy(1,1); setPK(1,2)=cauchy(1,2); setPK(1,3)=cauchy(1,3);
|
||||
setPK(2,1)=cauchy(2,1); setPK(2,2)=cauchy(2,2); setPK(2,3)=cauchy(2,3);
|
||||
setPK(3,1)=cauchy(3,1); setPK(3,2)=cauchy(3,2); setPK(3,3)=cauchy(3,3);
|
||||
}
|
||||
|
||||
initPK=0;
|
||||
}
|
||||
|
||||
//cauchystat:
|
||||
bigint step = update->ntimestep;
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
CauchyStat(step,F,Fi,Fdot,cauchy,setcauchy,setPK,vol,myvol0,deltat,alpha);
|
||||
|
||||
// use currentPK as new target:
|
||||
//p_target: 0 1 2 3 4 5
|
||||
// x y z yz xz xy
|
||||
|
||||
p_target[0]=-setPK(1,1);
|
||||
p_target[1]=-setPK(2,2);
|
||||
p_target[2]=-setPK(3,3);
|
||||
if (pstyle == TRICLINIC) {
|
||||
p_target[3]=-setPK(2,3);
|
||||
p_target[4]=-setPK(1,3);
|
||||
p_target[5]=-setPK(1,2);
|
||||
}
|
||||
|
||||
p_hydro = 0.0;
|
||||
for (int i = 0; i < 3; i++)
|
||||
if (p_flag[i]) {
|
||||
p_hydro += p_target[i];
|
||||
}
|
||||
p_hydro /= pdim;
|
||||
|
||||
// save information for Cauchystat restart
|
||||
setPKinit[0] = setcauchy(1,1);
|
||||
setPKinit[1] = setcauchy(1,2);
|
||||
setPKinit[2] = setcauchy(1,3);
|
||||
setPKinit[3] = setcauchy(2,2);
|
||||
setPKinit[4] = setcauchy(2,3);
|
||||
setPKinit[5] = setcauchy(3,3);
|
||||
restart_stored=1;
|
||||
|
||||
#undef H
|
||||
#undef Hdot
|
||||
#undef F
|
||||
#undef Fi
|
||||
#undef Fdot
|
||||
#undef invH0
|
||||
#undef cauchy
|
||||
#undef setcauchy
|
||||
#undef setPK
|
||||
#undef sigmahat
|
||||
|
||||
}
|
||||
|
||||
if((usePK==0) && (initRUN==1)) CauchyStat();
|
||||
if(initRUN==0){
|
||||
double* h = domain->h; // shape matrix in Voigt notation
|
||||
for (int i=0;i<6;i++) {
|
||||
h_old[i]=h[i];
|
||||
h_old[i]=domain->h[i];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
initRUN=1; // when run is initialized tensor[] not available (pressure on cell wall)
|
||||
|
||||
|
||||
// if deviatoric, recompute sigma each time p_target changes
|
||||
|
||||
if (deviatoric_flag) compute_sigma();
|
||||
}
|
||||
|
||||
|
|
@ -2622,7 +2409,222 @@ double FixNH::memory_usage()
|
|||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
CauchyStat
|
||||
initialize Cauchystat
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixNH::CauchyStat_init()
|
||||
{
|
||||
if (comm->me == 0) {
|
||||
if (screen) {
|
||||
fprintf(screen,"Using the Cauchystat fix with alpha=%f\n",alpha);
|
||||
if(restartPK==1) {
|
||||
fprintf(screen," (this is a continuation fix)\n");
|
||||
}
|
||||
else {
|
||||
fprintf(screen," (this is NOT a continuation fix)\n");
|
||||
}
|
||||
}
|
||||
if (logfile) {
|
||||
fprintf(logfile,"Using the Cauchystat with alpha=%f\n",alpha);
|
||||
if(restartPK==1) {
|
||||
fprintf(logfile," this is a continuation run\n");
|
||||
}
|
||||
else {
|
||||
fprintf(logfile," this is NOT a continuation run\n");
|
||||
}
|
||||
}
|
||||
}
|
||||
if(restartPK==1 && restart_stored==0) {
|
||||
error->all(FLERR,"Illegal cauchystat command. Continuation run"
|
||||
" must follow a previously equilibrated Cauchystat run");
|
||||
}
|
||||
if(alpha<=0.0) {
|
||||
error->all(FLERR,"Illegal cauchystat command. Alpha cannot be zero or negative");
|
||||
}
|
||||
|
||||
initRUN = 0;
|
||||
initPK = 1;
|
||||
usePK = 0;
|
||||
|
||||
#define H0(row,col) (H0[(row-1)][(col-1)])
|
||||
#define invH0(row,col) (invH0[(row-1)][(col-1)])
|
||||
|
||||
// initialize H0 to current box shape
|
||||
|
||||
int triclinic = domain->triclinic;
|
||||
double *h = domain->h;
|
||||
double *h_inv = domain->h_inv;
|
||||
double *boxhi = domain->boxhi;
|
||||
double *boxlo = domain->boxlo;
|
||||
double yz = domain->yz;
|
||||
double xz = domain->xz;
|
||||
double xy = domain->xy;
|
||||
|
||||
H0(1,1)=h[0]; H0(1,2)=h[5]; H0(1,3)=h[4];
|
||||
H0(2,1)=0.0; H0(2,2)=h[1]; H0(2,3)=h[3];
|
||||
H0(3,1)=0.0; H0(3,2)=0.0; H0(3,3)=h[2];
|
||||
|
||||
invH0(1,1)=h_inv[0]; invH0(1,2)=h_inv[5]; invH0(1,3)=h_inv[4];
|
||||
invH0(2,1)=0.0; invH0(2,2)=h_inv[1]; invH0(2,3)=h_inv[3];
|
||||
invH0(3,1)=0.0; invH0(3,2)=0.0; invH0(3,3)=h_inv[2];
|
||||
|
||||
CSvol0 = abs(MathExtra::det3(H0)); //find reference volume
|
||||
|
||||
#undef H0
|
||||
#undef invH0
|
||||
}
|
||||
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Cauchystat
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixNH::CauchyStat()
|
||||
{
|
||||
double* h = domain->h; // shape matrix in Voigt notation
|
||||
double* h_rate = domain->h_rate; // rate of box size/shape change in Voigt notation
|
||||
double H[3][3];
|
||||
double Hdot[3][3];
|
||||
|
||||
#define H(row,col) (H[(row-1)][(col-1)])
|
||||
#define Hdot(row,col) (Hdot[(row-1)][(col-1)])
|
||||
#define F(row,col) (F[(row-1)][(col-1)])
|
||||
#define Fi(row,col) (Fi[(row-1)][(col-1)])
|
||||
#define Fdot(row,col) (Fdot[(row-1)][(col-1)])
|
||||
#define invH0(row,col) (invH0[(row-1)][(col-1)])
|
||||
#define cauchy(row,col) (cauchy[(row-1)][(col-1)])
|
||||
#define setcauchy(row,col) (setcauchy[(row-1)][(col-1)])
|
||||
#define setPK(row,col) (setPK[(row-1)][(col-1)])
|
||||
#define sigmahat(row,col) (sigmahat[(row-1)][(col-1)])
|
||||
|
||||
H(1,1)=h[0]; H(1,2)=0.0; H(1,3)=0.0;
|
||||
H(2,1)=0.0; H(2,2)=h[1]; H(2,3)=0.0;
|
||||
H(3,1)=0.0; H(3,2)=0.0; H(3,3)=h[2];
|
||||
|
||||
for (int i=0;i<6;i++) {
|
||||
h_rate[i]=(h[i]-h_old[i])/update->dt;
|
||||
h_old[i]=h[i];
|
||||
}
|
||||
|
||||
Hdot(1,1)=h_rate[0]; Hdot(1,2)=0.0; Hdot(1,3)=0.0;
|
||||
Hdot(2,1)=0.0; Hdot(2,2)=h_rate[1]; Hdot(2,3)=0.0;
|
||||
Hdot(3,1)=0.0; Hdot(3,2)=0.0; Hdot(3,3)=h_rate[2];
|
||||
|
||||
if (domain->triclinic) {
|
||||
H(1,2)=h[5]; H(1,3)=h[4]; H(2,3)=h[3];
|
||||
Hdot(1,2)=h_rate[5]; Hdot(1,3)=h_rate[4]; Hdot(2,3)=h_rate[3];
|
||||
}
|
||||
|
||||
double F[3][3]={0.0};
|
||||
double FT[3][3]={0.0};
|
||||
double Fi[3][3]={0.0};
|
||||
double Fdot[3][3]={0.0};
|
||||
double J,vol;
|
||||
|
||||
MathExtra::times3(H,invH0,F);
|
||||
MathExtra::times3(Hdot,invH0,Fdot);
|
||||
MathExtra::invert3(F,Fi);
|
||||
J = MathExtra::det3(F);
|
||||
vol=CSvol0*J;
|
||||
|
||||
double deltat;
|
||||
deltat = update->dt; //increment of time step
|
||||
|
||||
// Current pressure on the cell boundaries:
|
||||
//tensor: 0 1 2 3 4 5
|
||||
// x y z xy xz yz
|
||||
|
||||
double *tensor = pressure->vector;
|
||||
|
||||
double cauchy[3][3];
|
||||
cauchy(1,1)=-tensor[0]; cauchy(1,2)=0.0; cauchy(1,3)=0.0;
|
||||
cauchy(2,1)=0.0; cauchy(2,2)=-tensor[1]; cauchy(2,3)=0.0;
|
||||
cauchy(3,1)=0.0; cauchy(3,2)=0.0; cauchy(3,3)=-tensor[2];
|
||||
|
||||
if (domain->triclinic) {
|
||||
cauchy(1,2)=-tensor[3]; cauchy(1,3)=-tensor[4];
|
||||
cauchy(2,1)=-tensor[3]; cauchy(2,3)=-tensor[5];
|
||||
cauchy(3,1)=-tensor[4]; cauchy(3,2)=-tensor[5];
|
||||
}
|
||||
|
||||
// target pressure on the cell boundaries:
|
||||
// p_target: 0 1 2 3 4 5
|
||||
// x y z yz xz xy
|
||||
|
||||
double setcauchy[3][3];
|
||||
setcauchy(1,1)=-p_target[0]; setcauchy(1,2)=0.0; setcauchy(1,3)=0.0;
|
||||
setcauchy(2,1)=0.0; setcauchy(2,2)=-p_target[1]; setcauchy(2,3)=0.0;
|
||||
setcauchy(3,1)=0.0; setcauchy(3,2)=0.0; setcauchy(3,3)=-p_target[2];
|
||||
|
||||
if (domain->triclinic) {
|
||||
setcauchy(1,2)=-p_target[5]; setcauchy(1,3)=-p_target[4];
|
||||
setcauchy(2,1)=-p_target[5]; setcauchy(2,3)=-p_target[3];
|
||||
setcauchy(3,1)=-p_target[4]; setcauchy(3,2)=-p_target[3];
|
||||
}
|
||||
|
||||
//initialize:
|
||||
if(initPK==1) {
|
||||
if(restartPK==1) {
|
||||
setPK(1,1)=setPKinit[0]; setPK(1,2)=setPKinit[1]; setPK(1,3)=setPKinit[2];
|
||||
setPK(2,1)=setPKinit[1]; setPK(2,2)=setPKinit[3]; setPK(2,3)=setPKinit[4];
|
||||
setPK(3,1)=setPKinit[2]; setPK(3,2)=setPKinit[4]; setPK(3,3)=setPKinit[5];
|
||||
}else {
|
||||
setPK(1,1)=cauchy(1,1); setPK(1,2)=cauchy(1,2); setPK(1,3)=cauchy(1,3);
|
||||
setPK(2,1)=cauchy(2,1); setPK(2,2)=cauchy(2,2); setPK(2,3)=cauchy(2,3);
|
||||
setPK(3,1)=cauchy(3,1); setPK(3,2)=cauchy(3,2); setPK(3,3)=cauchy(3,3);
|
||||
}
|
||||
initPK=0;
|
||||
}
|
||||
|
||||
CauchyStat_Step(update->ntimestep,F,Fi,Fdot,cauchy,setcauchy,setPK,vol,CSvol0,deltat,alpha);
|
||||
|
||||
// use currentPK as new target:
|
||||
//p_target: 0 1 2 3 4 5
|
||||
// x y z yz xz xy
|
||||
|
||||
p_target[0]=-setPK(1,1);
|
||||
p_target[1]=-setPK(2,2);
|
||||
p_target[2]=-setPK(3,3);
|
||||
if (pstyle == TRICLINIC) {
|
||||
p_target[3]=-setPK(2,3);
|
||||
p_target[4]=-setPK(1,3);
|
||||
p_target[5]=-setPK(1,2);
|
||||
}
|
||||
|
||||
p_hydro = 0.0;
|
||||
for (int i = 0; i < 3; i++)
|
||||
if (p_flag[i]) {
|
||||
p_hydro += p_target[i];
|
||||
}
|
||||
p_hydro /= pdim;
|
||||
|
||||
// save information for Cauchystat restart
|
||||
setPKinit[0] = setcauchy(1,1);
|
||||
setPKinit[1] = setcauchy(1,2);
|
||||
setPKinit[2] = setcauchy(1,3);
|
||||
setPKinit[3] = setcauchy(2,2);
|
||||
setPKinit[4] = setcauchy(2,3);
|
||||
setPKinit[5] = setcauchy(3,3);
|
||||
restart_stored=1;
|
||||
|
||||
#undef H
|
||||
#undef Hdot
|
||||
#undef F
|
||||
#undef Fi
|
||||
#undef Fdot
|
||||
#undef invH0
|
||||
#undef cauchy
|
||||
#undef setcauchy
|
||||
#undef setPK
|
||||
#undef sigmahat
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
CauchyStat_Step
|
||||
|
||||
Inputs:
|
||||
step : current step number
|
||||
|
|
@ -2632,7 +2634,6 @@ double FixNH::memory_usage()
|
|||
cauchy(3,3) : current cauchy stress tensor
|
||||
setcauchy(3,3) : target cauchy stress tensor
|
||||
alpha : parameter =0.01
|
||||
nstat : =1, flag for mod(step,nstat).ne.0
|
||||
setPK(3,3) : current PK stress tensor, at initialization (time=0) setPK=cauchy
|
||||
volume : current volume of the periodic box
|
||||
volume0 : initial volume of the periodic box
|
||||
|
|
@ -2641,28 +2642,26 @@ double FixNH::memory_usage()
|
|||
|
||||
Outputs:
|
||||
setPK(3,3) : PK stress tensor at the next time step
|
||||
------------------------------------------------------------------------- */
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixNH::CauchyStat(bigint step, double (&F)[3][3], double (&Fi)[3][3], double (&Fdot)[3][3], double (&cauchy)[3][3], double (&setcauchy)[3][3], double (&setPK)[3][3], double volume, double volume0, double deltat, double alpha)
|
||||
void FixNH::CauchyStat_Step(bigint step, double (&F)[3][3], double (&Fi)[3][3], double (&Fdot)[3][3], double (&cauchy)[3][3], double (&setcauchy)[3][3], double (&setPK)[3][3], double volume, double volume0, double deltat, double alpha)
|
||||
{
|
||||
|
||||
int nstat=1;
|
||||
|
||||
//macros to go from c to fortran style for arrays:
|
||||
#define F(row,col) (F[(row-1)][(col-1)])
|
||||
#define Fi(row,col) (Fi[(row-1)][(col-1)])
|
||||
#define Fdot(row,col) (Fdot[(row-1)][(col-1)])
|
||||
#define cauchy(row,col) (cauchy[(row-1)][(col-1)])
|
||||
#define setcauchy(row,col) (setcauchy[(row-1)][(col-1)])
|
||||
#define setPK(row,col) (setPK[(row-1)][(col-1)])
|
||||
|
||||
#define uv(row,col) (uv[(row-1)][(col-1)])
|
||||
#define deltastress(row) (deltastress[(row-1)])
|
||||
#define fdotvec(row) (fdotvec[(row-1)])
|
||||
#define dsdf(row,col) (dsdf[(row-1)][(col-1)])
|
||||
#define dsds(row,col) (dsds[(row-1)][(col-1)])
|
||||
#define deltaF(row) (deltaF[(row-1)])
|
||||
#define deltaPK(row) (deltaPK[(row-1)])
|
||||
#define F(row,col) (F[(row-1)][(col-1)])
|
||||
#define Fi(row,col) (Fi[(row-1)][(col-1)])
|
||||
#define Fdot(row,col) (Fdot[(row-1)][(col-1)])
|
||||
#define cauchy(row,col) (cauchy[(row-1)][(col-1)])
|
||||
#define setcauchy(row,col) (setcauchy[(row-1)][(col-1)])
|
||||
#define setPK(row,col) (setPK[(row-1)][(col-1)])
|
||||
#define uv(row,col) (uv[(row-1)][(col-1)])
|
||||
#define deltastress(row) (deltastress[(row-1)])
|
||||
#define fdotvec(row) (fdotvec[(row-1)])
|
||||
#define dsdf(row,col) (dsdf[(row-1)][(col-1)])
|
||||
#define dsds(row,col) (dsds[(row-1)][(col-1)])
|
||||
#define deltaF(row) (deltaF[(row-1)])
|
||||
#define deltaPK(row) (deltaPK[(row-1)])
|
||||
|
||||
//initialize local variables:
|
||||
int i,j,m,n;
|
||||
|
|
@ -2682,80 +2681,73 @@ void FixNH::CauchyStat(bigint step, double (&F)[3][3], double (&Fi)[3][3], doubl
|
|||
uv(5,1)=1; uv(5,2)=3;
|
||||
uv(6,1)=1; uv(6,2)=2;
|
||||
|
||||
if( (step % nstat) != 0) {
|
||||
//do nothing!
|
||||
}else {
|
||||
for(int ii = 1;ii <= 6;ii++) {
|
||||
i=uv(ii,1);
|
||||
j=uv(ii,2);
|
||||
deltastress(ii)=setcauchy(i,j)-cauchy(i,j);
|
||||
if(ii>3) deltastress(ii)=deltastress(ii)*2.0;
|
||||
fdotvec(ii)=Fdot(i,j)*deltat;
|
||||
}
|
||||
|
||||
for(int ii = 1;ii <= 6;ii++) {
|
||||
i=uv(ii,1);
|
||||
j=uv(ii,2);
|
||||
for(int jj = 1;jj <= 6;jj++) {
|
||||
m=uv(jj,1);
|
||||
n=uv(jj,2);
|
||||
dsds(ii,jj) = Fi(i,m)*Fi(j,n) + Fi(i,n)*Fi(j,m) + Fi(j,m)*Fi(i,n) + Fi(j,n)*Fi(i,m);
|
||||
for(int l = 1;l <= 3;l++) {
|
||||
for(int k = 1;k <= 3;k++) {
|
||||
dsdf(ii,jj) = dsdf(ii,jj) + cauchy(k,l)*( Fi(i,k)*Fi(j,l)*Fi(n,m) - Fi(i,m)*Fi(j,l)*Fi(n,k) - Fi(i,k)*Fi(j,m)*Fi(n,l) );
|
||||
}//k
|
||||
}//l
|
||||
}//jj
|
||||
}//ii
|
||||
|
||||
jac=volume/volume0;
|
||||
for(int ii = 1;ii <= 6;ii++) {
|
||||
for(int jj = 1;jj <= 6;jj++) {
|
||||
dsds(ii,jj)=dsds(ii,jj)*jac/4.0;
|
||||
dsdf(ii,jj)=dsdf(ii,jj)*jac;
|
||||
}//jj
|
||||
}//ii
|
||||
|
||||
for(int ii = 1;ii <= 6;ii++) {
|
||||
for(int jj = 1;jj <= 6;jj++) {
|
||||
deltaF(ii)=deltaF(ii)+dsdf(ii,jj)*fdotvec(jj); // deltaF=matmul(dsdf,fdotvec) in the fortran implementation
|
||||
}//jj
|
||||
}//ii
|
||||
|
||||
for(int ii = 1;ii <= 6;ii++) {
|
||||
for(int jj = 1;jj <= 6;jj++) {
|
||||
deltaPK(ii)=deltaPK(ii)+alpha*dsds(ii,jj)*deltastress(jj); // deltaPK=alpha*matmul(dsds,deltastress) + deltaF in the fortran implementation
|
||||
}//jj
|
||||
// this line applies alpha to both terms instead of just the delta sigma term.
|
||||
deltaPK(ii)=deltaPK(ii)+alpha*deltaF(ii);
|
||||
// this is the old version
|
||||
//deltaPK(ii)=deltaPK(ii)+deltaF(ii);
|
||||
}//ii
|
||||
|
||||
setPK(1,1)=setPK(1,1)+deltaPK(1); //equation (4) in SD-notes.pdf
|
||||
setPK(2,2)=setPK(2,2)+deltaPK(2);
|
||||
setPK(3,3)=setPK(3,3)+deltaPK(3);
|
||||
setPK(2,3)=setPK(2,3)+deltaPK(4);
|
||||
setPK(3,2)=setPK(3,2)+deltaPK(4);
|
||||
setPK(1,3)=setPK(1,3)+deltaPK(5);
|
||||
setPK(3,1)=setPK(3,1)+deltaPK(5);
|
||||
setPK(1,2)=setPK(1,2)+deltaPK(6);
|
||||
setPK(2,1)=setPK(2,1)+deltaPK(6);
|
||||
for(int ii = 1;ii <= 6;ii++) {
|
||||
i=uv(ii,1);
|
||||
j=uv(ii,2);
|
||||
deltastress(ii)=setcauchy(i,j)-cauchy(i,j);
|
||||
if(ii>3) deltastress(ii)=deltastress(ii)*2.0;
|
||||
fdotvec(ii)=Fdot(i,j)*deltat;
|
||||
}
|
||||
|
||||
//undefine macros:
|
||||
#undef F
|
||||
#undef Fi
|
||||
#undef Fdot
|
||||
#undef cauchy
|
||||
#undef setcauchy
|
||||
#undef setPK
|
||||
#undef uv
|
||||
#undef deltastress
|
||||
#undef fdotvec
|
||||
#undef dsdf
|
||||
#undef dsds
|
||||
#undef deltaF
|
||||
#undef deltaPK
|
||||
for(int ii = 1;ii <= 6;ii++) {
|
||||
i=uv(ii,1);
|
||||
j=uv(ii,2);
|
||||
for(int jj = 1;jj <= 6;jj++) {
|
||||
m=uv(jj,1);
|
||||
n=uv(jj,2);
|
||||
dsds(ii,jj) = Fi(i,m)*Fi(j,n) + Fi(i,n)*Fi(j,m) + Fi(j,m)*Fi(i,n) + Fi(j,n)*Fi(i,m);
|
||||
for(int l = 1;l <= 3;l++) {
|
||||
for(int k = 1;k <= 3;k++) {
|
||||
dsdf(ii,jj) = dsdf(ii,jj) + cauchy(k,l)*
|
||||
( Fi(i,k)*Fi(j,l)*Fi(n,m) - Fi(i,m)*Fi(j,l)*Fi(n,k) - Fi(i,k)*Fi(j,m)*Fi(n,l) );
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
jac=volume/volume0;
|
||||
for(int ii = 1;ii <= 6;ii++) {
|
||||
for(int jj = 1;jj <= 6;jj++) {
|
||||
dsds(ii,jj)=dsds(ii,jj)*jac/4.0;
|
||||
dsdf(ii,jj)=dsdf(ii,jj)*jac;
|
||||
}
|
||||
}
|
||||
|
||||
for(int ii = 1;ii <= 6;ii++) {
|
||||
for(int jj = 1;jj <= 6;jj++) {
|
||||
deltaF(ii)=deltaF(ii)+dsdf(ii,jj)*fdotvec(jj);
|
||||
}
|
||||
}
|
||||
|
||||
for(int ii = 1;ii <= 6;ii++) {
|
||||
for(int jj = 1;jj <= 6;jj++) {
|
||||
deltaPK(ii)=deltaPK(ii)+alpha*dsds(ii,jj)*deltastress(jj);
|
||||
}
|
||||
deltaPK(ii)=deltaPK(ii)+alpha*deltaF(ii);
|
||||
}
|
||||
|
||||
setPK(1,1)=setPK(1,1)+deltaPK(1);
|
||||
setPK(2,2)=setPK(2,2)+deltaPK(2);
|
||||
setPK(3,3)=setPK(3,3)+deltaPK(3);
|
||||
setPK(2,3)=setPK(2,3)+deltaPK(4);
|
||||
setPK(3,2)=setPK(3,2)+deltaPK(4);
|
||||
setPK(1,3)=setPK(1,3)+deltaPK(5);
|
||||
setPK(3,1)=setPK(3,1)+deltaPK(5);
|
||||
setPK(1,2)=setPK(1,2)+deltaPK(6);
|
||||
setPK(2,1)=setPK(2,1)+deltaPK(6);
|
||||
|
||||
#undef F
|
||||
#undef Fi
|
||||
#undef Fdot
|
||||
#undef cauchy
|
||||
#undef setcauchy
|
||||
#undef setPK
|
||||
#undef uv
|
||||
#undef deltastress
|
||||
#undef fdotvec
|
||||
#undef dsdf
|
||||
#undef dsds
|
||||
#undef deltaF
|
||||
#undef deltaPK
|
||||
|
||||
}
|
||||
|
|
|
|||
16
src/fix_nh.h
16
src/fix_nh.h
|
|
@ -144,17 +144,21 @@ class FixNH : public Fix {
|
|||
double H0[3][3]; //shape matrix for the undeformed cell
|
||||
double h_old[6]; //previous time step shape matrix for the undeformed cell
|
||||
double invH0[3][3]; //inverse of H0;
|
||||
double myvol0;
|
||||
double setPK[3][3];
|
||||
static double setPKinit[6];
|
||||
double alpha; //integration parameter cauchystat
|
||||
double CSvol0; //volume of undeformed cell
|
||||
double setPK[3][3]; //current set values of the PK stress (this is modified until the cauchy stress converges)
|
||||
static double setPKinit[6]; // initialization value of setPK for continuation runs
|
||||
double alpha; //integration parameter for the cauchystat
|
||||
int initPK; // 1 if setPK needs to be initialized either from cauchy or restart, else 0
|
||||
int usePK; // 0 if use CauchyStat else 1
|
||||
static int restartPK; // Read PK stress from the previous step
|
||||
static int restart_stored; // Read PK stress from the previous step
|
||||
static int restart_stored; // Read PK stress from the previous step
|
||||
int initRUN; // 0 if run not initialized (pressure->vector not computed yet), else 1 (pressure->vector available)
|
||||
|
||||
virtual void CauchyStat(bigint step, double (&F)[3][3], double (&Fi)[3][3], double (&Fdot)[3][3], double (&cauchy)[3][3], double (&setcauchy)[3][3], double (&setPK)[3][3], double volume, double volume0, double deltat, double alpha);
|
||||
virtual void CauchyStat_init();
|
||||
virtual void CauchyStat();
|
||||
virtual void CauchyStat_Step(bigint step, double (&F)[3][3], double (&Fi)[3][3], double (&Fdot)[3][3],
|
||||
double (&cauchy)[3][3], double (&setcauchy)[3][3], double (&setPK)[3][3],
|
||||
double volume, double volume0, double deltat, double alpha);
|
||||
|
||||
};
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue