Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78620210
test_coupling.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
Thu, Aug 22, 00:40
Size
6 KB
Mime Type
text/x-c++
Expires
Sat, Aug 24, 00:40 (2 d)
Engine
blob
Format
Raw Data
Handle
20059747
Attached To
rSPECMICP SpecMiCP / ReactMiCP
test_coupling.cpp
View Options
#include "catch.hpp"
#include "reactmicp/solver/staggers_base/staggers_base.hpp"
#include "reactmicp/solver/reactive_transport_solver.hpp"
#include <iostream>
namespace
test_internal
{
using
namespace
specmicp
;
using
namespace
specmicp
::
reactmicp
;
using
namespace
specmicp
::
reactmicp
::
solver
;
class
TestVariables
:
public
reactmicp
::
solver
::
VariablesBase
{
public
:
TestVariables
(
index_t
neq
)
:
m_neq
(
neq
),
predictor
(
Vector
::
Zero
(
neq
)),
displacement
(
Vector
::
Zero
(
neq
)),
velocity
(
Vector
::
Zero
(
neq
)),
transport_rate
(
Vector
::
Zero
(
neq
)),
chemistry_rate
(
Vector
::
Zero
(
neq
))
{}
index_t
get_neq
()
{
return
m_neq
;}
Vector
&
get_disp
()
{
return
displacement
;}
Vector
&
get_vel
()
{
return
velocity
;}
Vector
&
get_pred
()
{
return
predictor
;}
Vector
&
get_transport_rate
()
{
return
transport_rate
;}
Vector
&
get_chemistry_rate
()
{
return
chemistry_rate
;}
void
reset_main_variables
()
{
const
index_t
neq
=
get_neq
();
predictor
=
Vector
::
Zero
(
neq
);
displacement
=
Vector
::
Zero
(
neq
);
velocity
=
Vector
::
Zero
(
neq
);
transport_rate
=
Vector
::
Zero
(
neq
);
chemistry_rate
=
Vector
::
Zero
(
neq
);
}
private
:
index_t
m_neq
;
Vector
predictor
;
Vector
displacement
;
Vector
velocity
;
Vector
transport_rate
;
Vector
chemistry_rate
;
};
using
TestVariablesPtr
=
std
::
shared_ptr
<
TestVariables
>
;
inline
TestVariables
*
castvarptr
(
VariablesBase
*
varptr
)
{
return
static_cast
<
TestVariables
*
const
>
(
varptr
);
}
class
TransportStagger
:
public
reactmicp
::
solver
::
TransportStaggerBase
{
public
:
TransportStagger
(
Matrix
the_linear_system
)
:
linear_system
(
the_linear_system
)
{}
//! \brief Initialize the stagger at the beginning of the computation
virtual
void
initialize
(
VariablesBasePtr
var
)
{}
//! \brief Initialize the stagger at the beginning of an iteration
virtual
void
initialize_timestep
(
scalar_t
dt
,
VariablesBase
*
const
var
)
{
TestVariables
*
const
truevar
=
castvarptr
(
var
);
truevar
->
get_pred
()
=
truevar
->
get_disp
();
truevar
->
get_vel
().
setZero
();
m_dt
=
dt
;
}
//! \brief Solve the equation for the timestep
StaggerReturnCode
restart_timestep
(
VariablesBase
*
var
)
{
TestVariables
*
const
truevar
=
castvarptr
(
var
);
truevar
->
get_transport_rate
()
=
linear_system
*
truevar
->
get_disp
();
truevar
->
get_vel
()
=
truevar
->
get_transport_rate
()
+
truevar
->
get_chemistry_rate
();
truevar
->
get_disp
()
=
truevar
->
get_pred
()
+
m_dt
*
truevar
->
get_vel
();
return
StaggerReturnCode
::
ResidualMinimized
;
}
//! \brief Compute the residuals norm
scalar_t
get_residual
(
VariablesBase
*
var
)
{
TestVariables
*
const
truevar
=
castvarptr
(
var
);
Vector
residuals
=
truevar
->
get_vel
()
-
linear_system
*
truevar
->
get_disp
()
-
truevar
->
get_chemistry_rate
();
return
residuals
.
norm
();
}
//! \brief Compute the residuals norm
scalar_t
get_residual_0
(
VariablesBase
*
var
)
{
TestVariables
*
const
truevar
=
castvarptr
(
var
);
Vector
residuals
=
truevar
->
get_vel
()
-
linear_system
*
truevar
->
get_disp
()
-
truevar
->
get_chemistry_rate
();
return
residuals
.
norm
();
}
scalar_t
get_update
(
VariablesBase
*
var
)
{
TestVariables
*
const
truevar
=
castvarptr
(
var
);
return
truevar
->
get_vel
().
norm
();
}
private
:
scalar_t
m_dt
;
Matrix
linear_system
;
};
class
ChemistryStagger
:
public
reactmicp
::
solver
::
ChemistryStaggerBase
{
public
:
ChemistryStagger
(
Matrix
the_linear_system
)
:
linear_system
(
the_linear_system
)
{}
//! \brief Initialize the stagger at the beginning of the computation
virtual
void
initialize
(
VariablesBase
*
const
var
)
override
{}
//! \brief Initialize the stagger at the beginning of an iteration
virtual
void
initialize_timestep
(
scalar_t
dt
,
VariablesBase
*
const
var
)
override
{
m_dt
=
dt
;
}
//! \brief Solve the equation for the timestep
virtual
StaggerReturnCode
restart_timestep
(
VariablesBase
*
const
var
)
override
{
TestVariables
*
const
truevar
=
castvarptr
(
var
);
truevar
->
get_chemistry_rate
()
=
linear_system
*
truevar
->
get_disp
();
truevar
->
get_vel
()
=
truevar
->
get_transport_rate
()
+
truevar
->
get_chemistry_rate
();
truevar
->
get_disp
()
=
truevar
->
get_pred
()
+
m_dt
*
truevar
->
get_vel
();
return
StaggerReturnCode
::
ResidualMinimized
;
}
private
:
scalar_t
m_dt
;
Matrix
linear_system
;
};
class
VoidUpscalingStagger
:
public
UpscalingStaggerBase
{
//! \brief Initialize the stagger at the beginning of the computation
virtual
void
initialize
(
VariablesBase
*
const
var
)
override
{}
//! \brief Initialize the stagger at the beginning of an iteration
virtual
void
initialize_timestep
(
scalar_t
dt
,
VariablesBase
*
const
var
)
override
{}
//! \brief Solve the equation for the timestep
virtual
StaggerReturnCode
restart_timestep
(
VariablesBase
*
const
var
)
override
{
return
StaggerReturnCode
::
ResidualMinimized
;
}
};
}
//end namespace test_internal
TEST_CASE
(
"ReactiveTransport solver"
,
"[solver, fixed point]"
)
{
SECTION
(
"Test"
,
""
)
{
specmicp
::
index_t
neq
=
2
;
test_internal
::
TestVariablesPtr
var
=
std
::
make_shared
<
test_internal
::
TestVariables
>
(
neq
);
var
->
get_disp
()
<<
3
,
1
;
specmicp
::
Matrix
transport_system
(
neq
,
neq
);
transport_system
<<
3
,
-
1
,
2
,
-
2
;
specmicp
::
Matrix
chemistry_system
(
neq
,
neq
);
chemistry_system
<<
2
,
0
,
2
,
0
;
std
::
shared_ptr
<
specmicp
::
reactmicp
::
solver
::
TransportStaggerBase
>
tstag
=
std
::
make_shared
<
test_internal
::
TransportStagger
>
(
transport_system
);
std
::
shared_ptr
<
specmicp
::
reactmicp
::
solver
::
ChemistryStaggerBase
>
cstag
=
std
::
make_shared
<
test_internal
::
ChemistryStagger
>
(
chemistry_system
);
std
::
shared_ptr
<
specmicp
::
reactmicp
::
solver
::
UpscalingStaggerBase
>
ustag
=
std
::
make_shared
<
test_internal
::
VoidUpscalingStagger
>
();
specmicp
::
reactmicp
::
solver
::
ReactiveTransportSolver
solver
(
tstag
,
cstag
,
ustag
);
solver
.
get_options
().
maximum_iterations
=
100
;
specmicp
::
reactmicp
::
solver
::
ReactiveTransportReturnCode
retcode
=
solver
.
solve_timestep
(
0.1
,
var
);
REQUIRE
(
retcode
>
specmicp
::
reactmicp
::
solver
::
ReactiveTransportReturnCode
::
NotConvergedYet
);
std
::
cout
<<
"nb iterations 1 : "
<<
solver
.
get_perfs
().
nb_iterations
<<
" - residuals : "
<<
solver
.
get_perfs
().
residuals
<<
std
::
endl
;
retcode
=
solver
.
solve_timestep
(
0.1
,
var
);
REQUIRE
(
retcode
>
specmicp
::
reactmicp
::
solver
::
ReactiveTransportReturnCode
::
NotConvergedYet
);
std
::
cout
<<
"nb iterations 2 : "
<<
solver
.
get_perfs
().
nb_iterations
<<
" - residuals : "
<<
solver
.
get_perfs
().
residuals
<<
std
::
endl
;
}
}
Event Timeline
Log In to Comment