Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80956158
pair_smd_ulsph.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Sep 3, 11:06
Size
46 KB
Mime Type
text/x-c
Expires
Thu, Sep 5, 11:06 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20486832
Attached To
rLAMMPS lammps
pair_smd_ulsph.cpp
View Options
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include <math.h>
#include <float.h>
#include <stdlib.h>
#include <string.h>
#include "pair_smd_ulsph.h"
#include "atom.h"
#include "domain.h"
#include "force.h"
#include "update.h"
#include "modify.h"
#include "fix.h"
#include "comm.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "memory.h"
#include "error.h"
#include <stdio.h>
#include <iostream>
#include "smd_material_models.h"
#include "smd_math.h"
#include "smd_kernels.h"
using
namespace
SMD_Kernels
;
using
namespace
std
;
using
namespace
LAMMPS_NS
;
using
namespace
SMD_Math
;
#include <Eigen/Eigen>
using
namespace
Eigen
;
#define ARTIFICIAL_STRESS false
#define FORMAT1 "%60s : %g\n"
#define FORMAT2 "\n.............................. %s \n"
PairULSPH
::
PairULSPH
(
LAMMPS
*
lmp
)
:
Pair
(
lmp
)
{
// per-type arrays
Q1
=
NULL
;
eos
=
viscosity
=
strength
=
NULL
;
c0_type
=
NULL
;
c0
=
NULL
;
Lookup
=
NULL
;
artificial_stress
=
NULL
;
artificial_pressure
=
NULL
;
nmax
=
0
;
// make sure no atom on this proc such that initial memory allocation is correct
stressTensor
=
L
=
K
=
NULL
;
shepardWeight
=
NULL
;
smoothVel
=
NULL
;
numNeighs
=
NULL
;
F
=
NULL
;
rho
=
NULL
;
effm
=
NULL
;
velocity_gradient_required
=
false
;
// turn off computation of velocity gradient by default
density_summation
=
velocity_gradient
=
false
;
comm_forward
=
18
;
// this pair style communicates 18 doubles to ghost atoms
updateFlag
=
0
;
}
/* ---------------------------------------------------------------------- */
PairULSPH
::~
PairULSPH
()
{
if
(
allocated
)
{
//printf("... deallocating\n");
memory
->
destroy
(
Q1
);
memory
->
destroy
(
rho0
);
memory
->
destroy
(
eos
);
memory
->
destroy
(
viscosity
);
memory
->
destroy
(
strength
);
memory
->
destroy
(
c0_type
);
memory
->
destroy
(
Lookup
);
memory
->
destroy
(
artificial_pressure
);
memory
->
destroy
(
artificial_stress
);
delete
[]
onerad_dynamic
;
delete
[]
onerad_frozen
;
delete
[]
maxrad_dynamic
;
delete
[]
maxrad_frozen
;
delete
[]
K
;
delete
[]
shepardWeight
;
delete
[]
c0
;
delete
[]
smoothVel
;
delete
[]
stressTensor
;
delete
[]
L
;
delete
[]
numNeighs
;
delete
[]
F
;
delete
[]
rho
;
delete
[]
effm
;
}
}
/* ----------------------------------------------------------------------
*
* Re-compute mass density from scratch.
* Only used for plain fluid SPH with no physical viscosity models.
*
---------------------------------------------------------------------- */
void
PairULSPH
::
PreCompute_DensitySummation
()
{
double
*
radius
=
atom
->
radius
;
double
**
x
=
atom
->
x
;
double
*
rmass
=
atom
->
rmass
;
int
*
type
=
atom
->
type
;
int
*
ilist
,
*
jlist
,
*
numneigh
;
int
**
firstneigh
;
int
nlocal
=
atom
->
nlocal
;
int
inum
,
jnum
,
ii
,
jj
,
i
,
itype
,
jtype
,
j
;
double
h
,
irad
,
hsq
,
rSq
,
wf
;
Vector3d
dx
,
xi
,
xj
;
// set up neighbor list variables
inum
=
list
->
inum
;
ilist
=
list
->
ilist
;
numneigh
=
list
->
numneigh
;
firstneigh
=
list
->
firstneigh
;
// zero accumulators
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
rho
[
i
]
=
0.0
;
//shepardWeight[i] = 0.0;
}
/*
* only recompute mass density if density summation is used.
* otherwise, change in mass density is time-integrated
*/
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
itype
=
type
[
i
];
if
(
setflag
[
itype
][
itype
]
==
1
)
{
// initialize particle density with self-contribution.
h
=
2.0
*
radius
[
i
];
hsq
=
h
*
h
;
Poly6Kernel
(
hsq
,
h
,
0.0
,
domain
->
dimension
,
wf
);
rho
[
i
]
=
wf
*
rmass
[
i
];
// / shepardWeight[i];
//printf("SIC to rho is %f\n", rho[i]);
}
}
for
(
ii
=
0
;
ii
<
inum
;
ii
++
)
{
i
=
ilist
[
ii
];
itype
=
type
[
i
];
jlist
=
firstneigh
[
i
];
jnum
=
numneigh
[
i
];
irad
=
radius
[
i
];
xi
<<
x
[
i
][
0
],
x
[
i
][
1
],
x
[
i
][
2
];
for
(
jj
=
0
;
jj
<
jnum
;
jj
++
)
{
j
=
jlist
[
jj
];
j
&=
NEIGHMASK
;
xj
<<
x
[
j
][
0
],
x
[
j
][
1
],
x
[
j
][
2
];
dx
=
xj
-
xi
;
rSq
=
dx
.
squaredNorm
();
h
=
irad
+
radius
[
j
];
hsq
=
h
*
h
;
if
(
rSq
<
hsq
)
{
jtype
=
type
[
j
];
Poly6Kernel
(
hsq
,
h
,
rSq
,
domain
->
dimension
,
wf
);
if
(
setflag
[
itype
][
itype
]
==
1
)
{
rho
[
i
]
+=
wf
*
rmass
[
j
];
// / shepardWeight[i];
}
if
(
j
<
nlocal
)
{
if
(
setflag
[
jtype
][
jtype
]
==
1
)
{
rho
[
j
]
+=
wf
*
rmass
[
i
];
// / shepardWeight[j];
}
}
}
// end if check distance
}
// end loop over j
}
// end loop over i
}
/* ----------------------------------------------------------------------
*
* Compute shape matrix for kernel gradient correction and velocity gradient.
* This is used if material strength or viscosity models are employed.
*
---------------------------------------------------------------------- */
void
PairULSPH
::
PreCompute
()
{
double
**
atom_data9
=
atom
->
smd_data_9
;
double
*
radius
=
atom
->
radius
;
double
**
x
=
atom
->
x
;
double
**
x0
=
atom
->
x0
;
double
**
v
=
atom
->
vest
;
double
*
vfrac
=
atom
->
vfrac
;
int
*
type
=
atom
->
type
;
int
*
ilist
,
*
jlist
,
*
numneigh
;
int
**
firstneigh
;
int
nlocal
=
atom
->
nlocal
;
int
inum
,
jnum
,
ii
,
jj
,
i
,
itype
,
j
,
idim
;
double
wfd
,
h
,
irad
,
r
,
rSq
,
wf
,
ivol
,
jvol
;
Vector3d
dx
,
dv
,
g
,
du
;
Matrix3d
Ktmp
,
Ltmp
,
Ftmp
,
K3di
,
D
;
Vector3d
xi
,
xj
,
vi
,
vj
,
x0i
,
x0j
,
dx0
;
Matrix2d
K2di
,
K2d
;
// zero accumulators
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
itype
=
type
[
i
];
if
(
setflag
[
itype
][
itype
])
{
if
(
gradient_correction_flag
)
{
K
[
i
].
setZero
();
}
else
{
K
[
i
].
setIdentity
();
}
L
[
i
].
setZero
();
F
[
i
].
setZero
();
}
}
// set up neighbor list variables
inum
=
list
->
inum
;
ilist
=
list
->
ilist
;
numneigh
=
list
->
numneigh
;
firstneigh
=
list
->
firstneigh
;
for
(
ii
=
0
;
ii
<
inum
;
ii
++
)
{
i
=
ilist
[
ii
];
itype
=
type
[
i
];
jlist
=
firstneigh
[
i
];
jnum
=
numneigh
[
i
];
irad
=
radius
[
i
];
ivol
=
vfrac
[
i
];
// initialize Eigen data structures from LAMMPS data structures
for
(
idim
=
0
;
idim
<
3
;
idim
++
)
{
x0i
(
idim
)
=
x0
[
i
][
idim
];
xi
(
idim
)
=
x
[
i
][
idim
];
vi
(
idim
)
=
v
[
i
][
idim
];
}
for
(
jj
=
0
;
jj
<
jnum
;
jj
++
)
{
j
=
jlist
[
jj
];
j
&=
NEIGHMASK
;
for
(
idim
=
0
;
idim
<
3
;
idim
++
)
{
x0j
(
idim
)
=
x0
[
j
][
idim
];
xj
(
idim
)
=
x
[
j
][
idim
];
vj
(
idim
)
=
v
[
j
][
idim
];
}
dx
=
xj
-
xi
;
rSq
=
dx
.
squaredNorm
();
h
=
irad
+
radius
[
j
];
if
(
rSq
<
h
*
h
)
{
r
=
sqrt
(
rSq
);
jvol
=
vfrac
[
j
];
// distance vectors in current and reference configuration, velocity difference
dv
=
vj
-
vi
;
dx0
=
x0j
-
x0i
;
// kernel and derivative
spiky_kernel_and_derivative
(
h
,
r
,
domain
->
dimension
,
wf
,
wfd
);
//barbara_kernel_and_derivative(h, r, domain->dimension, wf, wfd);
// uncorrected kernel gradient
g
=
(
wfd
/
r
)
*
dx
;
/* build correction matrix for kernel derivatives */
if
(
gradient_correction_flag
)
{
Ktmp
=
-
g
*
dx
.
transpose
();
K
[
i
]
+=
jvol
*
Ktmp
;
}
// velocity gradient L
Ltmp
=
-
dv
*
g
.
transpose
();
L
[
i
]
+=
jvol
*
Ltmp
;
// deformation gradient F in Eulerian frame
du
=
dx
-
dx0
;
Ftmp
=
dv
*
g
.
transpose
();
F
[
i
]
+=
jvol
*
Ftmp
;
if
(
j
<
nlocal
)
{
if
(
gradient_correction_flag
)
{
K
[
j
]
+=
ivol
*
Ktmp
;
}
L
[
j
]
+=
ivol
*
Ltmp
;
F
[
j
]
+=
ivol
*
Ftmp
;
}
}
// end if check distance
}
// end loop over j
}
// end loop over i
/*
* invert shape matrix and compute corrected quantities
*/
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
itype
=
type
[
i
];
if
(
setflag
[
itype
][
itype
])
{
if
(
gradient_correction_flag
)
{
pseudo_inverse_SVD
(
K
[
i
]);
K
[
i
]
=
LimitEigenvalues
(
K
[
i
],
2.0
);
L
[
i
]
*=
K
[
i
];
F
[
i
]
*=
K
[
i
];
}
// end if (gradient_correction[itype]) {
/*
* accumulate strain increments
* we abuse the atom array "atom_data_9" for this purpose, which was originally designed to hold the deformation gradient.
*/
D
=
update
->
dt
*
0.5
*
(
L
[
i
]
+
L
[
i
].
transpose
());
atom_data9
[
i
][
0
]
+=
D
(
0
,
0
);
// xx
atom_data9
[
i
][
1
]
+=
D
(
1
,
1
);
// yy
atom_data9
[
i
][
2
]
+=
D
(
2
,
2
);
// zz
atom_data9
[
i
][
3
]
+=
D
(
0
,
1
);
// xy
atom_data9
[
i
][
4
]
+=
D
(
0
,
2
);
// xz
atom_data9
[
i
][
5
]
+=
D
(
1
,
2
);
// yz
}
// end if (setflag[itype][itype])
}
// end loop over i = 0 to nlocal
}
/* ---------------------------------------------------------------------- */
void
PairULSPH
::
compute
(
int
eflag
,
int
vflag
)
{
double
**
x
=
atom
->
x
;
double
**
v
=
atom
->
vest
;
double
**
vint
=
atom
->
v
;
// Velocity-Verlet algorithm velocities
double
**
f
=
atom
->
f
;
double
*
vfrac
=
atom
->
vfrac
;
double
*
de
=
atom
->
de
;
double
*
rmass
=
atom
->
rmass
;
double
*
radius
=
atom
->
radius
;
double
*
contact_radius
=
atom
->
contact_radius
;
double
**
atom_data9
=
atom
->
smd_data_9
;
int
*
type
=
atom
->
type
;
int
nlocal
=
atom
->
nlocal
;
int
i
,
j
,
ii
,
jj
,
jnum
,
itype
,
jtype
,
iDim
,
inum
;
double
r
,
wf
,
wfd
,
h
,
rSq
,
ivol
,
jvol
;
double
mu_ij
,
c_ij
,
rho_ij
;
double
delVdotDelR
,
visc_magnitude
,
deltaE
;
int
*
ilist
,
*
jlist
,
*
numneigh
;
int
**
firstneigh
;
Vector3d
fi
,
fj
,
dx
,
dv
,
f_stress
,
g
,
vinti
,
vintj
,
dvint
;
Vector3d
xi
,
xj
,
vi
,
vj
,
f_visc
,
sumForces
,
f_stress_new
;
Vector3d
gamma
,
f_hg
,
dx0
,
du_est
,
du
;
double
r_ref
,
weight
,
p
;
//int periodic = (domain->xperiodic || domain->yperiodic || domain->zperiodic);
double
ini_dist
;
Matrix3d
S
,
D
,
V
,
eye
;
eye
.
setIdentity
();
int
k
;
SelfAdjointEigenSolver
<
Matrix3d
>
es
;
if
(
eflag
||
vflag
)
ev_setup
(
eflag
,
vflag
);
else
evflag
=
vflag_fdotr
=
0
;
if
(
atom
->
nmax
>
nmax
)
{
//printf("... allocating in compute with nmax = %d\n", atom->nmax);
nmax
=
atom
->
nmax
;
delete
[]
K
;
K
=
new
Matrix3d
[
nmax
];
delete
[]
shepardWeight
;
shepardWeight
=
new
double
[
nmax
];
delete
[]
c0
;
c0
=
new
double
[
nmax
];
delete
[]
smoothVel
;
smoothVel
=
new
Vector3d
[
nmax
];
delete
[]
stressTensor
;
stressTensor
=
new
Matrix3d
[
nmax
];
delete
[]
L
;
L
=
new
Matrix3d
[
nmax
];
delete
[]
numNeighs
;
numNeighs
=
new
int
[
nmax
];
delete
[]
F
;
F
=
new
Matrix3d
[
nmax
];
delete
[]
rho
;
rho
=
new
double
[
nmax
];
delete
[]
effm
;
effm
=
new
double
[
nmax
];
}
// zero accumulators
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
shepardWeight
[
i
]
=
0.0
;
smoothVel
[
i
].
setZero
();
numNeighs
[
i
]
=
0
;
h
=
2.0
*
radius
[
i
];
r
=
0.0
;
spiky_kernel_and_derivative
(
h
,
r
,
domain
->
dimension
,
wf
,
wfd
);
}
/*
* if this is the very first step, zero the array which holds the accumulated strain
*/
if
(
update
->
ntimestep
==
0
)
{
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
itype
=
type
[
i
];
if
(
setflag
[
itype
][
itype
])
{
for
(
j
=
0
;
j
<
9
;
j
++
)
{
atom_data9
[
i
][
j
]
=
0.0
;
}
}
}
}
if
(
density_summation
)
{
//printf("dens summ\n");
PreCompute_DensitySummation
();
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
//compute volumes from rho
itype
=
type
[
i
];
if
(
setflag
[
itype
][
itype
])
{
vfrac
[
i
]
=
rmass
[
i
]
/
rho
[
i
];
}
}
}
if
(
velocity_gradient
)
{
PairULSPH
::
PreCompute
();
// get velocity gradient and kernel gradient correction
}
PairULSPH
::
AssembleStressTensor
();
/*
* QUANTITIES ABOVE HAVE ONLY BEEN CALCULATED FOR NLOCAL PARTICLES.
* NEED TO DO A FORWARD COMMUNICATION TO GHOST ATOMS NOW
*/
comm
->
forward_comm_pair
(
this
);
updateFlag
=
0
;
/*
* iterate over pairs of particles i, j and assign forces using pre-computed pressure
*/
// set up neighbor list variables
inum
=
list
->
inum
;
ilist
=
list
->
ilist
;
numneigh
=
list
->
numneigh
;
firstneigh
=
list
->
firstneigh
;
for
(
ii
=
0
;
ii
<
inum
;
ii
++
)
{
i
=
ilist
[
ii
];
itype
=
type
[
i
];
jlist
=
firstneigh
[
i
];
jnum
=
numneigh
[
i
];
ivol
=
vfrac
[
i
];
// initialize Eigen data structures from LAMMPS data structures
for
(
iDim
=
0
;
iDim
<
3
;
iDim
++
)
{
xi
(
iDim
)
=
x
[
i
][
iDim
];
vi
(
iDim
)
=
v
[
i
][
iDim
];
vinti
(
iDim
)
=
vint
[
i
][
iDim
];
}
for
(
jj
=
0
;
jj
<
jnum
;
jj
++
)
{
j
=
jlist
[
jj
];
j
&=
NEIGHMASK
;
xj
(
0
)
=
x
[
j
][
0
];
xj
(
1
)
=
x
[
j
][
1
];
xj
(
2
)
=
x
[
j
][
2
];
dx
=
xj
-
xi
;
rSq
=
dx
.
squaredNorm
();
h
=
radius
[
i
]
+
radius
[
j
];
if
(
rSq
<
h
*
h
)
{
// initialize Eigen data structures from LAMMPS data structures
for
(
iDim
=
0
;
iDim
<
3
;
iDim
++
)
{
vj
(
iDim
)
=
v
[
j
][
iDim
];
vintj
(
iDim
)
=
vint
[
j
][
iDim
];
}
r
=
sqrt
(
rSq
);
jtype
=
type
[
j
];
jvol
=
vfrac
[
j
];
// distance vectors in current and reference configuration, velocity difference
dv
=
vj
-
vi
;
dvint
=
vintj
-
vinti
;
// kernel and derivative
spiky_kernel_and_derivative
(
h
,
r
,
domain
->
dimension
,
wf
,
wfd
);
//barbara_kernel_and_derivative(h, r, domain->dimension, wf, wfd);
// uncorrected kernel gradient
g
=
(
wfd
/
r
)
*
dx
;
delVdotDelR
=
dx
.
dot
(
dv
)
/
(
r
+
0.1
*
h
);
// project relative velocity onto unit particle distance vector [m/s]
S
=
stressTensor
[
i
]
+
stressTensor
[
j
];
if
(
artificial_pressure
[
itype
][
jtype
]
>
0.0
)
{
p
=
S
.
trace
();
if
(
p
>
0.0
)
{
// we are in tension
r_ref
=
contact_radius
[
i
]
+
contact_radius
[
j
];
weight
=
Kernel_Cubic_Spline
(
r
,
h
)
/
Kernel_Cubic_Spline
(
r_ref
,
h
);
weight
=
pow
(
weight
,
4.0
);
S
-=
artificial_pressure
[
itype
][
jtype
]
*
weight
*
p
*
eye
;
}
}
/*
* artificial stress to control tensile instability
* Only works if particles are uniformly spaced initially.
*/
if
(
artificial_stress
[
itype
][
jtype
]
>
0.0
)
{
ini_dist
=
contact_radius
[
i
]
+
contact_radius
[
j
];
weight
=
Kernel_Cubic_Spline
(
r
,
h
)
/
Kernel_Cubic_Spline
(
ini_dist
,
h
);
weight
=
pow
(
weight
,
4.0
);
es
.
compute
(
S
);
D
=
es
.
eigenvalues
().
asDiagonal
();
for
(
k
=
0
;
k
<
3
;
k
++
)
{
if
(
D
(
k
,
k
)
>
0.0
)
{
D
(
k
,
k
)
-=
weight
*
artificial_stress
[
itype
][
jtype
]
*
D
(
k
,
k
);
}
}
V
=
es
.
eigenvectors
();
S
=
V
*
D
*
V
.
inverse
();
}
// compute forces
f_stress
=
-
ivol
*
jvol
*
S
*
g
;
// DO NOT TOUCH SIGN
/*
* artificial viscosity -- alpha is dimensionless
* MonaghanBalsara form of the artificial viscosity
*/
c_ij
=
0.5
*
(
c0
[
i
]
+
c0
[
j
]);
LimitDoubleMagnitude
(
delVdotDelR
,
1.1
*
c_ij
);
mu_ij
=
h
*
delVdotDelR
/
(
r
+
0.1
*
h
);
// units: [m * m/s / m = m/s]
rho_ij
=
0.5
*
(
rmass
[
i
]
/
ivol
+
rmass
[
j
]
/
jvol
);
visc_magnitude
=
0.5
*
(
Q1
[
itype
]
+
Q1
[
jtype
])
*
c_ij
*
mu_ij
/
rho_ij
;
f_visc
=
-
rmass
[
i
]
*
rmass
[
j
]
*
visc_magnitude
*
g
;
if
((
Lookup
[
HOURGLASS_CONTROL_AMPLITUDE
][
itype
]
>
0.0
)
&&
(
Lookup
[
HOURGLASS_CONTROL_AMPLITUDE
][
jtype
]
>
0.0
))
{
f_hg
=
ComputeHourglassForce
(
i
,
itype
,
j
,
jtype
,
dv
,
dx
,
g
,
c_ij
,
mu_ij
,
rho_ij
);
}
else
{
f_hg
.
setZero
();
}
sumForces
=
f_stress
+
f_visc
+
f_hg
;
// energy rate -- project velocity onto force vector
deltaE
=
sumForces
.
dot
(
dv
);
// apply forces to pair of particles
f
[
i
][
0
]
+=
sumForces
(
0
);
f
[
i
][
1
]
+=
sumForces
(
1
);
f
[
i
][
2
]
+=
sumForces
(
2
);
de
[
i
]
+=
deltaE
;
// accumulate smooth velocities
shepardWeight
[
i
]
+=
jvol
*
wf
;
smoothVel
[
i
]
+=
jvol
*
wf
*
dvint
;
numNeighs
[
i
]
+=
1
;
if
(
j
<
nlocal
)
{
f
[
j
][
0
]
-=
sumForces
(
0
);
f
[
j
][
1
]
-=
sumForces
(
1
);
f
[
j
][
2
]
-=
sumForces
(
2
);
de
[
j
]
+=
deltaE
;
shepardWeight
[
j
]
+=
ivol
*
wf
;
smoothVel
[
j
]
-=
ivol
*
wf
*
dvint
;
numNeighs
[
j
]
+=
1
;
}
// tally atomistic stress tensor
if
(
evflag
)
{
ev_tally_xyz
(
i
,
j
,
nlocal
,
0
,
0.0
,
0.0
,
sumForces
(
0
),
sumForces
(
1
),
sumForces
(
2
),
dx
(
0
),
dx
(
1
),
dx
(
2
));
}
}
}
}
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
itype
=
type
[
i
];
if
(
setflag
[
itype
][
itype
]
==
1
)
{
if
(
shepardWeight
[
i
]
!=
0.0
)
{
smoothVel
[
i
]
/=
shepardWeight
[
i
];
}
else
{
smoothVel
[
i
].
setZero
();
}
}
// end check if particle is SPH-type
}
// end loop over i = 0 to nlocal
if
(
vflag_fdotr
)
virial_fdotr_compute
();
}
/* ----------------------------------------------------------------------
Assemble total stress tensor with pressure, material sterength, and
viscosity contributions.
------------------------------------------------------------------------- */
void
PairULSPH
::
AssembleStressTensor
()
{
double
*
radius
=
atom
->
radius
;
double
*
vfrac
=
atom
->
vfrac
;
double
*
rmass
=
atom
->
rmass
;
double
*
eff_plastic_strain
=
atom
->
eff_plastic_strain
;
double
**
tlsph_stress
=
atom
->
smd_stress
;
double
*
e
=
atom
->
e
;
int
*
type
=
atom
->
type
;
int
i
,
itype
;
int
nlocal
=
atom
->
nlocal
;
Matrix3d
D
,
Ddev
,
W
,
V
,
sigma_diag
;
Matrix3d
eye
,
stressRate
,
StressRateDevJaumann
;
Matrix3d
sigmaInitial_dev
,
d_dev
,
sigmaFinal_dev
,
stressRateDev
,
oldStressDeviator
,
newStressDeviator
;
double
plastic_strain_increment
,
yieldStress
;
double
dt
=
update
->
dt
;
double
vol
,
newPressure
;
double
G_eff
=
0.0
;
// effective shear modulus
double
K_eff
;
// effective bulk modulus
double
M
,
p_wave_speed
;
double
rho
,
effectiveViscosity
;
Matrix3d
deltaStressDev
;
dtCFL
=
1.0e22
;
eye
.
setIdentity
();
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
{
itype
=
type
[
i
];
if
(
setflag
[
itype
][
itype
]
==
1
)
{
newStressDeviator
.
setZero
();
newPressure
=
0.0
;
stressTensor
[
i
].
setZero
();
vol
=
vfrac
[
i
];
rho
=
rmass
[
i
]
/
vfrac
[
i
];
effectiveViscosity
=
0.0
;
K_eff
=
0.0
;
G_eff
=
0.0
;
//printf("rho = %f\n", rho);
switch
(
eos
[
itype
])
{
default
:
error
->
one
(
FLERR
,
"unknown EOS."
);
break
;
case
NONE:
c0
[
i
]
=
1.0
;
break
;
case
EOS_TAIT:
TaitEOS_density
(
Lookup
[
EOS_TAIT_EXPONENT
][
itype
],
Lookup
[
REFERENCE_SOUNDSPEED
][
itype
],
Lookup
[
REFERENCE_DENSITY
][
itype
],
rho
,
newPressure
,
c0
[
i
]);
//printf("new pressure =%f\n", newPressure);
break
;
case
EOS_PERFECT_GAS:
PerfectGasEOS
(
Lookup
[
EOS_PERFECT_GAS_GAMMA
][
itype
],
vol
,
rmass
[
i
],
e
[
i
],
newPressure
,
c0
[
i
]);
break
;
case
EOS_LINEAR:
newPressure
=
Lookup
[
BULK_MODULUS
][
itype
]
*
(
rho
/
Lookup
[
REFERENCE_DENSITY
][
itype
]
-
1.0
);
//printf("p=%f, rho0=%f, rho=%f\n", newPressure, Lookup[REFERENCE_DENSITY][itype], rho);
c0
[
i
]
=
Lookup
[
REFERENCE_SOUNDSPEED
][
itype
];
break
;
}
K_eff
=
c0
[
i
]
*
c0
[
i
]
*
rho
;
// effective bulk modulus
/*
* ******************************* STRENGTH MODELS ************************************************
*/
if
(
strength
[
itype
]
!=
NONE
)
{
/*
* initial stress state: given by the unrotateted Cauchy stress.
* Assemble Eigen 3d matrix from stored stress state
*/
oldStressDeviator
(
0
,
0
)
=
tlsph_stress
[
i
][
0
];
oldStressDeviator
(
0
,
1
)
=
tlsph_stress
[
i
][
1
];
oldStressDeviator
(
0
,
2
)
=
tlsph_stress
[
i
][
2
];
oldStressDeviator
(
1
,
1
)
=
tlsph_stress
[
i
][
3
];
oldStressDeviator
(
1
,
2
)
=
tlsph_stress
[
i
][
4
];
oldStressDeviator
(
2
,
2
)
=
tlsph_stress
[
i
][
5
];
oldStressDeviator
(
1
,
0
)
=
oldStressDeviator
(
0
,
1
);
oldStressDeviator
(
2
,
0
)
=
oldStressDeviator
(
0
,
2
);
oldStressDeviator
(
2
,
1
)
=
oldStressDeviator
(
1
,
2
);
D
=
0.5
*
(
L
[
i
]
+
L
[
i
].
transpose
());
W
=
0.5
*
(
L
[
i
]
-
L
[
i
].
transpose
());
// spin tensor:: need this for Jaumann rate
d_dev
=
Deviator
(
D
);
switch
(
strength
[
itype
])
{
default
:
error
->
one
(
FLERR
,
"unknown strength model."
);
break
;
case
STRENGTH_LINEAR:
// here in a version with pressure part
// stressRateDev = Lookup[BULK_MODULUS][itype] * d_iso * eye + 2.0 * Lookup[SHEAR_MODULUS][itype] * d_dev;
// c0[i] = Lookup[REFERENCE_SOUNDSPEED][itype];
// newPressure = 0.0;
// here only stress deviator
stressRateDev
=
2.0
*
Lookup
[
SHEAR_MODULUS
][
itype
]
*
d_dev
;
//cout << "stress rate deviator is " << endl << stressRateDev << endl;
break
;
case
STRENGTH_LINEAR_PLASTIC:
yieldStress
=
Lookup
[
YIELD_STRENGTH
][
itype
]
+
Lookup
[
HARDENING_PARAMETER
][
itype
]
*
eff_plastic_strain
[
i
];
LinearPlasticStrength
(
Lookup
[
SHEAR_MODULUS
][
itype
],
yieldStress
,
oldStressDeviator
,
d_dev
,
dt
,
newStressDeviator
,
stressRateDev
,
plastic_strain_increment
);
eff_plastic_strain
[
i
]
+=
plastic_strain_increment
;
break
;
}
//double m = effective_longitudinal_modulus(itype, dt, d_iso, p_rate, d_dev, stressRate_dev, damage);
StressRateDevJaumann
=
stressRateDev
-
W
*
oldStressDeviator
+
oldStressDeviator
*
W
;
newStressDeviator
=
oldStressDeviator
+
dt
*
StressRateDevJaumann
;
tlsph_stress
[
i
][
0
]
=
newStressDeviator
(
0
,
0
);
tlsph_stress
[
i
][
1
]
=
newStressDeviator
(
0
,
1
);
tlsph_stress
[
i
][
2
]
=
newStressDeviator
(
0
,
2
);
tlsph_stress
[
i
][
3
]
=
newStressDeviator
(
1
,
1
);
tlsph_stress
[
i
][
4
]
=
newStressDeviator
(
1
,
2
);
tlsph_stress
[
i
][
5
]
=
newStressDeviator
(
2
,
2
);
// estimate effective shear modulus for time step stability
deltaStressDev
=
oldStressDeviator
-
newStressDeviator
;
G_eff
=
effective_shear_modulus
(
d_dev
,
deltaStressDev
,
dt
,
itype
);
}
// end if (strength[itype] != NONE)
if
(
viscosity
[
itype
]
!=
NONE
)
{
D
=
0.5
*
(
L
[
i
]
+
L
[
i
].
transpose
());
d_dev
=
Deviator
(
D
);
switch
(
viscosity
[
itype
])
{
default
:
error
->
one
(
FLERR
,
"unknown viscosity model."
);
break
;
case
VISCOSITY_NEWTON:
effectiveViscosity
=
Lookup
[
VISCOSITY_MU
][
itype
];
// double shear_rate = 2.0
// * sqrt(d_dev(0, 1) * d_dev(0, 1) + d_dev(0, 2) * d_dev(0, 2) + d_dev(1, 2) * d_dev(1, 2)); // 3d
//cout << "shear rate: " << shear_rate << endl;
//effectiveViscosity = PA6_270C(shear_rate);
//if (effectiveViscosity > 178.062e-6) {
// printf("effective visc is %f\n", effectiveViscosity);
//}
newStressDeviator
=
2.0
*
effectiveViscosity
*
d_dev
;
// newton original
//cout << "this is Ddev " << endl << d_dev << endl << endl;
break
;
}
}
// end if (viscosity[itype] != NONE)
/*
* assemble stress Tensor from pressure and deviatoric parts
*/
stressTensor
[
i
]
=
-
newPressure
*
eye
+
newStressDeviator
;
/*
* stable timestep based on speed-of-sound
*/
M
=
K_eff
+
4.0
*
G_eff
/
3.0
;
p_wave_speed
=
sqrt
(
M
/
rho
);
effm
[
i
]
=
G_eff
;
dtCFL
=
MIN
(
2
*
radius
[
i
]
/
p_wave_speed
,
dtCFL
);
/*
* stable timestep based on viscosity
*/
if
(
viscosity
[
itype
]
!=
NONE
)
{
dtCFL
=
MIN
(
4
*
radius
[
i
]
*
radius
[
i
]
*
rho
/
effectiveViscosity
,
dtCFL
);
}
/*
* kernel gradient correction
*/
if
(
gradient_correction_flag
)
{
stressTensor
[
i
]
*=
K
[
i
];
}
}
// end if (setflag[itype][itype] == 1)
}
// end loop over nlocal
//printf("stable timestep = %g\n", 0.1 * hMin * MaxBulkVelocity);
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void
PairULSPH
::
allocate
()
{
allocated
=
1
;
int
n
=
atom
->
ntypes
;
memory
->
create
(
setflag
,
n
+
1
,
n
+
1
,
"pair:setflag"
);
memory
->
create
(
Q1
,
n
+
1
,
"pair:Q1"
);
memory
->
create
(
rho0
,
n
+
1
,
"pair:Q2"
);
memory
->
create
(
c0_type
,
n
+
1
,
"pair:c0_type"
);
memory
->
create
(
eos
,
n
+
1
,
"pair:eosmodel"
);
memory
->
create
(
viscosity
,
n
+
1
,
"pair:viscositymodel"
);
memory
->
create
(
strength
,
n
+
1
,
"pair:strengthmodel"
);
memory
->
create
(
Lookup
,
MAX_KEY_VALUE
,
n
+
1
,
"pair:LookupTable"
);
memory
->
create
(
artificial_pressure
,
n
+
1
,
n
+
1
,
"pair:artificial_pressure"
);
memory
->
create
(
artificial_stress
,
n
+
1
,
n
+
1
,
"pair:artificial_stress"
);
memory
->
create
(
cutsq
,
n
+
1
,
n
+
1
,
"pair:cutsq"
);
// always needs to be allocated, even with granular neighborlist
/*
* initialize arrays to default values
*/
for
(
int
i
=
1
;
i
<=
n
;
i
++
)
{
for
(
int
j
=
i
;
j
<=
n
;
j
++
)
{
artificial_pressure
[
i
][
j
]
=
0.0
;
artificial_stress
[
i
][
j
]
=
0.0
;
setflag
[
i
][
j
]
=
0
;
}
}
onerad_dynamic
=
new
double
[
n
+
1
];
onerad_frozen
=
new
double
[
n
+
1
];
maxrad_dynamic
=
new
double
[
n
+
1
];
maxrad_frozen
=
new
double
[
n
+
1
];
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void
PairULSPH
::
settings
(
int
narg
,
char
**
arg
)
{
if
(
narg
!=
3
)
{
printf
(
"narg = %d
\n
"
,
narg
);
error
->
all
(
FLERR
,
"Illegal number of arguments for pair_style ulsph"
);
}
if
(
comm
->
me
==
0
)
{
printf
(
"
\n
>>========>>========>>========>>========>>========>>========>>========>>========
\n
"
);
printf
(
"... SMD / ULSPH PROPERTIES
\n\n
"
);
}
if
(
strcmp
(
arg
[
0
],
"*DENSITY_SUMMATION"
)
==
0
)
{
density_summation
=
true
;
density_continuity
=
false
;
if
(
comm
->
me
==
0
)
printf
(
"... density summation active
\n
"
);
}
else
if
(
strcmp
(
arg
[
0
],
"*DENSITY_CONTINUITY"
)
==
0
)
{
density_continuity
=
true
;
density_summation
=
false
;
if
(
comm
->
me
==
0
)
printf
(
"... density continuity active
\n
"
);
}
else
{
error
->
all
(
FLERR
,
"Illegal settings keyword for first keyword of pair style ulsph. Must be either *DENSITY_SUMMATION or *DENSITY_CONTINUITY"
);
}
if
(
strcmp
(
arg
[
1
],
"*VELOCITY_GRADIENT"
)
==
0
)
{
velocity_gradient
=
true
;
if
(
comm
->
me
==
0
)
printf
(
"... computation of velocity gradients active
\n
"
);
}
else
if
(
strcmp
(
arg
[
1
],
"*NO_VELOCITY_GRADIENT"
)
==
0
)
{
velocity_gradient
=
false
;
if
(
comm
->
me
==
0
)
printf
(
"... computation of velocity gradients NOT active
\n
"
);
}
else
{
error
->
all
(
FLERR
,
"Illegal settings keyword for first keyword of pair style ulsph. Must be either *VELOCITY_GRADIENT or *NO_VELOCITY_GRADIENT"
);
}
if
(
strcmp
(
arg
[
2
],
"*GRADIENT_CORRECTION"
)
==
0
)
{
gradient_correction_flag
=
true
;
if
(
comm
->
me
==
0
)
printf
(
"... first order correction of kernel gradients is active
\n
"
);
}
else
if
(
strcmp
(
arg
[
2
],
"*NO_GRADIENT_CORRECTION"
)
==
0
)
{
gradient_correction_flag
=
false
;
if
(
comm
->
me
==
0
)
printf
(
"... first order correction of kernel gradients is NOT active
\n
"
);
}
else
{
error
->
all
(
FLERR
,
"Illegal settings keyword for pair style ulsph"
);
}
// error check
//if ((gradient_correction_flag == true) && (density_summation)) {
// error->all(FLERR, "Cannot use *DENSITY_SUMMATION in combination with *YES_GRADIENT_CORRECTION");
//}
if
(
comm
->
me
==
0
)
printf
(
">>========>>========>>========>>========>>========>>========>>========>>========
\n
"
);
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void
PairULSPH
::
coeff
(
int
narg
,
char
**
arg
)
{
int
ioffset
,
iarg
,
iNextKwd
,
itype
,
jtype
;
char
str
[
128
];
std
::
string
s
,
t
;
if
(
narg
<
3
)
{
sprintf
(
str
,
"number of arguments for pair ulsph is too small!"
);
error
->
all
(
FLERR
,
str
);
}
if
(
!
allocated
)
allocate
();
/*
* if parameters are give in i,i form, i.e., no a cross interaction, set material parameters
*/
if
(
force
->
inumeric
(
FLERR
,
arg
[
0
])
==
force
->
inumeric
(
FLERR
,
arg
[
1
]))
{
itype
=
force
->
inumeric
(
FLERR
,
arg
[
0
]);
eos
[
itype
]
=
viscosity
[
itype
]
=
strength
[
itype
]
=
NONE
;
if
(
comm
->
me
==
0
)
{
printf
(
"
\n
>>========>>========>>========>>========>>========>>========>>========>>========
\n
"
);
printf
(
"...SMD / ULSPH PROPERTIES OF PARTICLE TYPE %d
\n\n
"
,
itype
);
}
/*
* read parameters which are common -- regardless of material / eos model
*/
ioffset
=
2
;
if
(
strcmp
(
arg
[
ioffset
],
"*COMMON"
)
!=
0
)
{
sprintf
(
str
,
"common keyword missing!"
);
error
->
all
(
FLERR
,
str
);
}
else
{
}
t
=
string
(
"*"
);
iNextKwd
=
-
1
;
for
(
iarg
=
ioffset
+
1
;
iarg
<
narg
;
iarg
++
)
{
s
=
string
(
arg
[
iarg
]);
if
(
s
.
compare
(
0
,
t
.
length
(),
t
)
==
0
)
{
iNextKwd
=
iarg
;
break
;
}
}
//printf("keyword following *COMMON is %s\n", arg[iNextKwd]);
if
(
iNextKwd
<
0
)
{
sprintf
(
str
,
"no *KEYWORD terminates *COMMON"
);
error
->
all
(
FLERR
,
str
);
}
if
(
iNextKwd
-
ioffset
!=
5
+
1
)
{
sprintf
(
str
,
"expected 5 arguments following *COMMON but got %d
\n
"
,
iNextKwd
-
ioffset
-
1
);
error
->
all
(
FLERR
,
str
);
}
Lookup
[
REFERENCE_DENSITY
][
itype
]
=
force
->
numeric
(
FLERR
,
arg
[
ioffset
+
1
]);
Lookup
[
REFERENCE_SOUNDSPEED
][
itype
]
=
force
->
numeric
(
FLERR
,
arg
[
ioffset
+
2
]);
Q1
[
itype
]
=
force
->
numeric
(
FLERR
,
arg
[
ioffset
+
3
]);
Lookup
[
HEAT_CAPACITY
][
itype
]
=
force
->
numeric
(
FLERR
,
arg
[
ioffset
+
4
]);
Lookup
[
HOURGLASS_CONTROL_AMPLITUDE
][
itype
]
=
force
->
numeric
(
FLERR
,
arg
[
ioffset
+
5
]);
Lookup
[
BULK_MODULUS
][
itype
]
=
Lookup
[
REFERENCE_SOUNDSPEED
][
itype
]
*
Lookup
[
REFERENCE_SOUNDSPEED
][
itype
]
*
Lookup
[
REFERENCE_DENSITY
][
itype
];
if
(
comm
->
me
==
0
)
{
printf
(
"material unspecific properties for SMD/ULSPH definition of particle type %d:
\n
"
,
itype
);
printf
(
FORMAT1
,
"reference density"
,
Lookup
[
REFERENCE_DENSITY
][
itype
]);
printf
(
FORMAT1
,
"reference speed of sound"
,
Lookup
[
REFERENCE_SOUNDSPEED
][
itype
]);
printf
(
FORMAT1
,
"linear viscosity coefficient"
,
Q1
[
itype
]);
printf
(
FORMAT1
,
"heat capacity [energy / (mass * temperature)]"
,
Lookup
[
HEAT_CAPACITY
][
itype
]);
printf
(
FORMAT1
,
"bulk modulus"
,
Lookup
[
BULK_MODULUS
][
itype
]);
printf
(
FORMAT1
,
"hourglass control amplitude"
,
Lookup
[
HOURGLASS_CONTROL_AMPLITUDE
][
itype
]);
}
/*
* read following material cards
*/
// if (comm->me == 0) {
// printf("next kwd is %s\n", arg[iNextKwd]);
// }
while
(
true
)
{
if
(
strcmp
(
arg
[
iNextKwd
],
"*END"
)
==
0
)
{
// if (comm->me == 0) {
// sprintf(str, "found *END");
// error->message(FLERR, str);
// }
break
;
}
ioffset
=
iNextKwd
;
if
(
strcmp
(
arg
[
ioffset
],
"*EOS_TAIT"
)
==
0
)
{
/*
* Tait EOS
*/
eos
[
itype
]
=
EOS_TAIT
;
//printf("reading *EOS_TAIT\n");
t
=
string
(
"*"
);
iNextKwd
=
-
1
;
for
(
iarg
=
ioffset
+
1
;
iarg
<
narg
;
iarg
++
)
{
s
=
string
(
arg
[
iarg
]);
if
(
s
.
compare
(
0
,
t
.
length
(),
t
)
==
0
)
{
iNextKwd
=
iarg
;
break
;
}
}
if
(
iNextKwd
<
0
)
{
sprintf
(
str
,
"no *KEYWORD terminates *EOS_TAIT"
);
error
->
all
(
FLERR
,
str
);
}
if
(
iNextKwd
-
ioffset
!=
1
+
1
)
{
sprintf
(
str
,
"expected 1 arguments following *EOS_TAIT but got %d
\n
"
,
iNextKwd
-
ioffset
-
1
);
error
->
all
(
FLERR
,
str
);
}
Lookup
[
EOS_TAIT_EXPONENT
][
itype
]
=
force
->
numeric
(
FLERR
,
arg
[
ioffset
+
1
]);
if
(
comm
->
me
==
0
)
{
printf
(
FORMAT2
,
"Tait EOS"
);
printf
(
FORMAT1
,
"Exponent"
,
Lookup
[
EOS_TAIT_EXPONENT
][
itype
]);
}
}
// end Tait EOS
else
if
(
strcmp
(
arg
[
ioffset
],
"*EOS_PERFECT_GAS"
)
==
0
)
{
/*
* Perfect Gas EOS
*/
eos
[
itype
]
=
EOS_PERFECT_GAS
;
//printf("reading *EOS_PERFECT_GAS\n");
t
=
string
(
"*"
);
iNextKwd
=
-
1
;
for
(
iarg
=
ioffset
+
1
;
iarg
<
narg
;
iarg
++
)
{
s
=
string
(
arg
[
iarg
]);
if
(
s
.
compare
(
0
,
t
.
length
(),
t
)
==
0
)
{
iNextKwd
=
iarg
;
break
;
}
}
if
(
iNextKwd
<
0
)
{
sprintf
(
str
,
"no *KEYWORD terminates *EOS_PERFECT_GAS"
);
error
->
all
(
FLERR
,
str
);
}
if
(
iNextKwd
-
ioffset
!=
1
+
1
)
{
sprintf
(
str
,
"expected 1 arguments following *EOS_PERFECT_GAS but got %d
\n
"
,
iNextKwd
-
ioffset
-
1
);
error
->
all
(
FLERR
,
str
);
}
Lookup
[
EOS_PERFECT_GAS_GAMMA
][
itype
]
=
force
->
numeric
(
FLERR
,
arg
[
ioffset
+
1
]);
if
(
comm
->
me
==
0
)
{
printf
(
FORMAT2
,
"Perfect Gas EOS"
);
printf
(
FORMAT1
,
"Heat Capacity Ratio Gamma"
,
Lookup
[
EOS_PERFECT_GAS_GAMMA
][
itype
]);
}
}
// end Perfect Gas EOS
else
if
(
strcmp
(
arg
[
ioffset
],
"*EOS_LINEAR"
)
==
0
)
{
/*
* Linear EOS
*/
eos
[
itype
]
=
EOS_LINEAR
;
t
=
string
(
"*"
);
iNextKwd
=
-
1
;
for
(
iarg
=
ioffset
+
1
;
iarg
<
narg
;
iarg
++
)
{
s
=
string
(
arg
[
iarg
]);
if
(
s
.
compare
(
0
,
t
.
length
(),
t
)
==
0
)
{
iNextKwd
=
iarg
;
break
;
}
}
if
(
iNextKwd
<
0
)
{
sprintf
(
str
,
"no *KEYWORD terminates *EOS_LINEAR"
);
error
->
all
(
FLERR
,
str
);
}
if
(
iNextKwd
-
ioffset
!=
0
+
1
)
{
sprintf
(
str
,
"expected 0 arguments following *EOS_LINEAR but got %d
\n
"
,
iNextKwd
-
ioffset
-
1
);
error
->
all
(
FLERR
,
str
);
}
if
(
comm
->
me
==
0
)
{
printf
(
FORMAT2
,
"Linear EOS"
);
printf
(
FORMAT1
,
"Bulk modulus"
,
Lookup
[
BULK_MODULUS
][
itype
]);
}
}
// end Linear EOS
else
if
(
strcmp
(
arg
[
ioffset
],
"*STRENGTH_LINEAR_PLASTIC"
)
==
0
)
{
if
(
velocity_gradient
!=
true
)
{
error
->
all
(
FLERR
,
"A strength model was requested but *VELOCITY_GRADIENT is not set"
);
}
/*
* linear elastic / ideal plastic material model with strength
*/
strength
[
itype
]
=
STRENGTH_LINEAR_PLASTIC
;
velocity_gradient_required
=
true
;
//printf("reading *LINEAR_PLASTIC\n");
t
=
string
(
"*"
);
iNextKwd
=
-
1
;
for
(
iarg
=
ioffset
+
1
;
iarg
<
narg
;
iarg
++
)
{
s
=
string
(
arg
[
iarg
]);
if
(
s
.
compare
(
0
,
t
.
length
(),
t
)
==
0
)
{
iNextKwd
=
iarg
;
break
;
}
}
if
(
iNextKwd
<
0
)
{
sprintf
(
str
,
"no *KEYWORD terminates *STRENGTH_LINEAR_PLASTIC"
);
error
->
all
(
FLERR
,
str
);
}
if
(
iNextKwd
-
ioffset
!=
3
+
1
)
{
sprintf
(
str
,
"expected 3 arguments following *STRENGTH_LINEAR_PLASTIC but got %d
\n
"
,
iNextKwd
-
ioffset
-
1
);
error
->
all
(
FLERR
,
str
);
}
Lookup
[
SHEAR_MODULUS
][
itype
]
=
force
->
numeric
(
FLERR
,
arg
[
ioffset
+
1
]);
Lookup
[
YIELD_STRENGTH
][
itype
]
=
force
->
numeric
(
FLERR
,
arg
[
ioffset
+
2
]);
Lookup
[
HARDENING_PARAMETER
][
itype
]
=
force
->
numeric
(
FLERR
,
arg
[
ioffset
+
3
]);
if
(
comm
->
me
==
0
)
{
printf
(
FORMAT2
,
"linear elastic / ideal plastic material mode"
);
printf
(
FORMAT1
,
"yield_strength"
,
Lookup
[
YIELD_STRENGTH
][
itype
]);
printf
(
FORMAT1
,
"constant hardening parameter"
,
Lookup
[
HARDENING_PARAMETER
][
itype
]);
printf
(
FORMAT1
,
"shear modulus"
,
Lookup
[
SHEAR_MODULUS
][
itype
]);
}
}
// end *STRENGTH_LINEAR_PLASTIC
else
if
(
strcmp
(
arg
[
ioffset
],
"*STRENGTH_LINEAR"
)
==
0
)
{
if
(
velocity_gradient
!=
true
)
{
error
->
all
(
FLERR
,
"A strength model was requested but *VELOCITY_GRADIENT is not set"
);
}
/*
* linear elastic / ideal plastic material model with strength
*/
strength
[
itype
]
=
STRENGTH_LINEAR
;
t
=
string
(
"*"
);
iNextKwd
=
-
1
;
for
(
iarg
=
ioffset
+
1
;
iarg
<
narg
;
iarg
++
)
{
s
=
string
(
arg
[
iarg
]);
if
(
s
.
compare
(
0
,
t
.
length
(),
t
)
==
0
)
{
iNextKwd
=
iarg
;
break
;
}
}
if
(
iNextKwd
<
0
)
{
sprintf
(
str
,
"no *KEYWORD terminates *STRENGTH_LINEAR"
);
error
->
all
(
FLERR
,
str
);
}
if
(
iNextKwd
-
ioffset
!=
1
+
1
)
{
sprintf
(
str
,
"expected 1 arguments following *STRENGTH_LINEAR but got %d
\n
"
,
iNextKwd
-
ioffset
-
1
);
error
->
all
(
FLERR
,
str
);
}
Lookup
[
SHEAR_MODULUS
][
itype
]
=
force
->
numeric
(
FLERR
,
arg
[
ioffset
+
1
]);
if
(
comm
->
me
==
0
)
{
printf
(
FORMAT2
,
"linear elastic strength model"
);
printf
(
FORMAT1
,
"shear modulus"
,
Lookup
[
SHEAR_MODULUS
][
itype
]);
}
}
// end *STRENGTH_LINEAR
else
if
(
strcmp
(
arg
[
ioffset
],
"*VISCOSITY_NEWTON"
)
==
0
)
{
if
(
velocity_gradient
!=
true
)
{
error
->
all
(
FLERR
,
"A viscosity model was requested but *VELOCITY_GRADIENT is not set"
);
}
/*
* linear elastic / ideal plastic material model with strength
*/
viscosity
[
itype
]
=
VISCOSITY_NEWTON
;
t
=
string
(
"*"
);
iNextKwd
=
-
1
;
for
(
iarg
=
ioffset
+
1
;
iarg
<
narg
;
iarg
++
)
{
s
=
string
(
arg
[
iarg
]);
if
(
s
.
compare
(
0
,
t
.
length
(),
t
)
==
0
)
{
iNextKwd
=
iarg
;
break
;
}
}
if
(
iNextKwd
<
0
)
{
sprintf
(
str
,
"no *KEYWORD terminates *VISCOSITY_NEWTON"
);
error
->
all
(
FLERR
,
str
);
}
if
(
iNextKwd
-
ioffset
!=
1
+
1
)
{
sprintf
(
str
,
"expected 1 arguments following *VISCOSITY_NEWTON but got %d
\n
"
,
iNextKwd
-
ioffset
-
1
);
error
->
all
(
FLERR
,
str
);
}
Lookup
[
VISCOSITY_MU
][
itype
]
=
force
->
numeric
(
FLERR
,
arg
[
ioffset
+
1
]);
if
(
comm
->
me
==
0
)
{
printf
(
FORMAT2
,
"Newton viscosity model"
);
printf
(
FORMAT1
,
"viscosity mu"
,
Lookup
[
VISCOSITY_MU
][
itype
]);
}
}
// end *STRENGTH_VISCOSITY_NEWTON
else
if
(
strcmp
(
arg
[
ioffset
],
"*ARTIFICIAL_PRESSURE"
)
==
0
)
{
/*
* use Monaghan's artificial pressure to prevent particle clumping
*/
t
=
string
(
"*"
);
iNextKwd
=
-
1
;
for
(
iarg
=
ioffset
+
1
;
iarg
<
narg
;
iarg
++
)
{
s
=
string
(
arg
[
iarg
]);
if
(
s
.
compare
(
0
,
t
.
length
(),
t
)
==
0
)
{
iNextKwd
=
iarg
;
break
;
}
}
if
(
iNextKwd
<
0
)
{
sprintf
(
str
,
"no *KEYWORD terminates *ARTIFICIAL_PRESSURE"
);
error
->
all
(
FLERR
,
str
);
}
if
(
iNextKwd
-
ioffset
!=
1
+
1
)
{
sprintf
(
str
,
"expected 1 arguments following *ARTIFICIAL_PRESSURE but got %d
\n
"
,
iNextKwd
-
ioffset
-
1
);
error
->
all
(
FLERR
,
str
);
}
artificial_pressure
[
itype
][
itype
]
=
force
->
numeric
(
FLERR
,
arg
[
ioffset
+
1
]);
if
(
comm
->
me
==
0
)
{
printf
(
FORMAT2
,
"Artificial Pressure is enabled."
);
printf
(
FORMAT1
,
"Artificial Pressure amplitude"
,
artificial_pressure
[
itype
][
itype
]);
}
}
// end *ARTIFICIAL_PRESSURE
else
if
(
strcmp
(
arg
[
ioffset
],
"*ARTIFICIAL_STRESS"
)
==
0
)
{
/*
* use Monaghan's artificial stress to prevent particle clumping
*/
t
=
string
(
"*"
);
iNextKwd
=
-
1
;
for
(
iarg
=
ioffset
+
1
;
iarg
<
narg
;
iarg
++
)
{
s
=
string
(
arg
[
iarg
]);
if
(
s
.
compare
(
0
,
t
.
length
(),
t
)
==
0
)
{
iNextKwd
=
iarg
;
break
;
}
}
if
(
iNextKwd
<
0
)
{
sprintf
(
str
,
"no *KEYWORD terminates *ARTIFICIAL_STRESS"
);
error
->
all
(
FLERR
,
str
);
}
if
(
iNextKwd
-
ioffset
!=
1
+
1
)
{
sprintf
(
str
,
"expected 1 arguments following *ARTIFICIAL_STRESS but got %d
\n
"
,
iNextKwd
-
ioffset
-
1
);
error
->
all
(
FLERR
,
str
);
}
artificial_stress
[
itype
][
itype
]
=
force
->
numeric
(
FLERR
,
arg
[
ioffset
+
1
]);
if
(
comm
->
me
==
0
)
{
printf
(
FORMAT2
,
"Artificial Stress is enabled."
);
printf
(
FORMAT1
,
"Artificial Stress amplitude"
,
artificial_stress
[
itype
][
itype
]);
}
}
// end *ARTIFICIAL_STRESS
else
{
sprintf
(
str
,
"unknown *KEYWORD: %s"
,
arg
[
ioffset
]);
error
->
all
(
FLERR
,
str
);
}
}
/*
* copy data which is looked up in inner pairwise loops from slow maps to fast arrays
*/
rho0
[
itype
]
=
Lookup
[
REFERENCE_DENSITY
][
itype
];
c0_type
[
itype
]
=
Lookup
[
REFERENCE_SOUNDSPEED
][
itype
];
setflag
[
itype
][
itype
]
=
1
;
/*
* error checks
*/
if
((
viscosity
[
itype
]
!=
NONE
)
&&
(
strength
[
itype
]
!=
NONE
))
{
sprintf
(
str
,
"cannot have both a strength and viscosity model for particle type %d"
,
itype
);
error
->
all
(
FLERR
,
str
);
}
if
(
eos
[
itype
]
==
NONE
)
{
sprintf
(
str
,
"must specify an EOS for particle type %d"
,
itype
);
error
->
all
(
FLERR
,
str
);
}
}
else
{
/*
* we are reading a cross-interaction line for particle types i, j
*/
itype
=
force
->
inumeric
(
FLERR
,
arg
[
0
]);
jtype
=
force
->
inumeric
(
FLERR
,
arg
[
1
]);
if
(
strcmp
(
arg
[
2
],
"*CROSS"
)
!=
0
)
{
sprintf
(
str
,
"ulsph cross interaction between particle type %d and %d requested, however, *CROSS keyword is missing"
,
itype
,
jtype
);
error
->
all
(
FLERR
,
str
);
}
if
(
setflag
[
itype
][
itype
]
!=
1
)
{
sprintf
(
str
,
"ulsph cross interaction between particle type %d and %d requested, however, properties of type %d have not yet been specified"
,
itype
,
jtype
,
itype
);
error
->
all
(
FLERR
,
str
);
}
if
(
setflag
[
jtype
][
jtype
]
!=
1
)
{
sprintf
(
str
,
"ulsph cross interaction between particle type %d and %d requested, however, properties of type %d have not yet been specified"
,
itype
,
jtype
,
jtype
);
error
->
all
(
FLERR
,
str
);
}
setflag
[
itype
][
jtype
]
=
1
;
setflag
[
jtype
][
itype
]
=
1
;
if
((
artificial_pressure
[
itype
][
itype
]
>
0.0
)
&&
(
artificial_pressure
[
jtype
][
jtype
]
>
0.0
))
{
artificial_pressure
[
itype
][
jtype
]
=
0.5
*
(
artificial_pressure
[
itype
][
itype
]
+
artificial_pressure
[
jtype
][
jtype
]);
artificial_pressure
[
jtype
][
itype
]
=
artificial_pressure
[
itype
][
jtype
];
}
else
{
artificial_pressure
[
itype
][
jtype
]
=
artificial_pressure
[
jtype
][
itype
]
=
0.0
;
}
if
((
artificial_stress
[
itype
][
itype
]
>
0.0
)
&&
(
artificial_stress
[
jtype
][
jtype
]
>
0.0
))
{
artificial_stress
[
itype
][
jtype
]
=
0.5
*
(
artificial_stress
[
itype
][
itype
]
+
artificial_stress
[
jtype
][
jtype
]);
artificial_stress
[
jtype
][
itype
]
=
artificial_stress
[
itype
][
jtype
];
}
else
{
artificial_stress
[
itype
][
jtype
]
=
artificial_stress
[
jtype
][
itype
]
=
0.0
;
}
if
(
comm
->
me
==
0
)
{
printf
(
">>========>>========>>========>>========>>========>>========>>========>>========
\n
"
);
}
}
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double
PairULSPH
::
init_one
(
int
i
,
int
j
)
{
if
(
!
allocated
)
allocate
();
if
(
setflag
[
i
][
j
]
==
0
)
error
->
all
(
FLERR
,
"All pair coeffs are not set"
);
// cutoff = sum of max I,J radii for
// dynamic/dynamic & dynamic/frozen interactions, but not frozen/frozen
double
cutoff
=
maxrad_dynamic
[
i
]
+
maxrad_dynamic
[
j
];
cutoff
=
MAX
(
cutoff
,
maxrad_frozen
[
i
]
+
maxrad_dynamic
[
j
]);
cutoff
=
MAX
(
cutoff
,
maxrad_dynamic
[
i
]
+
maxrad_frozen
[
j
]);
//printf("cutoff for pair sph/fluid = %f\n", cutoff);
return
cutoff
;
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void
PairULSPH
::
init_style
()
{
int
i
;
//printf(" in init style\n");
// request a granular neighbor list
int
irequest
=
neighbor
->
request
(
this
);
neighbor
->
requests
[
irequest
]
->
half
=
0
;
neighbor
->
requests
[
irequest
]
->
gran
=
1
;
// set maxrad_dynamic and maxrad_frozen for each type
// include future Fix pour particles as dynamic
for
(
i
=
1
;
i
<=
atom
->
ntypes
;
i
++
)
onerad_dynamic
[
i
]
=
onerad_frozen
[
i
]
=
0.0
;
double
*
radius
=
atom
->
radius
;
int
*
type
=
atom
->
type
;
int
nlocal
=
atom
->
nlocal
;
for
(
i
=
0
;
i
<
nlocal
;
i
++
)
onerad_dynamic
[
type
[
i
]]
=
MAX
(
onerad_dynamic
[
type
[
i
]],
radius
[
i
]);
MPI_Allreduce
(
&
onerad_dynamic
[
1
],
&
maxrad_dynamic
[
1
],
atom
->
ntypes
,
MPI_DOUBLE
,
MPI_MAX
,
world
);
MPI_Allreduce
(
&
onerad_frozen
[
1
],
&
maxrad_frozen
[
1
],
atom
->
ntypes
,
MPI_DOUBLE
,
MPI_MAX
,
world
);
}
/* ----------------------------------------------------------------------
neighbor callback to inform pair style of neighbor list to use
optional granular history list
------------------------------------------------------------------------- */
void
PairULSPH
::
init_list
(
int
id
,
NeighList
*
ptr
)
{
if
(
id
==
0
)
list
=
ptr
;
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double
PairULSPH
::
memory_usage
()
{
//printf("in memory usage\n");
return
11
*
nmax
*
sizeof
(
double
);
}
/* ---------------------------------------------------------------------- */
int
PairULSPH
::
pack_forward_comm
(
int
n
,
int
*
list
,
double
*
buf
,
int
pbc_flag
,
int
*
pbc
)
{
double
*
vfrac
=
atom
->
vfrac
;
double
*
eff_plastic_strain
=
atom
->
eff_plastic_strain
;
int
i
,
j
,
m
;
//printf("packing comm\n");
m
=
0
;
for
(
i
=
0
;
i
<
n
;
i
++
)
{
j
=
list
[
i
];
buf
[
m
++
]
=
vfrac
[
j
];
buf
[
m
++
]
=
c0
[
j
];
//2
buf
[
m
++
]
=
stressTensor
[
j
](
0
,
0
);
// pack symmetric stress tensor
buf
[
m
++
]
=
stressTensor
[
j
](
1
,
1
);
buf
[
m
++
]
=
stressTensor
[
j
](
2
,
2
);
buf
[
m
++
]
=
stressTensor
[
j
](
0
,
1
);
buf
[
m
++
]
=
stressTensor
[
j
](
0
,
2
);
buf
[
m
++
]
=
stressTensor
[
j
](
1
,
2
);
// 2 + 6 = 8
buf
[
m
++
]
=
F
[
j
](
0
,
0
);
// F is not symmetric
buf
[
m
++
]
=
F
[
j
](
0
,
1
);
buf
[
m
++
]
=
F
[
j
](
0
,
2
);
buf
[
m
++
]
=
F
[
j
](
1
,
0
);
buf
[
m
++
]
=
F
[
j
](
1
,
1
);
buf
[
m
++
]
=
F
[
j
](
1
,
2
);
buf
[
m
++
]
=
F
[
j
](
2
,
0
);
buf
[
m
++
]
=
F
[
j
](
2
,
1
);
buf
[
m
++
]
=
F
[
j
](
2
,
2
);
// 8 + 9 = 17
buf
[
m
++
]
=
eff_plastic_strain
[
j
];
// 18
}
return
m
;
}
/* ---------------------------------------------------------------------- */
void
PairULSPH
::
unpack_forward_comm
(
int
n
,
int
first
,
double
*
buf
)
{
double
*
vfrac
=
atom
->
vfrac
;
double
*
eff_plastic_strain
=
atom
->
eff_plastic_strain
;
int
i
,
m
,
last
;
m
=
0
;
last
=
first
+
n
;
for
(
i
=
first
;
i
<
last
;
i
++
)
{
vfrac
[
i
]
=
buf
[
m
++
];
c0
[
i
]
=
buf
[
m
++
];
// 2
stressTensor
[
i
](
0
,
0
)
=
buf
[
m
++
];
stressTensor
[
i
](
1
,
1
)
=
buf
[
m
++
];
stressTensor
[
i
](
2
,
2
)
=
buf
[
m
++
];
stressTensor
[
i
](
0
,
1
)
=
buf
[
m
++
];
stressTensor
[
i
](
0
,
2
)
=
buf
[
m
++
];
stressTensor
[
i
](
1
,
2
)
=
buf
[
m
++
];
// 2 + 6 = 8
stressTensor
[
i
](
1
,
0
)
=
stressTensor
[
i
](
0
,
1
);
stressTensor
[
i
](
2
,
0
)
=
stressTensor
[
i
](
0
,
2
);
stressTensor
[
i
](
2
,
1
)
=
stressTensor
[
i
](
1
,
2
);
F
[
i
](
0
,
0
)
=
buf
[
m
++
];
F
[
i
](
0
,
1
)
=
buf
[
m
++
];
F
[
i
](
0
,
2
)
=
buf
[
m
++
];
F
[
i
](
1
,
0
)
=
buf
[
m
++
];
F
[
i
](
1
,
1
)
=
buf
[
m
++
];
F
[
i
](
1
,
2
)
=
buf
[
m
++
];
F
[
i
](
2
,
0
)
=
buf
[
m
++
];
F
[
i
](
2
,
1
)
=
buf
[
m
++
];
F
[
i
](
2
,
2
)
=
buf
[
m
++
];
// 8 + 9 = 17
eff_plastic_strain
[
i
]
=
buf
[
m
++
];
// 18
}
}
/*
* EXTRACT
*/
void
*
PairULSPH
::
extract
(
const
char
*
str
,
int
&
i
)
{
//printf("in extract\n");
if
(
strcmp
(
str
,
"smd/ulsph/smoothVel_ptr"
)
==
0
)
{
return
(
void
*
)
smoothVel
;
}
else
if
(
strcmp
(
str
,
"smd/ulsph/stressTensor_ptr"
)
==
0
)
{
return
(
void
*
)
stressTensor
;
}
else
if
(
strcmp
(
str
,
"smd/ulsph/velocityGradient_ptr"
)
==
0
)
{
return
(
void
*
)
L
;
}
else
if
(
strcmp
(
str
,
"smd/ulsph/numNeighs_ptr"
)
==
0
)
{
return
(
void
*
)
numNeighs
;
}
else
if
(
strcmp
(
str
,
"smd/ulsph/dtCFL_ptr"
)
==
0
)
{
//printf("dtcfl = %f\n", dtCFL);
return
(
void
*
)
&
dtCFL
;
}
else
if
(
strcmp
(
str
,
"smd/ulsph/updateFlag_ptr"
)
==
0
)
{
return
(
void
*
)
&
updateFlag
;
}
else
if
(
strcmp
(
str
,
"smd/ulsph/effective_modulus_ptr"
)
==
0
)
{
return
(
void
*
)
effm
;
}
else
if
(
strcmp
(
str
,
"smd/ulsph/shape_matrix_ptr"
)
==
0
)
{
return
(
void
*
)
K
;
}
return
NULL
;
}
/* ----------------------------------------------------------------------
compute effective shear modulus by dividing rate of deviatoric stress with rate of shear deformation
------------------------------------------------------------------------- */
double
PairULSPH
::
effective_shear_modulus
(
const
Matrix3d
d_dev
,
const
Matrix3d
deltaStressDev
,
const
double
dt
,
const
int
itype
)
{
double
G_eff
;
// effective shear modulus, see Pronto 2d eq. 3.4.7
double
deltaStressDevSum
,
shearRateSq
,
strain_increment
;
if
(
domain
->
dimension
==
3
)
{
deltaStressDevSum
=
deltaStressDev
(
0
,
1
)
*
deltaStressDev
(
0
,
1
)
+
deltaStressDev
(
0
,
2
)
*
deltaStressDev
(
0
,
2
)
+
deltaStressDev
(
1
,
2
)
*
deltaStressDev
(
1
,
2
);
shearRateSq
=
d_dev
(
0
,
1
)
*
d_dev
(
0
,
1
)
+
d_dev
(
0
,
2
)
*
d_dev
(
0
,
2
)
+
d_dev
(
1
,
2
)
*
d_dev
(
1
,
2
);
}
else
{
deltaStressDevSum
=
deltaStressDev
(
0
,
1
)
*
deltaStressDev
(
0
,
1
);
shearRateSq
=
d_dev
(
0
,
1
)
*
d_dev
(
0
,
1
);
}
strain_increment
=
dt
*
dt
*
shearRateSq
;
if
(
strain_increment
>
1.0e-12
)
{
G_eff
=
0.5
*
sqrt
(
deltaStressDevSum
/
strain_increment
);
}
else
{
if
(
strength
[
itype
]
!=
NONE
)
{
G_eff
=
Lookup
[
SHEAR_MODULUS
][
itype
];
}
else
{
G_eff
=
0.0
;
}
}
return
G_eff
;
}
/* ----------------------------------------------------------------------
hourglass force for updated Lagrangian SPH
input: particles indices i, j, particle types ityep, jtype
------------------------------------------------------------------------- */
Vector3d
PairULSPH
::
ComputeHourglassForce
(
const
int
i
,
const
int
itype
,
const
int
j
,
const
int
jtype
,
const
Vector3d
dv
,
const
Vector3d
xij
,
const
Vector3d
g
,
const
double
c_ij
,
const
double
mu_ij
,
const
double
rho_ij
)
{
double
*
rmass
=
atom
->
rmass
;
Vector3d
dv_est
,
f_hg
;
double
visc_magnitude
;
dv_est
=
-
0.5
*
(
F
[
i
]
+
F
[
j
])
*
xij
;
double
hurz
=
dv_est
.
dot
(
dv
)
/
(
dv_est
.
norm
()
*
dv
.
norm
()
+
1.0e-16
);
if
(
hurz
<
0.0
)
{
visc_magnitude
=
0.5
*
(
Q1
[
itype
]
+
Q1
[
jtype
])
*
c_ij
*
mu_ij
/
rho_ij
;
f_hg
=
-
rmass
[
i
]
*
rmass
[
j
]
*
visc_magnitude
*
g
;
// printf(" f_hg = %f %f %f\n", f_hg(0), f_hg(1), f_hg(2));
// printf("\nnegative\n");
// printf(" dv_est = %f %f %f\n", dv_est(0), dv_est(1), dv_est(2));
// printf(" dv = %f %f %f\n", dv(0), dv(1), dv(2));
}
else
{
f_hg
.
setZero
();
}
return
f_hg
;
}
Event Timeline
Log In to Comment