Commit fe0d656b authored by Luca Pasquale's avatar Luca Pasquale
Browse files

Added ETA_ADR1DTest

parent 66e6a6dc
INCLUDE(TribitsAddExecutableAndTest)
INCLUDE(TribitsCopyFilesToBinaryDir)
INCLUDE_DIRECTORIES(${CMAKE_SOURCE_DIR})
TRIBITS_ADD_EXECUTABLE_AND_TEST(
ETA_ADR_1D
SOURCES main.cpp ETA_ADR1DTest.cpp
ARGS -c
NUM_MPI_PROCS 1
COMM serial mpi
)
//@HEADER
/*
*******************************************************************************
Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
This file is part of LifeV.
LifeV is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
LifeV is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with LifeV. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************
*/
//@HEADER
/**
@file ETA_ADR2DTest.cpp
@author L. Pasquale <luca.pasquale@mail.polimi.it>
@date 2012-11-20
*/
// ===================================================
//! Includes
// ===================================================
#include "ETA_ADR1DTest.hpp"
#include <lifev/eta/fem/ETFESpace.hpp>
#include <lifev/eta/expression/Integrate.hpp>
// ===================================================
//! Namespaces & define
// ===================================================
// ---------------------------------------------------------------
// We work in the LifeV namespace and define the mesh, matrix and
// vector types that we will need several times.
// ---------------------------------------------------------------
using namespace LifeV;
typedef RegionMesh<LinearLine> mesh_Type;
typedef MatrixEpetra<Real> matrix_Type;
typedef VectorEpetra vector_Type;
// ---------------------------------------------------------------
// We define then a function, which represents a velocity field,
// of magnitude 1 in the y direction (supposing an (x,y)
// reference frame).
// ---------------------------------------------------------------
// ===================================================
//! Functions
// ===================================================
Real betaFct ( const Real& /* t */, const Real& /* x */, const Real& /* y */, const Real& /* z */, const ID& /* i */ )
{
return 1;
}
Real fRhs ( const Real& /* t */, const Real& /*x*/, const Real& /*y*/, const Real& /* z */ , const ID& /* i */ )
{
return 2;
}
// ===================================================
//! Constructors
// ===================================================
ETA_ADR1DTest::ETA_ADR1DTest ()
{
#ifdef EPETRA_MPI
M_comm.reset ( new Epetra_MpiComm ( MPI_COMM_WORLD ) );
#else
M_comm.reset ( new Epetra_SerialComm() );
#endif
}
// ===================================================
//! Methods
// ===================================================
Real
ETA_ADR1DTest::run()
{
bool verbose (M_comm->MyPID() == 0);
// ---------------------------------------------------------------
// We define the mesh and parition it. We use the domain
// (-1,1) and a structured mesh.
// ---------------------------------------------------------------
if (verbose)
{
std::cout << " -- Building the mesh ... " << std::flush;
}
const UInt Nelements (10);
boost::shared_ptr< mesh_Type > fullMeshPtr (new mesh_Type);
regularMesh1D ( *fullMeshPtr, 0, Nelements, false,
1.0,
-1.0);
// Mesh partitioning doesn't work yet with 1D meshes
//MeshPartitioner< mesh_Type > meshPart (fullMeshPtr, M_comm);
//boost::shared_ptr< mesh_Type > meshPtr (meshPart.meshPartition() );
//fullMeshPtr.reset();
if (verbose)
{
std::cout << " done ! " << std::endl;
}
// ---------------------------------------------------------------
// We start by defining the finite element spaces. We use the type
// FESpace for the classical way and copy it in an ETFESpace for
// the ET assembly.
//
// In the end, both solution spaces display their respective
// number of degrees of freedom, which must be the same.
// ---------------------------------------------------------------
if (verbose)
{
std::cout << " -- Building FESpaces ... " << std::flush;
}
std::string uOrder ("P1");
std::string bOrder ("P1");
boost::shared_ptr<FESpace< mesh_Type, MapEpetra > > uSpace
( new FESpace< mesh_Type, MapEpetra > (fullMeshPtr, uOrder, 1, M_comm) );
boost::shared_ptr<FESpace< mesh_Type, MapEpetra > > betaSpace
( new FESpace< mesh_Type, MapEpetra > (fullMeshPtr, bOrder, 1, M_comm) );
if (verbose)
{
std::cout << " done ! " << std::endl;
}
if (verbose)
{
std::cout << " ---> Dofs: " << uSpace->dof().numTotalDof() << std::endl;
}
if (verbose)
{
std::cout << " -- Building ETFESpaces ... " << std::flush;
}
boost::shared_ptr<ETFESpace< mesh_Type, MapEpetra, 1, 1 > > ETuSpace
( new ETFESpace< mesh_Type, MapEpetra, 1, 1 > (fullMeshPtr, & (uSpace->refFE() ), & (uSpace->fe().geoMap() ), M_comm) );
boost::shared_ptr<ETFESpace< mesh_Type, MapEpetra, 1, 1 > > ETbetaSpace
( new ETFESpace< mesh_Type, MapEpetra, 1, 1 > (fullMeshPtr, & (betaSpace->refFE() ), & (betaSpace->fe().geoMap() ), M_comm) );
if (verbose)
{
std::cout << " done ! " << std::endl;
}
if (verbose)
{
std::cout << " ---> Dofs: " << ETuSpace->dof().numTotalDof() << std::endl;
}
// ---------------------------------------------------------------
// We interpolate then the advection function of the mesh at hand.
// This is performed with the classical FESpace only.
// ---------------------------------------------------------------
if (verbose)
{
std::cout << " -- Interpolating the advection field ... " << std::flush;
}
vector_Type beta (betaSpace->map(), Repeated);
betaSpace->interpolate (static_cast<FESpace< mesh_Type, MapEpetra >::function_Type> (betaFct), beta, 0.0);
if (verbose)
{
std::cout << " done! " << std::endl;
}
// ---------------------------------------------------------------
// We build now two matrices, one for each assembly procedure in
// order to compare them in the end. There is no difference in the
// construction of the matrices.
// ---------------------------------------------------------------
if (verbose)
{
std::cout << " -- Defining the matrices ... " << std::flush;
}
boost::shared_ptr<matrix_Type> systemMatrix (new matrix_Type ( uSpace->map() ) );
*systemMatrix *= 0.0;
boost::shared_ptr<matrix_Type> ETsystemMatrix (new matrix_Type ( ETuSpace->map() ) );
*ETsystemMatrix *= 0.0;
if (verbose)
{
std::cout << " done! " << std::endl;
}
// ---------------------------------------------------------------
// Definition of the RHS
if (verbose)
{
std::cout << " -- Defining and interpolating the RHS ... " << std::flush;
}
vector_Type rhs (uSpace->map(), Repeated);
rhs *= 0.0;
vector_Type ETrhs (uSpace->map(), Repeated);
ETrhs *= 0.0;
vector_Type fInterpolated (uSpace->map(), Repeated);
fInterpolated *= 0.0;
uSpace->interpolate ( static_cast<FESpace< mesh_Type, MapEpetra >::function_Type> (fRhs) , fInterpolated, 0.0 );
if (verbose)
{
std::cout << " done! " << std::endl;
}
// ---------------------------------------------------------------
// We come now to the assembly itself. We start with the classical
// way, which consists in instentiating an ADRAssembler and call
// the different methods from that class.
//
// Using a LifeChrono, we monitor the time required for the
// assembly.
//
// The terms assembled correspond, in the order of appearance, to
// - a laplacian term
// - an advection term
// - a reaction term (with coefficient 2.0), aka mass term.
//
// We also assemble the RHS
// ---------------------------------------------------------------
LifeChrono StdChrono;
StdChrono.start();
if (verbose)
{
std::cout << " -- Classical assembly ... " << std::flush;
}
ADRAssembler<mesh_Type, matrix_Type, vector_Type> adrAssembler;
adrAssembler.setup (uSpace, betaSpace);
adrAssembler.addDiffusion (systemMatrix, 1.0);
adrAssembler.addAdvection (systemMatrix, beta);
adrAssembler.addMass (systemMatrix, 2.0);
adrAssembler.addMassRhs (rhs, fInterpolated);
StdChrono.stop();
if (verbose)
{
std::cout << " done ! " << std::endl;
}
if (verbose)
{
std::cout << " Time : " << StdChrono.diff() << std::endl;
}
// ---------------------------------------------------------------
// We perform now the same assembly with the ET assembly and still
// monitor the timings required.
//
// Remark also that the Real constants and VectorSmall constants
// can be used directly in the expressions to integrate.
// ---------------------------------------------------------------
LifeChrono ETChrono;
ETChrono.start();
if (verbose)
{
std::cout << " -- ET assembly ... " << std::flush;
}
// Workaround for 1D case where grad(phi_j) is a VectorSmall<1>
// and not a double: we use a dot product with a unitary VectorSmall<1>
// to obtain its value as adouble
VectorSmall<1> oneVector;
oneVector[0]=1.0;
{
using namespace ExpressionAssembly;
integrate ( elements (ETuSpace->mesh() ),
uSpace->qr(),
ETuSpace,
ETuSpace,
dot( grad (phi_i) , grad (phi_j) )
+ dot ( grad (phi_j) , value(oneVector) ) * value (ETbetaSpace, beta)*phi_i
+ 2.0 * phi_i * phi_j
)
>> ETsystemMatrix;
integrate ( elements (ETuSpace->mesh() ),
uSpace->qr(),
ETuSpace,
value (ETuSpace, fInterpolated) *phi_i
)
>> ETrhs;
}
ETChrono.stop();
if (verbose)
{
std::cout << " done! " << std::endl;
}
if (verbose)
{
std::cout << " Time : " << ETChrono.diff() << std::endl;
}
// ---------------------------------------------------------------
// We finally need to check that both yield the same matrix. In
// that aim, we need to finalize both matrices.
// ---------------------------------------------------------------
if (verbose)
{
std::cout << " -- Closing the matrices and vectors... " << std::flush;
}
systemMatrix->globalAssemble();
ETsystemMatrix->globalAssemble();
rhs.globalAssemble();
ETrhs.globalAssemble();
if (verbose)
{
std::cout << " done ! " << std::endl;
}
// ---------------------------------------------------------------
// We compute now the matrix of the difference and finally the
// norm of the difference. This should be very low if the two
// matrices are identical.
// ---------------------------------------------------------------
if (verbose)
{
std::cout << " -- Computing the error ... " << std::flush;
}
boost::shared_ptr<matrix_Type> checkMatrix (new matrix_Type ( ETuSpace->map() ) );
*checkMatrix *= 0.0;
*checkMatrix += *systemMatrix;
*checkMatrix += (*ETsystemMatrix) * (-1);
checkMatrix->globalAssemble();
Real errorNorm ( checkMatrix->normInf() );
// ---------------------------------------------------------------
// Finally, we check that the two right hand sides are similar
// with the two assembly procedures. We also check that the two
// integrals correspond to their exact values.
// ---------------------------------------------------------------
vector_Type checkRhs (ETuSpace->map(), Repeated);
checkRhs = 0.0;
checkRhs += rhs;
checkRhs -= ETrhs;
checkRhs.globalAssemble();
vector_Type checkRhsUnique (checkRhs, Unique);
Real errorRhsNorm ( checkRhsUnique.normInf() );
if (verbose)
{
std::cout << " done ! " << std::endl;
}
if (verbose)
{
std::cout << " Matrix Error : " << errorNorm << std::endl;
}
if (verbose)
{
std::cout << " Rhs error : " << errorRhsNorm << std::endl;
}
// ---------------------------------------------------------------
// Finally we return the maximum of the two errors
// ---------------------------------------------------------------
return std::max (errorNorm, errorRhsNorm);
} // run
//@HEADER
/*
*******************************************************************************
Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
This file is part of LifeV.
LifeV is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
LifeV is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with LifeV. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************
*/
//@HEADER
/**
@file ETA_ADR2DTest.cpp
@author L. Pasquale <luca.pasquale@mail.polimi.it>
@date 2012-11-20
*/
#ifndef __ETA_ADR1DTEST_H
#define __ETA_ADR1DTEST_H 1
// Tell the compiler to ignore specific kind of warnings:
#pragma GCC diagnostic ignored "-Wunused-variable"
#pragma GCC diagnostic ignored "-Wunused-parameter"
#include <Epetra_ConfigDefs.h>
#ifdef EPETRA_MPI
#include <mpi.h>
#include <Epetra_MpiComm.h>
#else
#include <Epetra_SerialComm.h>
#endif
//Tell the compiler to restore the warning previously silented
#pragma GCC diagnostic warning "-Wunused-variable"
#pragma GCC diagnostic warning "-Wunused-parameter"
#include <lifev/core/LifeV.hpp>
#include <lifev/core/mesh/MeshPartitioner.hpp>
#include <lifev/core/mesh/RegionMesh1DStructured.hpp>
#include <lifev/core/mesh/RegionMesh.hpp>
#include <lifev/core/array/MatrixEpetra.hpp>
#include <lifev/core/array/VectorEpetra.hpp>
#include <boost/shared_ptr.hpp>
// ---------------------------------------------------------------
// We include two additional files which are useful to make the
// assembly using the "classical way". We also include a file to
// monitor the different timings.
// ---------------------------------------------------------------
#include <lifev/core/fem/FESpace.hpp>
#include <lifev/core/solver/ADRAssembler.hpp>
#include <lifev/core/util/LifeChrono.hpp>
// ---------------------------------------------------------------
// We work in the LifeV namespace and define the mesh, matrix and
// vector types that we will need several times.
// ---------------------------------------------------------------
/*!
@class ETA_ADR2DTest
@brief Simple ETA test to compare ETA to ADRAssembly in a 2D space
@author L. Pasquale <luca.pasquale@mail.polimi.it>
This test computes the matrix associated to the advection-diffusion-reaction problem: \n
\f$
-\alpha \delta u + \beta \cdot \nabla u +\gamma u = f
\f$ \n
with \f$\alpha=1\f$, \f$\beta=\left(\begin{array}{c}1\\0\end{array}\right)\f$, \f$\gamma=2\f$, \f$f=2\f$ and \f$u\f$ a scalar field \n
The domain is the square (-1,1)x(-1,1). \n
The matrix and the vector for the RHS are first assembled using the class ADRAssembler and then the ExpressionAssembly::integrate function of the ETA module. \n
if the two matrices (and vectors) are assembled correctly, they should be almost exactly the same. \n
*/
class ETA_ADR1DTest
// :
// public LifeV::Application
{
public:
/* @name Constructors and destructor
*/
//@{
//! Constructor
ETA_ADR1DTest ();
//! Destructor
~ETA_ADR1DTest () {}
//@}
/* @name Methods
*/
//@{
//! To lunch the simulation
LifeV::Real run();
//@}
private:
boost::shared_ptr<Epetra_Comm> M_comm;
bool M_verbose;
};
#endif /* __ETA_ADR2DTEST_H */
//@HEADER
/*
*******************************************************************************
Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
This file is part of LifeV.
LifeV is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
LifeV is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with LifeV. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************
*/
//@HEADER
/* ========================================================
Simple ETA test to compare ETA to ADRAssembly in a 2D space
*/