If you want to use symbolic calculations, you could also leverage SymEngine, and use Functions::SymbolicFunction
https://www.dealii.org/current/doxygen/deal.II/classFunctions_1_1SymbolicFunction.html In your code, you could create such object by making a std::unique_ptr<Functions::SymbolicFunction<dim>> fun; parse its expression from file, and then fun = std::make_unique<Functions::SymbolicFunction<dim>>(expression); In this way the generated function would compute the gradient symbolically. L. > On 9 Aug 2021, at 16:39, Marco Feder <marco.fede...@gmail.com> wrote: > > You're right, I totally missed the default parameter h in the constructor. > > Thanks, > Marco > Il giorno domenica 8 agosto 2021 alle 22:26:01 UTC+2 Timo Heister ha scritto: > FunctionParser uses a finite difference approximation for the gradient. I > think this is explained in the class documentation. > This explains the results you see... > > On Sun, Aug 8, 2021, 15:15 Marco Feder <marco....@gmail.com> wrote: > Ciao Luca! > > > Can you check if you get the same results with this order? > Yes, the rates are the same. > > > Are you calling `error_from_exact` with mapping as a first argument? > Yes > > Actually, I've fixed the issue after reading this: > > If the Solution<dim> class implements the Gradient, then > > ParsedConvergenceTable should use that. > > Indeed, I realised I wrote: > error_table.error_from_exact(mapping, dof_handler, solution, exact_solution); > > where exact_solution is a FunctionParser<dim> which will be initialised with > the expression of the analytical solution, so it doesn't give any info about > the Gradient. Instead, if I replace exact_solution in the last argument with > the "classical" Solution<dim>() (so that there's also an implementation of > Gradient) I have my expected rates. Apparently it wasn't using the Gradient > function. > > However, looking at the source code and in the documentation of > integrate_difference, I don't see how the H^1 norm was computed without > having an implementation of the Gradient. I would expect to see an exception > asking for its implementation. Was it using some FD-like approximation or am > I missing something? > > > Thanks, > Marco > > > Il giorno domenica 8 agosto 2021 alle 18:52:03 UTC+2 luca....@gmail.com ha > scritto: > The other option is that you are not calling the `error_from_exact` function > that takes a mapping argument, but the one without it, and your boundary > cells are introducing some error due to a worse mapping…. Are you calling > `error_from_exact` with mapping as a first argument? > > Luca > >> Il giorno 8 ago 2021, alle ore 6:49 PM, Luca Heltai <luca....@gmail.com> ha >> scritto: >> >> >> Ciao Marco, >> >>> I was expanding step-12 using a manufactured solution to check the order p >>> in H^1 (and p+1 in L^2) norm on a uniformly refined mesh for the DG upwind >>> method. I've used the ParsedConvergenceTable with a parameter file as >>> described in the docs. As Rate_key I'm using the DoFs, while as Rate_mode I >>> have reduction_rate_log2. >>> >>> With p=1 and p=2 everything is fine. But if I set the finite element degree >>> to 3, then the H^1 convergence rate decreases, as you can see in the >>> attached image. >>> <Screenshot 2021-08-07 at 17.56.06.png> >> >> Is it possible that you are using different quadrature rules in the two >> cases? Your image shows a deterioration of the error on the order of 1e-8 >> for H1, and 1e-12 for L2, which is very close to machine precision. >> >> Internally the parsed convergence table does the exact same thing you wrote >> explicitly. (If you check the source code, you’ll see that it calls >> integrate_difference and compute_global_error for each error type you >> specify in the parameter file). >> >>> This, however, doesn't happen if I use a classical ConvergenceTable. >>> Namely, I first compute the local error in each cell, and then the global >>> error in the classical way: >>> VectorTools::integrate_difference(mapping,dof_handler,solution, >>> Solution<dim>(),H1_error_per_cell, >>> QGauss<dim>(fe->tensor_degree()+1),VectorTools::H1_norm); >>> >>> const double H1_error = VectorTools::compute_global_error(triangulation, >>> H1_error_per_cell, VectorTools::H1_norm); //assuming I provided also the >>> gradient method for the Solution<dim> class >>> >>> >>> >>> Does anyone have any idea why this is happening? My guess was that while >>> computing the H^1 semi-norm the ParsedConvergenceTable class does some >>> approximation to compute the gradient from the exact solution expression >>> and hence that could be the source of the issue. Conversely, in the >>> "ConvergenceTable" way I do define explicitly the gradient of the exact >>> solution in the Solution<dim> class. >> >> If the Solution<dim> class implements the Gradient, then >> ParsedConvergenceTable should use that. >> >> You are calling “error_from_exact” of that class, right? >> >> This is where it is called: >> >> https://www.dealii.org/current/doxygen/deal.II/parsed__convergence__table_8h_source.html >> >> Line 621. As you see, the quadrature rule used is a Gauss formula of order >> (dh.get_fe().degree+1)*2. >> >> Can you check if you get the same results with this order? >> >> L. > > -- > 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+un...@googlegroups.com. > To view this discussion on the web visit > https://groups.google.com/d/msgid/dealii/dc72f870-1549-471c-b93b-854bfe3b67d9n%40googlegroups.com. > > -- > 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/bd3c5f62-86dc-4c39-b3da-a64418e7bc8cn%40googlegroups.com. -- 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/9EF78B23-6895-4AC6-9325-652BA2DC8C77%40gmail.com.
