Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F71804117
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
Sat, Jul 13, 05:26
Size
6 KB
Mime Type
text/x-c
Expires
Mon, Jul 15, 05:26 (2 d)
Engine
blob
Format
Raw Data
Handle
18996293
Attached To
rSPECMICP SpecMiCP / ReactMiCP
transport_program.cpp
View Options
#include "transport_program.hpp"
#include <iostream>
#define EPS_J 1e-8
namespace
specmicp
{
namespace
reactmicp
{
namespace
systems
{
namespace
siasaturated
{
void
SaturatedDiffusionProgram
::
number_equations
(
std
::
vector
<
index_t
>
list_fixed_nodes
)
{
m_ideq
=
Eigen
::
VectorXd
::
Zero
(
get_ndf
()
*
m_mesh
->
nnodes
());
for
(
auto
it
=
list_fixed_nodes
.
begin
();
it
<
list_fixed_nodes
.
end
();
++
it
)
{
for
(
index_t
component:
m_data
->
range_aqueous_component
())
{
m_ideq
(
get_dof_component
(
*
it
,
component
))
=
no_equation
;
}
}
index_t
neq
=
0
;
for
(
index_t
dof
=
0
;
dof
<
m_ideq
.
rows
();
++
dof
)
{
if
(
m_ideq
(
dof
)
!=
no_equation
)
{
m_ideq
(
dof
)
=
neq
;
++
neq
;
}
}
m_neq
=
neq
;
}
void
SaturatedDiffusionProgram
::
compute_residuals
(
const
Vector
&
displacement
,
const
Vector
&
velocity
,
Vector
&
residual
)
{
residual
=
Eigen
::
VectorXd
::
Zero
(
get_neq
());
for
(
index_t
element:
m_mesh
->
range_elements
())
{
element_residual_transport
(
element
,
displacement
,
velocity
,
residual
);
}
}
void
SaturatedDiffusionProgram
::
element_residual_transport
(
index_t
element
,
const
Eigen
::
VectorXd
&
displacement
,
const
Eigen
::
VectorXd
&
velocity
,
Eigen
::
VectorXd
&
residual
)
{
Eigen
::
VectorXd
element_residual
(
2
);
for
(
index_t
component:
m_data
->
range_aqueous_component
())
{
element_residual_transport_component
(
element
,
component
,
displacement
,
velocity
,
element_residual
);
for
(
index_t
en:
m_mesh
->
range_nen
())
{
const
index_t
node
=
m_mesh
->
get_node
(
element
,
en
);
const
index_t
id
=
m_ideq
(
get_dof_component
(
node
,
component
));
if
(
id
!=
no_equation
)
{
residual
(
id
)
+=
element_residual
(
en
);}
}
}
}
void
SaturatedDiffusionProgram
::
element_residual_transport_component
(
index_t
element
,
index_t
component
,
const
Vector
&
displacement
,
const
Vector
&
velocity
,
Vector
&
element_residual
)
{
Eigen
::
Matrix
<
scalar_t
,
2
,
2
>
mass
,
jacob
;
Eigen
::
Matrix
<
scalar_t
,
2
,
1
>
evelocity
,
edisplacement
;
scalar_t
mass_coeff
=
-
(
m_param
->
density_water
()
*
m_param
->
porosity
(
0
)
*
m_mesh
->
crosssection
()
*
m_mesh
->
get_dx
(
element
)
/
2
);
mass
<<
1
,
0
,
0
,
1
;
mass
*=
mass_coeff
;
scalar_t
flux_coeff
=
-
(
m_mesh
->
crosssection
()
/
m_mesh
->
get_dx
(
element
)
*
m_param
->
porosity
(
0
)
*
m_param
->
density_water
()
*
m_param
->
diffusion_coefficient
(
element
)
);
jacob
<<
1
,
-
1
,
-
1
,
1
;
jacob
*=
flux_coeff
;
evelocity
<<
velocity
(
get_dof_component
(
m_mesh
->
get_node
(
element
,
0
),
component
)),
velocity
(
get_dof_component
(
m_mesh
->
get_node
(
element
,
1
),
component
));
edisplacement
<<
displacement
(
get_dof_component
(
m_mesh
->
get_node
(
element
,
0
),
component
)),
displacement
(
get_dof_component
(
m_mesh
->
get_node
(
element
,
1
),
component
));
element_residual
=
mass
*
evelocity
+
jacob
*
edisplacement
;
for
(
index_t
en:
m_mesh
->
range_nen
())
{
element_residual
(
en
)
+=
(
m_mesh
->
crosssection
()
*
m_mesh
->
get_dx
(
element
)
/
2
*
external_flux
(
m_mesh
->
get_node
(
element
,
en
),
component
));
}
}
void
SaturatedDiffusionProgram
::
compute_jacobian
(
Vector
&
displacement
,
Vector
&
velocity
,
Eigen
::
SparseMatrix
<
scalar_t
>&
jacobian
,
scalar_t
alphadt
)
{
list_triplet_t
jacob
;
const
index_t
ncomp
=
m_data
->
nb_component
-
1
;
const
index_t
estimation
=
m_mesh
->
nnodes
()
*
(
ncomp
*
m_mesh
->
nen
);
jacob
.
reserve
(
estimation
);
for
(
index_t
element:
m_mesh
->
range_elements
())
{
element_jacobian_transport
(
element
,
displacement
,
velocity
,
jacob
,
alphadt
);
}
jacobian
=
Eigen
::
SparseMatrix
<
scalar_t
>
(
get_neq
(),
get_neq
());
jacobian
.
setFromTriplets
(
jacob
.
begin
(),
jacob
.
end
());
}
// Transport
// =========
void
SaturatedDiffusionProgram
::
element_jacobian_transport
(
index_t
element
,
Vector
&
displacement
,
Vector
&
velocity
,
list_triplet_t
&
jacobian
,
scalar_t
alphadt
)
{
for
(
index_t
component:
m_data
->
range_aqueous_component
())
{
Eigen
::
VectorXd
element_residual_orig
(
Eigen
::
VectorXd
::
Zero
(
2
));
element_residual_transport_component
(
element
,
component
,
displacement
,
velocity
,
element_residual_orig
);
for
(
index_t
en:
m_mesh
->
range_nen
())
{
Eigen
::
VectorXd
element_residual
(
Eigen
::
VectorXd
::
Zero
(
2
));
const
index_t
node
=
m_mesh
->
get_node
(
element
,
en
);
const
index_t
dof
=
get_dof_component
(
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_J
*
std
::
abs
(
tmp_v
);
if
(
h
==
0
)
h
=
EPS_J
;
velocity
(
dof
)
=
tmp_v
+
h
;
h
=
velocity
(
dof
)
-
tmp_v
;
displacement
(
dof
)
=
tmp_d
+
alphadt
*
h
;
element_residual_transport_component
(
element
,
component
,
displacement
,
velocity
,
element_residual
);
velocity
(
dof
)
=
tmp_v
;
displacement
(
dof
)
=
tmp_d
;
for
(
index_t
enr:
m_mesh
->
range_nen
())
{
const
index_t
noder
=
m_mesh
->
get_node
(
element
,
enr
);
const
index_t
idr
=
m_ideq
(
get_dof_component
(
noder
,
component
));
if
(
idr
==
no_equation
)
continue
;
jacobian
.
push_back
(
triplet_t
(
idr
,
idc
,
(
element_residual
(
enr
)
-
element_residual_orig
(
enr
))
/
h
));
}
}
}
}
void
SaturatedDiffusionProgram
::
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:
m_data
->
range_aqueous_component
())
{
const
index_t
dof
=
get_dof_component
(
node
,
component
);
const
index_t
id
=
m_ideq
(
dof
);
if
(
id
==
no_equation
)
continue
;
velocity
(
dof
)
+=
lambda
*
update
(
id
);
}
}
displacement
=
predictor
+
alpha_dt
*
velocity
;
}
}
// end namespace siasaturated
}
// end namespace systems
}
// end namespace reactmicp
}
// end namespace specmicp
Event Timeline
Log In to Comment