Not safe for unit system compatibility.
I tried to use PDELab with data types (DF
& RF
) that only allows operations that preserve the units (+
,-
) or that transform them appropriately (*
,/
). In this case, I used boost units.
For such a thing I created a new grid that "cast" any grid type to hold its units in a certain quantity, this because some grid types don't allow to modify ctype
directly. However, in order to illustrate the problem I will show it with YaspGrid
, which already allows changing its ctype
:
const int dim = 2;
using ctype = boost::units::quantity<boost::units::si::length>;
using Coordinates = Dune::EquidistantCoordinates<ctype,dim>;
using Grid = Dune::YaspGrid<dim,Coordinates>;
... // define L and N
std::shared_ptr<Grid> gridp = std::shared_ptr<Grid>(new Grid(L,N));
typedef Grid::LeafGridView GV;
GV gv = gridp->leafGridView();
So far, the domain is already setup and works perfectly. Now, the range part:
// choose some FEM with the desired domain and range
using DF = typename Grid::ctype;
using RF = boost::units::quantity<boost::units::si::pressure>;
typedef Dune::PDELab::P0LocalFiniteElementMap<DF,RF,dim> FEM;
FEM fem(Dune::GeometryTypes::cube(dim));
Now, everything as usual:
// Make grid function space used per component
using CON = Dune::PDELab::ConformingDirichletConstraints;
using VBE = Dune::PDELab::ISTL::VectorBackend<>;
using GFS = Dune::PDELab::GridFunctionSpace<GV,FEM,CON,VBE>;
GFS gfs(gv,fem);
gfs.name("p_w");
Here I notice that even GridFunction
transformations are working with the unit system:
// define the initial condition
auto ulambda = [](const auto& x){
Dune::FieldVector<RF,1> rv(0.0 * boost::units::si::pascal);
return rv;
};
auto u = Dune::PDELab::makeGridFunctionFromCallable(gv,ulambda);
// set up coefficient vector filled with initial condition
using Z = Dune::PDELab::Backend::Vector<GFS,RF>;
Z z(gfs);
Dune::PDELab::interpolate(u,gfs,z);
Setup of constraints also works.
auto blambda = [](const auto& x){return true;};
auto b = Dune::PDELab::
makeBoundaryConditionFromCallable(gv,blambda);
using CC = typename GFS::template ConstraintsContainer<RF>::Type;
CC cc;
Dune::PDELab::constraints(b,gfs,cc,true);
set_constrained_dofs(cc,0.0 * boost::units::si::pascal,z);
The VTKwriter
s of course don't work because they try to transform RF
to double
(see vtkexport.hh:78
), something that is usually not allowed by unit systems, however, shouldn't be hard to fix with some cast policy or an extra user GridFunction
that cast every unit quantity to double
.
Now, all of this is to be sure that all operations in local operators are at least correct in units. However, to reproduce the same problem I want to point out, any void local operator work.
using LOP = OnePhaseFlowFV; // in my test this operator is actually void.
LOP lop;
Finally, where the problem comes:
DF domain_value;
RF range_value;
auto jacobian_value = range_value/domain_value;
using JF = decltype(jacobian_value);
using MBE = Dune::PDELab::ISTL::BCRSMatrixBackend<>;
MBE mbe((int)pow(1+2*1,dim));
using GO0 = Dune::PDELab::GridOperator<GFS,GFS,LOP,MBE,DF,RF,JF,CC,CC>;
GO0 go0(gfs,cc,gfs,cc,lop,mbe);
This last line leads to compilation errors because of not safe unit operations. As far as I have seen it is coming from these specific points:
DefaultLocalAssambler
uses a weight
that is of the same type of the RangeField
, if I understand it well, weights must be dimensionless (e.g. dobule
). In particular, it fails because a unit value cannot be initialized without units (see defaultassambler.hh:96
and residualengine.hh:66
). Another further point of failure would be the definition of the type of time (see localassambler.hh:44
), which is generally different to the RangeField
.
Now, if we think about systems of PDE, DF
must be the same for every equation, but RF
is in general not. It seems to me that blocked vectors/matrices of the ISTL are able to manage that, however, I'm not really sure if it works now. Does it? (I'm talking specifically about the template arguments in GridOperator
which seems to only accept one type).
Do you see possible to fix it? or is there any deep PDELab concept that does not allow that?