Q: Question about extracting part of a vector

[Najwa Alshehri]

Dear group members,

I have one question, 

I am trying to calculate the L2 norm of the error of a solution. In particular, I have a vector of the solution u_omega and a FeSystem with two fes.

I am interested in computing the L2 norm of the error related to the first fe using integrate_difference  function. (Note here u_omega is not a blocked vector, however, it has two solutions stacked together in one vector). Can I extract the solution of the first part from u_omega?

Thank you in advance for your help.

Best,

Najwa

[Daniel Arndt]

Najwa,

The documentation of VectorTools::integrate_difference(https://www.dealii.org/current/doxygen/deal.II/namespaceVectorTools.html#a676190d2c897ac5da68a9c460fa95832) says

The additional argument weight allows to evaluate weighted norms. The weight function may be scalar, establishing a spatially variable weight in the domain for all components equally. This may be used, for instance, to only integrate over parts of the domain. The weight function may also be vector-valued, with as many components as the finite element: Then, different components get different weights. A typical application is when the error with respect to only one or a subset of the solution variables is to be computed, in which case the other components would have weight values equal to zero. The ComponentSelectFunction class is particularly useful for this purpose as it provides such a "mask" weight. The weight function is expected to be positive, but negative values are not filtered. The default value of this function, a null pointer, is interpreted as "no weighting function", i.e., weight=1 in the whole domain for all vector components uniformly. Best, Daniel

--
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/9741338e-2a31-418b-815d-277a5d7cb573n%40googlegroups.com.

[Najwa Alshehri]

Thank you Daniel for the clear quick answer. I will follow it.

Best,

Najwa

[Najwa Alshehri]

Hello again,

I have a follow-up question. Does this ComponentSelectFunction work also with vectors that are not written as blocked vectors? I have applied it before when I was dealing with blocked vectors and it worked perfectly. 

So, I did this, 

const ComponentSelectFunction<dim> primal_mask(0,2);

Later,

VectorTools::integrate_difference(omega2_dh,

u_omega2,

ExactSolution2<dim>(),

cellwise_errors2,

quadrature,

VectorTools::L2_norm

,&primal_mask);

const double u2_l2_error =

VectorTools::compute_global_error(triangulation_omega2,

cellwise_errors2,

VectorTools::L2_norm);

And I got the following error!!

An error occurred in line <455> of file <../include/deal.II/numerics/vector_tools_integrate_difference.templates.h> in function
    void dealii::VectorTools::internal::do_integrate_difference(const dealii::hp::MappingCollection<dim, spacedim>&, const dealii::DoFHandler<dim, spacedim>&, const InVector&, const dealii::Function<spacedim, typename InVector::value_type>&, OutVector&, const dealii::hp::QCollection<dim>&, const dealii::VectorTools::NormType&, const dealii::Function<spacedim>*, double) [with int dim = 2; int spacedim = 2; InVector = dealii::Vector<double>; OutVector = dealii::Vector<double>; typename InVector::value_type = double]
The violated condition was:
    ::dealii::deal_II_exceptions::internals::compare_for_equality(exact_solution.n_components, n_components)
Additional information:
    Dimension 1 not equal to 2.

 Obviously, we have dimensionality mismatching between u_omega2 and the exact solution. which means that the mask is not really picking up the first component of the solution u_omega2.  Do you have any suggestions to fix this issue?

Appreciate your help,

Najwa

[Daniel Arndt]

Najwa,

the relevant code reads something like

      const unsigned int n_components = dof.get_fe(0).n_components();
      AssertDimension(exact_solution.n_components, n_components);
      if (weight != nullptr)
        {
          Assert((weight->n_components == 1) ||
                   (weight->n_components == n_components),
                 ExcDimensionMismatch(weight->n_components, n_components));
        }

which means that your exact solution still needs to have as many components as the FiniteElement object used (and the number of components for ComponentSelectFunction must also match that).

ComponentSelectFunction then makes sure that all other components are ignored when computing the error.

Best,

Daniel

To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/c2544a5b-16b0-4d3d-8ccc-22e77ccf3038n%40googlegroups.com.

[Luca Heltai]

Dear Najwa,

how many components do the `ExactSolution2<dim>` function has?

The error complains about `exact_solution.n_components, n_components` not being equal, in particular it tells you that `ExactSolution2<dim>` only has one component.

If you built the ExactSolution class yourself, make sure at construction time you tell Function<dim> to allocate two components:

ExactSolution2<dim>::ExactSolution2<dim>( … )
: Functions::Function<dim>(2)

Best,
Luca.

> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/c2544a5b-16b0-4d3d-8ccc-22e77ccf3038n%40googlegroups.com.

[Najwa Alshehri]

Dear Daniel, your answer makes sense, Finally it worked.

Dear Luca, thank you for your answer. Indeed, the ExactSolution is a class I built myself. I followed your suggestion and it was the magic line that solved the issue. However, say that I have in another code instead of using a ExactSolution that I build myself, I used an FeFieldFunction for a reference solution, how would I make this function allocate two components?

Best,

Najwa

[Luca Heltai]

FEFieldFunction has the same number of components of the underlying dofhandler. If that has one component, then so does the FEFieldFunction. What you could do is wrap the FEFieldFunction into another function, e.g.

You could pass a lambda function as an argument to the constructor of the function, that uses your FEFieldFunction object, I.e.,

VectorFucntionFromScalarFuncObj<dim> my_fun([&](const auto &p){ return fe_field_func.value(p);}, 0, 2);

Luca

Il giorno 21 lug 2023, alle ore 22:05, Najwa Alshehri <najwaa...@gmail.com> ha scritto:

Dear Daniel, your answer makes sense, Finally it worked.

To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/697c9d76-0969-4c32-bd57-58371d02c21dn%40googlegroups.com.

[Najwa Alshehri]

Dear Luca,

This is interesting. I am exited to try this out. 

Thank you so much,

Najwa

[Najwa Alshehri]

Greetings everyone,

I want to express my gratitude for your assistance with my question, which has evolved into something more profound.

Regarding the function "VectorFucntionFromScalarFunctionObject," I have observed that it works for computing the L2 norm of the error. However, when I attempted to compute the H1_seminorm of the error, I noticed that the gradient was not implemented for this function. Please let me know if my understanding is incorrect, as I referred to the <deal.II/base/function.h> for this.

I am curious if there is a way to resolve this issue. My idea is to obtain the gradient of the "FeFeildFunction" since I know that the gradient is implemented for this function, as per <deal.II/numerics/fe_field_function.h>. Afterward, I can define a new function within a class, based on "VectorFucntionFromScalarFunctionObject," where the gradient of that function would be the gradient of the "FeFeildFunction."

Could you kindly advise me on the simplest way to accomplish this?

On a side note, although "VectorFucntionFromScalarFunctionObject" works with the L2 norm of the error, it is quite computationally intensive.

Best regards,
Najwa

[Wolfgang Bangerth]

On 7/26/23 12:17, Najwa Alshehri wrote:
>
> Regarding the function "VectorFucntionFromScalarFunctionObject," I have
> observed that it works for computing the L2 norm of the error. However, when I
> attempted to compute the H1_seminorm of the error, I noticed that the gradient
> was not implemented for this function. Please let me know if my understanding
> is incorrect, as I referred to the <deal.II/base/function.h> for this.

Yes, correct. Someone needs to implement
VectorFucntionFromScalarFunctionObject::gradient()
and/or
VectorFucntionFromScalarFunctionObject::gradient_list()
in exactly the same way as the value() and value_list() functions are implemented.

> Could you kindly advise me on the simplest way to accomplish this?

The easiest way is to create a class derived from
VectorFucntionFromScalarFunctionObject in your own project that implements the
missing function. The better way is to write a patch to the deal.II sources
that implements the missing function because then others can use your work in
the future as well :-)

> On a side note, although "VectorFucntionFromScalarFunctionObject" works with
> the L2 norm of the error, it is quite computationally intensive.

It isn't VectorFucntionFromScalarFunctionObject that's expensive, but the
FEFieldFunction. That's sort of inherently the case with this class. If you
can, it should be avoided and you should try an evaluate the FE field you have
at regular quadrature points, via FEValues.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Najwa Alshehri]

Dear Wolfgang,

Thank you for your answer. I will try to write a patch to the deal.II sources that implement the missing function. This might require some time. I will come back here if needed.

Thank you all for your quick answers,

Najwa

[Abbas Ballout]

Najwaa, 

I submitted a pull request recently about something similar. (I guess)
Maybe that would help. 

Abbas 
 

[Najwa Alshehri]

Hello Abbas,

Thank you for your response. I have actually written something, but I need to test it. It would be great to have a look at your work as well, so we can compare.

Best regards,

Najwa

--
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/DqH2auneWaY/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/dd5f93f8-fb58-45c6-9ab2-91ef6c5bdf56n%40googlegroups.com.

[Abbas Ballout]

Najwa, 
The files I edited are in the pull request here: https://github.com/dealii/dealii/pull/15805.
Let me know if I can help in a better way or of there is something I can do. 

Abbas. 

[Najwa Alshehri]

Dear team,

 

I want to express my gratitude for the support you have provided thus far. I would like to provide you with an update on the current situation and seek your assistance regarding a couple of issues.

 

Firstly, I have successfully implemented the missing function `VectorFunctionFromScalarFunctionObject::gradient()` in my project. Although it appears to be functioning correctly, the execution speed is extremely slow. In fact, it took until the next working day to obtain the results for only two cycles of refinement. Based on the response from @Wolfgang, it seems that the issue lies with the FeFieldFunction.

 

Allow me to provide some context. I have a reference solution serialized on a fine mesh. In my current code, my objective is to deserialize this solution and utilize it as the exact solution. To achieve this, I have employed the FeFieldFunction. Now, I have two scenarios in which I need to utilize the FeFieldFunction.

 

The first scenario involves directly using the FeFieldFunction when my computed solution consists of a single component. In this case, I have successfully calculated the L2 and H1 norms of the error using the `integrate_difference` and `compute_global_error` functions. The calculations proceeded smoothly, and the execution speed was reasonable.

 

However, the second scenario involves using the FeFieldFunction when the computed solution consists of two components. My objective is to calculate the error for the first component only. To accomplish this, I employed the new VectorFunctionFromScalarFunctionObject class locally. Unfortunately, the computation of the L2 norm in this case was considerably slow and the H1 norm was very slow.

 

Given the current situation, I have two questions that I hope you can assist me with, as you have done in the past:

 

1. Regarding my specific project, is there an alternative method, other than the FeFieldFunction, to save the deserialized solution as a function and utilize it as the exact solution? It is crucial that the solution be saved in the first component of a vector with two components.

 

2. In a more general sense, how can I test this function to determine if the sluggishness is solely attributable to the FeFieldFunction? Additionally, I would like to verify whether the function is functioning correctly or not, to be able to send it in a pull request to the dealii git repository.

 

Your guidance and expertise in resolving these issues would be greatly appreciated.

 

Thank you once again for your continued support.

 

Best regards,

Najwa

[Wolfgang Bangerth]

On 8/3/23 05:01, Najwa Alshehri wrote:
> The first scenario involves directly using the FeFieldFunction when my
> computed solution consists of a single component. In this case, I have
> successfully calculated the L2 and H1 norms of the error using the
> `integrate_difference` and `compute_global_error` functions. The calculations
> proceeded smoothly, and the execution speed was reasonable.
>
> However, the second scenario involves using the FeFieldFunction when the
> computed solution consists of two components. My objective is to calculate the
> error for the first component only. To accomplish this, I employed the new
> VectorFunctionFromScalarFunctionObject class locally. Unfortunately, the
> computation of the L2 norm in this case was considerably slow and the H1 norm
> was very slow.
>
> Given the current situation, I have two questions that I hope you can assist
> me with, as you have done in the past:
>
> 1. Regarding my specific project, is there an alternative method, other than
> the FeFieldFunction, to save the deserialized solution as a function and
> utilize it as the exact solution? It is crucial that the solution be saved in
> the first component of a vector with two components.
>
> 2. In a more general sense, how can I test this function to determine if the
> sluggishness is solely attributable to the FeFieldFunction? Additionally, I
> would like to verify whether the function is functioning correctly or not, to
> be able to send it in a pull request to the dealii git repository.

Najwa:
I don't know why the run time should be substantially slower for a system of
two components than for a single scalar problem. That might be worth
investigating. If you wanted to look into this, it would be useful to start
with as small a program that shows the issue.

In general, FEFieldFunction is slow because it doesn't know where the function
should be evaluated. If you give it a point, it first has to find the cell the
point is in, and then compute the reference coordinates of that point, etc.
You probably do this many times if you want to evaluate the solution on the
quadrature points of a whole triangulation. This will never be fast

There are perhaps better ways to deal with the situation. For example, instead
of doing the integration on the coarse mesh, you could take the coarse
solution, interpolate it onto the fine mesh, and then do the integration over
the fine mesh -- because there you can then just compute the values of both
the coarse and fine solution at the quadrature points of the fine mesh using
FEValues. This is fast, and the interpolation to the fine mesh should also be
reasonably fast.

There are also faster but more specialized alternatives to FEFieldFunction.
FEPointEvaluation is one option, though I don't know that class well enough.