Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F69645239
fix_rx.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, Jul 2, 20:01
Size
52 KB
Mime Type
text/x-c
Expires
Thu, Jul 4, 20:01 (1 d, 21 h)
Engine
blob
Format
Raw Data
Handle
18492203
Attached To
rLAMMPS lammps
fix_rx.cpp
View Options
/* ----------------------------------------------------------------------
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 <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "fix_rx.h"
#include "atom.h"
#include "error.h"
#include "group.h"
#include "modify.h"
#include "force.h"
#include "memory.h"
#include "comm.h"
#include "update.h"
#include "domain.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "math_special.h"
#include "pair_dpd_fdt_energy.h"
#include <float.h>
// DBL_EPSILON
#include <vector>
// std::vector<>
#include <algorithm>
// std::max
#include <cmath>
// std::fmod
using
namespace
LAMMPS_NS
;
using
namespace
FixConst
;
using
namespace
MathSpecial
;
enum
{
NONE
,
HARMONIC
};
enum
{
LUCY
};
#define MAXLINE 1024
#define DELTA 4
#define SparseKinetics_enableIntegralReactions (true)
#define SparseKinetics_invalidIndex (-1)
namespace
/* anonymous */
{
typedef
double
TimerType
;
TimerType
getTimeStamp
(
void
)
{
return
MPI_Wtime
();
}
double
getElapsedTime
(
const
TimerType
&
t0
,
const
TimerType
&
t1
)
{
return
t1
-
t0
;
}
}
// end namespace
/* ---------------------------------------------------------------------- */
FixRX
::
FixRX
(
LAMMPS
*
lmp
,
int
narg
,
char
**
arg
)
:
Fix
(
lmp
,
narg
,
arg
)
{
if
(
narg
<
7
||
narg
>
12
)
error
->
all
(
FLERR
,
"Illegal fix rx command"
);
restart_peratom
=
1
;
nevery
=
1
;
nreactions
=
maxparam
=
0
;
params
=
NULL
;
mol2param
=
NULL
;
pairDPDE
=
NULL
;
id_fix_species
=
NULL
;
id_fix_species_old
=
NULL
;
const
int
Verbosity
=
1
;
// Keep track of the argument list.
int
iarg
=
3
;
// Read the kinetic file in arg[3].
kineticsFile
=
arg
[
iarg
++
];
// Determine the local temperature averaging method in arg[4].
wtFlag
=
0
;
localTempFlag
=
NONE
;
{
char
*
word
=
arg
[
iarg
++
];
if
(
strcmp
(
word
,
"none"
)
==
0
){
wtFlag
=
0
;
localTempFlag
=
NONE
;
}
else
if
(
strcmp
(
word
,
"lucy"
)
==
0
){
wtFlag
=
LUCY
;
localTempFlag
=
HARMONIC
;
}
else
error
->
all
(
FLERR
,
"Illegal fix rx local temperature weighting technique"
);
}
// Select either sparse and dense matrix
// representations of the stoichiometric matrix.
useSparseKinetics
=
true
;
{
char
*
word
=
arg
[
iarg
++
];
if
(
strcmp
(
word
,
"sparse"
)
==
0
)
useSparseKinetics
=
true
;
else
if
(
strcmp
(
word
,
"dense"
)
==
0
)
useSparseKinetics
=
false
;
else
{
std
::
string
errmsg
=
"Illegal command "
+
std
::
string
(
word
)
+
" expected
\"
sparse
\"
or
\"
dense
\"\n
"
;
error
->
all
(
FLERR
,
errmsg
.
c_str
());
}
if
(
comm
->
me
==
0
and
Verbosity
>
1
){
std
::
string
msg
=
"FixRX: matrix format is "
;
if
(
useSparseKinetics
)
msg
+=
std
::
string
(
"sparse"
);
else
msg
+=
std
::
string
(
"dense"
);
error
->
message
(
FLERR
,
msg
.
c_str
());
}
}
// Determine the ODE solver/stepper strategy in arg[6].
odeIntegrationFlag
=
ODE_LAMMPS_RK4
;
{
char
*
word
=
arg
[
iarg
++
];
if
(
strcmp
(
word
,
"lammps_rk4"
)
==
0
||
strcmp
(
word
,
"rk4"
)
==
0
)
odeIntegrationFlag
=
ODE_LAMMPS_RK4
;
else
if
(
strcmp
(
word
,
"lammps_rkf45"
)
==
0
||
strcmp
(
word
,
"rkf45"
)
==
0
)
odeIntegrationFlag
=
ODE_LAMMPS_RKF45
;
else
{
std
::
string
errmsg
=
"Illegal ODE integration type: "
+
std
::
string
(
word
);
error
->
all
(
FLERR
,
errmsg
.
c_str
());
}
}
/// Set the default ODE parameters here. Modify with arg[].
/// 'minSteps' has a different meaning for RK4 and RKF45.
/// RK4: This is the # of steps that will be taken with h = dt_dpd / minSteps;
/// RKF45: This sets as h0 = dt_dpd / minSteps. If minSteps == 0, RKF45 will
/// estimate h0 internally. h will be adjusted as needed on subsequent steps.
minSteps
=
1
;
maxIters
=
100
;
relTol
=
1.0e-6
;
absTol
=
1.0e-8
;
diagnosticFrequency
=
0
;
for
(
int
i
=
0
;
i
<
numDiagnosticCounters
;
++
i
){
diagnosticCounter
[
i
]
=
0
;
diagnosticCounterPerODE
[
i
]
=
NULL
;
}
if
(
odeIntegrationFlag
==
ODE_LAMMPS_RK4
&&
narg
==
8
){
char
*
word
=
arg
[
iarg
++
];
minSteps
=
atoi
(
word
);
if
(
comm
->
me
==
0
and
Verbosity
>
1
){
char
msg
[
128
];
sprintf
(
msg
,
"FixRX: RK4 numSteps= %d"
,
minSteps
);
error
->
message
(
FLERR
,
msg
);
}
}
else
if
(
odeIntegrationFlag
==
ODE_LAMMPS_RK4
&&
narg
>
8
){
error
->
all
(
FLERR
,
"Illegal fix rx command. Too many arguments for RK4 solver."
);
}
else
if
(
odeIntegrationFlag
==
ODE_LAMMPS_RKF45
){
// Must have four options.
if
(
narg
<
11
)
error
->
all
(
FLERR
,
"Illegal fix rx command. Too few arguments for RKF45 solver."
);
minSteps
=
atoi
(
arg
[
iarg
++
]
);
maxIters
=
atoi
(
arg
[
iarg
++
]
);
relTol
=
strtod
(
arg
[
iarg
++
],
NULL
);
absTol
=
strtod
(
arg
[
iarg
++
],
NULL
);
if
(
iarg
<
narg
)
diagnosticFrequency
=
atoi
(
arg
[
iarg
++
]
);
// maxIters must be at least minSteps.
maxIters
=
std
::
max
(
minSteps
,
maxIters
);
if
(
comm
->
me
==
0
and
Verbosity
>
1
){
//printf("FixRX: RKF45 minSteps= %d maxIters= %d absTol= %e relTol= %e\n", minSteps, maxIters, absTol, relTol);
char
msg
[
128
];
sprintf
(
msg
,
"FixRX: RKF45 minSteps= %d maxIters= %d relTol= %.1e absTol= %.1e diagnosticFrequency= %d"
,
minSteps
,
maxIters
,
relTol
,
absTol
,
diagnosticFrequency
);
error
->
message
(
FLERR
,
msg
);
}
}
// Initialize/Create the sparse matrix database.
sparseKinetics_nu
=
NULL
;
sparseKinetics_nuk
=
NULL
;
sparseKinetics_inu
=
NULL
;
sparseKinetics_isIntegralReaction
=
NULL
;
sparseKinetics_maxReactants
=
0
;
sparseKinetics_maxProducts
=
0
;
sparseKinetics_maxSpecies
=
0
;
}
/* ---------------------------------------------------------------------- */
FixRX
::~
FixRX
()
{
// De-Allocate memory to prevent memory leak
for
(
int
ii
=
0
;
ii
<
nreactions
;
ii
++
){
delete
[]
stoich
[
ii
];
delete
[]
stoichReactants
[
ii
];
delete
[]
stoichProducts
[
ii
];
}
delete
[]
Arr
;
delete
[]
nArr
;
delete
[]
Ea
;
delete
[]
tempExp
;
delete
[]
stoich
;
delete
[]
stoichReactants
;
delete
[]
stoichProducts
;
delete
[]
kR
;
delete
[]
id_fix_species
;
delete
[]
id_fix_species_old
;
if
(
useSparseKinetics
){
memory
->
destroy
(
sparseKinetics_nu
);
memory
->
destroy
(
sparseKinetics_nuk
);
memory
->
destroy
(
sparseKinetics_inu
);
memory
->
destroy
(
sparseKinetics_isIntegralReaction
);
}
}
/* ---------------------------------------------------------------------- */
void
FixRX
::
post_constructor
()
{
int
maxspecies
=
1000
;
int
nUniqueSpecies
=
0
;
bool
match
;
for
(
int
i
=
0
;
i
<
modify
->
nfix
;
i
++
)
if
(
strncmp
(
modify
->
fix
[
i
]
->
style
,
"property/atom"
,
13
)
==
0
)
error
->
all
(
FLERR
,
"fix rx cannot be combined with fix property/atom"
);
char
**
tmpspecies
=
new
char
*
[
maxspecies
];
for
(
int
jj
=
0
;
jj
<
maxspecies
;
jj
++
)
tmpspecies
[
jj
]
=
NULL
;
// open file on proc 0
FILE
*
fp
;
fp
=
NULL
;
if
(
comm
->
me
==
0
)
{
fp
=
force
->
open_potential
(
kineticsFile
);
if
(
fp
==
NULL
)
{
char
str
[
128
];
sprintf
(
str
,
"Cannot open rx file %s"
,
kineticsFile
);
error
->
one
(
FLERR
,
str
);
}
}
// Assign species names to tmpspecies array and determine the number of unique species
int
n
,
nwords
;
char
line
[
MAXLINE
],
*
ptr
;
int
eof
=
0
;
char
*
word
;
while
(
1
)
{
if
(
comm
->
me
==
0
)
{
ptr
=
fgets
(
line
,
MAXLINE
,
fp
);
if
(
ptr
==
NULL
)
{
eof
=
1
;
fclose
(
fp
);
}
else
n
=
strlen
(
line
)
+
1
;
}
MPI_Bcast
(
&
eof
,
1
,
MPI_INT
,
0
,
world
);
if
(
eof
)
break
;
MPI_Bcast
(
&
n
,
1
,
MPI_INT
,
0
,
world
);
MPI_Bcast
(
line
,
n
,
MPI_CHAR
,
0
,
world
);
// strip comment, skip line if blank
if
((
ptr
=
strchr
(
line
,
'#'
)))
*
ptr
=
'\0'
;
nwords
=
atom
->
count_words
(
line
);
if
(
nwords
==
0
)
continue
;
// words = ptrs to all words in line
nwords
=
0
;
word
=
strtok
(
line
,
"
\t\n\r\f
"
);
while
(
word
!=
NULL
){
word
=
strtok
(
NULL
,
"
\t\n\r\f
"
);
match
=
false
;
for
(
int
jj
=
0
;
jj
<
nUniqueSpecies
;
jj
++
){
if
(
strcmp
(
word
,
tmpspecies
[
jj
])
==
0
){
match
=
true
;
break
;
}
}
if
(
!
match
){
if
(
nUniqueSpecies
+
1
>=
maxspecies
)
error
->
all
(
FLERR
,
"Exceeded the maximum number of species permitted in fix rx."
);
tmpspecies
[
nUniqueSpecies
]
=
new
char
[
strlen
(
word
)
+
1
];
strcpy
(
tmpspecies
[
nUniqueSpecies
],
word
);
nUniqueSpecies
++
;
}
word
=
strtok
(
NULL
,
"
\t\n\r\f
"
);
if
(
strcmp
(
word
,
"+"
)
!=
0
&&
strcmp
(
word
,
"="
)
!=
0
)
break
;
word
=
strtok
(
NULL
,
"
\t\n\r\f
"
);
}
}
atom
->
nspecies_dpd
=
nUniqueSpecies
;
nspecies
=
atom
->
nspecies_dpd
;
// new id = fix-ID + FIX_STORE_ATTRIBUTE
// new fix group = group for this fix
id_fix_species
=
NULL
;
id_fix_species_old
=
NULL
;
n
=
strlen
(
id
)
+
strlen
(
"_SPECIES"
)
+
1
;
id_fix_species
=
new
char
[
n
];
n
=
strlen
(
id
)
+
strlen
(
"_SPECIES_OLD"
)
+
1
;
id_fix_species_old
=
new
char
[
n
];
strcpy
(
id_fix_species
,
id
);
strcat
(
id_fix_species
,
"_SPECIES"
);
strcpy
(
id_fix_species_old
,
id
);
strcat
(
id_fix_species_old
,
"_SPECIES_OLD"
);
char
**
newarg
=
new
char
*
[
nspecies
+
5
];
char
**
newarg2
=
new
char
*
[
nspecies
+
5
];
newarg
[
0
]
=
id_fix_species
;
newarg
[
1
]
=
group
->
names
[
igroup
];
newarg
[
2
]
=
(
char
*
)
"property/atom"
;
newarg2
[
0
]
=
id_fix_species_old
;
newarg2
[
1
]
=
group
->
names
[
igroup
];
newarg2
[
2
]
=
(
char
*
)
"property/atom"
;
for
(
int
ii
=
0
;
ii
<
nspecies
;
ii
++
){
char
str1
[
2
+
strlen
(
tmpspecies
[
ii
])
+
1
];
char
str2
[
2
+
strlen
(
tmpspecies
[
ii
])
+
4
];
strcpy
(
str1
,
"d_"
);
strcpy
(
str2
,
"d_"
);
strncat
(
str1
,
tmpspecies
[
ii
],
strlen
(
tmpspecies
[
ii
]));
strncat
(
str2
,
tmpspecies
[
ii
],
strlen
(
tmpspecies
[
ii
]));
strncat
(
str2
,
"Old"
,
3
);
newarg
[
ii
+
3
]
=
new
char
[
strlen
(
str1
)
+
1
];
newarg2
[
ii
+
3
]
=
new
char
[
strlen
(
str2
)
+
1
];
strcpy
(
newarg
[
ii
+
3
],
str1
);
strcpy
(
newarg2
[
ii
+
3
],
str2
);
}
newarg
[
nspecies
+
3
]
=
(
char
*
)
"ghost"
;
newarg
[
nspecies
+
4
]
=
(
char
*
)
"yes"
;
newarg2
[
nspecies
+
3
]
=
(
char
*
)
"ghost"
;
newarg2
[
nspecies
+
4
]
=
(
char
*
)
"yes"
;
modify
->
add_fix
(
nspecies
+
5
,
newarg
);
fix_species
=
(
FixPropertyAtom
*
)
modify
->
fix
[
modify
->
nfix
-
1
];
modify
->
add_fix
(
nspecies
+
5
,
newarg2
);
fix_species_old
=
(
FixPropertyAtom
*
)
modify
->
fix
[
modify
->
nfix
-
1
];
if
(
nspecies
==
0
)
error
->
all
(
FLERR
,
"There are no rx species specified."
);
for
(
int
jj
=
0
;
jj
<
nspecies
;
jj
++
)
{
delete
[]
tmpspecies
[
jj
];
delete
[]
newarg
[
jj
+
3
];
delete
[]
newarg2
[
jj
+
3
];
}
delete
[]
newarg
;
delete
[]
newarg2
;
delete
[]
tmpspecies
;
read_file
(
kineticsFile
);
if
(
useSparseKinetics
)
this
->
initSparse
();
// set comm size needed by this Pair
comm_forward
=
nspecies
*
2
;
comm_reverse
=
2
;
}
/* ---------------------------------------------------------------------- */
void
FixRX
::
initSparse
()
{
const
int
Verbosity
=
1
;
if
(
comm
->
me
==
0
and
Verbosity
>
1
){
for
(
int
k
=
0
;
k
<
nspecies
;
++
k
)
printf
(
"atom->dname[%d]= %s
\n
"
,
k
,
atom
->
dname
[
k
]);
printf
(
"stoich[][]
\n
"
);
for
(
int
i
=
0
;
i
<
nreactions
;
++
i
){
int
nreac_i
=
0
,
nprod_i
=
0
;
printf
(
"%d: "
,
i
);
for
(
int
k
=
0
;
k
<
nspecies
;
++
k
){
printf
(
" %g"
,
stoich
[
i
][
k
]);
if
(
stoich
[
i
][
k
]
<
0.0
)
nreac_i
++
;
else
if
(
stoich
[
i
][
k
]
>
0.0
)
nprod_i
++
;
}
printf
(
" : %d %d
\n
"
,
nreac_i
,
nprod_i
);
}
printf
(
"stoichReactants[][]
\n
"
);
for
(
int
i
=
0
;
i
<
nreactions
;
++
i
){
int
nreac_i
=
0
;
printf
(
"%d: "
,
i
);
for
(
int
k
=
0
;
k
<
nspecies
;
++
k
){
printf
(
" %g"
,
stoichReactants
[
i
][
k
]);
if
(
stoichReactants
[
i
][
k
]
>
0.0
)
nreac_i
++
;
}
printf
(
" : %d
\n
"
,
nreac_i
);
}
printf
(
"stoichProducts[][]
\n
"
);
for
(
int
i
=
0
;
i
<
nreactions
;
++
i
){
int
nprod_i
=
0
;
printf
(
"%d: "
,
i
);
for
(
int
k
=
0
;
k
<
nspecies
;
++
k
){
printf
(
" %g"
,
stoichProducts
[
i
][
k
]);
if
(
stoichProducts
[
i
][
k
]
>
0.0
)
nprod_i
++
;
}
printf
(
" : %d
\n
"
,
nprod_i
);
}
}
// if (Verbose)
// 1) Measure the sparsity of stoich[][]
int
nzeros
=
0
;
int
mxprod
=
0
;
int
mxreac
=
0
;
int
mxspec
=
0
;
int
nIntegral
=
0
;
for
(
int
i
=
0
;
i
<
nreactions
;
++
i
){
int
nreac_i
=
0
,
nprod_i
=
0
;
std
::
string
pstr
,
rstr
;
bool
allAreIntegral
=
true
;
for
(
int
k
=
0
;
k
<
nspecies
;
++
k
){
if
(
stoichReactants
[
i
][
k
]
==
0
and
stoichProducts
[
i
][
k
]
==
0
)
nzeros
++
;
if
(
stoichReactants
[
i
][
k
]
>
0.0
){
allAreIntegral
&=
(
std
::
fmod
(
stoichReactants
[
i
][
k
],
1.0
)
==
0.0
);
nreac_i
++
;
if
(
rstr
.
length
()
>
0
)
rstr
+=
" + "
;
char
digit
[
6
];
sprintf
(
digit
,
"%4.1f "
,
stoichReactants
[
i
][
k
]);
rstr
+=
digit
;
rstr
+=
atom
->
dname
[
k
];
}
if
(
stoichProducts
[
i
][
k
]
>
0.0
){
allAreIntegral
&=
(
std
::
fmod
(
stoichProducts
[
i
][
k
],
1.0
)
==
0.0
);
nprod_i
++
;
if
(
pstr
.
length
()
>
0
)
pstr
+=
" + "
;
char
digit
[
6
];
sprintf
(
digit
,
"%4.1f "
,
stoichProducts
[
i
][
k
]);
pstr
+=
digit
;
pstr
+=
atom
->
dname
[
k
];
}
}
if
(
comm
->
me
==
0
and
Verbosity
>
1
)
printf
(
"rx%3d: %d %d %d ... %s %s %s
\n
"
,
i
,
nreac_i
,
nprod_i
,
allAreIntegral
,
rstr
.
c_str
(),
/*reversible[i]*/
(
false
)
?
"<=>"
:
"="
,
pstr
.
c_str
());
mxreac
=
std
::
max
(
mxreac
,
nreac_i
);
mxprod
=
std
::
max
(
mxprod
,
nprod_i
);
mxspec
=
std
::
max
(
mxspec
,
nreac_i
+
nprod_i
);
if
(
allAreIntegral
)
nIntegral
++
;
}
if
(
comm
->
me
==
0
and
Verbosity
>
1
){
char
msg
[
256
];
sprintf
(
msg
,
"FixRX: Sparsity of Stoichiometric Matrix= %.1f%% non-zeros= %d nspecies= %d nreactions= %d maxReactants= %d maxProducts= %d maxSpecies= %d integralReactions= %d"
,
100
*
(
double
(
nzeros
)
/
(
nspecies
*
nreactions
)),
nzeros
,
nspecies
,
nreactions
,
mxreac
,
mxprod
,
(
mxreac
+
mxprod
),
SparseKinetics_enableIntegralReactions
);
error
->
message
(
FLERR
,
msg
);
}
// Allocate the sparse matrix data.
{
sparseKinetics_maxSpecies
=
(
mxreac
+
mxprod
);
sparseKinetics_maxReactants
=
mxreac
;
sparseKinetics_maxProducts
=
mxprod
;
memory
->
create
(
sparseKinetics_nu
,
nreactions
,
sparseKinetics_maxSpecies
,
"sparseKinetics_nu"
);
memory
->
create
(
sparseKinetics_nuk
,
nreactions
,
sparseKinetics_maxSpecies
,
"sparseKinetics_nuk"
);
for
(
int
i
=
0
;
i
<
nreactions
;
++
i
)
for
(
int
k
=
0
;
k
<
sparseKinetics_maxSpecies
;
++
k
){
sparseKinetics_nu
[
i
][
k
]
=
0.0
;
sparseKinetics_nuk
[
i
][
k
]
=
SparseKinetics_invalidIndex
;
// Initialize with an invalid index.
}
if
(
SparseKinetics_enableIntegralReactions
){
memory
->
create
(
sparseKinetics_inu
,
nreactions
,
sparseKinetics_maxSpecies
,
"sparseKinetics_inu"
);
memory
->
create
(
sparseKinetics_isIntegralReaction
,
nreactions
,
"sparseKinetics_isIntegralReaction"
);
for
(
int
i
=
0
;
i
<
nreactions
;
++
i
){
sparseKinetics_isIntegralReaction
[
i
]
=
false
;
for
(
int
k
=
0
;
k
<
sparseKinetics_maxSpecies
;
++
k
)
sparseKinetics_inu
[
i
][
k
]
=
0
;
}
}
}
// Measure the distribution of the # of moles for the ::fastpowi function.
std
::
vector
<
int
>
nu_bin
(
10
);
for
(
int
i
=
0
;
i
<
nreactions
;
++
i
){
int
nreac_i
=
0
,
nprod_i
=
0
;
bool
isIntegral_i
=
true
;
for
(
int
k
=
0
;
k
<
nspecies
;
++
k
){
if
(
stoichReactants
[
i
][
k
]
>
0.0
){
const
int
idx
=
nreac_i
;
sparseKinetics_nu
[
i
][
idx
]
=
stoichReactants
[
i
][
k
];
sparseKinetics_nuk
[
i
][
idx
]
=
k
;
isIntegral_i
&=
(
std
::
fmod
(
stoichReactants
[
i
][
k
],
1.0
)
==
0.0
);
if
(
SparseKinetics_enableIntegralReactions
){
sparseKinetics_inu
[
i
][
idx
]
=
(
int
)
sparseKinetics_nu
[
i
][
idx
];
if
(
isIntegral_i
){
if
(
sparseKinetics_inu
[
i
][
idx
]
>=
nu_bin
.
size
())
nu_bin
.
resize
(
sparseKinetics_inu
[
i
][
idx
]
);
nu_bin
[
sparseKinetics_inu
[
i
][
idx
]
]
++
;
}
}
nreac_i
++
;
}
if
(
stoichProducts
[
i
][
k
]
>
0.0
){
const
int
idx
=
sparseKinetics_maxReactants
+
nprod_i
;
sparseKinetics_nu
[
i
][
idx
]
=
stoichProducts
[
i
][
k
];
sparseKinetics_nuk
[
i
][
idx
]
=
k
;
isIntegral_i
&=
(
std
::
fmod
(
sparseKinetics_nu
[
i
][
idx
],
1.0
)
==
0.0
);
if
(
SparseKinetics_enableIntegralReactions
){
sparseKinetics_inu
[
i
][
idx
]
=
(
int
)
sparseKinetics_nu
[
i
][
idx
];
if
(
isIntegral_i
){
if
(
sparseKinetics_inu
[
i
][
idx
]
>=
nu_bin
.
size
())
nu_bin
.
resize
(
sparseKinetics_inu
[
i
][
idx
]
);
nu_bin
[
sparseKinetics_inu
[
i
][
idx
]
]
++
;
}
}
nprod_i
++
;
}
}
if
(
SparseKinetics_enableIntegralReactions
)
sparseKinetics_isIntegralReaction
[
i
]
=
isIntegral_i
;
}
if
(
comm
->
me
==
0
and
Verbosity
>
1
){
for
(
int
i
=
1
;
i
<
nu_bin
.
size
();
++
i
)
if
(
nu_bin
[
i
]
>
0
)
printf
(
"nu_bin[%d] = %d
\n
"
,
i
,
nu_bin
[
i
]);
for
(
int
i
=
0
;
i
<
nreactions
;
++
i
){
std
::
string
pstr
,
rstr
;
for
(
int
kk
=
0
;
kk
<
sparseKinetics_maxReactants
;
kk
++
){
const
int
k
=
sparseKinetics_nuk
[
i
][
kk
];
if
(
k
!=
SparseKinetics_invalidIndex
){
if
(
rstr
.
length
()
>
0
)
rstr
+=
" + "
;
char
digit
[
6
];
if
(
SparseKinetics_enableIntegralReactions
and
sparseKinetics_isIntegralReaction
[
i
])
sprintf
(
digit
,
"%d "
,
sparseKinetics_inu
[
i
][
kk
]);
else
sprintf
(
digit
,
"%4.1f "
,
sparseKinetics_nu
[
i
][
kk
]);
rstr
+=
digit
;
rstr
+=
atom
->
dname
[
k
];
}
}
for
(
int
kk
=
sparseKinetics_maxReactants
;
kk
<
sparseKinetics_maxSpecies
;
kk
++
){
const
int
k
=
sparseKinetics_nuk
[
i
][
kk
];
if
(
k
!=
SparseKinetics_invalidIndex
){
if
(
pstr
.
length
()
>
0
)
pstr
+=
" + "
;
char
digit
[
6
];
if
(
SparseKinetics_enableIntegralReactions
and
sparseKinetics_isIntegralReaction
[
i
])
sprintf
(
digit
,
"%d "
,
sparseKinetics_inu
[
i
][
kk
]);
else
sprintf
(
digit
,
"%4.1f "
,
sparseKinetics_nu
[
i
][
kk
]);
pstr
+=
digit
;
pstr
+=
atom
->
dname
[
k
];
}
}
if
(
comm
->
me
==
0
and
Verbosity
>
1
)
printf
(
"rx%3d: %s %s %s
\n
"
,
i
,
rstr
.
c_str
(),
/*reversible[i]*/
(
false
)
?
"<=>"
:
"="
,
pstr
.
c_str
());
}
// end for nreactions
}
// end if Verbose
}
/* ---------------------------------------------------------------------- */
int
FixRX
::
setmask
()
{
int
mask
=
0
;
mask
|=
PRE_FORCE
;
return
mask
;
}
/* ---------------------------------------------------------------------- */
void
FixRX
::
init
()
{
pairDPDE
=
(
PairDPDfdtEnergy
*
)
force
->
pair_match
(
"dpd/fdt/energy"
,
1
);
if
(
pairDPDE
==
NULL
)
error
->
all
(
FLERR
,
"Must use pair_style dpd/fdt/energy with fix rx"
);
bool
eos_flag
=
false
;
for
(
int
i
=
0
;
i
<
modify
->
nfix
;
i
++
)
if
(
strcmp
(
modify
->
fix
[
i
]
->
style
,
"eos/table/rx"
)
==
0
)
eos_flag
=
true
;
if
(
!
eos_flag
)
error
->
all
(
FLERR
,
"fix rx requires fix eos/table/rx to be specified"
);
}
/* ---------------------------------------------------------------------- */
void
FixRX
::
setup_pre_force
(
int
vflag
)
{
int
nlocal
=
atom
->
nlocal
;
int
nghost
=
atom
->
nghost
;
int
*
mask
=
atom
->
mask
;
int
newton_pair
=
force
->
newton_pair
;
double
tmp
;
int
ii
;
if
(
localTempFlag
){
if
(
newton_pair
)
{
dpdThetaLocal
=
new
double
[
nlocal
+
nghost
];
for
(
ii
=
0
;
ii
<
nlocal
+
nghost
;
ii
++
)
dpdThetaLocal
[
ii
]
=
0.0
;
}
else
{
dpdThetaLocal
=
new
double
[
nlocal
];
for
(
ii
=
0
;
ii
<
nlocal
;
ii
++
)
dpdThetaLocal
[
ii
]
=
0.0
;
}
computeLocalTemperature
();
}
for
(
int
id
=
0
;
id
<
nlocal
;
id
++
)
for
(
int
ispecies
=
0
;
ispecies
<
nspecies
;
ispecies
++
){
tmp
=
atom
->
dvector
[
ispecies
][
id
];
atom
->
dvector
[
ispecies
+
nspecies
][
id
]
=
tmp
;
}
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
){
// Set the reaction rate constants to zero: no reactions occur at step 0
for
(
int
irxn
=
0
;
irxn
<
nreactions
;
irxn
++
)
kR
[
irxn
]
=
0.0
;
if
(
odeIntegrationFlag
==
ODE_LAMMPS_RK4
)
rk4
(
i
,
NULL
);
else
if
(
odeIntegrationFlag
==
ODE_LAMMPS_RKF45
)
rkf45
(
i
,
NULL
);
}
// Communicate the updated momenta and velocities to all nodes
comm
->
forward_comm_fix
(
this
);
if
(
localTempFlag
)
delete
[]
dpdThetaLocal
;
}
/* ---------------------------------------------------------------------- */
void
FixRX
::
pre_force
(
int
vflag
)
{
int
nlocal
=
atom
->
nlocal
;
int
nghost
=
atom
->
nghost
;
int
*
mask
=
atom
->
mask
;
double
*
dpdTheta
=
atom
->
dpdTheta
;
int
newton_pair
=
force
->
newton_pair
;
int
ii
;
double
theta
;
if
(
localTempFlag
){
if
(
newton_pair
)
{
dpdThetaLocal
=
new
double
[
nlocal
+
nghost
];
for
(
ii
=
0
;
ii
<
nlocal
+
nghost
;
ii
++
)
dpdThetaLocal
[
ii
]
=
0.0
;
}
else
{
dpdThetaLocal
=
new
double
[
nlocal
];
for
(
ii
=
0
;
ii
<
nlocal
;
ii
++
)
dpdThetaLocal
[
ii
]
=
0.0
;
}
computeLocalTemperature
();
}
TimerType
timer_localTemperature
=
getTimeStamp
();
// Zero the counters for the ODE solvers.
this
->
nSteps
=
this
->
nIters
=
this
->
nFuncs
=
this
->
nFails
=
0
;
if
(
odeIntegrationFlag
==
ODE_LAMMPS_RKF45
&&
diagnosticFrequency
==
1
)
{
memory
->
create
(
diagnosticCounterPerODE
[
StepSum
],
nlocal
,
"FixRX::diagnosticCounterPerODE"
);
memory
->
create
(
diagnosticCounterPerODE
[
FuncSum
],
nlocal
,
"FixRX::diagnosticCounterPerODE"
);
}
double
*
rwork
=
new
double
[
8
*
nspecies
+
nreactions
];
for
(
int
i
=
0
;
i
<
nlocal
;
i
++
)
if
(
mask
[
i
]
&
groupbit
){
if
(
localTempFlag
)
theta
=
dpdThetaLocal
[
i
];
else
theta
=
dpdTheta
[
i
];
//Compute the reaction rate constants
for
(
int
irxn
=
0
;
irxn
<
nreactions
;
irxn
++
)
kR
[
irxn
]
=
Arr
[
irxn
]
*
pow
(
theta
,
nArr
[
irxn
])
*
exp
(
-
Ea
[
irxn
]
/
force
->
boltz
/
theta
);
if
(
odeIntegrationFlag
==
ODE_LAMMPS_RK4
)
rk4
(
i
,
rwork
);
else
if
(
odeIntegrationFlag
==
ODE_LAMMPS_RKF45
)
rkf45
(
i
,
rwork
);
}
TimerType
timer_ODE
=
getTimeStamp
();
delete
[]
rwork
;
// Communicate the updated momenta and velocities to all nodes
comm
->
forward_comm_fix
(
this
);
if
(
localTempFlag
)
delete
[]
dpdThetaLocal
;
double
time_ODE
=
getElapsedTime
(
timer_localTemperature
,
timer_ODE
);
// Warn the user if a failure was detected in the ODE solver.
if
(
nFails
>
0
){
char
sbuf
[
128
];
sprintf
(
sbuf
,
"in FixRX::pre_force, ODE solver failed for %d atoms."
,
nFails
);
error
->
warning
(
FLERR
,
sbuf
);
}
// Compute and report ODE diagnostics, if requested.
if
(
odeIntegrationFlag
==
ODE_LAMMPS_RKF45
&&
diagnosticFrequency
!=
0
){
// Update the counters.
diagnosticCounter
[
StepSum
]
+=
nSteps
;
diagnosticCounter
[
FuncSum
]
+=
nFuncs
;
diagnosticCounter
[
TimeSum
]
+=
time_ODE
;
diagnosticCounter
[
AtomSum
]
+=
nlocal
;
diagnosticCounter
[
numDiagnosticCounters
-
1
]
++
;
if
(
(
diagnosticFrequency
>
0
&&
((
update
->
ntimestep
-
update
->
firststep
)
%
diagnosticFrequency
)
==
0
)
||
(
diagnosticFrequency
<
0
&&
update
->
ntimestep
==
update
->
laststep
)
)
this
->
odeDiagnostics
();
for
(
int
i
=
0
;
i
<
numDiagnosticCounters
;
++
i
)
if
(
diagnosticCounterPerODE
[
i
])
memory
->
destroy
(
diagnosticCounterPerODE
[
i
]
);
}
}
/* ---------------------------------------------------------------------- */
void
FixRX
::
read_file
(
char
*
file
)
{
nreactions
=
0
;
// open file on proc 0
FILE
*
fp
;
fp
=
NULL
;
if
(
comm
->
me
==
0
)
{
fp
=
force
->
open_potential
(
file
);
if
(
fp
==
NULL
)
{
char
str
[
128
];
sprintf
(
str
,
"Cannot open rx file %s"
,
file
);
error
->
one
(
FLERR
,
str
);
}
}
// Count the number of reactions from kinetics file
int
n
,
nwords
,
ispecies
;
char
line
[
MAXLINE
],
*
ptr
;
int
eof
=
0
;
while
(
1
)
{
if
(
comm
->
me
==
0
)
{
ptr
=
fgets
(
line
,
MAXLINE
,
fp
);
if
(
ptr
==
NULL
)
{
eof
=
1
;
fclose
(
fp
);
}
else
n
=
strlen
(
line
)
+
1
;
}
MPI_Bcast
(
&
eof
,
1
,
MPI_INT
,
0
,
world
);
if
(
eof
)
break
;
MPI_Bcast
(
&
n
,
1
,
MPI_INT
,
0
,
world
);
MPI_Bcast
(
line
,
n
,
MPI_CHAR
,
0
,
world
);
// strip comment, skip line if blank
if
((
ptr
=
strchr
(
line
,
'#'
)))
*
ptr
=
'\0'
;
nwords
=
atom
->
count_words
(
line
);
if
(
nwords
==
0
)
continue
;
nreactions
++
;
}
// open file on proc 0
if
(
comm
->
me
==
0
)
fp
=
force
->
open_potential
(
file
);
// read each reaction from kinetics file
eof
=
0
;
char
*
word
;
double
tmpStoich
;
double
sign
;
Arr
=
new
double
[
nreactions
];
nArr
=
new
double
[
nreactions
];
Ea
=
new
double
[
nreactions
];
tempExp
=
new
double
[
nreactions
];
stoich
=
new
double
*
[
nreactions
];
stoichReactants
=
new
double
*
[
nreactions
];
stoichProducts
=
new
double
*
[
nreactions
];
for
(
int
ii
=
0
;
ii
<
nreactions
;
ii
++
){
stoich
[
ii
]
=
new
double
[
nspecies
];
stoichReactants
[
ii
]
=
new
double
[
nspecies
];
stoichProducts
[
ii
]
=
new
double
[
nspecies
];
}
kR
=
new
double
[
nreactions
];
for
(
int
ii
=
0
;
ii
<
nreactions
;
ii
++
){
for
(
int
jj
=
0
;
jj
<
nspecies
;
jj
++
){
stoich
[
ii
][
jj
]
=
0.0
;
stoichReactants
[
ii
][
jj
]
=
0.0
;
stoichProducts
[
ii
][
jj
]
=
0.0
;
}
}
nreactions
=
0
;
sign
=
-
1.0
;
while
(
1
)
{
if
(
comm
->
me
==
0
)
{
ptr
=
fgets
(
line
,
MAXLINE
,
fp
);
if
(
ptr
==
NULL
)
{
eof
=
1
;
fclose
(
fp
);
}
else
n
=
strlen
(
line
)
+
1
;
}
MPI_Bcast
(
&
eof
,
1
,
MPI_INT
,
0
,
world
);
if
(
eof
)
break
;
MPI_Bcast
(
&
n
,
1
,
MPI_INT
,
0
,
world
);
MPI_Bcast
(
line
,
n
,
MPI_CHAR
,
0
,
world
);
// strip comment, skip line if blank
if
((
ptr
=
strchr
(
line
,
'#'
)))
*
ptr
=
'\0'
;
nwords
=
atom
->
count_words
(
line
);
if
(
nwords
==
0
)
continue
;
// words = ptrs to all words in line
nwords
=
0
;
word
=
strtok
(
line
,
"
\t\n\r\f
"
);
while
(
word
!=
NULL
){
tmpStoich
=
atof
(
word
);
word
=
strtok
(
NULL
,
"
\t\n\r\f
"
);
for
(
ispecies
=
0
;
ispecies
<
nspecies
;
ispecies
++
){
if
(
strcmp
(
word
,
&
atom
->
dname
[
ispecies
][
0
])
==
0
){
stoich
[
nreactions
][
ispecies
]
+=
sign
*
tmpStoich
;
if
(
sign
<
0.0
)
stoichReactants
[
nreactions
][
ispecies
]
+=
tmpStoich
;
else
stoichProducts
[
nreactions
][
ispecies
]
+=
tmpStoich
;
break
;
}
}
if
(
ispecies
==
nspecies
){
if
(
comm
->
me
)
{
fprintf
(
stderr
,
"%s mol fraction is not found in data file
\n
"
,
word
);
fprintf
(
stderr
,
"nspecies=%d ispecies=%d
\n
"
,
nspecies
,
ispecies
);
}
error
->
all
(
FLERR
,
"Illegal fix rx command"
);
}
word
=
strtok
(
NULL
,
"
\t\n\r\f
"
);
if
(
word
==
NULL
)
error
->
all
(
FLERR
,
"Missing parameters in reaction kinetic equation"
);
if
(
strcmp
(
word
,
"="
)
==
0
)
sign
=
1.0
;
if
(
strcmp
(
word
,
"+"
)
!=
0
&&
strcmp
(
word
,
"="
)
!=
0
){
if
(
word
==
NULL
)
error
->
all
(
FLERR
,
"Missing parameters in reaction kinetic equation"
);
Arr
[
nreactions
]
=
atof
(
word
);
word
=
strtok
(
NULL
,
"
\t\n\r\f
"
);
if
(
word
==
NULL
)
error
->
all
(
FLERR
,
"Missing parameters in reaction kinetic equation"
);
nArr
[
nreactions
]
=
atof
(
word
);
word
=
strtok
(
NULL
,
"
\t\n\r\f
"
);
if
(
word
==
NULL
)
error
->
all
(
FLERR
,
"Missing parameters in reaction kinetic equation"
);
Ea
[
nreactions
]
=
atof
(
word
);
sign
=
-
1.0
;
break
;
}
word
=
strtok
(
NULL
,
"
\t\n\r\f
"
);
}
nreactions
++
;
}
}
/* ---------------------------------------------------------------------- */
void
FixRX
::
setupParams
()
{
int
i
,
j
,
n
;
// set mol2param for all combinations
// must be a single exact match to lines read from file
memory
->
destroy
(
mol2param
);
memory
->
create
(
mol2param
,
nspecies
,
"pair:mol2param"
);
for
(
i
=
0
;
i
<
nspecies
;
i
++
)
{
n
=
-
1
;
for
(
j
=
0
;
j
<
nreactions
;
j
++
)
{
if
(
i
==
params
[
j
].
ispecies
)
{
if
(
n
>=
0
)
error
->
all
(
FLERR
,
"Potential file has duplicate entry"
);
n
=
j
;
}
}
mol2param
[
i
]
=
n
;
}
}
/* ---------------------------------------------------------------------- */
void
FixRX
::
rk4
(
int
id
,
double
*
rwork
)
{
double
*
k1
=
NULL
;
if
(
rwork
==
NULL
)
k1
=
new
double
[
6
*
nspecies
+
nreactions
];
else
k1
=
rwork
;
double
*
k2
=
k1
+
nspecies
;
double
*
k3
=
k2
+
nspecies
;
double
*
k4
=
k3
+
nspecies
;
double
*
y
=
k4
+
nspecies
;
double
*
yp
=
y
+
nspecies
;
double
*
dummyArray
=
yp
+
nspecies
;
// Passed to the rhs function.
const
int
numSteps
=
minSteps
;
const
double
h
=
update
->
dt
/
double
(
numSteps
);
// Update ConcOld
for
(
int
ispecies
=
0
;
ispecies
<
nspecies
;
ispecies
++
)
{
const
double
tmp
=
atom
->
dvector
[
ispecies
][
id
];
atom
->
dvector
[
ispecies
+
nspecies
][
id
]
=
tmp
;
y
[
ispecies
]
=
tmp
;
}
// Run the requested steps with h.
for
(
int
step
=
0
;
step
<
numSteps
;
step
++
)
{
// k1
rhs
(
0.0
,
y
,
k1
,
dummyArray
);
// k2
for
(
int
ispecies
=
0
;
ispecies
<
nspecies
;
ispecies
++
)
yp
[
ispecies
]
=
y
[
ispecies
]
+
0.5
*
h
*
k1
[
ispecies
];
rhs
(
0.0
,
yp
,
k2
,
dummyArray
);
// k3
for
(
int
ispecies
=
0
;
ispecies
<
nspecies
;
ispecies
++
)
yp
[
ispecies
]
=
y
[
ispecies
]
+
0.5
*
h
*
k2
[
ispecies
];
rhs
(
0.0
,
yp
,
k3
,
dummyArray
);
// k4
for
(
int
ispecies
=
0
;
ispecies
<
nspecies
;
ispecies
++
)
yp
[
ispecies
]
=
y
[
ispecies
]
+
h
*
k3
[
ispecies
];
rhs
(
0.0
,
yp
,
k4
,
dummyArray
);
for
(
int
ispecies
=
0
;
ispecies
<
nspecies
;
ispecies
++
)
y
[
ispecies
]
+=
h
*
(
k1
[
ispecies
]
/
6.0
+
k2
[
ispecies
]
/
3.0
+
k3
[
ispecies
]
/
3.0
+
k4
[
ispecies
]
/
6.0
);
}
// end for (int step...
// Store the solution back in atom->dvector.
for
(
int
ispecies
=
0
;
ispecies
<
nspecies
;
ispecies
++
){
if
(
y
[
ispecies
]
<
-
1.0e-10
)
error
->
one
(
FLERR
,
"Computed concentration in RK4 solver is < -1.0e-10"
);
else
if
(
y
[
ispecies
]
<
1e-15
)
y
[
ispecies
]
=
0.0
;
atom
->
dvector
[
ispecies
][
id
]
=
y
[
ispecies
];
}
if
(
rwork
==
NULL
)
delete
[]
k1
;
}
/* ---------------------------------------------------------------------- */
// f1 = dt*f(t,x)
// f2 = dt*f(t+ c20*dt,x + c21*f1)
// f3 = dt*f(t+ c30*dt,x + c31*f1 + c32*f2)
// f4 = dt*f(t+ c40*dt,x + c41*f1 + c42*f2 + c43*f3)
// f5 = dt*f(t+dt,x + c51*f1 + c52*f2 + c53*f3 + c54*f4)
// f6 = dt*f(t+ c60*dt,x + c61*f1 + c62*f2 + c63*f3 + c64*f4 + c65*f5)
//
// fifth-order runge-kutta integration
// x5 = x + b1*f1 + b3*f3 + b4*f4 + b5*f5 + b6*f6
// fourth-order runge-kutta integration
// x = x + a1*f1 + a3*f3 + a4*f4 + a5*f5
void
FixRX
::
rkf45_step
(
const
int
neq
,
const
double
h
,
double
y
[],
double
y_out
[],
double
rwk
[],
void
*
v_param
)
{
const
double
c21
=
0.25
;
const
double
c31
=
0.09375
;
const
double
c32
=
0.28125
;
const
double
c41
=
0.87938097405553
;
const
double
c42
=-
3.2771961766045
;
const
double
c43
=
3.3208921256258
;
const
double
c51
=
2.0324074074074
;
const
double
c52
=-
8.0
;
const
double
c53
=
7.1734892787524
;
const
double
c54
=-
0.20589668615984
;
const
double
c61
=-
0.2962962962963
;
const
double
c62
=
2.0
;
const
double
c63
=-
1.3816764132554
;
const
double
c64
=
0.45297270955166
;
const
double
c65
=-
0.275
;
const
double
a1
=
0.11574074074074
;
const
double
a3
=
0.54892787524366
;
const
double
a4
=
0.5353313840156
;
const
double
a5
=-
0.2
;
const
double
b1
=
0.11851851851852
;
const
double
b3
=
0.51898635477583
;
const
double
b4
=
0.50613149034201
;
const
double
b5
=-
0.18
;
const
double
b6
=
0.036363636363636
;
// local dependent variables (5 total)
double
*
f1
=
&
rwk
[
0
];
double
*
f2
=
&
rwk
[
neq
];
double
*
f3
=
&
rwk
[
2
*
neq
];
double
*
f4
=
&
rwk
[
3
*
neq
];
double
*
f5
=
&
rwk
[
4
*
neq
];
double
*
f6
=
&
rwk
[
5
*
neq
];
// scratch for the intermediate solution.
//double* ytmp = &rwk[6*neq];
double
*
ytmp
=
y_out
;
// 1)
rhs
(
0.0
,
y
,
f1
,
v_param
);
for
(
int
k
=
0
;
k
<
neq
;
k
++
){
f1
[
k
]
*=
h
;
ytmp
[
k
]
=
y
[
k
]
+
c21
*
f1
[
k
];
}
// 2)
rhs
(
0.0
,
ytmp
,
f2
,
v_param
);
for
(
int
k
=
0
;
k
<
neq
;
k
++
){
f2
[
k
]
*=
h
;
ytmp
[
k
]
=
y
[
k
]
+
c31
*
f1
[
k
]
+
c32
*
f2
[
k
];
}
// 3)
rhs
(
0.0
,
ytmp
,
f3
,
v_param
);
for
(
int
k
=
0
;
k
<
neq
;
k
++
)
{
f3
[
k
]
*=
h
;
ytmp
[
k
]
=
y
[
k
]
+
c41
*
f1
[
k
]
+
c42
*
f2
[
k
]
+
c43
*
f3
[
k
];
}
// 4)
rhs
(
0.0
,
ytmp
,
f4
,
v_param
);
for
(
int
k
=
0
;
k
<
neq
;
k
++
)
{
f4
[
k
]
*=
h
;
ytmp
[
k
]
=
y
[
k
]
+
c51
*
f1
[
k
]
+
c52
*
f2
[
k
]
+
c53
*
f3
[
k
]
+
c54
*
f4
[
k
];
}
// 5)
rhs
(
0.0
,
ytmp
,
f5
,
v_param
);
for
(
int
k
=
0
;
k
<
neq
;
k
++
)
{
f5
[
k
]
*=
h
;
ytmp
[
k
]
=
y
[
k
]
+
c61
*
f1
[
k
]
+
c62
*
f2
[
k
]
+
c63
*
f3
[
k
]
+
c64
*
f4
[
k
]
+
c65
*
f5
[
k
];
}
// 6)
rhs
(
0.0
,
ytmp
,
f6
,
v_param
);
for
(
int
k
=
0
;
k
<
neq
;
k
++
)
{
//const double f6 = h * ydot[k];
f6
[
k
]
*=
h
;
// 5th-order solution.
const
double
r5
=
b1
*
f1
[
k
]
+
b3
*
f3
[
k
]
+
b4
*
f4
[
k
]
+
b5
*
f5
[
k
]
+
b6
*
f6
[
k
];
// 4th-order solution.
const
double
r4
=
a1
*
f1
[
k
]
+
a3
*
f3
[
k
]
+
a4
*
f4
[
k
]
+
a5
*
f5
[
k
];
// Truncation error: difference between 4th and 5th-order solutions.
rwk
[
k
]
=
fabs
(
r5
-
r4
);
// Update solution.
//y_out[k] = y[k] + r5; // Local extrapolation
y_out
[
k
]
=
y
[
k
]
+
r4
;
}
return
;
}
int
FixRX
::
rkf45_h0
(
const
int
neq
,
const
double
t
,
const
double
t_stop
,
const
double
hmin
,
const
double
hmax
,
double
&
h0
,
double
y
[],
double
rwk
[],
void
*
v_params
)
{
// Set lower and upper bounds on h0, and take geometric mean as first trial value.
// Exit with this value if the bounds cross each other.
// Adjust upper bound based on ydot ...
double
hg
=
sqrt
(
hmin
*
hmax
);
//if (hmax < hmin)
//{
// h0 = hg;
// return;
//}
// Start iteration to find solution to ... {WRMS norm of (h0^2 y'' / 2)} = 1
double
*
ydot
=
rwk
;
double
*
y1
=
ydot
+
neq
;
double
*
ydot1
=
y1
+
neq
;
const
int
max_iters
=
10
;
bool
hnew_is_ok
=
false
;
double
hnew
=
hg
;
int
iter
=
0
;
// compute ydot at t=t0
rhs
(
t
,
y
,
ydot
,
v_params
);
while
(
1
)
{
// Estimate y'' with finite-difference ...
for
(
int
k
=
0
;
k
<
neq
;
k
++
)
y1
[
k
]
=
y
[
k
]
+
hg
*
ydot
[
k
];
// compute y' at t1
rhs
(
t
+
hg
,
y1
,
ydot1
,
v_params
);
// Compute WRMS norm of y''
double
yddnrm
=
0.0
;
for
(
int
k
=
0
;
k
<
neq
;
k
++
){
double
ydd
=
(
ydot1
[
k
]
-
ydot
[
k
])
/
hg
;
double
wterr
=
ydd
/
(
relTol
*
fabs
(
y
[
k
]
)
+
absTol
);
yddnrm
+=
wterr
*
wterr
;
}
yddnrm
=
sqrt
(
yddnrm
/
double
(
neq
)
);
//std::cout << "iter " << _iter << " hg " << hg << " y'' " << yddnrm << std::endl;
//std::cout << "ydot " << ydot[neq-1] << std::endl;
// should we accept this?
if
(
hnew_is_ok
||
iter
==
max_iters
){
hnew
=
hg
;
if
(
iter
==
max_iters
)
fprintf
(
stderr
,
"ERROR_HIN_MAX_ITERS
\n
"
);
break
;
}
// Get the new value of h ...
hnew
=
(
yddnrm
*
hmax
*
hmax
>
2.0
)
?
sqrt
(
2.0
/
yddnrm
)
:
sqrt
(
hg
*
hmax
);
// test the stopping conditions.
double
hrat
=
hnew
/
hg
;
// Accept this value ... the bias factor should bring it within range.
if
(
(
hrat
>
0.5
)
&&
(
hrat
<
2.0
)
)
hnew_is_ok
=
true
;
// If y'' is still bad after a few iterations, just accept h and give up.
if
(
(
iter
>
1
)
&&
hrat
>
2.0
)
{
hnew
=
hg
;
hnew_is_ok
=
true
;
}
//printf("iter=%d, yddnrw=%e, hnew=%e, hmin=%e, hmax=%e\n", iter, yddnrm, hnew, hmin, hmax);
hg
=
hnew
;
iter
++
;
}
// bound and bias estimate
h0
=
hnew
*
0.5
;
h0
=
fmax
(
h0
,
hmin
);
h0
=
fmin
(
h0
,
hmax
);
//printf("h0=%e, hmin=%e, hmax=%e\n", h0, hmin, hmax);
return
(
iter
+
1
);
}
void
FixRX
::
odeDiagnostics
(
void
)
{
TimerType
timer_start
=
getTimeStamp
();
// Compute:
// 1) Average # of ODE integrator steps and RHS evaluations per atom globally.
// 2) RMS # of ...
// 3) Average # of ODE steps and RHS evaluations per MPI task.
// 4) RMS # of ODE steps and RHS evaluations per MPI task.
// 5) MAX # of ODE steps and RHS evaluations per MPI task.
//
// ... 1,2 are for ODE control diagnostics.
// ... 3-5 are for load balancing diagnostics.
//
// To do this, we'll need to
// a) Allreduce (sum) the sum of nSteps / nFuncs. Dividing by atom->natoms
// gives the avg # of steps/funcs per atom globally.
// b) Reduce (sum) to root the sum of squares of the differences.
// i) Sum_i (steps_i - avg_steps_global)^2
// ii) Sum_i (funcs_i - avg_funcs_global)^2
// iii) (avg_steps_local - avg_steps_global)^2
// iv) (avg_funcs_local - avg_funcs_global)^2
const
int
numCounters
=
numDiagnosticCounters
-
1
;
// # of time-steps for averaging.
const
int
nTimes
=
this
->
diagnosticCounter
[
numDiagnosticCounters
-
1
];
// # of ODE's per time-step (on average).
//const int nODEs = this->diagnosticCounter[AtomSum] / nTimes;
// Sum up the sums from each task.
double
sums
[
numCounters
];
double
my_vals
[
numCounters
];
double
max_per_proc
[
numCounters
];
double
min_per_proc
[
numCounters
];
// Compute counters per dpd time-step.
for
(
int
i
=
0
;
i
<
numCounters
;
++
i
){
my_vals
[
i
]
=
this
->
diagnosticCounter
[
i
]
/
nTimes
;
//printf("my sum[%d] = %f %d\n", i, my_vals[i], comm->me);
}
MPI_Allreduce
(
my_vals
,
sums
,
numCounters
,
MPI_DOUBLE
,
MPI_SUM
,
world
);
MPI_Reduce
(
my_vals
,
max_per_proc
,
numCounters
,
MPI_DOUBLE
,
MPI_MAX
,
0
,
world
);
MPI_Reduce
(
my_vals
,
min_per_proc
,
numCounters
,
MPI_DOUBLE
,
MPI_MIN
,
0
,
world
);
const
double
nODEs
=
sums
[
numCounters
-
1
];
double
avg_per_atom
[
numCounters
],
avg_per_proc
[
numCounters
];
// Averages per-ODE and per-proc per time-step.
for
(
int
i
=
0
;
i
<
numCounters
;
++
i
){
avg_per_atom
[
i
]
=
sums
[
i
]
/
nODEs
;
avg_per_proc
[
i
]
=
sums
[
i
]
/
comm
->
nprocs
;
}
// Sum up the differences from each task.
double
sum_sq
[
2
*
numCounters
];
double
my_sum_sq
[
2
*
numCounters
];
for
(
int
i
=
0
;
i
<
numCounters
;
++
i
){
double
diff_i
=
my_vals
[
i
]
-
avg_per_proc
[
i
];
my_sum_sq
[
i
]
=
diff_i
*
diff_i
;
}
double
max_per_ODE
[
numCounters
],
min_per_ODE
[
numCounters
];
// Process the per-ODE RMS of the # of steps/funcs
if
(
diagnosticFrequency
==
1
){
double
my_max
[
numCounters
],
my_min
[
numCounters
];
const
int
nlocal
=
atom
->
nlocal
;
const
int
*
mask
=
atom
->
mask
;
for
(
int
i
=
0
;
i
<
numCounters
;
++
i
){
my_sum_sq
[
i
+
numCounters
]
=
0
;
my_max
[
i
]
=
0
;
my_min
[
i
]
=
DBL_MAX
;
if
(
diagnosticCounterPerODE
[
i
]
!=
NULL
){
for
(
int
j
=
0
;
j
<
nlocal
;
++
j
)
if
(
mask
[
j
]
&
groupbit
){
double
diff
=
double
(
diagnosticCounterPerODE
[
i
][
j
])
-
avg_per_atom
[
i
];
my_sum_sq
[
i
+
numCounters
]
+=
diff
*
diff
;
my_max
[
i
]
=
std
::
max
(
my_max
[
i
],
(
double
)
diagnosticCounterPerODE
[
i
][
j
]
);
my_min
[
i
]
=
std
::
min
(
my_min
[
i
],
(
double
)
diagnosticCounterPerODE
[
i
][
j
]
);
}
}
}
MPI_Reduce
(
my_sum_sq
,
sum_sq
,
2
*
numCounters
,
MPI_DOUBLE
,
MPI_SUM
,
0
,
world
);
MPI_Reduce
(
my_max
,
max_per_ODE
,
numCounters
,
MPI_DOUBLE
,
MPI_MAX
,
0
,
world
);
MPI_Reduce
(
my_min
,
min_per_ODE
,
numCounters
,
MPI_DOUBLE
,
MPI_MIN
,
0
,
world
);
}
else
MPI_Reduce
(
my_sum_sq
,
sum_sq
,
numCounters
,
MPI_DOUBLE
,
MPI_SUM
,
0
,
world
);
TimerType
timer_stop
=
getTimeStamp
();
double
time_local
=
getElapsedTime
(
timer_start
,
timer_stop
);
if
(
comm
->
me
==
0
){
char
smesg
[
128
];
#define print_mesg(smesg) {\
if (screen) fprintf(screen,"%s\n", smesg); \
if (logfile) fprintf(logfile,"%s\n", smesg); }
sprintf
(
smesg
,
"FixRX::ODE Diagnostics: # of steps |# of rhs evals| run-time (sec)"
);
print_mesg
(
smesg
);
sprintf
(
smesg
,
" AVG per ODE : %-12.5g | %-12.5g | %-12.5g"
,
avg_per_atom
[
0
],
avg_per_atom
[
1
],
avg_per_atom
[
2
]);
print_mesg
(
smesg
);
// only valid for single time-step!
if
(
diagnosticFrequency
==
1
){
double
rms_per_ODE
[
numCounters
];
for
(
int
i
=
0
;
i
<
numCounters
;
++
i
)
rms_per_ODE
[
i
]
=
sqrt
(
sum_sq
[
i
+
numCounters
]
/
nODEs
);
sprintf
(
smesg
,
" RMS per ODE : %-12.5g | %-12.5g "
,
rms_per_ODE
[
0
],
rms_per_ODE
[
1
]);
print_mesg
(
smesg
);
sprintf
(
smesg
,
" MAX per ODE : %-12.5g | %-12.5g "
,
max_per_ODE
[
0
],
max_per_ODE
[
1
]);
print_mesg
(
smesg
);
sprintf
(
smesg
,
" MIN per ODE : %-12.5g | %-12.5g "
,
min_per_ODE
[
0
],
min_per_ODE
[
1
]);
print_mesg
(
smesg
);
}
sprintf
(
smesg
,
" AVG per Proc : %-12.5g | %-12.5g | %-12.5g"
,
avg_per_proc
[
0
],
avg_per_proc
[
1
],
avg_per_proc
[
2
]);
print_mesg
(
smesg
);
if
(
comm
->
nprocs
>
1
){
double
rms_per_proc
[
numCounters
];
for
(
int
i
=
0
;
i
<
numCounters
;
++
i
)
rms_per_proc
[
i
]
=
sqrt
(
sum_sq
[
i
]
/
comm
->
nprocs
);
sprintf
(
smesg
,
" RMS per Proc : %-12.5g | %-12.5g | %-12.5g"
,
rms_per_proc
[
0
],
rms_per_proc
[
1
],
rms_per_proc
[
2
]);
print_mesg
(
smesg
);
sprintf
(
smesg
,
" MAX per Proc : %-12.5g | %-12.5g | %-12.5g"
,
max_per_proc
[
0
],
max_per_proc
[
1
],
max_per_proc
[
2
]);
print_mesg
(
smesg
);
sprintf
(
smesg
,
" MIN per Proc : %-12.5g | %-12.5g | %-12.5g"
,
min_per_proc
[
0
],
min_per_proc
[
1
],
min_per_proc
[
2
]);
print_mesg
(
smesg
);
}
sprintf
(
smesg
,
" AVG'd over %d time-steps"
,
nTimes
);
print_mesg
(
smesg
);
sprintf
(
smesg
,
" AVG'ing took %g sec"
,
time_local
);
print_mesg
(
smesg
);
#undef print_mesg
}
// Reset the counters.
for
(
int
i
=
0
;
i
<
numDiagnosticCounters
;
++
i
)
diagnosticCounter
[
i
]
=
0
;
return
;
}
void
FixRX
::
rkf45
(
int
id
,
double
*
rwork
)
{
// Rounding coefficient.
const
double
uround
=
DBL_EPSILON
;
// Adaption limit (shrink or grow)
const
double
adaption_limit
=
4.0
;
//double *y = new double[8*nspecies + nreactions];
double
*
y
=
NULL
;
if
(
rwork
==
NULL
)
y
=
new
double
[
8
*
nspecies
+
nreactions
];
else
y
=
rwork
;
double
*
rhstmp
=
y
+
8
*
nspecies
;
const
int
neq
=
nspecies
;
// Update ConcOld and initialize the ODE solution vector y[].
for
(
int
ispecies
=
0
;
ispecies
<
nspecies
;
ispecies
++
){
const
double
tmp
=
atom
->
dvector
[
ispecies
][
id
];
atom
->
dvector
[
ispecies
+
nspecies
][
id
]
=
tmp
;
y
[
ispecies
]
=
tmp
;
}
// Integration length.
const
double
t_stop
=
update
->
dt
;
// DPD time-step.
// Safety factor on the adaption. very specific but not necessary .. 0.9 is common.
const
double
hsafe
=
0.840896415
;
// Time rounding factor.
const
double
tround
=
t_stop
*
uround
;
// Counters for diagnostics.
int
nst
=
0
;
// # of steps (accepted)
int
nit
=
0
;
// # of iterations total
int
nfe
=
0
;
// # of RHS evaluations
// Min/Max step-size limits.
const
double
h_min
=
100.0
*
tround
;
const
double
h_max
=
(
minSteps
>
0
)
?
t_stop
/
double
(
minSteps
)
:
t_stop
;
// Set the initial step-size. 0 forces an internal estimate ... stable Euler step size.
double
h
=
(
minSteps
>
0
)
?
t_stop
/
double
(
minSteps
)
:
0.0
;
double
t
=
0.0
;
if
(
h
<
h_min
){
//fprintf(stderr,"hin not implemented yet\n");
//exit(-1);
nfe
=
rkf45_h0
(
neq
,
t
,
t_stop
,
h_min
,
h_max
,
h
,
y
,
y
+
neq
,
rhstmp
);
}
//printf("t= %e t_stop= %e h= %e\n", t, t_stop, h);
// Integrate until we reach the end time.
while
(
fabs
(
t
-
t_stop
)
>
tround
){
double
*
yout
=
y
+
neq
;
double
*
eout
=
yout
+
neq
;
// Take a trial step.
rkf45_step
(
neq
,
h
,
y
,
yout
,
eout
,
rhstmp
);
// Estimate the solution error.
// ... weighted 2-norm of the error.
double
err2
=
0.0
;
for
(
int
k
=
0
;
k
<
neq
;
k
++
){
const
double
wterr
=
eout
[
k
]
/
(
relTol
*
fabs
(
y
[
k
]
)
+
absTol
);
err2
+=
wterr
*
wterr
;
}
double
err
=
fmax
(
uround
,
sqrt
(
err2
/
double
(
nspecies
)
));
// Accept the solution?
if
(
err
<=
1.0
||
h
<=
h_min
){
t
+=
h
;
nst
++
;
for
(
int
k
=
0
;
k
<
neq
;
k
++
)
y
[
k
]
=
yout
[
k
];
}
// Adjust h for the next step.
double
hfac
=
hsafe
*
sqrt
(
sqrt
(
1.0
/
err
)
);
// Limit the adaption.
hfac
=
fmax
(
hfac
,
1.0
/
adaption_limit
);
hfac
=
fmin
(
hfac
,
adaption_limit
);
// Apply the adaption factor...
h
*=
hfac
;
// Limit h.
h
=
fmin
(
h
,
h_max
);
h
=
fmax
(
h
,
h_min
);
// Stretch h if we're within 5% ... and we didn't just fail.
if
(
err
<=
1.0
&&
(
t
+
1.05
*
h
)
>
t_stop
)
h
=
t_stop
-
t
;
// And don't overshoot the end.
if
(
t
+
h
>
t_stop
)
h
=
t_stop
-
t
;
nit
++
;
nfe
+=
6
;
if
(
maxIters
&&
nit
>
maxIters
){
//fprintf(stderr,"atom[%d] took too many iterations in rkf45 %d %e %e\n", id, nit, t, t_stop);
nFails
++
;
break
;
// We should set an error here so that the solution is not used!
}
}
// end while
nSteps
+=
nst
;
nIters
+=
nit
;
nFuncs
+=
nfe
;
//if (diagnosticFrequency == 1 && diagnosticCounterPerODE[StepSum] != NULL)
if
(
diagnosticCounterPerODE
[
StepSum
]
!=
NULL
){
diagnosticCounterPerODE
[
StepSum
][
id
]
=
nst
;
diagnosticCounterPerODE
[
FuncSum
][
id
]
=
nfe
;
}
//printf("id= %d nst= %d nit= %d\n", id, nst, nit);
// Store the solution back in atom->dvector.
for
(
int
ispecies
=
0
;
ispecies
<
nspecies
;
ispecies
++
){
if
(
y
[
ispecies
]
<
-
1.0e-10
)
error
->
one
(
FLERR
,
"Computed concentration in RKF45 solver is < -1.0e-10"
);
else
if
(
y
[
ispecies
]
<
1e-20
)
y
[
ispecies
]
=
0.0
;
atom
->
dvector
[
ispecies
][
id
]
=
y
[
ispecies
];
}
if
(
rwork
==
NULL
)
delete
[]
y
;
}
/* ---------------------------------------------------------------------- */
int
FixRX
::
rhs
(
double
t
,
const
double
*
y
,
double
*
dydt
,
void
*
params
)
{
// Use the sparse format instead.
if
(
useSparseKinetics
)
return
this
->
rhs_sparse
(
t
,
y
,
dydt
,
params
);
else
return
this
->
rhs_dense
(
t
,
y
,
dydt
,
params
);
}
/* ---------------------------------------------------------------------- */
int
FixRX
::
rhs_dense
(
double
t
,
const
double
*
y
,
double
*
dydt
,
void
*
params
)
{
double
rxnRateLawForward
;
double
*
rxnRateLaw
=
(
double
*
)
params
;
double
VDPD
=
domain
->
xprd
*
domain
->
yprd
*
domain
->
zprd
/
atom
->
natoms
;
double
concentration
;
int
nspecies
=
atom
->
nspecies_dpd
;
for
(
int
ispecies
=
0
;
ispecies
<
nspecies
;
ispecies
++
)
dydt
[
ispecies
]
=
0.0
;
// Construct the reaction rate laws
for
(
int
jrxn
=
0
;
jrxn
<
nreactions
;
jrxn
++
){
rxnRateLawForward
=
kR
[
jrxn
];
for
(
int
ispecies
=
0
;
ispecies
<
nspecies
;
ispecies
++
){
concentration
=
y
[
ispecies
]
/
VDPD
;
rxnRateLawForward
*=
pow
(
concentration
,
stoichReactants
[
jrxn
][
ispecies
]);
}
rxnRateLaw
[
jrxn
]
=
rxnRateLawForward
;
}
// Construct the reaction rates for each species
for
(
int
ispecies
=
0
;
ispecies
<
nspecies
;
ispecies
++
)
for
(
int
jrxn
=
0
;
jrxn
<
nreactions
;
jrxn
++
)
dydt
[
ispecies
]
+=
stoich
[
jrxn
][
ispecies
]
*
VDPD
*
rxnRateLaw
[
jrxn
];
return
0
;
}
/* ---------------------------------------------------------------------- */
int
FixRX
::
rhs_sparse
(
double
t
,
const
double
*
y
,
double
*
dydt
,
void
*
v_params
)
const
{
double
*
_rxnRateLaw
=
(
double
*
)
v_params
;
const
double
VDPD
=
domain
->
xprd
*
domain
->
yprd
*
domain
->
zprd
/
atom
->
natoms
;
#define kFor (this->kR)
#define kRev (NULL)
#define rxnRateLaw (_rxnRateLaw)
#define conc (dydt)
#define maxReactants (this->sparseKinetics_maxReactants)
#define maxSpecies (this->sparseKinetics_maxSpecies)
#define nuk (this->sparseKinetics_nuk)
#define nu (this->sparseKinetics_nu)
#define inu (this->sparseKinetics_inu)
#define isIntegral(idx) (SparseKinetics_enableIntegralReactions \
&& this->sparseKinetics_isIntegralReaction[idx])
for
(
int
k
=
0
;
k
<
nspecies
;
++
k
)
conc
[
k
]
=
y
[
k
]
/
VDPD
;
// Construct the reaction rate laws
for
(
int
i
=
0
;
i
<
nreactions
;
++
i
)
{
double
rxnRateLawForward
;
if
(
isIntegral
(
i
)){
rxnRateLawForward
=
kFor
[
i
]
*
powint
(
conc
[
nuk
[
i
][
0
]
],
inu
[
i
][
0
]);
for
(
int
kk
=
1
;
kk
<
maxReactants
;
++
kk
){
const
int
k
=
nuk
[
i
][
kk
];
if
(
k
==
SparseKinetics_invalidIndex
)
break
;
//if (k != SparseKinetics_invalidIndex)
rxnRateLawForward
*=
powint
(
conc
[
k
],
inu
[
i
][
kk
]
);
}
}
else
{
rxnRateLawForward
=
kFor
[
i
]
*
pow
(
conc
[
nuk
[
i
][
0
]
],
nu
[
i
][
0
]);
for
(
int
kk
=
1
;
kk
<
maxReactants
;
++
kk
){
const
int
k
=
nuk
[
i
][
kk
];
if
(
k
==
SparseKinetics_invalidIndex
)
break
;
//if (k != SparseKinetics_invalidIndex)
rxnRateLawForward
*=
pow
(
conc
[
k
],
nu
[
i
][
kk
]
);
}
}
rxnRateLaw
[
i
]
=
rxnRateLawForward
;
}
// Construct the reaction rates for each species from the
// Stoichiometric matrix and ROP vector.
for
(
int
k
=
0
;
k
<
nspecies
;
++
k
)
dydt
[
k
]
=
0.0
;
for
(
int
i
=
0
;
i
<
nreactions
;
++
i
){
// Reactants ...
dydt
[
nuk
[
i
][
0
]
]
-=
nu
[
i
][
0
]
*
rxnRateLaw
[
i
];
for
(
int
kk
=
1
;
kk
<
maxReactants
;
++
kk
){
const
int
k
=
nuk
[
i
][
kk
];
if
(
k
==
SparseKinetics_invalidIndex
)
break
;
//if (k != SparseKinetics_invalidIndex)
dydt
[
k
]
-=
nu
[
i
][
kk
]
*
rxnRateLaw
[
i
];
}
// Products ...
dydt
[
nuk
[
i
][
maxReactants
]
]
+=
nu
[
i
][
maxReactants
]
*
rxnRateLaw
[
i
];
for
(
int
kk
=
maxReactants
+
1
;
kk
<
maxSpecies
;
++
kk
){
const
int
k
=
nuk
[
i
][
kk
];
if
(
k
==
SparseKinetics_invalidIndex
)
break
;
//if (k != SparseKinetics_invalidIndex)
dydt
[
k
]
+=
nu
[
i
][
kk
]
*
rxnRateLaw
[
i
];
}
}
// Add in the volume factor to convert to the proper units.
for
(
int
k
=
0
;
k
<
nspecies
;
++
k
)
dydt
[
k
]
*=
VDPD
;
#undef kFor
#undef kRev
#undef rxnRateLaw
#undef conc
#undef maxReactants
#undef maxSpecies
#undef nuk
#undef nu
#undef inu
#undef isIntegral
//#undef invalidIndex
return
0
;
}
/* ---------------------------------------------------------------------- */
void
FixRX
::
computeLocalTemperature
()
{
int
i
,
j
,
ii
,
jj
,
inum
,
jnum
,
itype
,
jtype
;
double
xtmp
,
ytmp
,
ztmp
,
delx
,
dely
,
delz
;
double
rsq
;
int
*
ilist
,
*
jlist
,
*
numneigh
,
**
firstneigh
;
double
**
x
=
atom
->
x
;
int
*
type
=
atom
->
type
;
int
nlocal
=
atom
->
nlocal
;
int
nghost
=
atom
->
nghost
;
int
newton_pair
=
force
->
newton_pair
;
// local temperature variables
double
wij
=
0.0
;
double
*
dpdTheta
=
atom
->
dpdTheta
;
// Initialize the local density and local temperature arrays
if
(
newton_pair
)
{
sumWeights
=
new
double
[
nlocal
+
nghost
];
for
(
ii
=
0
;
ii
<
nlocal
+
nghost
;
ii
++
)
sumWeights
[
ii
]
=
0.0
;
}
else
{
sumWeights
=
new
double
[
nlocal
];
for
(
ii
=
0
;
ii
<
nlocal
;
ii
++
)
dpdThetaLocal
[
ii
]
=
0.0
;
}
inum
=
pairDPDE
->
list
->
inum
;
ilist
=
pairDPDE
->
list
->
ilist
;
numneigh
=
pairDPDE
->
list
->
numneigh
;
firstneigh
=
pairDPDE
->
list
->
firstneigh
;
// loop over neighbors of my atoms
for
(
ii
=
0
;
ii
<
inum
;
ii
++
)
{
i
=
ilist
[
ii
];
xtmp
=
x
[
i
][
0
];
ytmp
=
x
[
i
][
1
];
ztmp
=
x
[
i
][
2
];
itype
=
type
[
i
];
jlist
=
firstneigh
[
i
];
jnum
=
numneigh
[
i
];
for
(
jj
=
0
;
jj
<
jnum
;
jj
++
)
{
j
=
jlist
[
jj
];
j
&=
NEIGHMASK
;
jtype
=
type
[
j
];
delx
=
xtmp
-
x
[
j
][
0
];
dely
=
ytmp
-
x
[
j
][
1
];
delz
=
ztmp
-
x
[
j
][
2
];
rsq
=
delx
*
delx
+
dely
*
dely
+
delz
*
delz
;
if
(
rsq
<
pairDPDE
->
cutsq
[
itype
][
jtype
])
{
double
rcut
=
sqrt
(
pairDPDE
->
cutsq
[
itype
][
jtype
]);
double
rij
=
sqrt
(
rsq
);
double
ratio
=
rij
/
rcut
;
// Lucy's Weight Function
if
(
wtFlag
==
LUCY
){
wij
=
(
1.0
+
3.0
*
ratio
)
*
(
1.0
-
ratio
)
*
(
1.0
-
ratio
)
*
(
1.0
-
ratio
);
dpdThetaLocal
[
i
]
+=
wij
/
dpdTheta
[
j
];
if
(
newton_pair
||
j
<
nlocal
)
dpdThetaLocal
[
j
]
+=
wij
/
dpdTheta
[
i
];
}
sumWeights
[
i
]
+=
wij
;
if
(
newton_pair
||
j
<
nlocal
)
sumWeights
[
j
]
+=
wij
;
}
}
}
if
(
newton_pair
)
comm
->
reverse_comm_fix
(
this
);
// self-interaction for local temperature
for
(
i
=
0
;
i
<
nlocal
;
i
++
){
// Lucy Weight Function
if
(
wtFlag
==
LUCY
){
wij
=
1.0
;
dpdThetaLocal
[
i
]
+=
wij
/
dpdTheta
[
i
];
}
sumWeights
[
i
]
+=
wij
;
// Normalized local temperature
dpdThetaLocal
[
i
]
=
dpdThetaLocal
[
i
]
/
sumWeights
[
i
];
if
(
localTempFlag
==
HARMONIC
)
dpdThetaLocal
[
i
]
=
1.0
/
dpdThetaLocal
[
i
];
}
delete
[]
sumWeights
;
}
/* ---------------------------------------------------------------------- */
int
FixRX
::
pack_forward_comm
(
int
n
,
int
*
list
,
double
*
buf
,
int
pbc_flag
,
int
*
pbc
)
{
int
ii
,
jj
,
m
;
double
tmp
;
m
=
0
;
for
(
ii
=
0
;
ii
<
n
;
ii
++
)
{
jj
=
list
[
ii
];
for
(
int
ispecies
=
0
;
ispecies
<
nspecies
;
ispecies
++
){
tmp
=
atom
->
dvector
[
ispecies
][
jj
];
buf
[
m
++
]
=
tmp
;
tmp
=
atom
->
dvector
[
ispecies
+
nspecies
][
jj
];
buf
[
m
++
]
=
tmp
;
}
}
return
m
;
}
/* ---------------------------------------------------------------------- */
void
FixRX
::
unpack_forward_comm
(
int
n
,
int
first
,
double
*
buf
)
{
int
ii
,
m
,
last
;
double
tmp
;
m
=
0
;
last
=
first
+
n
;
for
(
ii
=
first
;
ii
<
last
;
ii
++
){
for
(
int
ispecies
=
0
;
ispecies
<
nspecies
;
ispecies
++
){
tmp
=
buf
[
m
++
];
atom
->
dvector
[
ispecies
][
ii
]
=
tmp
;
tmp
=
buf
[
m
++
];
atom
->
dvector
[
ispecies
+
nspecies
][
ii
]
=
tmp
;
}
}
}
/* ---------------------------------------------------------------------- */
int
FixRX
::
pack_reverse_comm
(
int
n
,
int
first
,
double
*
buf
)
{
int
i
,
m
,
last
;
m
=
0
;
last
=
first
+
n
;
for
(
i
=
first
;
i
<
last
;
i
++
)
{
buf
[
m
++
]
=
dpdThetaLocal
[
i
];
buf
[
m
++
]
=
sumWeights
[
i
];
}
return
m
;
}
/* ---------------------------------------------------------------------- */
void
FixRX
::
unpack_reverse_comm
(
int
n
,
int
*
list
,
double
*
buf
)
{
int
i
,
j
,
m
;
m
=
0
;
for
(
i
=
0
;
i
<
n
;
i
++
)
{
j
=
list
[
i
];
dpdThetaLocal
[
j
]
+=
buf
[
m
++
];
sumWeights
[
j
]
+=
buf
[
m
++
];
}
}
Event Timeline
Log In to Comment