Skip to content

GitLab

Commit 50f91d96 authored by Timo Koch's avatar Timo Koch
Browse files

[column] Introduce epsilons. Shorten test. 

Reduce the end time of the simulation when testing.
*ALWAYS* use epsilons for floating point comparisons!
parent eecfcdde
Expand all Hide whitespace changes
Inline Side-by-side
lecture/mm/columnxylene/CMakeLists.txt
View file @ 50f91d96
......@@ -4,8 +4,8 @@ add_dumux_test(columnxyleneexercise columnxyleneexercise columnxyleneexercise.cc
python ${dumux_INCLUDE_DIRS}/bin/runtest.py
--script fuzzy
--files ${CMAKE_SOURCE_DIR}/lecture/references/columnxylol-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/columnxylol-00039.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/columnxyleneexercise -ParameterFile ${CMAKE_CURRENT_SOURCE_DIR}/columnxyleneexercise.input -Grid.File ${CMAKE_CURRENT_SOURCE_DIR}/grids/column.dgf -TimeManager.OutputInterval 100")
${CMAKE_CURRENT_BINARY_DIR}/columnxylol-00008.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/columnxyleneexercise -Grid.File ${CMAKE_CURRENT_SOURCE_DIR}/grids/column.dgf -TimeManager.OutputInterval 100 -TimeManager.TEnd 600")
# headers for installation and headercheck
install(FILES
......
lecture/mm/columnxylene/columnxyleneexercise.input
View file @ 50f91d96
......@@ -2,8 +2,9 @@
[TimeManager]
DtInitial = 1 # [s]
TEnd = 3600 # [s]
TEnd = 1800 # [s]
OutputInterval = 1
MaxTimeStepSize = 1
[Grid]
File = ./grids/column.dgf
......
lecture/mm/columnxylene/columnxyleneproblem.hh
View file @ 50f91d96
......@@ -31,8 +31,6 @@
#include <dumux/porousmediumflow/3p3c/implicit/model.hh>
#include <dumux/porousmediumflow/implicit/problem.hh>
//#include <dumux/common/deprecated.hh>
#include "columnxylenespatialparams.hh"
#define ISOTHERMAL 0
......@@ -60,20 +58,8 @@ SET_TYPE_PROP(ColumnProblem,
// Enable gravity
SET_BOOL_PROP(ColumnProblem, ProblemEnableGravity, true);
// Use forward differences instead of central differences
SET_INT_PROP(ColumnProblem, ImplicitNumericDifferenceMethod, 1);
// Maximum tolerated relative error in the Newton method
SET_SCALAR_PROP(ColumnProblem, NewtonMaxRelativeShift, 1e-4);
// Write newton convergence
SET_BOOL_PROP(ColumnProblem, NewtonWriteConvergence, false);
// Use line search
SET_BOOL_PROP(ColumnProblem, NewtonUseLineSearch, true);
// Set the maximum time step
SET_SCALAR_PROP(ColumnProblem, TimeManagerMaxTimeStepSize, 1.);
}
......@@ -296,39 +282,39 @@ private:
void initial_(PrimaryVariables &values,
const GlobalPosition &globalPos) const
{
Scalar y = globalPos[1];
const Scalar y = globalPos[1];
values[temperatureIdx] = 296.15;
values[pressureIdx] = 1.e5;
if(y > 1.2-eps_){
if(y > 1.2 - eps_){
values[switch2Idx] = 0.112; // almost no contaminant component
values[switch1Idx] = 0.005;
} else if(y < 1.2-0.3){ // extended domain
} else if(y < 1.2-0.3 + eps_){ // extended domain
values[switch2Idx] = 1.e-4; // almost no contaminant component
values[switch1Idx] = 0.005;
} else {
values[switch1Idx] = 0.005;
if((y<=1.2-0.001)&&(y>=1.2-0.0148)) values[switch2Idx] = 0+((1.2-y)/0.0148)*0.112;
if((y<1.2-0.0148)&&(y>=1.2-0.0296)) values[switch2Idx] = 0.112+(((1.2-y)-0.0148)/0.0148)*(0.120-0.112);
if((y<1.2-0.0296)&&(y>=1.2-0.0444)) values[switch2Idx] = 0.120+(((1.2-y)-0.0296)/0.0148)*(0.125-0.120);
if((y<1.2-0.0444)&&(y>=1.2-0.0592)) values[switch2Idx] = 0.125+(((1.2-y)-0.0444)/0.0148)*(0.137-0.125);
if((y<1.2-0.0592)&&(y>=1.2-0.0740)) values[switch2Idx] = 0.137+(((1.2-y)-0.0592)/0.0148)*(0.150-0.137);
if((y<1.2-0.0740)&&(y>=1.2-0.0888)) values[switch2Idx] = 0.150+(((1.2-y)-0.0740)/0.0148)*(0.165-0.150);
if((y<1.2-0.0888)&&(y>=1.2-0.1036)) values[switch2Idx] = 0.165+(((1.2-y)-0.0888)/0.0148)*(0.182-0.165);
if((y<1.2-0.1036)&&(y>=1.2-0.1184)) values[switch2Idx] = 0.182+(((1.2-y)-0.1036)/0.0148)*(0.202-0.182);
if((y<1.2-0.1184)&&(y>=1.2-0.1332)) values[switch2Idx] = 0.202+(((1.2-y)-0.1184)/0.0148)*(0.226-0.202);
if((y<1.2-0.1332)&&(y>=1.2-0.1480)) values[switch2Idx] = 0.226+(((1.2-y)-0.1332)/0.0148)*(0.257-0.226);
if((y<1.2-0.1480)&&(y>=1.2-0.1628)) values[switch2Idx] = 0.257+(((1.2-y)-0.1480)/0.0148)*(0.297-0.257);
if((y<1.2-0.1628)&&(y>=1.2-0.1776)) values[switch2Idx] = 0.297+(((1.2-y)-0.1628)/0.0148)*(0.352-0.297);
if((y<1.2-0.1776)&&(y>=1.2-0.1924)) values[switch2Idx] = 0.352+(((1.2-y)-0.1776)/0.0148)*(0.426-0.352);
if((y<1.2-0.1924)&&(y>=1.2-0.2072)) values[switch2Idx] = 0.426+(((1.2-y)-0.1924)/0.0148)*(0.522-0.426);
if((y<1.2-0.2072)&&(y>=1.2-0.2220)) values[switch2Idx] = 0.522+(((1.2-y)-0.2072)/0.0148)*(0.640-0.522);
if((y<1.2-0.2220)&&(y>=1.2-0.2368)) values[switch2Idx] = 0.640+(((1.2-y)-0.2220)/0.0148)*(0.767-0.640);
if((y<1.2-0.2368)&&(y>=1.2-0.2516)) values[switch2Idx] = 0.767+(((1.2-y)-0.2368)/0.0148)*(0.878-0.767);
if((y<1.2-0.2516)&&(y>=1.2-0.2664)) values[switch2Idx] = 0.878+(((1.2-y)-0.2516)/0.0148)*(0.953-0.878);
if((y<1.2-0.2664)&&(y>=1.2-0.2812)) values[switch2Idx] = 0.953+(((1.2-y)-0.2664)/0.0148)*(0.988-0.953);
if((y<1.2-0.2812)&&(y>=1.2-0.3000)) values[switch2Idx] = 0.988;
if((y<1.2-0.001+eps_)&&(y>1.2-0.0148-eps_)) values[switch2Idx] = 0+((1.2-y)/0.0148)*0.112;
if((y<1.2-0.0148-eps_)&&(y>1.2-0.0296-eps_)) values[switch2Idx] = 0.112+(((1.2-y)-0.0148)/0.0148)*(0.120-0.112);
if((y<1.2-0.0296-eps_)&&(y>1.2-0.0444-eps_)) values[switch2Idx] = 0.120+(((1.2-y)-0.0296)/0.0148)*(0.125-0.120);
if((y<1.2-0.0444-eps_)&&(y>1.2-0.0592-eps_)) values[switch2Idx] = 0.125+(((1.2-y)-0.0444)/0.0148)*(0.137-0.125);
if((y<1.2-0.0592-eps_)&&(y>1.2-0.0740-eps_)) values[switch2Idx] = 0.137+(((1.2-y)-0.0592)/0.0148)*(0.150-0.137);
if((y<1.2-0.0740-eps_)&&(y>1.2-0.0888-eps_)) values[switch2Idx] = 0.150+(((1.2-y)-0.0740)/0.0148)*(0.165-0.150);
if((y<1.2-0.0888-eps_)&&(y>1.2-0.1036-eps_)) values[switch2Idx] = 0.165+(((1.2-y)-0.0888)/0.0148)*(0.182-0.165);
if((y<1.2-0.1036-eps_)&&(y>1.2-0.1184-eps_)) values[switch2Idx] = 0.182+(((1.2-y)-0.1036)/0.0148)*(0.202-0.182);
if((y<1.2-0.1184-eps_)&&(y>1.2-0.1332-eps_)) values[switch2Idx] = 0.202+(((1.2-y)-0.1184)/0.0148)*(0.226-0.202);
if((y<1.2-0.1332-eps_)&&(y>1.2-0.1480-eps_)) values[switch2Idx] = 0.226+(((1.2-y)-0.1332)/0.0148)*(0.257-0.226);
if((y<1.2-0.1480-eps_)&&(y>1.2-0.1628-eps_)) values[switch2Idx] = 0.257+(((1.2-y)-0.1480)/0.0148)*(0.297-0.257);
if((y<1.2-0.1628-eps_)&&(y>1.2-0.1776-eps_)) values[switch2Idx] = 0.297+(((1.2-y)-0.1628)/0.0148)*(0.352-0.297);
if((y<1.2-0.1776-eps_)&&(y>1.2-0.1924-eps_)) values[switch2Idx] = 0.352+(((1.2-y)-0.1776)/0.0148)*(0.426-0.352);
if((y<1.2-0.1924-eps_)&&(y>1.2-0.2072-eps_)) values[switch2Idx] = 0.426+(((1.2-y)-0.1924)/0.0148)*(0.522-0.426);
if((y<1.2-0.2072-eps_)&&(y>1.2-0.2220-eps_)) values[switch2Idx] = 0.522+(((1.2-y)-0.2072)/0.0148)*(0.640-0.522);
if((y<1.2-0.2220-eps_)&&(y>1.2-0.2368-eps_)) values[switch2Idx] = 0.640+(((1.2-y)-0.2220)/0.0148)*(0.767-0.640);
if((y<1.2-0.2368-eps_)&&(y>1.2-0.2516-eps_)) values[switch2Idx] = 0.767+(((1.2-y)-0.2368)/0.0148)*(0.878-0.767);
if((y<1.2-0.2516-eps_)&&(y>1.2-0.2664-eps_)) values[switch2Idx] = 0.878+(((1.2-y)-0.2516)/0.0148)*(0.953-0.878);
if((y<1.2-0.2664-eps_)&&(y>1.2-0.2812-eps_)) values[switch2Idx] = 0.953+(((1.2-y)-0.2664)/0.0148)*(0.988-0.953);
if((y<1.2-0.2812-eps_)&&(y>1.2-0.3000-eps_)) values[switch2Idx] = 0.988;
}
}
......
lecture/mm/columnxylene/columnxyleneproblem2.hh
View file @ 50f91d96
......@@ -63,18 +63,9 @@ SET_TYPE_PROP(ColumnProblem,
// Enable gravity
SET_BOOL_PROP(ColumnProblem, ProblemEnableGravity, true);
// Use forward differences instead of central differences
SET_INT_PROP(ColumnProblem, ImplicitNumericDifferenceMethod, 1);
// Maximum tolerated relative error in the Newton method
SET_SCALAR_PROP(ColumnProblem, NewtonMaxRelativeShift, 1e-4);
// Write newton convergence
SET_BOOL_PROP(ColumnProblem, NewtonWriteConvergence, false);
// Use line search
SET_BOOL_PROP(ColumnProblem, NewtonUseLineSearch, true);
// Set the maximum time step
SET_SCALAR_PROP(ColumnProblem, TimeManagerMaxTimeStepSize, 1.);
}
......@@ -280,52 +271,9 @@ public:
}
bool shouldWriteRestartFile() const
{
return 0;
}
/*
bool shouldWriteOutput() const
{
return
this->timeManager().timeStepIndex() == 0 ||
this->timeManager().timeStepIndex() == 10 ||
this->timeManager().timeStepIndex() == 20 ||
this->timeManager().timeStepIndex() == 50 ||
this->timeManager().timeStepIndex() == 100 ||
this->timeManager().timeStepIndex() == 200 ||
this->timeManager().timeStepIndex() == 300 ||
this->timeManager().timeStepIndex() == 400 ||
this->timeManager().timeStepIndex() == 500 ||
this->timeManager().timeStepIndex() == 600 ||
this->timeManager().timeStepIndex() == 700 ||
this->timeManager().timeStepIndex() == 800 ||
this->timeManager().timeStepIndex() == 900 ||
this->timeManager().timeStepIndex() == 1000 ||
this->timeManager().timeStepIndex() == 1100 ||
this->timeManager().timeStepIndex() == 1200 ||
this->timeManager().timeStepIndex() == 1300 ||
this->timeManager().timeStepIndex() == 1400 ||
this->timeManager().timeStepIndex() == 1500 ||
this->timeManager().timeStepIndex() == 1600 ||
this->timeManager().timeStepIndex() == 1700 ||
this->timeManager().timeStepIndex() == 1800 ||
this->timeManager().timeStepIndex() == 1900 ||
this->timeManager().timeStepIndex() == 2000 ||
this->timeManager().timeStepIndex() == 2100 ||
this->timeManager().timeStepIndex() == 2200 ||
this->timeManager().timeStepIndex() == 2300 ||
this->timeManager().timeStepIndex() == 2400 ||
this->timeManager().timeStepIndex() == 2500 ||
this->timeManager().timeStepIndex() == 2600 ||
this->timeManager().timeStepIndex() == 2700 ||
this->timeManager().timeStepIndex() == 2800 ||
this->timeManager().timeStepIndex() == 2900 ||
this->timeManager().willBeFinished();
}
*/
{
return 0;
}
private:
// internal method for the initial condition (reused for the
......@@ -343,26 +291,26 @@ private:
values[switch1Idx] = 0.005;
} else {
values[switch1Idx] = 0.005;
if((y<=0.3-0.001)&&(y>=0.3-0.0148)) values[switch2Idx] = 0+((0.3-y)/0.0148)*0.112;
if((y<0.3-0.0148)&&(y>=0.3-0.0296)) values[switch2Idx] = 0.112+(((0.3-y)-0.0148)/0.0148)*(0.120-0.112);
if((y<0.3-0.0296)&&(y>=0.3-0.0444)) values[switch2Idx] = 0.120+(((0.3-y)-0.0296)/0.0148)*(0.125-0.120);
if((y<0.3-0.0444)&&(y>=0.3-0.0592)) values[switch2Idx] = 0.125+(((0.3-y)-0.0444)/0.0148)*(0.137-0.125);
if((y<0.3-0.0592)&&(y>=0.3-0.0740)) values[switch2Idx] = 0.137+(((0.3-y)-0.0592)/0.0148)*(0.150-0.137);
if((y<0.3-0.0740)&&(y>=0.3-0.0888)) values[switch2Idx] = 0.150+(((0.3-y)-0.0740)/0.0148)*(0.165-0.150);
if((y<0.3-0.0888)&&(y>=0.3-0.1036)) values[switch2Idx] = 0.165+(((0.3-y)-0.0888)/0.0148)*(0.182-0.165);
if((y<0.3-0.1036)&&(y>=0.3-0.1184)) values[switch2Idx] = 0.182+(((0.3-y)-0.1036)/0.0148)*(0.202-0.182);
if((y<0.3-0.1184)&&(y>=0.3-0.1332)) values[switch2Idx] = 0.202+(((0.3-y)-0.1184)/0.0148)*(0.226-0.202);
if((y<0.3-0.1332)&&(y>=0.3-0.1480)) values[switch2Idx] = 0.226+(((0.3-y)-0.1332)/0.0148)*(0.257-0.226);
if((y<0.3-0.1480)&&(y>=0.3-0.1628)) values[switch2Idx] = 0.257+(((0.3-y)-0.1480)/0.0148)*(0.297-0.257);
if((y<0.3-0.1628)&&(y>=0.3-0.1776)) values[switch2Idx] = 0.297+(((0.3-y)-0.1628)/0.0148)*(0.352-0.297);
if((y<0.3-0.1776)&&(y>=0.3-0.1924)) values[switch2Idx] = 0.352+(((0.3-y)-0.1776)/0.0148)*(0.426-0.352);
if((y<0.3-0.1924)&&(y>=0.3-0.2072)) values[switch2Idx] = 0.426+(((0.3-y)-0.1924)/0.0148)*(0.522-0.426);
if((y<0.3-0.2072)&&(y>=0.3-0.2220)) values[switch2Idx] = 0.522+(((0.3-y)-0.2072)/0.0148)*(0.640-0.522);
if((y<0.3-0.2220)&&(y>=0.3-0.2368)) values[switch2Idx] = 0.640+(((0.3-y)-0.2220)/0.0148)*(0.767-0.640);
if((y<0.3-0.2368)&&(y>=0.3-0.2516)) values[switch2Idx] = 0.767+(((0.3-y)-0.2368)/0.0148)*(0.878-0.767);
if((y<0.3-0.2516)&&(y>=0.3-0.2664)) values[switch2Idx] = 0.878+(((0.3-y)-0.2516)/0.0148)*(0.953-0.878);
if((y<0.3-0.2664)&&(y>=0.3-0.2812)) values[switch2Idx] = 0.953+(((0.3-y)-0.2664)/0.0148)*(0.988-0.953);
if((y<0.3-0.2812)&&(y>=0.3-0.3000)) values[switch2Idx] = 0.988;
if((y<=0.3-0.001+eps_)&&(y>0.3-0.0148-eps_)) values[switch2Idx] = 0+((0.3-y)/0.0148)*0.112;
if((y<0.3-0.0148-eps_)&&(y>0.3-0.0296-eps_)) values[switch2Idx] = 0.112+(((0.3-y)-0.0148)/0.0148)*(0.120-0.112);
if((y<0.3-0.0296-eps_)&&(y>0.3-0.0444-eps_)) values[switch2Idx] = 0.120+(((0.3-y)-0.0296)/0.0148)*(0.125-0.120);
if((y<0.3-0.0444-eps_)&&(y>0.3-0.0592-eps_)) values[switch2Idx] = 0.125+(((0.3-y)-0.0444)/0.0148)*(0.137-0.125);
if((y<0.3-0.0592-eps_)&&(y>0.3-0.0740-eps_)) values[switch2Idx] = 0.137+(((0.3-y)-0.0592)/0.0148)*(0.150-0.137);
if((y<0.3-0.0740-eps_)&&(y>0.3-0.0888-eps_)) values[switch2Idx] = 0.150+(((0.3-y)-0.0740)/0.0148)*(0.165-0.150);
if((y<0.3-0.0888-eps_)&&(y>0.3-0.1036-eps_)) values[switch2Idx] = 0.165+(((0.3-y)-0.0888)/0.0148)*(0.182-0.165);
if((y<0.3-0.1036-eps_)&&(y>0.3-0.1184-eps_)) values[switch2Idx] = 0.182+(((0.3-y)-0.1036)/0.0148)*(0.202-0.182);
if((y<0.3-0.1184-eps_)&&(y>0.3-0.1332-eps_)) values[switch2Idx] = 0.202+(((0.3-y)-0.1184)/0.0148)*(0.226-0.202);
if((y<0.3-0.1332-eps_)&&(y>0.3-0.1480-eps_)) values[switch2Idx] = 0.226+(((0.3-y)-0.1332)/0.0148)*(0.257-0.226);
if((y<0.3-0.1480-eps_)&&(y>0.3-0.1628-eps_)) values[switch2Idx] = 0.257+(((0.3-y)-0.1480)/0.0148)*(0.297-0.257);
if((y<0.3-0.1628-eps_)&&(y>0.3-0.1776-eps_)) values[switch2Idx] = 0.297+(((0.3-y)-0.1628)/0.0148)*(0.352-0.297);
if((y<0.3-0.1776-eps_)&&(y>0.3-0.1924-eps_)) values[switch2Idx] = 0.352+(((0.3-y)-0.1776)/0.0148)*(0.426-0.352);
if((y<0.3-0.1924-eps_)&&(y>0.3-0.2072-eps_)) values[switch2Idx] = 0.426+(((0.3-y)-0.1924)/0.0148)*(0.522-0.426);
if((y<0.3-0.2072-eps_)&&(y>0.3-0.2220-eps_)) values[switch2Idx] = 0.522+(((0.3-y)-0.2072)/0.0148)*(0.640-0.522);
if((y<0.3-0.2220-eps_)&&(y>0.3-0.2368-eps_)) values[switch2Idx] = 0.640+(((0.3-y)-0.2220)/0.0148)*(0.767-0.640);
if((y<0.3-0.2368-eps_)&&(y>0.3-0.2516-eps_)) values[switch2Idx] = 0.767+(((0.3-y)-0.2368)/0.0148)*(0.878-0.767);
if((y<0.3-0.2516-eps_)&&(y>0.3-0.2664-eps_)) values[switch2Idx] = 0.878+(((0.3-y)-0.2516)/0.0148)*(0.953-0.878);
if((y<0.3-0.2664-eps_)&&(y>0.3-0.2812-eps_)) values[switch2Idx] = 0.953+(((0.3-y)-0.2664)/0.0148)*(0.988-0.953);
if((y<0.3-0.2812-eps_)&&(y>0.3-0.3000-eps_)) values[switch2Idx] = 0.988;
}
}
......
lecture/mm/columnxylene/columnxylenespatialparams.hh
View file @ 50f91d96
......@@ -154,20 +154,6 @@ public:
fineMaterialParams_.setRhoBulk(1500.);
}
~ColumnSpatialParams()
{}
/*!
* \brief Update the spatial parameters with the flow solution
* after a timestep.
*
* \param globalSolution The global solution vector
*/
void update(const SolutionVector &globalSolution)
{
};
/*!
* \brief Apply the intrinsic permeability tensor to a pressure
* potential gradient.
......@@ -278,9 +264,14 @@ public:
private:
bool isFineMaterial_(const GlobalPosition &pos) const
{
if (0.90 <= pos[1])
if (0.90 < pos[1] + eps_)
{
return true;
else return false;
}
else
{
return false;
}
};
Scalar fineK_;
......@@ -296,6 +287,8 @@ private:
MaterialLawParams coarseMaterialParams_;
Scalar lambdaSolid_;
static constexpr Scalar eps_ = 1.5e-7;
};
}
......
Supports Markdown
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