Skip to content

GitLab

Commit 9d10212a authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

some corrections of warnings

parent 442888a5
Hide whitespace changes
Inline Side-by-side
.gitlab-ci.yml
View file @ 9d10212a
......@@ -9,6 +9,8 @@ cache:
#- ./contrib/ci-setup
#- dunecontrol --opts=/duneci/opts.clang --current all
#- dunecontrol --current make test
#only:
#- master
dune:git--gcc:
image: duneci/dune-fufem:git
......@@ -16,3 +18,5 @@ dune:git--gcc:
- ./contrib/ci-setup
- dunecontrol --opts=/duneci/opts.gcc --current all
- dunecontrol --current make test
only:
- master
CMakeLists.txt
View file @ 9d10212a
cmake_minimum_required(VERSION 3.0)
project(dune-amdis CXX)
# set(ALBERTA_ROOT /opt/software/alberta)
# set(ALBERTA_INCLUDE_DIR /opt/software/alberta/include)
set(ALBERTA_ROOT /opt/software/alberta)
set(ALBERTA_INCLUDE_DIR /opt/software/alberta/include)
# set(CMAKE_PREFIX_PATH /opt/software/dune/lib/cmake)
# set(UG_DIR /opt/software/dune/lib/cmake/ug)
# set(Vc_DIR /opt/software/dune/lib/cmake/Vc)
......
dune/amdis/Initfile.hpp
View file @ 9d10212a
......@@ -102,7 +102,7 @@ namespace AMDiS
boost::char_separator<char> sep(",; ");
Tokenizer tokens(valStr, sep);
int i = 0;
size_t i = 0;
for (auto token : tokens)
{
AMDIS_TEST_EXIT(i < dim, "Vector data exceeds array dimension!");
......
dune/amdis/Operator.hpp
View file @ 9d10212a
......@@ -69,6 +69,13 @@ namespace AMDiS
return addSOTImpl(toTerm(A), _i, _j);
}
static std::shared_ptr<Self> create()
{
return std::make_shared<Self>();
}
/// Extract the polynomial degree from \p rowFeSpace and \p colFeSpace.
template <class RowFeSpace, class ColFeSpace>
void init(RowFeSpace const& rowFeSpace,
......
dune/amdis/linear_algebra/mtl/DOFVector.hpp
View file @ 9d10212a
......@@ -110,7 +110,7 @@ namespace AMDiS
{
// NOTE: This needs a change in mtl dense_vector implementation, i.e. we
// have to add an alias for the function change_dim with the name resize.
Dune::Functions::interpolate(feSpace, *vector, f);
Dune::Functions::interpolate(feSpace, wrapper(*vector), f);
}
/// Calls the copy assignment operator of the BaseVector \ref vector
......
dune/amdis/test/CMakeLists.txt
View file @ 9d10212a
......@@ -3,7 +3,7 @@ dune_add_test(NAME test_duneamdis
TARGET duneamdis
COMPILE_ONLY)
find_package(GTest REQUIRED)
# find_package(GTest REQUIRED)
set(projects "test1" "test2")
foreach(project ${projects})
......
dune/amdis/test/test1.cc
View file @ 9d10212a
......@@ -2,7 +2,7 @@
#include <dune/amdis/ProblemStat.hpp>
#include <dune/amdis/Literals.hpp>
#include <gtest/gtest.h>
#define TEST_EQ(expr1, expr2) AMDIS_TEST_EXIT( (expr1) == (expr2), "'" << #expr1 << "' != '" << #expr2 << "'" )
#ifndef DIM
#define DIM 2
......@@ -21,18 +21,6 @@ namespace {
using Geometry = typename MeshView::template Codim<0>::Geometry;
using RefElements = Dune::ReferenceElements<typename Geometry::ctype, Geometry::mydimension>;
class AMDiSTester : public ::testing::Test
{
protected:
// create a problemStat
virtual void SetUp() override
{
prob = std::make_shared<TestProblem>("test");
prob->initialize(INIT_ALL);
}
std::shared_ptr<TestProblem> prob;
};
template <class FeSpace>
......@@ -76,65 +64,80 @@ namespace {
};
TEST_F(AMDiSTester, CreateDOFVector)
class AMDiSTester
{
using FeSpace = typename TestProblem::template FeSpace<0>;
// construction of DOFVectors
DOFVector<FeSpace, double> vec1(prob->getFeSpace(0_c), "vec1");
auto vec2 = make_dofvector(prob->getFeSpace(0_c), "vec2");
// retrieving DOFVectors from problemStat
auto vec3 = prob->getSolution(0_c);
auto vec5 = prob->getSolution<0>();
auto vec6 = prob->getSolution(index_<0>());
auto vec7 = prob->getSolution()->getDOFVector(0_c);
auto vec8 = prob->getSolution()->getDOFVector<0>();
auto vec9 = prob->getSolution()->getDOFVector(index_<0>());
auto& sys_vec = *prob->getSolution();
auto vec10 = sys_vec[0_c];
auto vec11 = sys_vec[index_<0>()];
// construction of SystemVector
using FeSpaces = typename TestProblem::FeSpaces;
SystemVector<FeSpaces, double> sys_vec1(*prob->getFeSpaces(), prob->getComponentNames());
auto sys_vec2 = make_systemvector(*prob->getFeSpaces(), prob->getComponentNames());
auto sys_vec3 = make_systemvector(*prob->getFeSpaces(), "sys_vec");
public:
// create a problemStat
AMDiSTester(std::string name)
: prob(std::make_shared<TestProblem>(name))
{
prob->initialize(INIT_ALL);
}
// retrieving systemVector from problemStat
auto sys_vec4 = *prob->getSolution();
}
void test1() // CreateDOFVector
{
using FeSpace = typename TestProblem::template FeSpace<0>;
// construction of DOFVectors
DOFVector<FeSpace, double> vec1(prob->getFeSpace(0_c), "vec1");
auto vec2 = make_dofvector(prob->getFeSpace(0_c), "vec2");
// retrieving DOFVectors from problemStat
auto vec3 = prob->getSolution(0_c);
auto vec5 = prob->getSolution<0>();
auto vec6 = prob->getSolution(index_<0>());
auto vec7 = prob->getSolution()->getDOFVector(0_c);
auto vec8 = prob->getSolution()->getDOFVector<0>();
auto vec9 = prob->getSolution()->getDOFVector(index_<0>());
auto& sys_vec = *prob->getSolution();
auto vec10 = sys_vec[0_c];
auto vec11 = sys_vec[index_<0>()];
// construction of SystemVector
using FeSpaces = typename TestProblem::FeSpaces;
SystemVector<FeSpaces, double> sys_vec1(*prob->getFeSpaces(), prob->getComponentNames());
auto sys_vec2 = make_systemvector(*prob->getFeSpaces(), prob->getComponentNames());
auto sys_vec3 = make_systemvector(*prob->getFeSpaces(), "sys_vec");
// retrieving systemVector from problemStat
auto sys_vec4 = *prob->getSolution();
}
TEST_F(AMDiSTester, FillDOFVector)
{
auto& vec = prob->getSolution(0_c);
auto const& feSpace = vec.getFeSpace();
// interpolate function to dofvector
vec.interpol([](auto const& x) { return x*x; });
// test whether interpolation is fine
using FeSpace = std::decay_t< decltype(feSpace) >;
LocalEvaluation<FeSpace> evaluator(feSpace);
for (auto const& element : elements(feSpace.gridView())) {
evaluator.bind(element);
void test2() // FillDOFVector
{
auto& vec = prob->getSolution(0_c);
auto const& feSpace = vec.getFeSpace();
auto const& refElement = RefElements::general(element.type());
// interpolate function to dofvector
vec.interpol([](auto const& x) { return x*x; });
size_t nVertices = element.subEntities(DIM);
for (size_t i = 0; i < nVertices; ++i) {
auto pos = refElement.position(i, DIM);
auto data = evaluator.eval(vec, pos);
// test whether interpolation is fine
using FeSpace = std::decay_t< decltype(feSpace) >;
LocalEvaluation<FeSpace> evaluator(feSpace);
for (auto const& element : elements(feSpace.gridView())) {
evaluator.bind(element);
auto const& refElement = RefElements::general(element.type());
auto x = element.geometry().global(pos);
EXPECT_EQ(data, x*x);
size_t nVertices = element.subEntities(DIM);
for (size_t i = 0; i < nVertices; ++i) {
auto pos = refElement.position(i, DIM);
auto data = evaluator.eval(vec, pos);
auto x = element.geometry().global(pos);
TEST_EQ(data, x*x);
}
}
}
}
protected:
std::shared_ptr<TestProblem> prob;
};
} // end namespace
......@@ -142,9 +145,10 @@ int main(int argc, char** argv)
{
AMDiS::init(argc, argv);
::testing::InitGoogleTest(&argc, argv);
int error_code = RUN_ALL_TESTS();
auto tester = std::make_shared<AMDiSTester>("test");
tester->test1();
tester->test2();
AMDiS::finalize();
return error_code;
return 0;
}
\ No newline at end of file
dune/amdis/test/test2.cc
View file @ 9d10212a
......@@ -2,8 +2,6 @@
#include <dune/amdis/ProblemStat.hpp>
#include <dune/amdis/Literals.hpp>
#include <gtest/gtest.h>
#ifndef DIM
#define DIM 2
#endif
......@@ -12,92 +10,96 @@
#endif
using namespace AMDiS;
namespace {
using TestParam = ProblemStatTraits<DIM, DOW, 1>;
using TestProblem = ProblemStat<TestParam>;
using Op = TestProblem::OperatorType;
class AMDiSTester : public ::testing::Test
class AMDiSTester
{
protected:
public:
// create a problemStat
virtual void SetUp() override
AMDiSTester(std::string name)
: prob(std::make_shared<TestProblem>(name))
{
prob = std::make_shared<TestProblem>("test");
prob->initialize(INIT_ALL);
}
std::shared_ptr<TestProblem> prob;
};
TEST_F(AMDiSTester, AddOperators)
{
auto const& u = prob->getSolution<0>();
// define operators
Op op1;
op1.addZOT( 1.0 );
op1.addZOT( constant(1.0) );
// value of DOFVector at QPs
op1.addZOT( valueOf(u) );
op1.addZOT( valueOfFunc(u, [](auto const& U) { return U; }) );
op1.addZOT( valueOfFunc(u, [](auto const& U) { return U; }, 1) );
// partial derivative of DOFVector
op1.addZOT( derivativeOf(u, 0) );
op1.addZOT( derivativeOf(u, 1) );
// op1.addFOT( gradientOf(u), GRD_PSI );
// evaluate functor at coordinates
op1.addZOT( eval([](auto const& x) { return x[0]; }) );
op1.addZOT( [](auto const& x) { return 1.0; } );
// add operator to system
prob->addMatrixOperator(op1, 0, 0);
prob->addVectorOperator(op1, 0);
// define operators with factor argument
Op op2;
op2.addZOT( valueOf(u, 1.0) );
op2.addZOT( derivativeOf(u, 0, 1.0) );
op2.addZOT( derivativeOf(u, 1, 1.0) );
// add operator to system, with optional factor as pointer
double* factor = new double(1.0);
prob->addMatrixOperator(op2, 0, 0, factor);
prob->addVectorOperator(op2, 0, factor);
}
void test1() // AddOperators
{
auto const& u = prob->getSolution<0>();
// define operators
Op op1;
op1.addZOT( 1.0 );
op1.addZOT( constant(1.0) );
// value of DOFVector at QPs
op1.addZOT( valueOf(u) );
op1.addZOT( valueOfFunc(u, [](auto const& U) { return U; }) );
op1.addZOT( valueOfFunc(u, [](auto const& U) { return U; }, 1) );
// partial derivative of DOFVector
op1.addZOT( derivativeOf(u, 0) );
op1.addZOT( derivativeOf(u, 1) );
// op1.addFOT( gradientOf(u), GRD_PSI );
// evaluate functor at coordinates
op1.addZOT( eval([](auto const& x) { return x[0]; }) );
op1.addZOT( [](auto const& x) { return 1.0; } );
// add operator to system
prob->addMatrixOperator(op1, 0, 0);
prob->addVectorOperator(op1, 0);
// define operators with factor argument
Op op2;
op2.addZOT( valueOf(u, 1.0) );
op2.addZOT( derivativeOf(u, 0, 1.0) );
op2.addZOT( derivativeOf(u, 1, 1.0) );
// add operator to system, with optional factor as pointer
double* factor = new double(1.0);
prob->addMatrixOperator(op2, 0, 0, factor);
prob->addVectorOperator(op2, 0, factor);
}
TEST_F(AMDiSTester, MatrixVector)
{
auto const& u = prob->getSolution<0>();
void test2() // MatrixVector
{
auto const& u = prob->getSolution<0>();
using GradOp = std::decay_t< decltype(gradientOf(u)) >;
using DerivOp = std::decay_t< decltype(derivativeOf(u,0)) >;
using Jacobian = typename GradOp::value_type;
using Derivative = typename DerivOp::value_type;
Jacobian v, w;
double erg = v[0].dot(w[0]);
}
using GradOp = std::decay_t< decltype(gradientOf(u)) >;
using DerivOp = std::decay_t< decltype(derivativeOf(u,0)) >;
using Jacobian = typename GradOp::value_type;
using Derivative = typename DerivOp::value_type;
Jacobian v, w;
double erg = v[0].dot(w[0]);
}
protected:
std::shared_ptr<TestProblem> prob;
};
} // end namespace
int main(int argc, char** argv)
{
int main(int argc, char* argv[])
{
AMDiS::init(argc, argv);
::testing::InitGoogleTest(&argc, argv);
int error_code = RUN_ALL_TESTS();
auto tester = std::make_shared<AMDiSTester>("test");
tester->test1();
tester->test2();
AMDiS::finalize();
return error_code;
return 0;
}
\ No newline at end of file
init/heat.json.2d
View file @ 9d10212a
......@@ -33,7 +33,7 @@
},
"adapt": {
"timestep": 0.1,
"timestep": 0.001,
"start time": 0.0,
"end time": 1.0
}
......
init/navier_stokes.json.2d
View file @ 9d10212a
......@@ -10,11 +10,11 @@
"mesh": "stokesMesh",
"solver" : {
"name": "gmres",
"max iteration": 500,
"name": "bicgstab_ell",
"ell": 3,
"max iteration": 5000,
"reduction": 1e-6,
"info": 1,
"left precon": "diag"
"info": 1
},
"output": {
......@@ -29,7 +29,7 @@
},
"adapt": {
"timestep": 0.001,
"timestep": 1.0e-3,
"start time": 0.0,
"end time": 1.0
}
......
src/navier_stokes.cc
View file @ 9d10212a
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <iostream>
#include <ctime>
#include <cmath>
#include <dune/amdis/AMDiS.hpp>
#include <dune/amdis/AdaptInstationary.hpp>
#include <dune/amdis/ProblemStat.hpp>
#include <dune/amdis/ProblemInstat.hpp>
#ifndef DIM
#define DIM 2
......@@ -21,6 +19,7 @@ using namespace AMDiS;
// 3 components: velocity with polynomial degree 2 and pressure with polynomial degree 1
using StokesParam = ProblemStatTraits<DIM, DOW, 2, 2, 1>;
using StokesProblem = ProblemStat<StokesParam>;
using StokesProblemInstat = ProblemInstat<StokesParam>;
int main(int argc, char** argv)
{
......@@ -28,6 +27,9 @@ int main(int argc, char** argv)
StokesProblem prob("stokes");
prob.initialize(INIT_ALL);
StokesProblemInstat probInstat("stokes", prob);
probInstat.initialize(INIT_UH_OLD);
double viscosity = 1.0;
double density = 1.0;
......@@ -36,43 +38,55 @@ int main(int argc, char** argv)
Parameters::get("stokes->density", density);
Parameters::get("stokes->boundary velocity", vel);
using DOFVectorU = std::decay_t< decltype(prob.getSolution<0>()) >;
using DOFVectorP = std::decay_t< decltype(prob.getSolution<DOW>()) >;
std::array<DOFVectorU*, DOW> velocity;
For<0,DOW>::loop([&](auto const _i) { velocity[_i] = &prob.getSolution(_i); });
AdaptInfo adaptInfo("adapt");
double tau = adaptInfo.getTimestep();
// define the differential operators
using Op = StokesProblem::OperatorType;
For<0,DOW>::loop([&](auto const _i)
{
// <1/tau * u_i, v_i> + <viscosity*grad(u_i), grad(v_i)>
Op* opL = new Op;
opL->addZOT( constant(density/tau) );
opL->addSOT( constant(viscosity) );
prob.addMatrixOperator(*opL, _i, _i);
// <1/tau * u_i^old, v_i>
Op* opRhs = new Op;
opRhs->addZOT( valueOf(prob.getSolution(_i), density/tau) );
prob.addVectorOperator(*opRhs, _i);
// <1/tau * u_i, v_i>
auto opTime = Op::create(), opTimeOld = Op::create();
opTime->addZOT( density );
opTimeOld->addZOT( valueOf(prob.getSolution(_i), density) );
prob.addMatrixOperator(opTime, _i, _i, probInstat.getInvTau());
prob.addVectorOperator(opTimeOld, _i, probInstat.getInvTau());
#if 0
// <(u * nabla)u_i^old, v_i>
For<0, DOW>::loop([&](auto const _j)
{
auto opNonlin = Op::create();
opNonlin->addZOT( derivativeOf(prob.getSolution(_i), _j, density) );
prob.addMatrixOperator(opNonlin, _i, _j);
});
#else
// <(u^old * nabla)u_i, v_i>
For<0, DOW>::loop([&](auto const _j)
{
auto opNonlin = Op::create();
opNonlin->addFOT( valueOf(prob.getSolution(_j), density), _j, GRD_PHI );
prob.addMatrixOperator(opNonlin, _i, _i);
});
#endif
// <viscosity*grad(u_i), grad(v_i)>
auto opL = Op::create();
opL->addSOT( viscosity );
prob.addMatrixOperator(opL, _i, _i);
// <p, d_i(v_i)>
Op* opB = new Op;
opB->addFOT( constant(1.0), _i, GRD_PSI );
prob.addMatrixOperator(*opB, _i, DOW);
auto opP = Op::create();
opP->addFOT( 1.0, _i, GRD_PSI );
prob.addMatrixOperator(opP, _i, DOW);
// <d_i(u_i), q>
Op* opDiv = new Op;
opDiv->addFOT( constant(1.0), _i, GRD_PHI );
prob.addMatrixOperator(*opDiv, DOW, _i);
// <(u^old * nabla)u_i, v_i>
// NOTE: only simple linearization of nonlin-term, since gradientOf() not yet implemented
For<0, DOW>::loop([&](auto const _j)
{
Op* opNonlin = new Op;
opNonlin->addFOT( valueOf(prob.getSolution(_j), density), _j, GRD_PHI );
prob.addMatrixOperator(*opNonlin, _i, _i);
});
auto opDiv = Op::create();
opDiv->addFOT( 1.0, _i, GRD_PHI );
prob.addMatrixOperator(opDiv, DOW, _i);
});
// define boundary regions
......@@ -90,23 +104,10 @@ int main(int argc, char** argv)
prob.addDirichletBC(left, 0, 0, zero);
prob.addDirichletBC(left, 1, 1, parabolic);
// set initial values for the solver (maybe not necessary)
*prob.getSolution() = 0.0;
double t = 0.0;
for (t = adaptInfo.getStartTime();
t < adaptInfo.getEndTime();
t+= adaptInfo.getTimestep())
{
adaptInfo.setTime(t);
AMDIS_MSG("Timestep t = " << t);
prob.writeFiles(adaptInfo);
// assemble and solve system
prob.buildAfterCoarsen(adaptInfo, Flag(0));
prob.solve(adaptInfo);
}
AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo);
adapt.adapt();
// output solution
prob.writeFiles(adaptInfo);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment