Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F65450623
diffusion.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
Mon, Jun 3, 21:41
Size
3 KB
Mime Type
text/x-c
Expires
Wed, Jun 5, 21:41 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
17980007
Attached To
rSPECMICP SpecMiCP / ReactMiCP
diffusion.cpp
View Options
#include "dfpmsolver/parabolic_driver.hpp"
#include "dfpm/meshes/uniform_mesh1d.hpp"
#include "dfpm/meshes/axisymmetric_uniform_mesh1d.hpp"
#include "reactmicp/systems/saturated_diffusion/transport_program.hpp"
#include "reactmicp/systems/saturated_diffusion/variables.hpp"
#include "reactmicp/systems/saturated_diffusion/transport_parameters.hpp"
#include "specmicp/reaction_path.hpp"
void
diffusion_test
()
{
specmicp
::
database
::
Database
dbhandler
(
"../data/cemdata_specmicp.js"
);
std
::
vector
<
std
::
string
>
to_keep
=
{
"H2O"
,
"H[+]"
,
"HCO3[-]"
};
dbhandler
.
keep_only_components
(
to_keep
);
specmicp
::
RawDatabasePtr
thedatabase
=
dbhandler
.
get_database
();
int
nb_nodes
=
100
+
1
;
double
dx
=
0.0025
;
double
radius
=
3.5
+
dx
*
1
;
double
height
=
8
;
//specmicp::reactmicp::mesh::Mesh1DPtr themesh = specmicp::reactmicp::mesh::uniform_mesh1d(nb_nodes, 0.0050, 5);
specmicp
::
reactmicp
::
mesh
::
Mesh1DPtr
themesh
=
specmicp
::
reactmicp
::
mesh
::
axisymmetric_uniform_mesh1d
(
nb_nodes
,
radius
,
dx
,
height
);
auto
parameters
=
std
::
make_shared
<
specmicp
::
reactmicp
::
systems
::
siasaturated
::
SaturatedDiffusionTransportParameters
>
(
nb_nodes
,
1e-8
,
0.2
);
specmicp
::
reactmicp
::
systems
::
siasaturated
::
SIASaturatedVariables
variables
(
thedatabase
,
themesh
,
parameters
);
specmicp
::
reactmicp
::
systems
::
siasaturated
::
SaturatedDiffusionProgram
the_program
(
themesh
,
thedatabase
,
parameters
);
specmicp
::
Vector
&
concentrations
=
variables
.
total_mobile_concentrations
();
concentrations
.
setOnes
();
concentrations
.
block
(
0
,
0
,
the_program
.
get_ndf
(),
1
).
setZero
();
for
(
int
node
=
0
;
node
<
50
+
1
;
++
node
)
{
concentrations
(
variables
.
get_dof_total_concentration
(
node
,
1
))
=
0.0
;
concentrations
(
variables
.
get_dof_total_concentration
(
node
,
2
))
=
0.0
;
parameters
->
diffusion_coefficient
(
node
)
=
1e-5
;
}
the_program
.
number_equations
({
0
,});
specmicp
::
dfpmsolver
::
ParabolicDriver
<
specmicp
::
reactmicp
::
systems
::
siasaturated
::
SaturatedDiffusionProgram
>
solver
(
the_program
);
solver
.
get_options
().
sparse_solver
=
specmicp
::
SparseSolver
::
GMRES
;
solver
.
get_options
().
alpha
=
1
;
double
initial_sum
=
0
;
for
(
specmicp
::
index_t
node:
themesh
->
range_nodes
())
{
initial_sum
+=
concentrations
(
variables
.
get_dof_total_concentration
(
node
,
2
))
*
themesh
->
get_volume_cell
(
node
);
}
double
explicit_dt
=
std
::
pow
(
themesh
->
get_dx
(
0
),
2
)
/
parameters
->
diffusion_coefficient
(
0
);
std
::
cout
<<
"Explicit dt : "
<<
explicit_dt
<<
std
::
endl
;
double
dt
=
2000
;
//0.1*explicit_dt;
double
total
=
0
;
double
target
=
10
*
24
*
3600
;
int
k
=
0
;
std
::
cout
<<
0.0
<<
"
\t
"
;
std
::
cout
<<
initial_sum
<<
"
\t
"
<<
0.0
<<
"
\t
"
<<
std
::
endl
;
while
(
total
<
target
)
{
//std::cout << total << std::endl;
specmicp
::
dfpmsolver
::
ParabolicDriverReturnCode
retcode
=
solver
.
solve_timestep
(
dt
,
concentrations
);
//specmicp_assert(retcode == specmicp::dfpmsolver::ParabolicDriverReturnCode::ResidualMinimized);
total
+=
dt
;
if
(
k
%
10
==
0
)
{
std
::
cout
<<
total
/
24.0
/
3600.0
<<
"
\t
"
;
double
sum
=
0
;
for
(
specmicp
::
index_t
node:
themesh
->
range_nodes
())
{
sum
+=
concentrations
(
variables
.
get_dof_total_concentration
(
node
,
2
))
*
themesh
->
get_volume_cell
(
node
);
}
std
::
cout
<<
sum
<<
"
\t
"
<<
(
initial_sum
-
sum
)
/
(
2
*
3.5
*
M_PI
*
height
*
1e-4
)
<<
"
\t
"
<<
std
::
endl
;
}
++
k
;
}
}
int
main
()
{
diffusion_test
();
}
Event Timeline
Log In to Comment