I revised the appendix from my last message a little bit and attache now a
minimal working example (just 140 lines) along with a CMakeLists.txt.
After checking the profiling results from valgrind, the combination of
MappingQ with FE_Q takes *not* the fast path.
For info: I use dealii version 9.3.2
Best,
Simon
Am Do., 20. Okt. 2022 um 18:11 Uhr schrieb Simon Wiesheier <
simon.wieshe...@gmail.com>:
> " When you use FEPointEvaluation, you should construct it only once and
> re-use the same object for different points. Furthermore, you should also
> avoid to create "p_dofs" and the "std::vector" near the I was not clear
> with my original message. Anyway, the problem is the FEValues object that
> gets used. I am confused by your other message that you use FE_Q together
> with MappingQ - that combination should be supported and if it is not, we
> should take a look at a (reduced) code from you. "
>
> I added a snippet of my code (see appendix) in which I describe the logic
> as to what I am doing with FEPointEvaluation.
> In fact, constructing FEPointEvaluation (and the vector p_dofs) once and
> re-using them brings only minor changes as the overall costs are dominated
> by the call to reinit().
> But, of course, it helps at least.
>
> I am surprised too that the fast path is not used. Maybe you can identify
> a problem in my code.
> Thank you!
>
> Best,
> Simon
>
> Am Do., 20. Okt. 2022 um 17:02 Uhr schrieb Martin Kronbichler <
> kronbichler.mar...@gmail.com>:
>
>> Dear Simon,
>>
>> When you use FEPointEvaluation, you should construct it only once and
>> re-use the same object for different points. Furthermore, you should also
>> avoid to create "p_dofs" and the "std::vector" near the I was not clear
>> with my original message. Anyway, the problem is the FEValues object that
>> gets used. I am confused by your other message that you use FE_Q together
>> with MappingQ - that combination should be supported and if it is not, we
>> should take a look at a (reduced) code from you.
>>
>> Regarding the high timings: There is some parallelization by tasks that
>> gets done inside the constructor of FEValues. This has good intents for the
>> case that we are in 3D and have a reasonable amount of work to do. However,
>> you are in 1D (if I read your code correctly), and then it is having
>> adverse effects. The reason is that the constructor of FEValues is very
>> likely completely dominated by memory allocation. When we have 1 thread,
>> everything is fine, but when we have multiple threads working they will
>> start to interfere with each other when the request memory through
>> malloc(), which has to be coordinated by the operating system (and thus
>> gets slower). In fact, the big gap between compute time and wall time shows
>> that there is a lot of time wasted by "system time" that does not do actual
>> work on the cores.
>>
>> I guess the library could have a better measure of when to spawn tasks in
>> FEValues in similar context, but it is a lot of work to get this right.
>> (This is why I keep avoiding it in critical functions.)
>>
>> Best,
>> Martin
>>
>>
>> On 20.10.22 16:47, Simon Wiesheier wrote:
>>
>> Update:
>>
>> I profiled my program with valgrind --tool=callgrind and could figure out
>> that
>> FEPointEvaluation creates an FEValues object along with a quadrature
>> object under the hood.
>> Closer inspection revealed that all constructors, destructors,...
>> associated with FEPointEvaluation
>> need roughly 5000 instructions more (per call!).
>> That said, FEValues is indeed the faster approach, at least for FE_Q
>> elements.
>>
>> export DEAL_II_NUM_THREADS=1
>> eliminated the gap between cpu and wall time.
>> Using FEValues directly, I get cpu time 19.8 seconds
>> and in the case of FEPointEvaluation cpu time = 21.9 seconds;
>> Wall times are in the same ballpark.
>> Out of curiosity, why produces multi-threading such high wall times (200
>> seconds) in my case?.
>>
>> These times are far too big given that the solution of the linear system
>> takes only about 13 seconds.
>> But based on what all of you have said, there is probably no other to way
>> to implement my problem.
>>
>> Best
>> Simon
>>
>> Am Do., 20. Okt. 2022 um 11:55 Uhr schrieb Simon Wiesheier <
>> simon.wieshe...@gmail.com>:
>>
>>> Dear Martin and Wolfgang,
>>>
>>> " You seem to be looking for FEPointEvaluation. That class is shown in
>>> step-19 and provides, for simple FiniteElement types, a much faster way to
>>> evaluate solutions at arbitrary points within a cell. Do you want to give
>>> it a try? "
>>>
>>> I implemented the FEPointEvaluation approach like this:
>>>
>>> FEPointEvaluation<1,1> fe_eval(mapping,
>>> FE_Q<1>(1),
>>> update_gradients |
>>> update_values);
>>> fe_eval.reinit(cell,
>>> make_array_view(std::vector<Point<1>>{ref_point_energy_vol}));
>>> Vector<double> p_dofs(2);
>>> cell->get_dof_values(solution_global, p_dofs);
>>> fe_eval.evaluate(make_array_view(p_dofs),
>>> EvaluationFlags::values |
>>> EvaluationFlags::gradients);
>>> double val = fe_eval.get_value(0);
>>> Tensor<1,1> grad = fe_eval.get_gradient(0);
>>>
>>> I am using FE_Q elements of degree one and a MappingQ object also of
>>> degree one.
>>>
>>> Frankly, I do not really understand the measured computation times.
>>> My program has several loadsteps with nested Newton iterations:
>>> Loadstep 1:
>>> Assembly 1: cpu time 12.8 sec wall time 268.7 sec
>>> Assembly 2: cpu time 17.7 sec wall time 275.2 sec
>>> Assembly 3: cpu time 22.3 sec wall time 272.6 sec
>>> Assembly 4: cpu time 23.8 sec wall time 271.3sec
>>> Loadstep 2:
>>> Assembly 1: cpu time 14.3 sec wall time 260.0 sec
>>> Assembly 2: cpu time 16.9 sec wall time 262.1 sec
>>> Assembly 3: cpu time 18.5 sec wall time 270.6 sec
>>> Assembly 4: cpu time 17.1 sec wall time 262.2 sec
>>> ...
>>>
>>> Using FEValues instead of FEPointEvaluation, the results are:
>>> Loadstep 1:
>>> Assembly 1: cpu time 23.9 sec wall time 171.0 sec
>>> Assembly 2: cpu time 32.5 sec wall time 168.9 sec
>>> Assembly 3: cpu time 33.2 sec wall time 168.0 sec
>>> Assembly 4: cpu time 32.7 sec wall time 166.9 sec
>>> Loadstep 2:
>>> Assembly 1: cpu time 24.9 sec wall time 168.0 sec
>>> Assembly 2: cpu time 34.7 sec wall time 167.3 sec
>>> Assembly 3: cpu time 33.9 sec wall time 167.8 sec
>>> Assembly 4: cpu time 34.3 sec wall time 167.7 sec
>>> ...
>>>
>>> Clearly, the fluctuations using FEValues are smaller than in case of
>>> FEPointEvaluation.
>>> Anyway, using FEPointEvaluation the cpu time is smaller but the wall
>>> time substantially bigger.
>>> If I am not mistaken, the values cpu time 34.3 sec and wall time 167.7
>>> sec mean that
>>> the cpu needs 34.3 sec to execute my assembly routine and has to wait in
>>> the
>>> remaining 167.7-34.3 seconds.
>>> This huge gap between cpu and wall time has to be related to what I do
>>> with FEValues or FEPointEvaluation
>>> as cpu and wall time are nearly balanced if I use either neither of them.
>>> What might be the problem?
>>>
>>> Best
>>> Simon
>>>
>>>
>>>
>>>
>>>
>>> Am Mi., 19. Okt. 2022 um 22:34 Uhr schrieb Wolfgang Bangerth <
>>> bange...@colostate.edu>:
>>>
>>>> On 10/19/22 08:45, Simon Wiesheier wrote:
>>>> >
>>>> > What I want to do boils down to the following:
>>>> > Given the reference co-ordinates of a point 'p', along with the cell
>>>> on
>>>> > which 'p' lives,
>>>> > give me the value and gradient of a finite element function evaluated
>>>> at
>>>> > 'p'.
>>>> >
>>>> > My idea was to create a quadrature object with 'p' being the only
>>>> > quadrature point and pass this
>>>> > quadrature object to the FEValues object and finally do the
>>>> > .reinit(cell) call (then, of course, get_function_values()...)
>>>> > 'p' is different for all (2.5 million) quadrature points, which is
>>>> why I
>>>> > create the FEValues object so many times.
>>>>
>>>> It's worth pointing out that is exactly what
>>>> VectorTools::point_values()
>>>> does.
>>>>
>>>> (As others have already mentioned, if you want to do that many many
>>>> times over, this is too expensive and you should be using
>>>> FEPointEvaluation instead.)
>>>>
>>>> Best
>>>> W.
>>>>
>>>> --
>>>> ------------------------------------------------------------------------
>>>> Wolfgang Bangerth email:
>>>> bange...@colostate.edu
>>>> www:
>>>> http://www.math.colostate.edu/~bangerth/
>>>>
>>>> --
>>>> The deal.II project is located at http://www.dealii.org/
>>>> For mailing list/forum options, see
>>>> https://groups.google.com/d/forum/dealii?hl=en
>>>> ---
>>>> You received this message because you are subscribed to the Google
>>>> Groups "deal.II User Group" group.
>>>> To unsubscribe from this group and stop receiving emails from it, send
>>>> an email to dealii+unsubscr...@googlegroups.com.
>>>> To view this discussion on the web visit
>>>> https://groups.google.com/d/msgid/dealii/cd1c8fa0-443d-b7bf-b433-f5ab033a247c%40colostate.edu
>>>> .
>>>>
>>> --
>> The deal.II project is located at http://www.dealii.org/
>> For mailing list/forum options, see
>> https://groups.google.com/d/forum/dealii?hl=en
>> ---
>> You received this message because you are subscribed to the Google Groups
>> "deal.II User Group" group.
>> To unsubscribe from this group and stop receiving emails from it, send an
>> email to dealii+unsubscr...@googlegroups.com.
>> To view this discussion on the web visit
>> https://groups.google.com/d/msgid/dealii/CAM50jEt3LjNS%2BF%3D2pGUSHRNbaHy9EzRfogpyBxb51-7M%3DALxrA%40mail.gmail.com
>> <https://groups.google.com/d/msgid/dealii/CAM50jEt3LjNS%2BF%3D2pGUSHRNbaHy9EzRfogpyBxb51-7M%3DALxrA%40mail.gmail.com?utm_medium=email&utm_source=footer>
>> .
>>
>> --
>> The deal.II project is located at http://www.dealii.org/
>> For mailing list/forum options, see
>> https://groups.google.com/d/forum/dealii?hl=en
>> ---
>> You received this message because you are subscribed to a topic in the
>> Google Groups "deal.II User Group" group.
>> To unsubscribe from this topic, visit
>> https://groups.google.com/d/topic/dealii/uAplhH99yg4/unsubscribe.
>> To unsubscribe from this group and all its topics, send an email to
>> dealii+unsubscr...@googlegroups.com.
>> To view this discussion on the web visit
>> https://groups.google.com/d/msgid/dealii/24ba8ef3-39f4-af6c-7296-02677827d6d3%40gmail.com
>> <https://groups.google.com/d/msgid/dealii/24ba8ef3-39f4-af6c-7296-02677827d6d3%40gmail.com?utm_medium=email&utm_source=footer>
>> .
>>
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/CAM50jEt1Y2rYwq_2gW1vXsUmVYdYnvh4hP5taKuCw7PATYOGAA%40mail.gmail.com.
##
# CMake script for the step-1 tutorial program:
##
# Set the name of the project and target:
SET(TARGET "minimal_working_example")
# Declare all source files the target consists of. Here, this is only
# the one step-X.cc file, but as you expand your project you may wish
# to add other source files as well. If your project becomes much larger,
# you may want to either replace the following statement by something like
# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
SET(TARGET_SRC
${TARGET}.cc
)
# Usually, you will not need to modify anything beyond this point...
CMAKE_MINIMUM_REQUIRED(VERSION 3.1.0)
FIND_PACKAGE(deal.II 9.3.0
HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
)
IF(NOT ${deal.II_FOUND})
MESSAGE(FATAL_ERROR "\n"
"*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
"You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to
cmake\n"
"or set an environment variable \"DEAL_II_DIR\" that contains this path."
)
ENDIF()
DEAL_II_INITIALIZE_CACHED_VARIABLES()
PROJECT(${TARGET})
DEAL_II_INVOKE_AUTOPILOT()
#include <deal.II/base/point.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/timer.h>
#include <deal.II/matrix_free/fe_point_evaluation.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/mapping_q_eulerian.h>
#include <deal.II/lac/vector.h>
#include <deal.II/numerics/vector_tools.h>
using namespace dealii;
class EnergyVol
{
public:
EnergyVol();
protected:
Triangulation<1> triangulation_vol;
DoFHandler<1> dof_handler_vol;
FE_Q<1>
fe_vol;
QGauss<1>
qf_cell_vol;
MappingQ<1>
mapping_vol;
FEPointEvaluation<1,1> fe_eval_vol;
Vector<double> solution_vol;
//some member functions like make_grid(), system_setup(),...
};
EnergyVol::EnergyVol()
:
dof_handler_vol(triangulation_vol),
fe_vol(1),
qf_cell_vol(2),
mapping_vol(1),
fe_eval_vol(mapping_vol,
fe_vol,
update_gradients | update_values)
{
GridGenerator::hyper_cube(triangulation_vol,
1.0, 2.0);
triangulation_vol.refine_global(2);
dof_handler_vol.distribute_dofs(fe_vol);
solution_vol.reinit(dof_handler_vol.n_dofs());
VectorTools::interpolate<1>( dof_handler_vol,
Functions::ZeroFunction<1>(),
solution_vol
);
}
class EnergyDev
{
public:
EnergyDev();
protected:
Triangulation<1>
triangulation_dev;
DoFHandler<1> dof_handler_dev;
FE_Q<1> fe_dev;
QGauss<1>
qf_cell_dev;
MappingQ<1> mapping_dev;
FEPointEvaluation<1,1> fe_eval_dev;
Vector<double> solution_dev;
};
EnergyDev::EnergyDev()
:
dof_handler_dev(triangulation_dev),
fe_dev(1),
qf_cell_dev(2),
mapping_dev(1),
fe_eval_dev(mapping_dev,
fe_dev,
update_gradients | update_values)
{
GridGenerator::hyper_cube(triangulation_dev,
1.0, 2.0);
triangulation_dev.refine_global(4);
dof_handler_dev.distribute_dofs(fe_dev);
DoFRenumbering::Cuthill_McKee(dof_handler_dev);
solution_dev.reinit(dof_handler_dev.n_dofs());
VectorTools::interpolate<1>( dof_handler_dev,
Functions::ZeroFunction<1>(),
solution_dev
);
}
class Energy : public EnergyVol, public EnergyDev
{
public:
Energy();
~Energy();
void compute_stress();
};
Energy::Energy()
{}
Energy::~Energy()
{
dof_handler_vol.clear();
this->dof_handler_dev.clear();
}
//call this function from main
void Energy::compute_stress()
{
TimerOutput timer(std::cout, TimerOutput::every_call_and_summary,
TimerOutput::cpu_and_wall_times);
for(unsigned int i=0; i<1000; i++)
{
typename DoFHandler<1>::active_cell_iterator cell_vol =
dof_handler_vol.begin_active(); // of course, different in every call
Point<1> point_vol(0.5); // of course, different in every call
fe_eval_vol.reinit(cell_vol,
make_array_view(std::vector<Point<1>>{point_vol}));
Vector<double> p_dofs(2);
cell_vol->get_dof_values(solution_vol, p_dofs);
fe_eval_vol.evaluate(make_array_view(p_dofs),
EvaluationFlags::values |
EvaluationFlags::gradients);
//
typename DoFHandler<1>::active_cell_iterator cell_dev =
this->dof_handler_dev.begin_active(); // of course, different in every call
Point<1> point_dev(0.5); // of course, different in every call
this->fe_eval_dev.reinit(cell_dev,
make_array_view(std::vector<Point<1>>{point_dev}));
Vector<double> dev_dofs(2);
cell_dev->get_dof_values(this->solution_dev, dev_dofs);
this->fe_eval_dev.evaluate(make_array_view(dev_dofs),
EvaluationFlags::values |
EvaluationFlags::gradients);
double p = fe_eval_vol.get_value(0);
Tensor<1,1> dp_dJ = fe_eval_vol.get_gradient(0);
double s = fe_eval_dev.get_value(0);
Tensor<1,1> s_wrt_IC = fe_eval_dev.get_gradient(0);
// do something with the values and gradients ...
}
}
int main ()
{
using namespace dealii;
Energy energy;
energy.compute_stress();
}
