Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
Marco De Lucia
RedModRphree
Commits
57e6d53e
Commit
57e6d53e
authored
May 06, 2018
by
Marco De Lucia
Browse files
preparing demos and included files (db, dataset, ...)
parent
131808cd
Changes
5
Hide whitespace changes
Inline
Side-by-side
demo/00Index
0 → 100644
View file @
57e6d53e
demo-equilibrium 1D reactive transport with equilibrium chemistry
demo-kinetics 1D reactive transport with kinetic chemistry
demo/demo-equilibrium.R
0 → 100644
View file @
57e6d53e
### Licence: LGPL version 2.1
### Time-stamp: "Last modified 2016-04-15 17:17:41 delucia"
library
(
RedModRphree
)
Matplot
<-
function
(
sample
,
sur
,
sim
,
head
=
FALSE
,
...
)
{
rs
<-
sur
[[
sample
]]
$
C
rn
<-
sim
[[
sample
]]
$
C
if
(
head
){
cat
(
"Surrogate:\n"
)
print
(
head
(
rs
))
cat
(
"Simulation:\n"
)
print
(
head
(
rn
))
}
matplot
(
rn
[,
c
(
"Ca"
,
"Mg"
,
"Calcite"
,
"Dolomite"
)],
type
=
"l"
,
lwd
=
2
,
lty
=
"solid"
,
col
=
c
(
"red"
,
"black"
,
"blue"
,
"olivedrab3"
),
...
)
matplot
(
rs
[,
c
(
"Ca"
,
"Mg"
,
"Calcite"
,
"Dolomite"
)],
type
=
"l"
,
lwd
=
2
,
lty
=
"dashed"
,
col
=
c
(
"red"
,
"black"
,
"blue"
,
"olivedrab3"
),
add
=
TRUE
)
legend
(
"topright"
,
c
(
"Ca"
,
"Mg"
,
"Calcite"
,
"Dolomite"
),
text.col
=
c
(
"red"
,
"black"
,
"blue"
,
"olivedrab3"
),
bty
=
"n"
)
}
base
<-
c
(
"SOLUTION 1"
,
"units mol/kgw"
,
"temp 25.0"
,
"water 1"
,
"pH 9.91 charge"
,
"pe 4.0"
,
"C 1.2279E-04"
,
"Ca 1.2279E-04"
,
"Mg 0.001"
,
"Cl 0.002"
,
"PURE 1"
,
"O2g -0.1675 10"
,
"KINETICS 1"
,
"-steps 100"
,
"-bad_step_max 1000"
,
"Calcite"
,
"-m 0.000207"
,
"-parms 3.20"
,
"Dolomite"
,
"-m 0.0"
,
"-parms 0.32"
,
"END"
)
initsim
<-
c
(
"SELECTED_OUTPUT"
,
"-high_precision true"
,
"-reset false"
,
"USER_PUNCH"
,
"-head C Ca Cl Mg H e O2g Calcite Dolomite"
,
"10 PUNCH TOT(\"C\"), TOT(\"Ca\"), TOT(\"Cl\"), TOT(\"Mg\"), ACT(\"H+\"), ACT(\"e-\"), EQUI(\"O2g\"), EQUI(\"Calcite\"), EQUI(\"Dolomite\")"
,
"SOLUTION 1"
,
"units mol/kgw"
,
"temp 25.0"
,
"water 1"
,
"pH 9.91 charge"
,
"pe 4.0"
,
"C 1.2279E-04"
,
"Ca 1.2279E-04"
,
"Cl 1E-12"
,
"Mg 1E-12"
,
"PURE 1"
,
"O2g -0.6788 10.0"
,
"Calcite 0.0 2.07E-3"
,
"Dolomite 0.0 0.0"
,
"END"
)
db
<-
RPhreeFile
(
"db/pqc_strip_kin.dat"
,
is.db
=
TRUE
)
phrLoadDatabaseString
(
db
)
selout
<-
c
(
"SELECTED_OUTPUT"
,
"-high_precision true"
,
"-reset false"
,
"-time"
,
"-soln"
,
"-temperature true"
,
"-water true"
,
"-pH"
,
"-pe"
,
"-totals C Ca Cl Mg"
,
"-kinetic_reactants Calcite Dolomite"
,
"-equilibrium O2g"
)
prop
<-
c
(
"C"
,
"Ca"
,
"Cl"
,
"Mg"
,
"pH"
,
"pe"
,
"O2g"
,
"Calcite"
,
"Dolomite"
)
maxsim
<-
20
## minimum number of single chemical simulations in each time step to spawn parallelization
grid
<-
50
## grid discretization (number of elements)
L
<-
0.5
## total system length [m]
U
<-
9.375e-6
## imposed (constant) Darcy velocity [m3/s]
Cf
<-
0.9
## Courant number to add some numerical dispersion
## boundary conditions (concentration of the injected fluid entering into the 1D column)
bound
<-
c
(
H
=
1E-7
,
e
=
1E-4
,
Mg
=
0.001
,
Cl
=
0.002
)
setup
<-
list
(
n
=
50
,
len
=
L
,
U
=
U
,
base
=
base
,
first
=
selout
,
Cf
=
Cf
,
bound
=
bound
,
dt
=
210
,
prop
=
prop
,
immobile
=
c
(
7
,
8
,
9
),
kin
=
c
(
8
,
9
),
ann
=
list
(
O2g
=
-0.1675
),
initsim
=
initsim
)
oo
<-
ReactTranspBalanceKin
(
setup
=
setup
,
db
=
db
,
step
=
"fix_dt"
,
maxtime
=
2100
,
procs
=
4
,
maxsim
=
maxsim
,
writeout
=
FALSE
,
reduce
=
TRUE
,
surrogate
=
FALSE
)
of
<-
ReactTranspBalanceKin
(
setup
=
setup
,
db
=
db
,
step
=
"fix_dt"
,
maxtime
=
2100
,
procs
=
4
,
maxsim
=
maxsim
,
writeout
=
FALSE
,
reduce
=
FALSE
,
surrogate
=
FALSE
)
Matplot
(
101
,
sur
=
oo
,
sim
=
of
,
xlab
=
"Domain element"
,
ylab
=
"Ca, Mg [molal], Calcite, Dolomite [mol]"
,
main
=
"RandomForest surrogate"
)
load
(
file
=
"Rdata/20180411_CalcKin_Surr_EGU.Rdata"
)
MySurr
<-
function
(
state
,
model
){
require
(
caret
)
excl
<-
which
(
colnames
(
state
)
==
"O2g"
)
remember_state_names
<-
colnames
(
state
)
des
<-
state
[,
-
excl
]
if
(
NA
%in%
des
)
{
print
(
des
)
stop
(
"Mysurr: NA in data"
)
}
if
(
NaN
%in%
des
)
{
print
(
des
)
stop
(
"Mysurr: NaN in data"
)
}
order_for_surrogate
<-
c
(
"C"
,
"Ca"
,
"Cl"
,
"Mg"
,
"pH"
,
"pe"
,
"Calcite"
,
"Dolomite"
)
rem_attr
<-
attr
(
state
,
"immobile"
)
ss
<-
des
[,
order_for_surrogate
]
pred
<-
lapply
(
names
(
model
),
function
(
x
)
predict.train
(
model
[[
x
]],
ss
))
out
<-
cbind
(
as.data.frame
(
pred
),
state
[,
excl
])
colnames
(
out
)
<-
c
(
names
(
model
),
"O2g"
)
## order back to how it is expected in "state"
out
<-
out
[,
remember_state_names
]
attr
(
out
,
"immobile"
)
<-
rem_attr
return
(
out
)
}
so
<-
ReactTranspBalanceKin
(
setup
=
setup
,
db
=
db
,
step
=
"fix_dt"
,
maxtime
=
21000
,
procs
=
4
,
maxsim
=
maxsim
,
writeout
=
FALSE
,
reduce
=
TRUE
,
surrogate
=
TRUE
,
surrogate.FUN
=
MySurr
,
model
=
mods
,
baleq
=
baleq
,
tol
=
1E-5
,
call_pqc
=
TRUE
)
sf
<-
ReactTranspBalanceKin
(
setup
=
setup
,
db
=
db
,
step
=
"fix_dt"
,
maxtime
=
21000
,
procs
=
4
,
maxsim
=
maxsim
,
writeout
=
FALSE
,
reduce
=
FALSE
,
surrogate
=
TRUE
,
surrogate.FUN
=
MySurr
,
model
=
mods
,
baleq
=
baleq
,
tol
=
1E-5
,
call_pqc
=
TRUE
)
par
(
mfrow
=
c
(
2
,
1
))
Matplot
(
101
,
sur
=
so
,
sim
=
oo
,
xlab
=
"Domain element"
,
ylab
=
"Ca, Mg [molal], Calcite, Dolomite [mol]"
,
main
=
"RandomForest surrogate"
)
Matplot
(
101
,
sur
=
sf
,
sim
=
of
,
xlab
=
"Domain element"
,
ylab
=
"Ca, Mg [molal], Calcite, Dolomite [mol]"
,
main
=
"RandomForest surrogate"
)
##### test equilibrium
prop
<-
c
(
"C(4)"
,
"Ca"
,
"Cl"
,
"Mg"
,
"pH"
,
"pe"
,
"Calcite"
,
"Dolomite"
)
## PHREEQC mojo to output selection
sel
<-
RPhreeCheckSel
(
c
(
punch
=
TRUE
))
selec
<-
c
(
"SELECTED_OUTPUT"
,
"-high_precision true"
,
"-reset false"
,
"USER_PUNCH"
,
"-head C Ca Cl Mg pH pe Calcite Dolomite"
,
"10 PUNCH TOT(\"C(4)\"), TOT(\"Ca\"), TOT(\"Cl\"), TOT(\"Mg\"), -LA(\"H+\"), -LA(\"e-\"), EQUI(\"Calcite\"), EQUI(\"Dolomite\")"
)
## load the initial phreeqc script, which also defines initial conditions
geoinp
<-
RPhreeFile
(
"Surrogates/calcite_pqc2.pqc"
)
db2
<-
RPhreeFile
(
"Peter/phreeqc.dat"
,
is.db
=
TRUE
)
phrLoadDatabaseString
(
db2
)
maxsim
<-
50
## minimum number of single chemical simulations in each time step to spawn parallelization
grid
<-
100
## grid discretization (number of elements)
L
<-
0.5
## total system length [m]
U
<-
9.375e-6
## imposed (constant) Darcy velocity [m3/s]
Cf
<-
0.8
## Courant number to add some numerical dispersion
## boundary conditions (concentration of the injected fluid entering into the 1D column)
bound
<-
c
(
H
=
1E-7
,
e
=
1E-4
,
Mg
=
0.001
,
Cl
=
0.002
)
#, "O(0)"=0.04,
grid
<-
50
## grid discretization (number of elements)
setup
<-
list
(
n
=
grid
,
len
=
L
,
U
=
U
,
base
=
geoinp
,
first
=
selec
,
Cf
=
Cf
,
bound
=
bound
,
prop
=
prop
,
immobile
=
c
(
7
,
8
),
initsim
=
c
(
selec
,
geoinp
))
source
(
"RedModRphree/R/Rphree_ReactTrans.R"
)
source
(
"RedModRphree/R/Rphree_Kinetics.R"
)
eo
<-
ReactTranspBalanceEq
(
setup
=
setup
,
db
=
db2
,
step
=
"time"
,
maxtime
=
21000
,
procs
=
4
,
maxsim
=
maxsim
,
writeout
=
FALSE
,
reduce
=
TRUE
,
surrogate
=
FALSE
)
ef
<-
ReactTranspBalanceEq
(
setup
=
setup
,
db
=
db2
,
step
=
"time"
,
maxtime
=
21000
,
procs
=
4
,
maxsim
=
maxsim
,
writeout
=
FALSE
,
reduce
=
FALSE
,
surrogate
=
FALSE
)
load
(
file
=
"Rdata/20180314_calcite_pqc_4.Rdata"
)
load
(
file
=
"Surrogates/20180316_calcite_pqc_parRF.Rdata"
)
MySurr
<-
function
(
state
,
model
){
order_for_surrogate
<-
c
(
"C"
,
"Ca"
,
"Cl"
,
"Mg"
,
"Calcite"
,
"Dolomite"
,
"pH"
,
"pe"
)
remember_names
<-
colnames
(
state
)
rem_attr
<-
attr
(
state
,
"immobile"
)
ss
<-
state
[,
order_for_surrogate
]
pred
<-
lapply
(
names
(
model
),
function
(
x
)
as.numeric
(
predict
(
model
[[
x
]],
ss
)))
names
(
pred
)
<-
names
(
model
)
out
<-
data.matrix
(
as.data.frame
(
pred
))
## order back to how it is expected in "state"
out
<-
out
[,
remember_names
]
attr
(
out
,
"immobile"
)
<-
rem_attr
return
(
out
)
}
ro
<-
ReactTranspBalanceEq
(
setup
=
setup
,
db
=
db2
,
step
=
"time"
,
maxtime
=
21000
,
procs
=
4
,
maxsim
=
maxsim
,
writeout
=
FALSE
,
reduce
=
TRUE
,
surrogate
=
TRUE
,
surrogate.FUN
=
MySurr
,
model
=
mods
,
baleq
=
baleq
,
tol
=
1E-5
,
call_pqc
=
TRUE
)
rf
<-
ReactTranspBalanceEq
(
setup
=
setup
,
db
=
db2
,
step
=
"time"
,
maxtime
=
21000
,
procs
=
4
,
maxsim
=
maxsim
,
writeout
=
FALSE
,
reduce
=
FALSE
,
surrogate
=
TRUE
,
surrogate.FUN
=
MySurr
,
model
=
mods
,
baleq
=
baleq
,
tol
=
1E-5
,
call_pqc
=
TRUE
)
Matplot
(
25
,
sur
=
otr
,
sim
=
sr
,
xlab
=
"Domain element"
,
ylab
=
"Ca, Mg [molal], Calcite, Dolomite [mol]"
,
main
=
"RandomForest surrogate"
)
s1
<-
AdvectionPQC
(
state_T
,
inflow
=
bound
,
dx
=
L
/
50
,
dt
=
853.3333
)
prop
<-
c
(
"C"
,
"Ca"
,
"Cl"
,
"Mg"
,
"pH"
,
"pe"
,
"O2g"
,
"Calcite"
,
"Dolomite"
)
distprop
<-
c
(
"C"
,
"Ca"
,
"Cl"
,
"Mg"
,
"pH"
,
"pe"
,
"O2g"
,
"Calcite"
,
"Dolomite"
)
maxsim
<-
20
## minimum number of single chemical simulations in each time step to spawn parallelization
grid
<-
50
## grid discretization (number of elements)
L
<-
0.5
## total system length [m]
U
<-
9.375e-6
## imposed (constant) Darcy velocity [m3/s]
Cf
<-
0.9
## Courant number to add some numerical dispersion
## boundary conditions (concentration of the injected fluid entering into the 1D column)
bound
<-
c
(
H
=
1E-7
,
e
=
1E-4
,
Mg
=
0.001
,
Cl
=
0.002
)
#, "O(0)"=0.04,
exper
<-
expand.grid
(
grid
=
c
(
50
,
100
,
200
),
tol
=
c
(
1E-6
,
1E-5
,
1E-4
))
simsur
<-
vector
(
mode
=
"list"
,
length
=
nrow
(
exper
))
timesur
<-
vector
(
mode
=
"list"
,
length
=
nrow
(
exper
))
outbalsur
<-
vector
(
mode
=
"list"
,
length
=
nrow
(
exper
))
source
(
"Rphree_SurrogateKin.R"
)
for
(
i
in
seq
(
nrow
(
exper
)))
{
setup
<-
list
(
n
=
exper
$
grid
[
i
],
len
=
L
,
U
=
U
,
base
=
base
,
first
=
selout
,
Cf
=
Cf
,
bound
=
bound
,
dt
=
500
,
prop
=
prop
,
distprop
=
distprop
,
immobile
=
c
(
7
,
8
,
9
),
kin
=
c
(
8
,
9
),
sel
=
sel
,
ann
=
list
(
O2g
=
-0.1675
),
initsim
=
initsim
,
init_npunch
=
9
)
timesur
[[
i
]]
<-
system.time
(
simsur
[[
i
]]
<-
ReactTranspCompleteOutputBalKin
(
setup
=
setup
,
db
=
db
,
step
=
"fix_dt"
,
maxtime
=
2100
,
procs
=
4
,
maxsim
=
maxsim
,
writeout
=
FALSE
,
reduce
=
TRUE
,
surrogate
=
TRUE
,
surrogate.FUN
=
MySurr
,
model
=
mods
,
baleq
=
baleq
,
tol
=
exper
$
tol
[
i
],
call_pqc
=
TRUE
)
)
outbalsur
[[
i
]]
<-
out_bal
}
save
(
list
=
c
(
"des"
,
"res"
,
"samp"
,
"design"
,
"result"
,
"timesur"
,
"outbalsur"
,
"simsur"
,
"mods"
,
"MySurr"
,
"baleq"
),
file
=
"Rdata/20180411_CalcKin_Surr_EGU.Rdata"
)
load
(
file
=
"Rdata/20180411_CalcKin_Surr_EGU.Rdata"
)
data.table
::
fwrite
(
design
,
file
=
"TestCalcitePQC_Kinetic_design.csv"
,
append
=
FALSE
,
quote
=
"auto"
,
sep
=
","
,
eol
=
"\n"
,
na
=
"NA"
,
dec
=
"."
,
row.names
=
FALSE
,
col.names
=
TRUE
)
data.table
::
fwrite
(
result
,
file
=
"TestCalcitePQC_Kinetic_result.csv"
,
append
=
FALSE
,
quote
=
"auto"
,
sep
=
","
,
eol
=
"\n"
,
na
=
"NA"
,
dec
=
"."
,
row.names
=
FALSE
,
col.names
=
TRUE
)
## data.table::fwrite(design, file = "TestCalcitePQC_3_design.csv", append = FALSE, quote = "auto",
## sep = ",", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = TRUE)
data.table
::
fwrite
(
result
,
file
=
"TestCalcitePQC_3_results.csv"
,
append
=
FALSE
,
quote
=
"auto"
,
sep
=
","
,
eol
=
"\n"
,
na
=
"NA"
,
dec
=
"."
,
row.names
=
FALSE
,
col.names
=
TRUE
)
subexp1
<-
subset
(
exper
,
grid
==
50
|
grid
==
100
|
grid
==
200
)
toto
<-
rep
(
c
(
1
,
2
,
3
),
3
)
par
(
mfrow
=
c
(
3
,
3
))
for
(
i
in
seq
(
nrow
(
subexp1
)))
{
Matplot
(
11
,
sur
=
simsur
[[
i
]],
sim
=
expsim
[[
toto
[
i
]
]],
xlab
=
"Domain element"
,
ylab
=
"Ca, Mg [molal], Calcite, Dolomite [mol]"
,
main
=
paste
(
"RandomForest surrogate, tol="
,
subexp1
$
tol
[
i
]))
}
###################################### for debug and development
load
(
file
=
"Rdata/20180411_CalcKin_Surr_EGU.Rdata"
)
base
<-
RPhreeFile
(
"Surrogates/stub_kin.pqc"
)
initsim
<-
RPhreeFile
(
"Surrogates/stub_kin_init.pqc"
)
sel
<-
RPhreeCheckSel
(
c
(
punch
=
TRUE
))
db
<-
RPhreeFile
(
"db/pqc_strip_kin.dat"
,
is.db
=
TRUE
)
selout
<-
c
(
"SELECTED_OUTPUT"
,
"-high_precision true"
,
"-reset false"
,
"-time"
,
"-soln"
,
"-temperature true"
,
"-water true"
,
"-pH"
,
"-pe"
,
"-totals C Ca Cl Mg"
,
"-kinetic_reactants Calcite Dolomite"
,
"-equilibrium O2g"
)
prop
<-
c
(
"C"
,
"Ca"
,
"Cl"
,
"Mg"
,
"pH"
,
"pe"
,
"O2g"
,
"Calcite"
,
"Dolomite"
)
distprop
<-
c
(
"C"
,
"Ca"
,
"Cl"
,
"Mg"
,
"pH"
,
"pe"
,
"O2g"
,
"Calcite"
,
"Dolomite"
)
maxsim
<-
10
## minimum number of single chemical simulations in each time step to spawn parallelization
grid
<-
50
## grid discretization (number of elements)
L
<-
0.5
## total system length [m]
U
<-
9.375e-6
## imposed (constant) Darcy velocity [m3/s]
Cf
<-
0.9
## Courant number to add some numerical dispersion
## boundary conditions (concentration of the injected fluid entering into the 1D column)
bound
<-
c
(
H
=
1E-7
,
e
=
1E-4
,
Mg
=
0.001
,
Cl
=
0.002
)
#, "O(0)"=0.04,
source
(
"Rphree_SurrogateKin.R"
)
setup
<-
list
(
n
=
grid
,
len
=
L
,
U
=
U
,
base
=
base
,
first
=
selout
,
Cf
=
Cf
,
bound
=
bound
,
prop
=
prop
,
distprop
=
distprop
,
immobile
=
c
(
7
,
8
,
9
),
kin
=
c
(
8
,
9
),
sel
=
sel
,
ann
=
list
(
O2g
=
-0.1675
),
initsim
=
initsim
,
init_npunch
=
9
)
library
(
profvis
)
profvis
({
gc
()
tracemem
(
Res_C
)
sim
<-
ReactTranspCompleteOutputBalKin
(
setup
=
setup
,
db
=
db
,
step
=
"time"
,
maxtime
=
21000
,
procs
=
1
,
maxsim
=
maxsim
,
writeout
=
FALSE
,
reduce
=
FALSE
,
surrogate
=
FALSE
)
})
library
(
caret
)
MySurr
<-
function
(
state
,
model
){
excl
<-
which
(
colnames
(
state
)
==
"O2g"
)
remember_state_names
<-
colnames
(
state
)
des
<-
state
[,
-
excl
]
if
(
NA
%in%
des
)
{
print
(
des
)
stop
(
"Mysurr: NA in data"
)
}
if
(
NaN
%in%
des
)
{
print
(
des
)
stop
(
"Mysurr: NaN in data"
)
}
order_for_surrogate
<-
c
(
"C"
,
"Ca"
,
"Cl"
,
"Mg"
,
"pH"
,
"pe"
,
"Calcite"
,
"Dolomite"
)
rem_attr
<-
attr
(
state
,
"immobile"
)
ss
<-
des
[,
order_for_surrogate
]
pred
<-
parallel
::
mclapply
(
names
(
model
),
function
(
x
)
predict.train
(
model
[[
x
]],
ss
),
mc.cores
=
4
)
out
<-
cbind
(
as.data.frame
(
pred
),
state
[,
excl
])
colnames
(
out
)
<-
c
(
names
(
model
),
"O2g"
)
out
<-
out
[,
remember_state_names
]
attr
(
out
,
"immobile"
)
<-
rem_attr
return
(
out
)
}
setup
<-
list
(
n
=
20
,
len
=
L
,
U
=
U
,
base
=
base
,
first
=
selout
,
Cf
=
Cf
,
bound
=
bound
,
prop
=
prop
,
distprop
=
distprop
,
immobile
=
c
(
7
,
8
,
9
),
kin
=
c
(
8
,
9
),
sel
=
sel
,
ann
=
list
(
O2g
=
-0.1675
),
initsim
=
initsim
,
init_npunch
=
9
)
Rprof
(
tmp
<-
tempfile
())
##, memory.profiling = TRUE, gc.profiling = TRUE)
sim
<-
ReactTranspCompleteOutputBalKin
(
setup
=
setup
,
db
=
db
,
step
=
"time"
,
maxtime
=
4200
,
procs
=
1
,
maxsim
=
maxsim
,
writeout
=
FALSE
,
tol
=
1E-5
,
reduce
=
FALSE
,
surrogate
=
TRUE
,
surrogate.FUN
=
MySurr
,
model
=
mods
,
baleq
=
baleq
,
call_pqc
=
TRUE
)
Rprof
()
proftools
::
plotProfileCallGraph
(
readProfileData
(
tmp
),
score
=
"total"
)
proftools
::
plotProfileCallGraph
(
readProfileData
(
tmp
),
style
=
plain.style
,
score
=
"total"
)
library
(
profvis
)
profvis
({
sim
<-
ReactTranspCompleteOutputBalKin
(
setup
=
setup
,
db
=
db
,
step
=
"time"
,
maxtime
=
4200
,
procs
=
1
,
maxsim
=
maxsim
,
writeout
=
FALSE
,
tol
=
1E-5
,
reduce
=
FALSE
,
surrogate
=
TRUE
,
surrogate.FUN
=
MySurr
,
model
=
mods
,
baleq
=
baleq
,
call_pqc
=
TRUE
)
})
demo/demo-kinetics.R
0 → 100644
View file @
57e6d53e
### Licence: LGPL version 2.1
### Time-stamp: "Last modified 2016-04-15 17:17:41 delucia"
library
(
RedModRphree
)
db
<-
RPhreeFile
(
system.file
(
"extdata"
,
"phreeqc_kin.dat"
,
package
=
"RedModRphree"
),
is.db
=
TRUE
)
Matplot
<-
function
(
sample
,
sur
,
sim
,
head
=
FALSE
,
...
)
{
rs
<-
sur
[[
sample
]]
$
C
rn
<-
sim
[[
sample
]]
$
C
if
(
head
){
cat
(
"Surrogate:\n"
)
print
(
head
(
rs
))
cat
(
"Simulation:\n"
)
print
(
head
(
rn
))
}
matplot
(
rn
[,
c
(
"Ca"
,
"Mg"
,
"Calcite"
,
"Dolomite"
)],
type
=
"l"
,
lwd
=
2
,
lty
=
"solid"
,
col
=
c
(
"red"
,
"black"
,
"blue"
,
"olivedrab3"
),
...
)
matplot
(
rs
[,
c
(
"Ca"
,
"Mg"
,
"Calcite"
,
"Dolomite"
)],
type
=
"l"
,
lwd
=
2
,
lty
=
"dashed"
,
col
=
c
(
"red"
,
"black"
,
"blue"
,
"olivedrab3"
),
add
=
TRUE
)
legend
(
"topright"
,
c
(
"Ca"
,
"Mg"
,
"Calcite"
,
"Dolomite"
),
text.col
=
c
(
"red"
,
"black"
,
"blue"
,
"olivedrab3"
),
bty
=
"n"
)
}
base
<-
c
(
"SOLUTION 1"
,
"units mol/kgw"
,
"temp 25.0"
,
"water 1"
,
"pH 9.91 charge"
,
"pe 4.0"
,
"C 1.2279E-04"
,
"Ca 1.2279E-04"
,
"Mg 0.001"
,
"Cl 0.002"
,
"PURE 1"
,
"O2g -0.1675 10"
,
"KINETICS 1"
,
"-steps 100"
,
"-bad_step_max 1000"
,
"Calcite"
,
"-m 0.000207"
,
"-parms 3.20"
,
"Dolomite"
,
"-m 0.0"
,
"-parms 0.32"
,
"END"
)
initsim
<-
c
(
"SELECTED_OUTPUT"
,
"-high_precision true"
,
"-reset false"
,
"USER_PUNCH"
,
"-head C Ca Cl Mg H e O2g Calcite Dolomite"
,
"10 PUNCH TOT(\"C\"), TOT(\"Ca\"), TOT(\"Cl\"), TOT(\"Mg\"), ACT(\"H+\"), ACT(\"e-\"), EQUI(\"O2g\"), EQUI(\"Calcite\"), EQUI(\"Dolomite\")"
,
"SOLUTION 1"
,
"units mol/kgw"
,
"temp 25.0"
,
"water 1"
,
"pH 9.91 charge"
,
"pe 4.0"
,
"C 1.2279E-04"
,
"Ca 1.2279E-04"
,
"Cl 1E-12"
,
"Mg 1E-12"
,
"PURE 1"
,
"O2g -0.6788 10.0"
,
"Calcite 0.0 2.07E-3"
,
"Dolomite 0.0 0.0"
,
"END"
)
db
<-
RPhreeFile
(
"db/pqc_strip_kin.dat"
,
is.db
=
TRUE
)
phrLoadDatabaseString
(
db
)
selout
<-
c
(
"SELECTED_OUTPUT"
,
"-high_precision true"
,
"-reset false"
,
"-time"
,
"-soln"
,
"-temperature true"
,
"-water true"
,
"-pH"
,
"-pe"
,
"-totals C Ca Cl Mg"
,
"-kinetic_reactants Calcite Dolomite"
,
"-equilibrium O2g"
)
prop
<-
c
(
"C"
,
"Ca"
,
"Cl"
,
"Mg"
,
"pH"
,
"pe"
,
"O2g"
,
"Calcite"
,
"Dolomite"
)
maxsim
<-
20
## minimum number of single chemical simulations in each time step to spawn parallelization
grid
<-
50
## grid discretization (number of elements)
L
<-
0.5
## total system length [m]
U
<-
9.375e-6
## imposed (constant) Darcy velocity [m3/s]
Cf
<-
0.9
## Courant number to add some numerical dispersion
## boundary conditions (concentration of the injected fluid entering into the 1D column)
bound
<-
c
(
H
=
1E-7
,
e
=
1E-4
,
Mg
=
0.001
,
Cl
=
0.002
)
## pack all together in a "setup" list
setup
<-
list
(
n
=
50
,
len
=
L
,
U
=
U
,
base
=
base
,
first
=
selout
,
Cf
=
Cf
,
bound
=
bound
,
dt
=
210
,
prop
=
prop
,
immobile
=
c
(
7
,
8
,
9
),
kin
=
c
(
8
,
9
),
ann
=
list
(
O2g
=
-0.1675
),
initsim
=
initsim
)
oo
<-
ReactTranspBalanceKin
(
setup
=
setup
,
db
=
db
,
step
=
"fix_dt"
,
maxtime
=
2100
,
procs
=
4
,
maxsim
=
maxsim
,
writeout
=
FALSE
,
reduce
=
TRUE
,
surrogate
=
FALSE
)
of
<-
ReactTranspBalanceKin
(
setup
=
setup
,
db
=
db
,
step
=
"fix_dt"
,
maxtime
=
2100
,
procs
=
4
,
maxsim
=
maxsim
,
writeout
=
FALSE
,
reduce
=
FALSE
,
surrogate
=
FALSE
)
Matplot
(
101
,
sur
=
oo
,
sim
=
of
,
xlab
=
"Domain element"
,
ylab
=
"Ca, Mg [molal], Calcite, Dolomite [mol]"
,
main
=
"RandomForest surrogate"
)
load
(
file
=
"Rdata/20180411_CalcKin_Surr_EGU.Rdata"
)
MySurr
<-
function
(
state
,
model
){
require
(
caret
)
excl
<-
which
(
colnames
(
state
)
==
"O2g"
)
remember_state_names
<-
colnames
(
state
)
des
<-
state
[,
-
excl
]
if
(
NA
%in%
des
)
{
print
(
des
)
stop
(
"Mysurr: NA in data"
)
}
if
(
NaN
%in%
des
)
{
print
(
des
)
stop
(
"Mysurr: NaN in data"
)
}
order_for_surrogate
<-
c
(
"C"
,
"Ca"
,
"Cl"
,
"Mg"
,
"pH"
,
"pe"
,
"Calcite"
,
"Dolomite"
)
rem_attr
<-
attr
(
state
,
"immobile"
)
ss
<-
des
[,
order_for_surrogate
]
pred
<-
lapply
(
names
(
model
),
function
(
x
)
predict.train
(
model
[[
x
]],
ss
))
out
<-
cbind
(
as.data.frame
(
pred
),
state
[,
excl
])
colnames
(
out
)
<-
c
(
names
(
model
),
"O2g"
)
## order back to how it is expected in "state"
out
<-
out
[,
remember_state_names
]
attr
(
out
,
"immobile"
)
<-
rem_attr
return
(
out
)
}
so
<-
ReactTranspBalanceKin
(
setup
=
setup
,
db
=
db
,
step
=
"fix_dt"
,
maxtime
=
21000
,
procs
=
4
,
maxsim
=
maxsim
,
writeout
=
FALSE
,
reduce
=
TRUE
,
surrogate
=
TRUE
,
surrogate.FUN
=
MySurr
,
model
=
mods
,
baleq
=
baleq
,
tol
=
1E-5
,
call_pqc
=
TRUE
)
sf
<-
ReactTranspBalanceKin
(
setup
=
setup
,
db
=
db
,
step
=
"fix_dt"
,
maxtime
=
21000
,
procs
=
4
,
maxsim
=
maxsim
,
writeout
=
FALSE
,
reduce
=
FALSE
,
surrogate
=
TRUE
,
surrogate.FUN
=
MySurr
,
model
=
mods
,
baleq
=
baleq
,
tol
=
1E-5
,
call_pqc
=
TRUE
)
par
(
mfrow
=
c
(
2
,
1
))
Matplot
(
101
,
sur
=
so
,
sim
=
oo
,
xlab
=
"Domain element"
,
ylab
=
"Ca, Mg [molal], Calcite, Dolomite [mol]"
,
main
=
"RandomForest surrogate"
)
Matplot
(
101
,
sur
=
sf
,
sim
=
of
,
xlab
=
"Domain element"
,
ylab
=
"Ca, Mg [molal], Calcite, Dolomite [mol]"
,
main
=
"RandomForest surrogate"
)
inst/extdata/phreeqc_kin.dat
0 → 100644
View file @
57e6d53e
### This is the standard "phreeqc.dat" stripped of EXCHANGE and
### SURFACE and with the RATES for Calcite and Dolomite to use with
### RedModRphree
### Time-stamp: "Last modified 2018-05-06 14:36:23 delucia"
# PHREEQC.DAT for calculating pressure dependence of reactions, with
# molal volumina of aqueous species and of minerals, and
# critical temperatures and pressures of gases used in Peng-Robinson's EOS.
# Details are given at the end of this file.
SOLUTION_MASTER_SPECIES
#
#element species alk gfw_formula element_gfw
#
H H+ -1.0 H 1.008
H(0) H2 0 H
H(1) H+ -1.0 0
E e- 0 0.0 0
O H2O 0 O 16.0
O(0) O2 0 O
O(-2) H2O 0 0
Ca Ca+2 0 Ca 40.08
Mg Mg+2 0 Mg 24.312
Na Na+ 0 Na 22.9898
K K+ 0 K 39.102
Fe Fe+2 0 Fe 55.847
Fe(+2) Fe+2 0 Fe
Fe(+3) Fe+3 -2.0 Fe
Mn Mn+2 0 Mn 54.938
Mn(+2) Mn+2 0 Mn
Mn(+3) Mn+3 0 Mn
Al Al+3 0 Al 26.9815
Ba Ba+2 0 Ba 137.34
Sr Sr+2 0 Sr 87.62
Si H4SiO4 0 SiO2 28.0843
Cl Cl- 0 Cl 35.453
C CO3-2 2.0 HCO3 12.0111
C(+4) CO3-2 2.0 HCO3
C(-4) CH4 0 CH4
Alkalinity CO3-2 1.0 Ca0.5(CO3)0.5 50.05
S SO4-2 0 SO4 32.064
S(6) SO4-2 0 SO4
S(-2) HS- 1.0 S
N NO3- 0 N 14.0067
N(+5) NO3- 0 N
N(+3) NO2- 0 N
N(0) N2 0 N
N(-3) NH4+ 0 N 14.0067
#Amm AmmH+ 0 AmmH 17.0
B H3BO3 0 B 10.81
P PO4-3 2.0 P 30.9738