Skip to content
GitLab
Projects
Groups
Snippets
Help
Loading...
Help
What's new
7
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in / Register
Toggle navigation
Open sidebar
dorie
dorie
Commits
19288ec4
Commit
19288ec4
authored
Oct 04, 2018
by
Santiago Ospina
Browse files
Options
Browse Files
Download
Email Patches
Plain Diff
clean up transport implementation
Signed-off-by:
Santiago Ospina
<
santiago.ospina@iup.uni-heidelberg.de
>
parent
e4145e9e
Changes
9
Pipelines
1
Hide whitespace changes
Inline
Side-by-side
Showing
9 changed files
with
274 additions
and
269 deletions
+274
-269
doc/default_files/parameters-transport.xml
doc/default_files/parameters-transport.xml
+2
-2
dune/dorie/interface/transport_simulation.cc
dune/dorie/interface/transport_simulation.cc
+19
-19
dune/dorie/interface/transport_simulation.hh
dune/dorie/interface/transport_simulation.hh
+21
-21
dune/dorie/interface/util.hh
dune/dorie/interface/util.hh
+4
-12
dune/dorie/solver/transport_boundary.hh
dune/dorie/solver/transport_boundary.hh
+86
-86
dune/dorie/solver/transport_initial.hh
dune/dorie/solver/transport_initial.hh
+1
-1
dune/dorie/solver/transport_operator_FV.hh
dune/dorie/solver/transport_operator_FV.hh
+129
-116
dune/dorie/test/test-simulation-dummy-transport.cc
dune/dorie/test/test-simulation-dummy-transport.cc
+10
-10
dune/dorie/test/test-simulation-dummy-transport.mini.in
dune/dorie/test/test-simulation-dummy-transport.mini.in
+2
-2
No files found.
doc/default_files/parameters-transport.xml
View file @
19288ec4
...
...
@@ -174,8 +174,8 @@ adding an empty line, make text **bold** or ``monospaced``.
</category>
<category
name=
"transport.parameters"
>
<parameter
name=
"molecularDifu
s
sion"
>
<definition>
Molecular difu
s
sion of the fluid phase
</definition>
<parameter
name=
"molecularDif
f
usion"
>
<definition>
Molecular dif
f
usion of the fluid phase
</definition>
<values>
float
>
0
</values>
<comment>
2E-9
</comment>
</parameter>
...
...
dune/dorie/interface/transport_simulation.cc
View file @
19288ec4
...
...
@@ -79,16 +79,16 @@ TransportSimulation<Traits>::TransportSimulation(
template
<
typename
Traits
>
void
TransportSimulation
<
Traits
>::
operator_setup
()
{
if
(
!
_
gf
W
ater
F
lux
)
if
(
!
gf
_w
ater
_f
lux
)
DUNE_THROW
(
Dune
::
Exception
,
" "
);
if
(
!
_
gf
S
aturation
)
if
(
!
gf
_s
aturation
)
DUNE_THROW
(
Dune
::
Exception
,
" "
);
// --- Local Operators ---
double
D_
eff
=
inifile
.
get
<
double
>
(
"parameters.molecularDifu
s
sion"
);
slop
=
std
::
make_unique
<
SLOP
>
(
*
sboundary
,
_
gf
W
ater
F
lux
,
_
gf
S
aturation
,
D_
eff
/*inifile,*param,gv,*fboundary,*fsource*/
);
double
diff_co
eff
=
inifile
.
get
<
double
>
(
"parameters.molecularDif
f
usion"
);
slop
=
std
::
make_unique
<
SLOP
>
(
*
sboundary
,
gf
_w
ater
_f
lux
,
gf
_s
aturation
,
diff_co
eff
/*inifile,*param,gv,*fboundary,*fsource*/
);
tlop
=
std
::
make_unique
<
TLOP
>
(
_
gf
S
aturation
/*inifile,*param,gv*/
);
tlop
=
std
::
make_unique
<
TLOP
>
(
gf
_s
aturation
/*inifile,*param,gv*/
);
Dune
::
PDELab
::
constraints
(
*
this
->
gfs
,
*
this
->
cc
,
false
);
// --- Grid Operators ---
...
...
@@ -145,14 +145,14 @@ void TransportSimulation<Traits>::step ()
this
->
suggest_timestep
(
controller
->
getDT
());
// Check whether time intervals of grid functions are valid
if
(
Dune
::
FloatCmp
::
gt
(
_
water
F
lux
I
nterval
.
begin
()
,
t
+
_
dt
))
if
(
Dune
::
FloatCmp
::
gt
(
water
_f
lux
_i
nterval
.
begin
,
t
+
dt
))
DUNE_THROW
(
Dune
::
Exception
,
" "
);
if
(
Dune
::
FloatCmp
::
gt
(
_
saturation
I
nterval
.
begin
()
,
t
+
_
dt
))
if
(
Dune
::
FloatCmp
::
gt
(
saturation
_i
nterval
.
begin
,
t
+
dt
))
DUNE_THROW
(
Dune
::
Exception
,
" "
);
if
(
Dune
::
FloatCmp
::
lt
(
_
water
F
lux
I
nterval
.
end
()
,
t
+
_
dt
))
DUNE_THROW
(
Dune
::
Exception
,
_
water
F
lux
I
nterval
.
end
()
<<
" < "
<<
t
+
_
dt
);
if
(
Dune
::
FloatCmp
::
lt
(
_
saturation
I
nterval
.
end
()
,
t
+
_
dt
))
if
(
Dune
::
FloatCmp
::
lt
(
water
_f
lux
_i
nterval
.
end
,
t
+
dt
))
DUNE_THROW
(
Dune
::
Exception
,
water
_f
lux
_i
nterval
.
end
<<
" < "
<<
t
+
dt
);
if
(
Dune
::
FloatCmp
::
lt
(
saturation
_i
nterval
.
end
,
t
+
dt
))
DUNE_THROW
(
Dune
::
Exception
,
" "
);
const
bool
solver_warnings
=
verbose
>
0
&&
helper
.
rank
()
==
0
?
...
...
@@ -164,13 +164,13 @@ void TransportSimulation<Traits>::step ()
if
(
not
solver_warnings
)
dwarn
.
push
(
false
);
osm
->
apply
(
t
,
_
dt
,
*
u
,
*
unext
);
osm
->
apply
(
t
,
dt
,
*
u
,
*
unext
);
if
(
not
solver_warnings
)
dwarn
.
pop
();
u
=
unext
;
if
(
not
ts_param
->
implicit
())
this
->
suggest_timestep
(
slop
->
suggestTimestep
(
_
dt
));
this
->
suggest_timestep
(
slop
->
suggestTimestep
(
dt
));
}
catch
(
Dune
::
ISTLError
&
e
){
exception
=
true
;
...
...
@@ -192,9 +192,9 @@ void TransportSimulation<Traits>::step ()
template
<
typename
Traits
>
void
TransportSimulation
<
Traits
>::
update_adapters
()
const
{
C
w
=
std
::
make_shared
<
GFSolute
>
(
gfs
,
u
);
assert
(
_
gf
S
aturation
);
C
total
=
std
::
make_shared
<
GFTotalSolute
>
(
*
C
w
,
*
_
gf
S
aturation
);
c_
w
=
std
::
make_shared
<
GFSolute
>
(
gfs
,
u
);
assert
(
gf
_s
aturation
);
c_
total
=
std
::
make_shared
<
GFTotalSolute
>
(
*
c_
w
,
*
gf
_s
aturation
);
}
template
<
typename
Traits
>
...
...
@@ -206,11 +206,11 @@ void TransportSimulation<Traits>::write_data () const
update_adapters
();
if
(
inifile
.
get
<
bool
>
(
"output.vertexData"
))
{
vtkwriter
->
addVertexData
(
C
w
,
"solute"
);
vtkwriter
->
addVertexData
(
C
total
,
"solute (total)"
);
vtkwriter
->
addVertexData
(
c_
w
,
"solute"
);
vtkwriter
->
addVertexData
(
c_
total
,
"solute (total)"
);
}
else
{
vtkwriter
->
addCellData
(
C
w
,
"solute"
);
vtkwriter
->
addCellData
(
C
total
,
"solute (total)"
);
vtkwriter
->
addCellData
(
c_
w
,
"solute"
);
vtkwriter
->
addCellData
(
c_
total
,
"solute (total)"
);
}
try
{
...
...
dune/dorie/interface/transport_simulation.hh
View file @
19288ec4
...
...
@@ -224,22 +224,22 @@ protected:
std
::
shared_ptr
<
U
>
u
;
mutable
std
::
shared_ptr
<
GFSolute
>
C
w
;
mutable
std
::
shared_ptr
<
GFTotalSolute
>
C
total
;
std
::
shared_ptr
<
GFWaterFlux
>
_
gf
W
ater
F
lux
;
std
::
shared_ptr
<
GFSaturation
>
_
gf
S
aturation
;
mutable
std
::
shared_ptr
<
GFSolute
>
c_
w
;
mutable
std
::
shared_ptr
<
GFTotalSolute
>
c_
total
;
std
::
shared_ptr
<
GFWaterFlux
>
gf
_w
ater
_f
lux
;
std
::
shared_ptr
<
GFSaturation
>
gf
_s
aturation
;
TimeInterval
_
water
F
lux
I
nterval
=
TimeInterval
(
std
::
numeric_limits
<
TimeField
>::
min
(),
std
::
numeric_limits
<
TimeField
>::
min
()
)
;
TimeInterval
_
saturation
I
nterval
=
TimeInterval
(
std
::
numeric_limits
<
TimeField
>::
min
(),
std
::
numeric_limits
<
TimeField
>::
min
()
)
;
TimeInterval
water
_f
lux
_i
nterval
=
TimeInterval
{
std
::
numeric_limits
<
TimeField
>::
min
(),
std
::
numeric_limits
<
TimeField
>::
min
()
}
;
TimeInterval
saturation
_i
nterval
=
TimeInterval
{
std
::
numeric_limits
<
TimeField
>::
min
(),
std
::
numeric_limits
<
TimeField
>::
min
()
}
;
std
::
unique_ptr
<
Writer
>
vtkwriter
;
const
int
verbose
;
const
Dune
::
VTK
::
OutputType
output_type
;
TimeField
_
dt
;
TimeField
dt
;
bool
operator_setup_flag
=
true
;
public:
...
...
@@ -257,23 +257,23 @@ public:
);
// TODO: Define data exchange with richards
void
set_waterflux
(
TimeInterval
time_interval
,
std
::
shared_ptr
<
GFWaterFlux
>
gf
W
ater
F
lux
)
void
set_waterflux
(
TimeInterval
_
time_interval
,
std
::
shared_ptr
<
GFWaterFlux
>
_
gf
_w
ater
_f
lux
)
{
if
(
gf
W
ater
F
lux
)
if
(
_
gf
_w
ater
_f
lux
)
{
_
water
F
lux
I
nterval
=
time_interval
;
_
gf
W
ater
F
lux
=
gf
W
ater
F
lux
;
water
_f
lux
_i
nterval
=
_
time_interval
;
gf
_w
ater
_f
lux
=
_
gf
_w
ater
_f
lux
;
}
else
DUNE_THROW
(
Dune
::
Exception
,
"Pointer to gfWaterFlux is invalid!"
);
}
// TODO: Define data exchange with richards
void
set_saturation
(
TimeInterval
time_interval
,
std
::
shared_ptr
<
GFSaturation
>
gf
S
aturation
)
void
set_saturation
(
TimeInterval
_
time_interval
,
std
::
shared_ptr
<
GFSaturation
>
_
gf
_s
aturation
)
{
if
(
gf
S
aturation
)
if
(
_
gf
_s
aturation
)
{
_
saturation
I
nterval
=
time_interval
;
_
gf
S
aturation
=
gf
S
aturation
;
saturation
_i
nterval
=
_
time_interval
;
gf
_s
aturation
=
_
gf
_s
aturation
;
}
else
DUNE_THROW
(
Dune
::
Exception
,
"Pointer to gfWaterFlux is invalid!"
);
}
...
...
@@ -295,17 +295,17 @@ public:
return
controller
->
getTime
();
}
void
suggest_timestep
(
TimeField
dt
)
override
void
suggest_timestep
(
TimeField
_
dt
)
override
{
_
dt
=
std
::
min
(
dt
,
controller
->
getDT
());
dt
=
std
::
min
(
_
dt
,
controller
->
getDT
());
}
// TODO: Define data exchange with richards
std
::
shared_ptr
<
GFSolute
>
get_solute
()
{
if
(
!
C
w
)
if
(
!
c_
w
)
DUNE_THROW
(
Dune
::
Exception
,
"Solute in invalid state!"
);
return
C
w
;
return
c_
w
;
}
protected:
...
...
dune/dorie/interface/util.hh
View file @
19288ec4
...
...
@@ -194,18 +194,10 @@ public:
}
};
template
<
class
TF
>
class
TimeInterval
{
public:
using
TimeField
=
TF
;
TimeInterval
(
TimeField
begin
,
TimeField
end
)
:
_begin
(
begin
),
_end
(
end
)
{}
const
TimeField
&
begin
()
const
{
return
_begin
;}
TimeField
&
begin
()
{
return
_begin
;}
const
TimeField
&
end
()
const
{
return
_end
;}
TimeField
&
end
()
{
return
_end
;}
private:
TimeField
_begin
,
_end
;
template
<
typename
TF
>
struct
TimeInterval
{
TF
begin
;
TF
end
;
};
...
...
dune/dorie/solver/transport_boundary.hh
View file @
19288ec4
...
...
@@ -9,106 +9,106 @@
#include "util_boundary_condition.hh"
namespace
Dune
{
namespace
Dorie
{
namespace
Dorie
{
/// Boundary type and condition value queries
/** This class containts functions that return the type of boundary conditions
* and the values of the boundary condition parameters. Both types of queries
* can be time dependent.
*/
template
<
typename
Traits
>
class
SoluteBoundary
{
private:
/// Boundary type and condition value queries
/** This class containts functions that return the type of boundary conditions
* and the values of the boundary condition parameters. Both types of queries
* can be time dependent.
*/
template
<
typename
Traits
>
class
SoluteBoundary
{
private:
typedef
typename
Traits
::
RangeField
RF
;
typedef
typename
Traits
::
Domain
Domain
;
typedef
typename
Traits
::
IntersectionDomain
ID
;
typedef
typename
Traits
::
Intersection
Intersection
;
typedef
typename
Traits
::
RangeField
RF
;
typedef
typename
Traits
::
Domain
Domain
;
typedef
typename
Traits
::
IntersectionDomain
ID
;
typedef
typename
Traits
::
Intersection
Intersection
;
enum
{
dim
=
Traits
::
dim
};
enum
{
dim
=
Traits
::
dim
};
//! Object for handling the BC data file and function queries
std
::
unique_ptr
<
BCReadoutInterface
<
Traits
>>
bcDataHandler
;
//! Object for handling the BC data file and function queries
std
::
unique_ptr
<
BCReadoutInterface
<
Traits
>>
bcDataHandler
;
public:
public:
/// Read BC file type and create data handling object
/** \see Dune::Dorie::BCReadoutInterface
* \throw Dune::IOError BC File Type not supported
*/
SoluteBoundary
(
const
Dune
::
ParameterTree
&
config
)
{
std
::
string
bcFileType
=
config
.
get
<
std
::
string
>
(
"boundary.fileType"
);
if
(
bcFileType
==
"rectangularGrid"
)
{
bcDataHandler
=
std
::
make_unique
<
RectangularGrid
<
Traits
>>
(
config
);
}
else
DUNE_THROW
(
IOError
,
"unknown bcFileType specified!"
);
/// Read BC file type and create data handling object
/** \see Dune::Dorie::BCReadoutInterface
* \throw Dune::IOError BC File Type not supported
*/
SoluteBoundary
(
const
Dune
::
ParameterTree
&
config
)
{
std
::
string
bcFileType
=
config
.
get
<
std
::
string
>
(
"boundary.fileType"
);
if
(
bcFileType
==
"rectangularGrid"
)
{
bcDataHandler
=
std
::
make_unique
<
RectangularGrid
<
Traits
>>
(
config
);
}
else
DUNE_THROW
(
IOError
,
"unknown bcFileType specified!"
);
}
/// Return next time at which any boundary condition changes. Return -1.0 by default
/** \param time Current time
*/
RF
getNextTimeStamp
(
const
RF
&
time
)
const
{
if
(
bcDataHandler
)
return
bcDataHandler
->
getNextTimeStamp
(
time
);
return
-
1.0
;
}
/// Return next time at which any boundary condition changes. Return -1.0 by default
/** \param time Current time
*/
RF
getNextTimeStamp
(
const
RF
&
time
)
const
{
if
(
bcDataHandler
)
return
bcDataHandler
->
getNextTimeStamp
(
time
);
return
-
1.0
;
}
/// Return Boundary Condition Type at certain position and time
/** \param is Intersection enitity
* \param x Position in local entity coordinates
* \param time Time value
* \param xGlobal Global coordinate
* \return Boundary condition type enumerator
*/
BoundaryCondition
::
Type
bc
(
const
Intersection
&
is
,
const
ID
&
x
,
const
RF
&
time
)
const
{
const
Domain
xGlobal
=
is
.
geometry
().
global
(
x
);
if
(
is
.
boundary
())
{
return
bcDataHandler
->
getBCtype
(
xGlobal
,
time
);
}
return
BoundaryCondition
::
Other
;
// unknown
/// Return Boundary Condition Type at certain position and time
/** \param is Intersection enitity
* \param x Position in local entity coordinates
* \param time Time value
* \param xGlobal Global coordinate
* \return Boundary condition type enumerator
*/
BoundaryCondition
::
Type
bc
(
const
Intersection
&
is
,
const
ID
&
x
,
const
RF
&
time
)
const
{
const
Domain
xGlobal
=
is
.
geometry
().
global
(
x
);
if
(
is
.
boundary
())
{
return
bcDataHandler
->
getBCtype
(
xGlobal
,
time
);
}
return
BoundaryCondition
::
Other
;
// unknown
}
/**
* @brief Dirichlet boundary condition at certain position and time
* \param is Intersection enitity
* \param x Position in local entity coordinates
* \param time Time value
* \param xglobal Global coordinate
* \return Value of matric head
*/
RF
g
(
const
Intersection
&
is
,
const
ID
&
x
,
const
RF
&
time
)
const
{
const
Domain
xGlobal
=
is
.
geometry
().
global
(
x
);
if
(
is
.
boundary
())
{
if
(
bcDataHandler
)
return
bcDataHandler
->
getDirichletValue
(
xGlobal
,
time
);
}
return
0.0
;
/**
* @brief Dirichlet boundary condition at certain position and time
* \param is Intersection enitity
* \param x Position in local entity coordinates
* \param time Time value
* \param xglobal Global coordinate
* \return Value of matric head
*/
RF
g
(
const
Intersection
&
is
,
const
ID
&
x
,
const
RF
&
time
)
const
{
const
Domain
xGlobal
=
is
.
geometry
().
global
(
x
);
if
(
is
.
boundary
())
{
if
(
bcDataHandler
)
return
bcDataHandler
->
getDirichletValue
(
xGlobal
,
time
);
}
return
0.0
;
}
/**
* @brief Neumann boundary condition at certain position and time
* \param is Intersection enitity
* \param x Position in local entity coordinates
* \param time Time value
* \param xglobal Global coordinate
* \return Value of flux
*/
RF
j
(
const
Intersection
&
is
,
const
ID
&
x
,
const
RF
&
time
)
const
{
const
Domain
xGlobal
=
is
.
geometry
().
global
(
x
);
if
(
is
.
boundary
())
{
if
(
bcDataHandler
)
return
bcDataHandler
->
getNeumannValue
(
xGlobal
,
time
);
}
return
0.0
;
/**
* @brief Neumann boundary condition at certain position and time
* \param is Intersection enitity
* \param x Position in local entity coordinates
* \param time Time value
* \param xglobal Global coordinate
* \return Value of flux
*/
RF
j
(
const
Intersection
&
is
,
const
ID
&
x
,
const
RF
&
time
)
const
{
const
Domain
xGlobal
=
is
.
geometry
().
global
(
x
);
if
(
is
.
boundary
())
{
if
(
bcDataHandler
)
return
bcDataHandler
->
getNeumannValue
(
xGlobal
,
time
);
}
};
}
return
0.0
;
}
};
}
}
#endif
dune/dorie/solver/transport_initial.hh
View file @
19288ec4
...
...
@@ -4,7 +4,7 @@
#include <dune/pdelab/function/const.hh>
namespace
Dune
{
namespace
Dorie
{
namespace
Dorie
{
template
<
typename
T
>
using
SoluteInitial
=
Dune
::
PDELab
::
ConstGridFunction
<
typename
T
::
GridView
,
typename
T
::
RangeField
>
;
...
...
dune/dorie/solver/transport_operator_FV.hh
View file @
19288ec4
...
...
@@ -80,20 +80,22 @@ public:
enum
{
doPatternSkeleton
=
true
};
// Residual assembly flags
enum
{
doAlphaVolume
=
true
};
//
//
enum { doAlphaVolume = true };
enum
{
doAlphaBoundary
=
true
};
enum
{
doAlphaSkeleton
=
true
};
enum
{
doLambdaVolume
=
true
};
//
// enum { doLambdaVolume = true };
// enum { doLambdaBoundary = true };
private:
static_assert
(
std
::
is_same
<
typename
GridFunctionWaterFlux
::
Traits
::
GridViewType
,
typename
GridFunctionSaturation
::
Traits
::
GridViewType
>::
value
,
"TransportFVSpatialOperator: GridFunctionWaterFlux and GridFunctionSaturation has to use the same grid view."
);
static_assert
(
std
::
is_same
<
typename
GridFunctionWaterFlux
::
Traits
::
GridViewType
,
typename
GridFunctionSaturation
::
Traits
::
GridViewType
>::
value
,
"TransportFVSpatialOperator: GridFunctionWaterFlux and"
"GridFunctionSaturation has to use the same grid view."
);
// Define domain field
using
DF
=
typename
GridFunctionSaturation
::
Traits
::
DomainFieldType
;
using
Flux
=
typename
GridFunctionWaterFlux
::
Traits
::
RangeType
;
using
Water
Flux
=
typename
GridFunctionWaterFlux
::
Traits
::
RangeType
;
using
Saturation
=
typename
GridFunctionSaturation
::
Traits
::
RangeType
;
// Extract world dimension
...
...
@@ -101,15 +103,23 @@ private:
public:
/**
* @brief Constructor of TransportFVSpatialOperator
*
* @param boundary The boundary terms
* @param[in] gf_flux The grid function of water flux
* @param[in] gf_sat The grid function of water saturation
* @param[in] diff_coeff The diffusion coefficient
*/
TransportFVSpatialOperator
(
Boundary
&
boundary
,
std
::
shared_ptr
<
GridFunctionWaterFlux
>
gf
F
lux
,
std
::
shared_ptr
<
GridFunctionSaturation
>
gf
Saturation
,
double
D_
eff
std
::
shared_ptr
<
GridFunctionWaterFlux
>
gf
_f
lux
,
std
::
shared_ptr
<
GridFunctionSaturation
>
gf
_sat
,
double
diff_co
eff
)
:
_boundary
(
boundary
)
,
_gf
F
lux
(
gf
F
lux
)
,
_gf
Saturation
(
gfSaturation
)
,
_
D_eff
(
D_
eff
)
,
_gf
_f
lux
(
gf
_f
lux
)
,
_gf
_sat
(
gf_sat
)
,
_
diff_coeff
(
diff_co
eff
)
{}
/*---------------------------------------------------------------------*//**
...
...
@@ -164,73 +174,73 @@ public:
Traits
::
LocalBasisType
::
Traits
::
RangeFieldType
;
// Get entity entities from both sides of the intersection
auto
entity_i
nside
=
ig
.
inside
();
auto
entity_o
utside
=
ig
.
outside
();
const
auto
&
entity_i
=
ig
.
inside
();
const
auto
&
entity_o
=
ig
.
outside
();
// Get geometries
auto
geo
=
ig
.
geometry
();
auto
geo_i
nside
=
entity_i
nside
.
geometry
();
auto
geo_o
utside
=
entity_o
utside
.
geometry
();
auto
geo_i
=
entity_i
.
geometry
();
auto
geo_o
=
entity_o
.
geometry
();
// Get geometry of intersection in local coordinates of neighbor entities
auto
geo_in_inside
=
ig
.
geometryInInside
();
// Center in face's reference element
auto
ref_el
=
referenceElement
(
geo
);
auto
face_local
=
ref_el
.
position
(
0
,
0
);
auto
geo_in_i
=
ig
.
geometryInInside
();
// Face volume for integration
auto
face_volume
=
geo
.
volume
();
// Entity centers in references elements
auto
ref_el_inside
=
referenceElement
(
geo_inside
);
auto
ref_el_outside
=
referenceElement
(
geo_outside
);
auto
local_inside
=
ref_el_inside
.
position
(
0
,
0
);
auto
local_outside
=
ref_el_outside
.
position
(
0
,
0
);
// Evaluate diffusion coeff. from either side and take harmonic average
auto
D_eff_inside
=
_D_eff
;
auto
D_eff_outside
=
_D_eff
;
// Evaluate diffusion coeff. from both sides and take harmonic average
auto
diff_coeff_i
=
_diff_coeff
;
// FIXME
auto
diff_coeff_o
=
_diff_coeff
;
// FIXME
auto
normal_face
=
ig
.
centerUnitOuterNormal
();
auto
D_eff_avg
=
2.0
/
(
1.0
/
(
D_eff_inside
+
1E-30
)
+
1.0
/
(
D_eff_outside
+
1E-30
));
auto
diff_coeff_avg
=
2.0
/
(
1.0
/
(
diff_coeff_i
+
1E-30
)
+
1.0
/
(
diff_coeff_o
+
1E-30
));
auto
entity_inside_local
=
geo_in_inside
.
global
(
face_local
);
// Get center position of the face w.r.t inside entity
auto
p_center_face_i
=
geo_in_i
.
center
(
);
// Evaluate convective term (i.e water flux)
Flux
flux_i
nside
;
_gf
F
lux
->
evaluate
(
entity_i
nside
,
entity_inside_local
,
flux_i
nside
);
auto
vn
=
flux_i
nside
*
normal_face
;
WaterFlux
water_
flux_i
;
_gf
_f
lux
->
evaluate
(
entity_i
,
p_center_face_i
,
water_
flux_i
);
auto
water_flux_n
=
water_
flux_i
*
normal_face
;
// Inside
unknown
value
// Inside
solute
value
RF
u
=
x_i
(
lfsu_i
,
0
);
#if 1
// Upwinding
u
=
(
v
n
>=
0
)
?
x_i
(
lfsu_i
,
0
)
:
x_o
(
lfsu_o
,
0
);
u
=
(
water_flux_
n
>=
0
)
?
x_i
(
lfsu_i
,
0
)
:
x_o
(
lfsu_o
,
0
);
#endif
// Entity centers in references elements
auto
ref_el_i
=
referenceElement
(
geo_i
);
auto
ref_el_o
=
referenceElement
(
geo_o
);
auto
p_center_i
=
ref_el_i
.
position
(
0
,
0
);
auto
p_center_o
=
ref_el_o
.
position
(
0
,
0
);
// Entity centers in global coordinates
auto
global_inside
=
geo_i
nside
.
global
(
local_inside
);
auto
global_outside
=
geo_o
utside
.
global
(
local_outside
);
auto
p_center_i_global
=
geo_i
.
global
(
p_center_i
);
auto
p_center_o_global
=
geo_o
.
global
(
p_center_o
);
// Distance between the two entity centers
global_inside
-=
global_outside
;
auto
distance
=
global_inside
.
two_norm
();
p_center_i_global
-=
p_center_o_global
;
auto
distance
=
p_center_i_global
.
two_norm
();
// Finite difference of u between the two entities
RF
dudn
=
(
x_o
(
lfsu_o
,
0
)
-
x_i
(
lfsu_i
,
0
))
/
distance
;
// Evaluate saturation
typename
GridFunctionSaturation
::
Traits
::
RangeType
sat_inside
;
_gfSaturation
->
evaluate
(
ig
.
inside
(),
entity_inside_local
,
sat_inside
);
typename
GridFunctionSaturation
::
Traits
::
RangeType
sat_i
;
_gf_sat
->
evaluate
(
entity_i
,
p_center_i
,
sat_i
);
// Solute flux in normal direction w.r.t the intersection
auto
normal_flux
=
u
*
water_flux_n
-
diff_coeff_avg
*
dudn
*
sat_i
;
// Symmetric contribution to residual on inside element
r_i
.
accumulate
(
lfs
u
_i
,
0
,
(
u
*
vn
-
D_eff_avg
*
dudn
*
sat_inside
)
*
face_volume
);
r_o
.
accumulate
(
lfs
u
_o
,
0
,
-
(
u
*
vn
-
D_eff_avg
*
dudn
*
sat_inside
)
*
face_volume
);
r_i
.
accumulate
(
lfs
v
_i
,
0
,
normal_flux
*
face_volume
);
r_o
.
accumulate
(
lfs
v
_o
,
0
,
-
normal_flux
*
face_volume
);
if
(
_first_stage
)
_dtmin
=
std
::
min
(
_dtmin
,
suggest
T
imestep
I
ntersection
(
ig
));
_dt
_
min
=
std
::
min
(
_dt
_
min
,
suggest
_t
imestep
_i
ntersection
(
ig
));
}
...
...
@@ -254,12 +264,16 @@ public:
void
alpha_boundary
(
const
IG
&
ig
,
const
LFSU
&
lfsu_i
,
const
X
&
x_i
,
const
LFSV
&
lfsv_i
,
R
&
r_i
)
const
{
// Get entities
const
auto
&
entitiy_face
=
ig
.
intersection
();
const
auto
&
entity_i
=
ig
.
inside
();
// Get geometries
auto
geo
=
ig
.
geometry
();
auto
geo
=
entitiy_face
.
geometry
();
auto
ref_el
=
referenceElement
(
geo
);
auto
face_local
=
ref_el
.
position
(
0
,
0
);
auto
p_center_face
=
ref_el
.
position
(
0
,
0
);