Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F81105286
diffusion.hpp
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
Wed, Sep 4, 05:44
Size
10 KB
Mime Type
text/x-c++
Expires
Fri, Sep 6, 05:44 (2 d)
Engine
blob
Format
Raw Data
Handle
20495946
Attached To
rSPECMICP SpecMiCP / ReactMiCP
diffusion.hpp
View Options
/*-------------------------------------------------------
- Module : reactmicp/systems/diffusion
- File : diffusion.hpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget <fabieng@princeton.edu>, Princeton University
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name of the Princeton University nor the
names of its contributors may be used to endorse or promote products
derived from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
---------------------------------------------------------*/
#ifndef SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_HPP
#define SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_HPP
#include <vector>
#include <memory>
#include "database/data_container.hpp"
#include "reactmicp/common.hpp"
#include "reactmicp/micpsolvers/parabolicprogram.hpp"
#include "reactmicp/meshes/mesh1d.hpp"
#include "diffusion_numbering.hpp"
#include "diffusion_secondary.hpp"
namespace
specmicp
{
namespace
reactmicp
{
namespace
systems
{
//! \brief Parameters for the diffusion system
class
DiffusionParameter
{
public
:
DiffusionParameter
(
double
diffusion_coefficient
,
double
porosity
)
:
m_diff
(
diffusion_coefficient
),
m_poro
(
porosity
)
{}
//! \brief Density of water (kg/m^3)
double
density_water
()
{
return
1e3
;}
//! \brief Diffusion coefficient (element by element) (m^2/s)
double
diffusion_coefficient
(
ind_t
element
)
{
return
m_diff
;}
//! \brief Porosity (at a node (function of composition))
double
porosity
(
ind_t
node
)
{
return
m_poro
;}
private
:
double
m_diff
;
double
m_poro
;
};
//! \brief Saturated diffusion in reactive porous media
class
DiffusionProgram
:
public
solvers
::
ParabolicProgram
<
DiffusionProgram
>
{
public
:
using
Variables
=
ParabolicVariables
;
const
int
no_equation
=
-
1
;
DiffusionProgram
(
std
::
shared_ptr
<
mesh
::
Mesh1D
>
themesh
,
std
::
shared_ptr
<
database
::
DataContainer
>
thedatabase
,
std
::
shared_ptr
<
DiffusionParameter
>
parameters
,
diffusion
::
ListBoundaryCondition
&
list_bcs
,
Variables
&
var
)
:
m_mesh
(
themesh
),
m_data
(
thedatabase
),
m_param
(
parameters
),
m_ideq
(
themesh
->
nnodes
(),
thedatabase
,
list_bcs
,
var
),
m_second
(
themesh
->
nnodes
(),
thedatabase
,
m_ideq
)
{
}
DiffusionProgram
(
std
::
shared_ptr
<
mesh
::
Mesh1D
>
themesh
,
std
::
shared_ptr
<
database
::
DataContainer
>
thedatabase
,
std
::
shared_ptr
<
DiffusionParameter
>
parameters
,
diffusion
::
ListBoundaryCondition
&
list_bcs
,
const
EquilibriumState
&
initialstate
,
Variables
&
var
);
//! \brief Return the number of degree of freedoms
ind_t
get_neq
()
const
{
return
m_ideq
.
get_neq
();}
ind_t
get_ndf
()
const
{
return
m_ideq
.
get_ndf
();}
//! \brief Return the residuals
void
get_residuals
(
const
ParabolicVariables
&
variable
,
Vector
&
residual
);
//! brief Return the jacobian
void
get_jacobian
(
ParabolicVariables
&
variable
,
Vector
&
residual
,
SparseMatrix
&
jacobian
,
double
alphadt
);
//! \brief Return the index of the MiCP equations
const
std
::
vector
<
bool
>&
micp_index
()
const
{
return
m_is_micp
;}
//! \brief Update the variables
void
update_variables
(
ParabolicVariables
&
x
,
const
Eigen
::
VectorXd
&
predictor
,
const
Eigen
::
VectorXd
&
update
,
double
factor
,
double
alpha_dt
);
//! \brief set the predictor
void
set_predictor
(
ParabolicVariables
&
x
,
Eigen
::
VectorXd
&
predictor
,
double
alpha
,
double
dt
);
//! \brief compute secondary concentrations
bool
hook_start_iteration
(
const
ParabolicVariables
&
x
,
double
norm_residual
);
//! \brief Return the maximum lambda allowed for the linesearch
double
maximum_lambda
(
const
ParabolicVariables
&
x
,
const
Eigen
::
VectorXd
&
update
,
double
alphadt
);
//! \brief copy mineral amount from variables
void
copy_cp_variables
(
const
ParabolicVariables
&
x
,
Eigen
::
VectorXd
&
moles_mineral
);
//! \brief projection step
void
projection
(
ParabolicVariables
&
x
);
// todo or not........
// double update_equation(int id_eq,
// double update,
// double predictor,
// double alphadt
// )
// {
// // find _id eq ....
// }
private
:
// Residuals
// =========
// transport
// ---------
//! \brief Compute the residuals for the diffusion equations
void
residual_transport
(
const
Variables
&
variable
,
Vector
&
residual
);
//! \brief Contribution of the diffusion operator to the residuals
void
element_residual_transport
(
ind_t
element
,
const
Variables
&
variable
,
Vector
&
residual
);
//! \brief Contribution of the minerals to the residuals of the diffusion equations
double
nodal_mineral_transient_term_transport
(
ind_t
node
,
ind_t
component
,
const
Variables
&
variable
);
void
element_residual_transport_component
(
ind_t
element
,
int
component
,
const
Variables
&
variable
,
Vector
&
element_residual
);
// speciation
// ----------
//! \brief Compute the residual for the speciation problem
void
residual_speciation
(
const
Variables
&
variable
,
Vector
&
residual
);
//! \brief Compute the residuals of the aqueous mass balances
void
nodal_residual_massbalance
(
ind_t
node
,
const
Variables
&
variable
,
Vector
&
residual
);
//! \brief compute the residual for the aqueous mass balance of 'component' at 'node'
double
nodal_component_residual_massbalance
(
ind_t
node
,
int
component
,
const
Variables
&
variable
);
//! \brief Compute the residuals of the mineral speciation
void
nodal_residual_mineral
(
ind_t
node
,
const
Variables
&
variable
,
Vector
&
residual
);
//! \brief compute the residual for the stability equation of 'mineral' at 'node'
double
nodal_mineral_residual_mineral
(
ind_t
node
,
int
mineral
,
const
Variables
&
variable
);
// Jacobian
// ========
// transport
// ---------
//! \brief Contribution of the diffusion equation to the jacobian
void
jacobian_transport
(
Variables
&
variable
,
Vector
&
residual
,
list_triplet_t
&
jacobian
,
double
alphadt
);
//! \brief Element by element assembly procedure
void
element_jacobian_transport
(
ind_t
element
,
Variables
&
variable
,
Vector
&
residual
,
list_triplet_t
&
jacobian
,
double
alphadt
);
// speciation
// ----------
//! \brief Contribution of the speciation problem to the mass balance
void
jacobian_speciation
(
Variables
&
variable
,
list_triplet_t
&
jacobian
,
double
alphadt
);
//! \brief Contribution of the mass balance at 'node' to the jacobian
void
nodal_jacobian_massbalance
(
ind_t
node
,
const
Variables
&
variable
,
list_triplet_t
&
jacobian
,
double
alphadt
);
//! \brief Contribution of the mineral at equilibrium problem to the jacobian
void
nodal_jacobian_mineral
(
ind_t
node
,
const
Variables
&
variable
,
list_triplet_t
&
jacobian
,
double
alphadt
);
//! \brief Finite difference jacobian for speciation
void
nodal_jacobian_speciation_fd
(
ind_t
node
,
Variables
&
variable
,
list_triplet_t
&
jacobian
,
double
alphadt
);
// Initialization
// ==============
void
initialize_no_bc_with_init
(
const
EquilibriumState
&
initialstate
,
ParabolicVariables
&
var
);
// Update
// ======
double
update_diffusion
(
double
update
,
double
predictor
,
double
alphadt
)
{
return
predictor
+
alphadt
*
update
;
}
double
update_massbalance
(
double
update
,
double
predictor
,
double
alphadt
)
{
return
predictor
+
alphadt
*
update
;
}
double
update_mineral
(
double
update
,
double
predictor
,
double
alphadt
)
{
return
predictor
+
alphadt
*
update
;
}
// Member
// =======
std
::
shared_ptr
<
mesh
::
Mesh1D
>
m_mesh
;
std
::
shared_ptr
<
database
::
DataContainer
>
m_data
;
std
::
shared_ptr
<
DiffusionParameter
>
m_param
;
diffusion
::
DiffusionNumbering
m_ideq
;
diffusion
::
DiffusionSecondaryVariables
m_second
;
std
::
vector
<
bool
>
m_is_micp
;
//SparseVector m_boundary_condition;
};
}
// end namespace systems
}
// end namespace reactmicp
}
// end namespace specmicp
#endif
// SPECMICP_REACTMICP_SYSTEMS_DIFFUSION_HPP
Event Timeline
Log In to Comment