Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F82121707
transport_program.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, Sep 9, 17:15
Size
10 KB
Mime Type
text/x-c
Expires
Wed, Sep 11, 17:15 (2 d)
Engine
blob
Format
Raw Data
Handle
20656141
Attached To
rSPECMICP SpecMiCP / ReactMiCP
transport_program.cpp
View Options
#include "transport_program.hpp"
#include "dfpm/meshes/mesh1d.hpp"
#include "variables.hpp"
namespace
specmicp
{
namespace
reactmicp
{
namespace
systems
{
namespace
satdiff
{
//class SaturatedDiffusion::
SaturatedDiffusion
::
SaturatedDiffusion
(
SaturatedVariablesPtr
variables
,
std
::
vector
<
index_t
>
list_fixed_nodes
)
:
m_mesh
(
variables
->
get_mesh
()),
m_variables
(
variables
),
m_is_in_residual_computation
(
false
),
m_ndf
(
2
*
variables
->
nb_component
()),
m_tot_ndf
(
2
*
variables
->
nb_component
()
*
variables
->
get_mesh
()
->
nb_nodes
())
{
number_equations
(
list_fixed_nodes
,
{},
{
0
,});
}
SaturatedDiffusion
::
SaturatedDiffusion
(
SaturatedVariablesPtr
variables
,
std
::
vector
<
index_t
>
list_fixed_nodes
,
std
::
map
<
index_t
,
index_t
>
list_slave_nodes
,
std
::
vector
<
index_t
>
list_immobile_components
)
:
m_mesh
(
variables
->
get_mesh
()),
m_variables
(
variables
),
m_is_in_residual_computation
(
false
),
m_ndf
(
2
*
variables
->
nb_component
()),
m_tot_ndf
(
2
*
variables
->
nb_component
()
*
variables
->
get_mesh
()
->
nb_nodes
())
{
number_equations
(
list_fixed_nodes
,
list_slave_nodes
,
list_immobile_components
);
}
void
SaturatedDiffusion
::
number_equations
(
std
::
vector
<
index_t
>
list_fixed_nodes
,
std
::
map
<
index_t
,
index_t
>
list_slave_nodes
,
std
::
vector
<
index_t
>
list_immobile_components
)
{
m_ideq
.
resizeLike
(
m_variables
->
displacement
());
m_ideq
.
setZero
();
// flag fixed nodes
for
(
index_t
node:
list_fixed_nodes
)
{
for
(
index_t
component
=
1
;
component
<
m_variables
->
nb_component
();
++
component
)
{
m_ideq
(
m_variables
->
dof_aqueous_concentration
(
node
,
component
))
=
no_equation
;
}
}
// flag slaves nodes
// we flag them by making their ideq more negative than no_equation
for
(
auto
slave_pair:
list_slave_nodes
)
{
for
(
index_t
component
=
1
;
component
<
m_variables
->
nb_component
();
++
component
)
{
m_ideq
(
m_variables
->
dof_aqueous_concentration
(
slave_pair
.
first
,
component
))
=
no_equation
-
1
;
}
}
// set equation numbers
index_t
neq
=
0
;
for
(
index_t
node
=
0
;
node
<
m_mesh
->
nb_nodes
();
++
node
)
{
for
(
index_t
component:
list_immobile_components
)
{
m_ideq
(
m_variables
->
dof_aqueous_concentration
(
node
,
component
))
=
no_equation
;
m_ideq
(
m_variables
->
dof_aqueous_concentration
(
node
,
component
))
=
no_equation
;
}
for
(
index_t
component
=
0
;
component
<
m_variables
->
nb_component
();
++
component
)
{
index_t
dof
=
m_variables
->
dof_aqueous_concentration
(
node
,
component
);
if
(
m_ideq
(
dof
)
>
no_equation
)
// don't attribute an equation number if it is a slave or a fixed node
{
m_ideq
(
dof
)
=
neq
;
++
neq
;
}
dof
=
m_variables
->
dof_solid_concentration
(
node
,
component
);
m_ideq
(
dof
)
=
no_equation
;
}
}
// slave nodes
// attribute the correct equation number
for
(
auto
slave_pair:
list_slave_nodes
)
{
for
(
index_t
component
=
1
;
component
<
m_variables
->
nb_component
();
++
component
)
{
m_ideq
(
m_variables
->
dof_aqueous_concentration
(
slave_pair
.
first
,
component
))
=
m_ideq
(
m_variables
->
dof_aqueous_concentration
(
slave_pair
.
second
,
component
));
}
}
m_neq
=
neq
;
}
void
SaturatedDiffusion
::
compute_residuals
(
const
Vector
&
displacement
,
const
Vector
&
velocity
,
Vector
&
residual
)
{
residual
=
Vector
::
Zero
(
get_neq
());
m_is_in_residual_computation
=
true
;
for
(
index_t
element:
m_mesh
->
range_elements
())
{
for
(
index_t
component
=
1
;
component
<
m_variables
->
nb_component
();
++
component
)
{
Eigen
::
Vector2d
element_residual
;
element_residual
.
setZero
();
residuals_element_component
(
element
,
component
,
displacement
,
velocity
,
element_residual
);
for
(
index_t
en
=
0
;
en
<
2
;
++
en
)
{
const
index_t
node
=
m_mesh
->
get_node
(
element
,
en
);
const
index_t
id
=
m_ideq
(
m_variables
->
dof_aqueous_concentration
(
node
,
component
));
if
(
id
!=
no_equation
)
{
residual
(
id
)
+=
element_residual
(
en
);}
}
}
}
m_is_in_residual_computation
=
false
;
}
void
SaturatedDiffusion
::
residuals_element_component
(
index_t
element
,
index_t
component
,
const
Vector
&
displacement
,
const
Vector
&
velocity
,
Eigen
::
Vector2d
&
element_residual
)
{
const
scalar_t
mass_coeff_0
=
-
m_mesh
->
get_volume_cell_element
(
element
,
0
);
const
scalar_t
mass_coeff_1
=
-
m_mesh
->
get_volume_cell_element
(
element
,
1
);
const
index_t
node_0
=
m_mesh
->
get_node
(
element
,
0
);
const
index_t
node_1
=
m_mesh
->
get_node
(
element
,
1
);
const
scalar_t
diff_coeff
=
1.0
/
(
0.5
/
m_variables
->
diffusion_coefficient
(
node_0
)
+
0.5
/
m_variables
->
diffusion_coefficient
(
node_1
));
scalar_t
flux_coeff
=
-
(
m_mesh
->
get_face_area
(
element
)
/
m_mesh
->
get_dx
(
element
)
*
diff_coeff
)
;
const
index_t
dof_0
=
m_variables
->
dof_aqueous_concentration
(
node_0
,
component
);
const
index_t
dof_1
=
m_variables
->
dof_aqueous_concentration
(
node_1
,
component
);
// diffusion
scalar_t
flux_diffusion
=
flux_coeff
*
(
displacement
(
dof_0
)
-
displacement
(
dof_1
));
element_residual
(
0
)
=
flux_diffusion
;
element_residual
(
1
)
=
-
flux_diffusion
;
// advection
if
(
m_variables
->
fluid_velocity
(
element
)
!=
0.0
)
{
scalar_t
flux_advection
=
(
m_mesh
->
get_face_area
(
element
))
*
m_variables
->
fluid_velocity
(
element
);
if
(
m_variables
->
fluid_velocity
(
element
)
>
0
)
{
flux_advection
*=
(
displacement
(
dof_0
)
-
displacement
(
dof_1
));
element_residual
(
1
)
+=
flux_advection
;
}
else
{
flux_advection
*=
(
displacement
(
dof_1
)
-
displacement
(
dof_0
));
element_residual
(
0
)
-=
flux_advection
;
}
}
if
(
m_is_in_residual_computation
)
{
m_variables
->
aqueous_concentration
(
node_0
,
component
,
m_variables
->
transport_rate
())
+=
element_residual
(
0
);
m_variables
->
aqueous_concentration
(
node_1
,
component
,
m_variables
->
transport_rate
())
+=
element_residual
(
1
);
}
// velocity
element_residual
(
0
)
+=
mass_coeff_0
*
(
velocity
(
dof_0
)
*
m_variables
->
porosity
(
node_0
)
+
m_variables
->
vel_porosity
(
node_0
)
*
displacement
(
dof_0
));
element_residual
(
1
)
+=
mass_coeff_1
*
(
velocity
(
dof_1
)
*
m_variables
->
porosity
(
node_1
)
+
m_variables
->
vel_porosity
(
node_1
)
*
displacement
(
dof_1
));
// external rate
element_residual
(
0
)
+=
mass_coeff_0
*
m_variables
->
solid_concentration
(
node_0
,
component
,
m_variables
->
chemistry_rate
());
element_residual
(
1
)
+=
mass_coeff_1
*
m_variables
->
solid_concentration
(
node_1
,
component
,
m_variables
->
chemistry_rate
());
}
void
SaturatedDiffusion
::
compute_jacobian
(
Vector
&
displacement
,
Vector
&
velocity
,
Eigen
::
SparseMatrix
<
scalar_t
>&
jacobian
,
scalar_t
alphadt
)
{
dfpm
::
list_triplet_t
jacob
;
const
index_t
ncomp
=
m_variables
->
nb_component
();
const
index_t
estimation
=
m_mesh
->
nb_nodes
()
*
(
ncomp
*
m_mesh
->
nen
);
jacob
.
reserve
(
estimation
);
for
(
index_t
element:
m_mesh
->
range_elements
())
{
jacobian_element
(
element
,
displacement
,
velocity
,
jacob
,
alphadt
);
}
jacobian
=
Eigen
::
SparseMatrix
<
scalar_t
>
(
get_neq
(),
get_neq
());
jacobian
.
setFromTriplets
(
jacob
.
begin
(),
jacob
.
end
());
}
void
SaturatedDiffusion
::
jacobian_element
(
index_t
element
,
Vector
&
displacement
,
Vector
&
velocity
,
dfpm
::
list_triplet_t
&
jacobian
,
scalar_t
alphadt
)
{
for
(
index_t
component
=
1
;
component
<
m_variables
->
nb_component
();
++
component
)
{
Eigen
::
Vector2d
element_residual_orig
(
Eigen
::
Vector2d
::
Zero
());
residuals_element_component
(
element
,
component
,
displacement
,
velocity
,
element_residual_orig
);
for
(
index_t
en
=
0
;
en
<
2
;
++
en
)
{
Eigen
::
Vector2d
element_residual
(
Eigen
::
Vector2d
::
Zero
());
const
index_t
node
=
m_mesh
->
get_node
(
element
,
en
);
const
index_t
dof
=
m_variables
->
dof_aqueous_concentration
(
node
,
component
);
const
index_t
idc
=
m_ideq
(
dof
);
if
(
idc
==
no_equation
)
continue
;
const
scalar_t
tmp_v
=
velocity
(
dof
);
const
scalar_t
tmp_d
=
displacement
(
dof
);
scalar_t
h
=
eps_jacobian
*
std
::
abs
(
tmp_v
);
if
(
h
<
1e-4
*
eps_jacobian
)
h
=
eps_jacobian
;
velocity
(
dof
)
=
tmp_v
+
h
;
h
=
velocity
(
dof
)
-
tmp_v
;
displacement
(
dof
)
=
tmp_d
+
alphadt
*
h
;
residuals_element_component
(
element
,
component
,
displacement
,
velocity
,
element_residual
);
velocity
(
dof
)
=
tmp_v
;
displacement
(
dof
)
=
tmp_d
;
for
(
index_t
enr
=
0
;
enr
<
2
;
++
enr
)
{
const
index_t
noder
=
m_mesh
->
get_node
(
element
,
enr
);
const
index_t
idr
=
m_ideq
(
m_variables
->
dof_aqueous_concentration
(
noder
,
component
));
if
(
idr
==
no_equation
)
continue
;
jacobian
.
push_back
(
dfpm
::
triplet_t
(
idr
,
idc
,
(
element_residual
(
enr
)
-
element_residual_orig
(
enr
))
/
h
));
}
}
}
}
//! \brief Update the solutions
void
SaturatedDiffusion
::
update_solution
(
const
Vector
&
update
,
scalar_t
lambda
,
scalar_t
alpha_dt
,
Vector
&
predictor
,
Vector
&
displacement
,
Vector
&
velocity
)
{
for
(
index_t
node:
m_mesh
->
range_nodes
())
{
for
(
index_t
component
=
1
;
component
<
m_variables
->
nb_component
();
++
component
)
{
const
index_t
dof
=
m_variables
->
dof_aqueous_concentration
(
node
,
component
);
const
index_t
id
=
m_ideq
(
dof
);
if
(
id
==
no_equation
)
continue
;
velocity
(
dof
)
+=
lambda
*
update
(
id
);
}
}
//displacement = m_variables->predictor() + alpha_dt*velocity;
displacement
=
predictor
+
alpha_dt
*
velocity
;
}
}
// end namespace satdiff
}
// end namespace systems
}
// end namespace reactmicp
}
// end namespace specmicp
Event Timeline
Log In to Comment