Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
30 commits
Select commit Hold shift + click to select a range
8358586
[PL] added Jacobian assembler interface and implementations
chleh Aug 19, 2016
2f0be12
[PL] re-introduced VectorMatrixAssembler
chleh Aug 19, 2016
1175265
[MaL] tools for mapping std::vector -> Eigen::Matrix
chleh Aug 19, 2016
9b47442
[PL] made Process Jacobian-assembly-ready
chleh Aug 19, 2016
0043119
[PL] made specific processes Jacobian-assembly-ready
chleh Aug 19, 2016
fff563f
[AppL] ProjectData constructs Jacobian assembler
chleh Aug 19, 2016
469267f
[NL] interface changes
chleh Aug 19, 2016
527a7da
[T] adapted to time integration interface changes
chleh Aug 19, 2016
caa581c
[T] added Jacobian assembler tests
chleh Aug 19, 2016
e0db088
[App] added Newton-Raphson ctests
chleh Aug 19, 2016
e9efbc5
[T] fixed tests for PETSc
chleh Aug 22, 2016
cf469b5
fixed compilation errors
chleh Aug 24, 2016
feb6b82
[MaL] templated EigenMapTools
chleh Aug 24, 2016
1eb5719
[PL] adapted to changed EigenMapTools
chleh Aug 24, 2016
557a341
[T] adapted to changed EigenMapTools
chleh Aug 24, 2016
ade940f
[PL] use EigenMapTools in local assemblers.
chleh Aug 24, 2016
a3ba818
[PL] added Jacobian assembler docu
chleh Aug 24, 2016
e5e6d5d
[PL] documented VectorMatrixAssembler
chleh Aug 24, 2016
8d0437b
renamed toZeroedMatrix -> createZeroedMatrix
chleh Sep 1, 2016
7ade23e
[NL] improved documentation
chleh Sep 1, 2016
dfc4437
[PL] improved error message
chleh Sep 1, 2016
063d004
fixup App/ProjectData
endJunction Sep 11, 2016
5d482e0
fixup GWF
endJunction Sep 11, 2016
358a355
fixup HeatConduction
endJunction Sep 11, 2016
d9cd80d
fixup SD
endJunction Sep 11, 2016
2446e04
fixup TES
endJunction Sep 11, 2016
ec2f050
[PL] fix compilation errors
chleh Sep 19, 2016
398cac6
[Data] added Newton ctests
chleh Sep 19, 2016
96dee82
[App] marked TES Newton test as large
chleh Sep 20, 2016
07bc9ef
[PL] 0 --> 0.0
chleh Sep 20, 2016
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
21 changes: 14 additions & 7 deletions Applications/ApplicationsLib/ProjectData.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@
#include "MeshLib/Mesh.h"

#include "NumLib/ODESolver/ConvergenceCriterion.h"
#include "ProcessLib/CreateJacobianAssembler.h"

// FileIO
#include "GeoLib/IO/XmlIO/Boost/BoostXmlGmlInterface.h"
Expand Down Expand Up @@ -253,6 +254,9 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config)

std::unique_ptr<ProcessLib::Process> process;

auto jacobian_assembler = ProcessLib::createJacobianAssembler(
process_config.getConfigSubtreeOptional("jacobian_assembler"));

if (type == "GROUNDWATER_FLOW")
{
// The existence check of the in the configuration referenced
Expand All @@ -261,17 +265,20 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config)
// several meshes. Then we have to assign the referenced mesh
// here.
process = ProcessLib::GroundwaterFlow::createGroundwaterFlowProcess(
*_mesh_vec[0], _process_variables, _parameters, process_config);
*_mesh_vec[0], std::move(jacobian_assembler),
_process_variables, _parameters, process_config);
}
else if (type == "TES")
{
process = ProcessLib::TES::createTESProcess(
*_mesh_vec[0], _process_variables, _parameters, process_config);
*_mesh_vec[0], std::move(jacobian_assembler),
_process_variables, _parameters, process_config);
}
else if (type == "HEAT_CONDUCTION")
{
process = ProcessLib::HeatConduction::createHeatConductionProcess(
*_mesh_vec[0], _process_variables, _parameters, process_config);
*_mesh_vec[0], std::move(jacobian_assembler),
_process_variables, _parameters, process_config);
}
else if (type == "SMALL_DEFORMATION")
{
Expand All @@ -280,14 +287,14 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config)
case 2:
process = ProcessLib::SmallDeformation::
createSmallDeformationProcess<2>(
*_mesh_vec[0], _process_variables, _parameters,
process_config);
*_mesh_vec[0], std::move(jacobian_assembler),
_process_variables, _parameters, process_config);
break;
case 3:
process = ProcessLib::SmallDeformation::
createSmallDeformationProcess<3>(
*_mesh_vec[0], _process_variables, _parameters,
process_config);
*_mesh_vec[0], std::move(jacobian_assembler),
_process_variables, _parameters, process_config);
break;
default:
OGS_FATAL(
Expand Down
54 changes: 52 additions & 2 deletions Applications/CLI/Tests.cmake
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,18 @@ if(NOT OGS_USE_MPI)
cube_1x1x1_hex_${mesh_size}.vtu cube_${mesh_size}_pcs_0_ts_1_t_1.000000.vtu Linear_1_to_minus1 pressure
)

AddTest(
NAME GroundWaterFlowProcess_cube_1x1x1_${mesh_size}_Newton
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For discussion: This problem is linear, why newton? what does it test?

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The other tests of the linear GWFlow use Picard, which also doesn't make sense in a way.
It is tested that the central differences assembly works.

PATH Elliptic/cube_1x1x1_GroundWaterFlow
EXECUTABLE ogs
EXECUTABLE_ARGS cube_${mesh_size}_newton.prj
WRAPPER time
TESTER vtkdiff
ABSTOL 1e-15 RELTOL 1e-15
DIFF_DATA
cube_1x1x1_hex_${mesh_size}.vtu cube_${mesh_size}_newton_pcs_0_ts_1_t_1.000000.vtu Linear_1_to_minus1 pressure
)

AddTest(
NAME GroundWaterFlowProcess_cube_1x1x1_Neumann_${mesh_size}
PATH Elliptic/cube_1x1x1_GroundWaterFlow
Expand Down Expand Up @@ -198,7 +210,7 @@ if(NOT OGS_USE_MPI)
TESTER vtkdiff
ABSTOL 1e-14 RELTOL 1e-14
DIFF_DATA
cube_1e3_top_neumann_pcs_0_ts_1_t_1.000000.vtu cube_1e3_top_neumann_pcs_0_ts_1_t_1.000000.vtu pressure pressure
cube_1e3_top_neumann.vtu cube_1e3_top_neumann_pcs_0_ts_1_t_1.000000.vtu pressure pressure
)
AddTest(
NAME GroundWaterFlowProcess_cube_bottom
Expand All @@ -209,7 +221,30 @@ if(NOT OGS_USE_MPI)
TESTER vtkdiff
ABSTOL 1e-14 RELTOL 1e-14
DIFF_DATA
cube_1e3_bottom_neumann_pcs_0_ts_1_t_1.000000.vtu cube_1e3_bottom_neumann_pcs_0_ts_1_t_1.000000.vtu pressure pressure
cube_1e3_bottom_neumann.vtu cube_1e3_bottom_neumann_pcs_0_ts_1_t_1.000000.vtu pressure pressure
)
# Some Neumann BC tests -- Newton
AddTest(
NAME GroundWaterFlowProcess_cube_top_Newton
PATH Elliptic/cube_1x1x1_GroundWaterFlow
EXECUTABLE ogs
EXECUTABLE_ARGS cube_1e3_top_neumann_newton.prj
WRAPPER time
TESTER vtkdiff
ABSTOL 1e-14 RELTOL 1e-14
DIFF_DATA
cube_1e3_top_neumann.vtu cube_1e3_top_neumann_newton_pcs_0_ts_1_t_1.000000.vtu pressure pressure
)
AddTest(
NAME GroundWaterFlowProcess_cube_bottom_Newton
PATH Elliptic/cube_1x1x1_GroundWaterFlow
EXECUTABLE ogs
EXECUTABLE_ARGS cube_1e3_bottom_neumann_newton.prj
WRAPPER time
TESTER vtkdiff
ABSTOL 1e-14 RELTOL 1e-14
DIFF_DATA
cube_1e3_bottom_neumann.vtu cube_1e3_bottom_neumann_newton_pcs_0_ts_1_t_1.000000.vtu pressure pressure
)

# TES tests
Expand Down Expand Up @@ -245,6 +280,21 @@ if(NOT OGS_USE_MPI)
tes_zeolite_discharge_large_pcs_0_ts_28_t_1_000000.vtu tes_zeolite_discharge_large_pcs_0_ts_28_t_1.000000.vtu solid_density fct_solid_density
)

AddTest(
NAME LARGE_TES_zeolite_discharge_Newton
PATH Parabolic/TES/1D
EXECUTABLE ogs
EXECUTABLE_ARGS tes-1D-zeolite-discharge-small-newton.prj
WRAPPER time
TESTER vtkdiff
ABSTOL 1.5e-3 RELTOL 1.5e-3
DIFF_DATA
tes_zeolite_discharge_small_ts_19_t_0_100000.vtu tes_zeolite_discharge_small_newton_pcs_0_ts_32_t_0.100000.vtu pressure pressure
tes_zeolite_discharge_small_ts_19_t_0_100000.vtu tes_zeolite_discharge_small_newton_pcs_0_ts_32_t_0.100000.vtu temperature temperature
# tes_zeolite_discharge_small_ts_19_t_0_100000.vtu tes_zeolite_discharge_small_newton_pcs_0_ts_32_t_0.100000.vtu vapour_partial_pressure vapour_partial_pressure
tes_zeolite_discharge_small_ts_19_t_0_100000.vtu tes_zeolite_discharge_small_newton_pcs_0_ts_32_t_0.100000.vtu solid_density solid_density
)

AddTest(
NAME 1D_HeatConduction_dirichlet
PATH Parabolic/T/1D_dirichlet
Expand Down
193 changes: 193 additions & 0 deletions MathLib/LinAlg/Eigen/EigenMapTools.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,193 @@
/**
* \copyright
* Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
*/

#ifndef MATHLIB_EIGENMAPTOOLS_H
#define MATHLIB_EIGENMAPTOOLS_H

#include <cassert>
#include <vector>
#include <Eigen/Core>

namespace MathLib
{
/*! Creates an Eigen mapped matrix having its entries stored in the given
* \c data vector.
*
* \return An Eigen mapped matrix of the given size. All values of the matrix
* are set to zero.
*
* \pre The passed \c data vector must have zero size.
* \post The \c data has size \c rows * \c cols.
*
* \note The data vector will have the same storage order (row or column major)
* as the requested matrix type.
*/
template <typename Matrix>
Eigen::Map<Matrix> createZeroedMatrix(std::vector<double>& data,
Eigen::MatrixXd::Index rows,
Eigen::MatrixXd::Index cols)
{
static_assert(Matrix::IsRowMajor || Matrix::IsVectorAtCompileTime,
"The default storage order in OGS is row major storage for "
"dense matrices.");
assert(Matrix::RowsAtCompileTime == Eigen::Dynamic ||
Matrix::RowsAtCompileTime == rows);
assert(Matrix::ColsAtCompileTime == Eigen::Dynamic ||
Matrix::ColsAtCompileTime == cols);
assert(data.empty()); // in order that resize fills the vector with zeros.

data.resize(rows * cols);
return {data.data(), rows, cols};
}

/*! Creates an Eigen mapped matrix having its entries stored in the given
* \c data vector.
*
* This is a convienence method which makes the specification of dynamically
* allocated Eigen matrices as return type easier.
*/
inline Eigen::Map<
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
createZeroedMatrix(std::vector<double>& data,
Eigen::MatrixXd::Index rows,
Eigen::MatrixXd::Index cols)
{
return createZeroedMatrix<
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>(
data, rows, cols);
}

/*! Creates an Eigen mapped matrix from the given data vector.
*
* \attention The data vector must have the same storage order (row or column
* major) as the requested matrix type.
*/
template <typename Matrix>
Eigen::Map<const Matrix> toMatrix(std::vector<double> const& data,
Eigen::MatrixXd::Index rows,
Eigen::MatrixXd::Index cols)
{
static_assert(Matrix::IsRowMajor || Matrix::IsVectorAtCompileTime,
"The default storage order in OGS is row major storage for "
"dense matrices.");
assert(Matrix::RowsAtCompileTime == Eigen::Dynamic ||
Matrix::RowsAtCompileTime == rows);
assert(Matrix::ColsAtCompileTime == Eigen::Dynamic ||
Matrix::ColsAtCompileTime == cols);
assert(static_cast<Eigen::MatrixXd::Index>(data.size()) == rows * cols);

return {data.data(), rows, cols};
}

/*! Creates an Eigen mapped matrix from the given data vector.
*
* This is a convienence method which makes the specification of dynamically
* allocated Eigen matrices as return type easier.
*/
inline Eigen::Map<
const Eigen::
Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
toMatrix(std::vector<double> const& data,
Eigen::MatrixXd::Index rows,
Eigen::MatrixXd::Index cols)
{
return toMatrix<
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>(
data, rows, cols);
}

/*! Creates an Eigen mapped matrix from the given data vector.
*
* \attention The data vector must have the same storage order (row or column
* major) as the requested matrix type.
*/
template <typename Matrix>
Eigen::Map<Matrix> toMatrix(std::vector<double>& data,
Eigen::MatrixXd::Index rows,
Eigen::MatrixXd::Index cols)
{
static_assert(Matrix::IsRowMajor || Matrix::IsVectorAtCompileTime,
"The default storage order in OGS is row major storage for "
"dense matrices.");
assert(Matrix::RowsAtCompileTime == Eigen::Dynamic ||
Matrix::RowsAtCompileTime == rows);
assert(Matrix::ColsAtCompileTime == Eigen::Dynamic ||
Matrix::ColsAtCompileTime == cols);
assert(static_cast<Eigen::MatrixXd::Index>(data.size()) == rows * cols);

return {data.data(), rows, cols};
}

/*! Creates an Eigen mapped matrix from the given data vector.
*
* This is a convienence method which makes the specification of dynamically
* allocated Eigen matrices as return type easier.
*/
inline Eigen::Map<
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
toMatrix(std::vector<double>& data,
Eigen::MatrixXd::Index rows,
Eigen::MatrixXd::Index cols)
{
return toMatrix<
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>(
data, rows, cols);
}

/*! Creates an Eigen mapped vector having its entries stored in the given
* \c data vector.
*
* \return An Eigen mapped vector of the given size. All values of the vector
* are set to zero.
*
* \pre The passed \c data vector must have zero size.
* \post The \c data has size \c size.
*/
template <typename Vector>
Eigen::Map<Vector> createZeroedVector(std::vector<double>& data,
Eigen::VectorXd::Index size)
{
static_assert(Vector::IsVectorAtCompileTime, "A vector type is required.");
assert(Vector::SizeAtCompileTime == Eigen::Dynamic ||
Vector::SizeAtCompileTime == size);
assert(data.empty()); // in order that resize fills the vector with zeros.

data.resize(size);
return {data.data(), size};
}

//! Creates an Eigen mapped vector from the given data vector.
template <typename Vector>
Eigen::Map<const Vector> toVector(std::vector<double> const& data,
Eigen::VectorXd::Index size)
{
static_assert(Vector::IsVectorAtCompileTime, "A vector type is required.");
assert(Vector::SizeAtCompileTime == Eigen::Dynamic ||
Vector::SizeAtCompileTime == size);
assert(static_cast<Eigen::VectorXd::Index>(data.size()) == size);

return {data.data(), size};
}

//! Creates an Eigen mapped vector from the given data vector.
template <typename Vector>
Eigen::Map<Vector> toVector(std::vector<double>& data,
Eigen::VectorXd::Index size)
{
static_assert(Vector::IsVectorAtCompileTime, "A vector type is required.");
assert(Vector::SizeAtCompileTime == Eigen::Dynamic ||
Vector::SizeAtCompileTime == size);
assert(static_cast<Eigen::VectorXd::Index>(data.size()) == size);

return {data.data(), size};
}

} // MathLib

#endif // MATHLIB_EIGENMAPTOOLS_H
9 changes: 4 additions & 5 deletions NumLib/ODESolver/NonlinearSolver.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ namespace NumLib
void NonlinearSolver<NonlinearSolverTag::Picard>::assemble(
GlobalVector const& x) const
{
_equation_system->assembleMatricesPicard(x);
_equation_system->assemble(x);
}

bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
Expand Down Expand Up @@ -61,7 +61,7 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(

BaseLib::RunTime time_assembly;
time_assembly.start();
sys.assembleMatricesPicard(x);
sys.assemble(x);
sys.getA(A);
sys.getRhs(rhs);
INFO("[time] Assembly took %g s.", time_assembly.elapsed());
Expand Down Expand Up @@ -164,7 +164,7 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
void NonlinearSolver<NonlinearSolverTag::Newton>::assemble(
GlobalVector const& x) const
{
_equation_system->assembleResidualNewton(x);
_equation_system->assemble(x);
// TODO if the equation system would be reset to nullptr after each
// assemble() or solve() call, the user would be forced to set the
// equation every time and could not forget it.
Expand Down Expand Up @@ -205,9 +205,8 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(

BaseLib::RunTime time_assembly;
time_assembly.start();
sys.assembleResidualNewton(x);
sys.assemble(x);
sys.getResidual(x, res);
sys.assembleJacobian(x);
sys.getJacobian(J);
INFO("[time] Assembly took %g s.", time_assembly.elapsed());

Expand Down
Loading