Skip to content
Snippets Groups Projects
Commit 15aac202 authored by Sander, Oliver's avatar Sander, Oliver
Browse files

Refactor code to make it more flexible

parent 886a9260
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/compute-disc-error.cc
+ 144
− 87
View file @ 15aac202
......@@ -39,9 +39,9 @@ const int blocksize = TargetSpace::TangentVector::dimension;
using namespace Dune;
template <class GridView, int order>
void measureEOC(const GridView gridView,
const GridView referenceGridView,
const ParameterTree& parameterSet)
void measureDiscreteEOC(const GridView gridView,
const GridView referenceGridView,
const ParameterTree& parameterSet)
{
typedef std::vector<TargetSpace> SolutionType;
......@@ -51,6 +51,11 @@ void measureEOC(const GridView gridView,
typedef Dune::Functions::PQkNodalBasis<GridView, order> FEBasis;
FEBasis feBasis(gridView);
FEBasis referenceFEBasis(referenceGridView);
using FufemFEBasis = DuneFunctionsBasis<FEBasis>;
FufemFEBasis fufemReferenceFEBasis(referenceFEBasis);
FufemFEBasis fufemFEBasis(feBasis);
//////////////////////////////////////////////////////////////////////////////////
// Read the data whose error is to be measured
......@@ -70,108 +75,132 @@ void measureEOC(const GridView gridView,
for (size_t i=0; i<x.size(); i++)
x[i] = TargetSpace(embeddedX[i]);
// The numerical solution, as a grid function
GFE::EmbeddedGlobalGFEFunction<FufemFEBasis, TargetSpace> numericalSolution(fufemFEBasis, x);
///////////////////////////////////////////////////////////////////////////
// Read the reference configuration
///////////////////////////////////////////////////////////////////////////
EmbeddedVectorType embeddedReferenceX(referenceFEBasis.size());
inFile.open(parameterSet.get<std::string>("referenceData"), std::ios_base::binary);
if (not inFile)
DUNE_THROW(IOError, "File " << parameterSet.get<std::string>("referenceData") << " could not be opened.");
GenericVector::readBinary(inFile, embeddedReferenceX);
SolutionType referenceX(embeddedReferenceX.size());
for (size_t i=0; i<referenceX.size(); i++)
referenceX[i] = TargetSpace(embeddedReferenceX[i]);
// The reference solution, as a grid function
GFE::EmbeddedGlobalGFEFunction<FufemFEBasis, TargetSpace> referenceSolution(fufemReferenceFEBasis, referenceX);
/////////////////////////////////////////////////////////////////
// Measure the discretization error
/////////////////////////////////////////////////////////////////
if (parameterSet.get<std::string>("discretizationErrorMode")=="analytical")
double l2ErrorSquared = 0;
double h1ErrorSquared = 0;
HierarchicSearch<typename GridView::Grid,typename GridView::IndexSet> hierarchicSearch(gridView.grid(), gridView.indexSet());
for (const auto& rElement : elements(referenceGridView))
{
using FufemFEBasis = DuneFunctionsBasis<FEBasis>;
FufemFEBasis fufemFEBasis(feBasis);
const auto& quadRule = QuadratureRules<double, dim>::rule(rElement.type(), 6);
for (const auto& qp : quadRule)
{
auto integrationElement = rElement.geometry().integrationElement(qp.position());
// Read reference solution and its derivative into a PythonFunction
typedef VirtualDifferentiableFunction<FieldVector<double, dim>, TargetSpace::CoordinateType> FBase;
auto globalPos = rElement.geometry().global(qp.position());
Python::Module module = Python::import(parameterSet.get<std::string>("referenceSolution"));
auto referenceSolution = module.get("fdf").toC<std::shared_ptr<FBase>>();
auto element = hierarchicSearch.findEntity(globalPos);
auto localPos = element.geometry().local(globalPos);
// The numerical solution, as a grid function
GFE::EmbeddedGlobalGFEFunction<FufemFEBasis, TargetSpace> numericalSolution(fufemFEBasis, x);
auto diff = referenceSolution(rElement, qp.position()) - numericalSolution(element, localPos);
// QuadratureRule for the integral of the L^2 error
QuadratureRuleKey quadKey(dim,6);
l2ErrorSquared += integrationElement * qp.weight() * diff.two_norm2();
// Compute the embedded L^2 error
double l2Error = DiscretizationError<GridView>::computeL2Error(&numericalSolution,
referenceSolution.get(),
quadKey);
auto derDiff = referenceSolution.derivative(rElement, qp.position()) - numericalSolution.derivative(element, localPos);
// Compute the embedded H^1 error
double h1Error = DiscretizationError<GridView>::computeH1HalfNormDifferenceSquared(gridView,
&numericalSolution,
referenceSolution.get(),
quadKey);
h1ErrorSquared += integrationElement * qp.weight() * derDiff.frobenius_norm2();
std::cout << "elements: " << gridView.size(0)
<< " "
<< "L^2 error: " << l2Error
<< " ";
std::cout << "H^1 error: " << std::sqrt(l2Error*l2Error + h1Error) << std::endl;
}
}
if (parameterSet.get<std::string>("discretizationErrorMode")=="discrete")
{
FEBasis referenceFEBasis(referenceGridView);
std::cout << "levels: " << gridView.grid().maxLevel()+1
<< " "
<< "L^2 error: " << std::sqrt(l2ErrorSquared)
<< " "
<< "H^1 error: " << std::sqrt(l2ErrorSquared + h1ErrorSquared)
<< std::endl;
}
// Reference configuration
EmbeddedVectorType embeddedReferenceX(referenceFEBasis.size());
inFile.open(parameterSet.get<std::string>("referenceData"), std::ios_base::binary);
if (not inFile)
DUNE_THROW(IOError, "File " << parameterSet.get<std::string>("referenceData") << " could not be opened.");
GenericVector::readBinary(inFile, embeddedReferenceX);
template <class GridView, int order>
void measureAnalyticalEOC(const GridView gridView,
const ParameterTree& parameterSet)
{
typedef std::vector<TargetSpace> SolutionType;
SolutionType referenceX(embeddedReferenceX.size());
for (size_t i=0; i<referenceX.size(); i++)
referenceX[i] = TargetSpace(embeddedReferenceX[i]);
//////////////////////////////////////////////////////////////////////////////////
// Construct the scalar function space bases corresponding to the GFE space
//////////////////////////////////////////////////////////////////////////////////
using FufemFEBasis = DuneFunctionsBasis<FEBasis>;
FufemFEBasis fufemReferenceFEBasis(referenceFEBasis);
FufemFEBasis fufemFEBasis(feBasis);
typedef Dune::Functions::PQkNodalBasis<GridView, order> FEBasis;
FEBasis feBasis(gridView);
// The numerical solution, as a grid function
GFE::EmbeddedGlobalGFEFunction<FufemFEBasis, TargetSpace> referenceSolution(fufemReferenceFEBasis, referenceX);
//////////////////////////////////////////////////////////////////////////////////
// Read the data whose error is to be measured
//////////////////////////////////////////////////////////////////////////////////
// The numerical solution, as a grid function
GFE::EmbeddedGlobalGFEFunction<FufemFEBasis, TargetSpace> numericalSolution(fufemFEBasis, x);
// Input data
typedef BlockVector<TargetSpace::CoordinateType> EmbeddedVectorType;
double l2ErrorSquared = 0;
double h1ErrorSquared = 0;
EmbeddedVectorType embeddedX(feBasis.size());
std::ifstream inFile(parameterSet.get<std::string>("simulationData"), std::ios_base::binary);
if (not inFile)
DUNE_THROW(IOError, "File " << parameterSet.get<std::string>("simulationData") << " could not be opened.");
GenericVector::readBinary(inFile, embeddedX);
inFile.close();
HierarchicSearch<typename GridView::Grid,typename GridView::IndexSet> hierarchicSearch(gridView.grid(), gridView.indexSet());
SolutionType x(embeddedX.size());
for (size_t i=0; i<x.size(); i++)
x[i] = TargetSpace(embeddedX[i]);
std::cout << "l2: " << l2ErrorSquared << std::endl;
for (const auto& rElement : elements(referenceGridView))
{
const auto& quadRule = QuadratureRules<double, dim>::rule(rElement.type(), 6);
/////////////////////////////////////////////////////////////////
// Measure the discretization error
/////////////////////////////////////////////////////////////////
for (const auto& qp : quadRule)
{
auto integrationElement = rElement.geometry().integrationElement(qp.position());
using FufemFEBasis = DuneFunctionsBasis<FEBasis>;
FufemFEBasis fufemFEBasis(feBasis);
auto globalPos = rElement.geometry().global(qp.position());
// Read reference solution and its derivative into a PythonFunction
typedef VirtualDifferentiableFunction<FieldVector<double, dim>, TargetSpace::CoordinateType> FBase;
auto element = hierarchicSearch.findEntity(globalPos);
auto localPos = element.geometry().local(globalPos);
Python::Module module = Python::import(parameterSet.get<std::string>("referenceSolution"));
auto referenceSolution = module.get("fdf").toC<std::shared_ptr<FBase>>();
auto diff = referenceSolution(rElement, qp.position()) - numericalSolution(element, localPos);
// The numerical solution, as a grid function
GFE::EmbeddedGlobalGFEFunction<FufemFEBasis, TargetSpace> numericalSolution(fufemFEBasis, x);
l2ErrorSquared += integrationElement * qp.weight() * diff.two_norm2();
// QuadratureRule for the integral of the L^2 error
QuadratureRuleKey quadKey(dim,6);
auto derDiff = referenceSolution.derivative(rElement, qp.position()) - numericalSolution.derivative(element, localPos);
// Compute the embedded L^2 error
double l2Error = DiscretizationError<GridView>::computeL2Error(&numericalSolution,
referenceSolution.get(),
quadKey);
h1ErrorSquared += integrationElement * qp.weight() * derDiff.frobenius_norm2();
// Compute the embedded H^1 error
double h1Error = DiscretizationError<GridView>::computeH1HalfNormDifferenceSquared(gridView,
&numericalSolution,
referenceSolution.get(),
quadKey);
}
}
std::cout << "l2: " << l2ErrorSquared << std::endl;
std::cout << "levels: " << gridView.grid().maxLevel()+1
<< " "
<< "L^2 error: " << std::sqrt(l2ErrorSquared)
<< " "
<< "H^1 error: " << std::sqrt(l2ErrorSquared + h1ErrorSquared)
<< std::endl;
}
std::cout << "elements: " << gridView.size(0)
<< " "
<< "L^2 error: " << l2Error
<< " ";
std::cout << "H^1 error: " << std::sqrt(l2Error*l2Error + h1Error) << std::endl;
}
......@@ -231,25 +260,53 @@ int main (int argc, char *argv[]) try
referenceGrid->globalRefine(parameterSet.get<int>("numReferenceLevels")-1);
// Do the actual measurement
const int order = parameterSet.get<int>("order");
switch (order)
if (parameterSet.get<std::string>("discretizationErrorMode")=="discrete")
{
switch (order)
{
case 1:
measureDiscreteEOC<GridType::LeafGridView,1>(grid->leafGridView(), referenceGrid->leafGridView(), parameterSet);
break;
case 2:
measureDiscreteEOC<GridType::LeafGridView,2>(grid->leafGridView(), referenceGrid->leafGridView(), parameterSet);
break;
case 3:
measureDiscreteEOC<GridType::LeafGridView,3>(grid->leafGridView(), referenceGrid->leafGridView(), parameterSet);
break;
default:
DUNE_THROW(NotImplemented, "Order '" << order << "' is not implemented");
}
}
if (parameterSet.get<std::string>("discretizationErrorMode")=="analytical")
{
case 1:
measureEOC<GridType::LeafGridView,1>(grid->leafGridView(), referenceGrid->leafGridView(), parameterSet);
break;
switch (order)
{
case 1:
measureAnalyticalEOC<GridType::LeafGridView,1>(grid->leafGridView(), parameterSet);
break;
case 2:
measureEOC<GridType::LeafGridView,2>(grid->leafGridView(), referenceGrid->leafGridView(), parameterSet);
break;
case 2:
measureAnalyticalEOC<GridType::LeafGridView,2>(grid->leafGridView(), parameterSet);
break;
case 3:
measureEOC<GridType::LeafGridView,3>(grid->leafGridView(), referenceGrid->leafGridView(), parameterSet);
break;
case 3:
measureAnalyticalEOC<GridType::LeafGridView,3>(grid->leafGridView(), parameterSet);
break;
default:
DUNE_THROW(NotImplemented, "Order '" << order << "' is not implemented");
default:
DUNE_THROW(NotImplemented, "Order '" << order << "' is not implemented");
}
}
return 0;
}
catch (Exception e)
......
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