Coefficient field (material properties) reconstruction through MFEM- Wave Equation

[aaelmeli]

Hi

I am working on two types of the inverse problem (reconstruction of the material properties) that involves solving the wave equation in the frequency domain for multiple right-hand sides. This problem is well known in geophysics (seismic inversion) and in biomedical imaging (Ultrasound imaging). Typically, this problem involves thousands of forward solves (FE simulation) as a result of sweeping overall frequency content and all sources (right-hand sides) at each iteration to the optimization (reconstruction).

Although I am just starting learning MFEM and trying to switch to this new mindset, I can imagine how powerful MFEM is especially for multiphysics simulation. However, I am still have some issues in deciding whether MFEM is a suitable choice for my problems.

To be specific, the optimum framework that I am looking for is preferably having the following:

1- Can handle complex arithmetics (needed for solving the wave equation in the frequency domain)
2- Ability to modify the integration rule for specific region/regions of the domain (e.g. mid-point integration), this is needed to model a variant of the perfectly matched layer (PMDL) that acts as wave absorber.
3- Ability to update material properties after each optimization iteration in runtime without the need for updating the input.i file after each iteration.

So, if you could please give some pointers on the previous 3 points and some pointer on how to model an arbitrary, completely heterogeneous 3D domain (i.e each finite element in the mesh has different material property).

I really appreciate your response.

[psocratis]

Hi @aaelmeli,

Here are some answers to your questions.

1. Yes, MFEM can handle complex formulations. It does not use complex arithmetic but rather forms the 2x2 system for real and imaginary parts. Fortunately, this is hidden from the user., i.e, the user provides the 2 bilinear (linear) forms (1 for the real and 1 for the imaginary part of the differential operator) and everything else is done under the hood. I suggest you take a look at ex22 and ex25. In fact ex25 is using a PML formulation for the time-harmonic Maxwell equations.

2. Yes but it might need a minor workaround. I think this is not possible out of the box. However, one possible way I can think of how to do this is to use element attributes. More precisely, you can set the integration rule for a specific integrator and you can "restrict" an integrator to be valid in a specific region of the domain, using element attributes.

3. Yes, I believe you can do this using "time-dependent" coefficients.

You can define the heterogeneous domain using a Coefficient. For example, a PWConstCoefficient can be used to assign a different constant value to each attribute (and you can assign a different attribute to each element). You can also build more complicated (non-constant) PWcoefficients in a similar way as the ones found in coefficient.[ch]pp

I hope this helps

Best,

Socratis

[rcarson3]

@aaelmeli it sounds like for 2 you might be talking about wanting something like a B-bar formulation that's commonly used in solid mechanics type problems.

You could easily extend the NonlinearFormIntegrator, BilinearFormIntegrator, or LinearFormIntegrator classes to do this. I for example have implemented this in my quasi-statics solid mechanics code, ExaConstit, here: https://github.com/LLNL/ExaConstit/blob/exaconstit-dev/src/mechanics_integrators.cpp#L991-L1154 . So, if you want something like that it should just require a few extra lines of code.