Code fails when mesh adaptation is incorperated into a code formulated with the 'varf' keyword

Hi all, I’m currently trying to solve a toy problem using the varf keyword. So far, I have been able to solve the problem without mesh adaptation using both the problem keyword and varf keyword, and I have also been able to solve the problem with mesh adaptation using the problem keyword. However, issues arise when I try to introduce the same mesh adaptation to the code formulated using the varf keyword; the exact error message that I recieve is dependent upon which solver/preconditioner I am using. If I use UMFPACK64 with sparsesolver, I recieve no error message at all and the code simply terminates itself upon reaching the first mesh adaptation. Alternatively, if I use PETSc with "-pc_type gamg" , the code teminates in the same place with the following error code in the terminal:

[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see https://petsc.org/release/faq/#valgrind and https://petsc.org/release/faq/
[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run
[0]PETSC ERROR: to get more information on the crash.
[0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.

I don’t understand how this could be a memory related issue for such a simple PDE system. Can anyone help me with this problem? Perhaps I’ve just written my code incorrectly. Many thanks!

My code formulated using the problem keyword (with mesh adaptation):
SimplerProblemLinearCoupled_problem.edp (2.2 KB)

My code formulated using the varf keyword (without mesh adaptation):
SimplerProblemLinearCoupled_varf.edp (2.6 KB)

My code formulated using the varf keyword (with mesh adaptation):
SimplerProblemLinearCoupled_varf_adaptive.edp (2.9 KB)

If you don’t know how to do mesh adaptation with PETSc, you can have a look at this example: FreeFem-sources/examples/hpddm/laplace-adapt-3d-PETSc.edp at develop · FreeFem/FreeFem-sources · GitHub.

My guess is that you are not properly resizing the Mat after mesh adaptation, thus the error. See the above example for how to fix this.

Hi, thanks for the reply. The example you’ve provided looks quite different to my own code, I’m not particularly famililar with the methods presented but I will take a look. Can you see any obvious mistakes in my code or is it just not possible to use the adaptmesh( ) function in a code formulated with the varf keyword? Also, I use matrix, not Mat; is this a problem?

One additional point, for some reason it occasionally gives the following error: Assertion fail : (mu>0 && y.N()%mu==0). Not sure if this is a clue at all.

Sorry, I thought you had a question specific to PETSc, but it appears you have the error even in sequential. Have you run the code with debugging turned on?

The immediate problem is a size problem.
If you don’t set the size of aux and just define it on first assignment
the code will run, see the changes here, but seg faults later.
IIRC I had some issue earlier running out of memory using regions or
something but there are probably quirks in the mesh adaptation
that may not work indefinitely.

SimplerProblemLinearCoupled_varf_adaptive.edp (3.1 KB)

Hi prj, yes exactly, it doesn’t even work in sequential unforrtunately. I’m using VScode with a FreeFem extension, which sadly does not feature a debugger. Do you know of a way to debug FreeFem code?

Hi Mike, thanks for the reply! Your adjustment almost does the trick, however it appears that it is also necessary to resize sol and solold. When I print out the sizes of sol and solold, they each remain at 6200 while the aux, b and Xh all change with the mesh adaptation. Any ideas of how to resize sol and solold? It doesn’t seem to work when I try to resize them in the same in which aux, b and Xh are resized.

Just sol = 0; solold = 0; will resize both functions.

IIRC when I ran into a problem it was because adapting a mesh with regions
had a memory leak but i just used indicator variable instead.
That probably is not relevant here.

I hesitated to go further with size issues since I couldn’t remember what does what …
I think there may have also been an issue with number of dofs and array asize
depending on the problem. Some values of tgv are supposed to eliminate
fixed equations so I would imagine the A and rhs shrink.

Hi Mike, you might’ve given me enough info to work with. I’ll update this thread if I get anywhere. Many thanks to you and prj for your input! And if anything occurs to you, do let me know.

Hi, I seem to have sorted the issue; now the code works with mesh adaptation while using the varf keyword. Here’s the code that seems to work:
SimplerProblemLinearCoupled_varf_adaptive_fixed.edp (3.1 KB)

However, I’ve unearthed an unexpected problem and that is the question of how to resize the value of the unknowns at the previous timestep (uu1 and uu2) to the newly adapted mesh. I’ve never given this much thought and always assumed that the proper technique was to simply include the following lines of code just after calling the adaptmesh( ) function:

  uu1=uu1;     
  uu2=uu2;

But, what I’ve noticed is that the inclusion of this line of code causes signifcant changes to uu1 and uu2 in subsequent mesh adaptations during a given timestep, which doesn’t seem right. Also, the solution obtained when these lines of codes are included is quite different to the solution obtained when they are absent. My question simply is, what is the correct procedure here? I’ve attached the code formulated with the problem keyword below. On lines 111-120 I’ve inlcuded several different possible lines of code (options A to D), one of which may be the correct procedure. Maybe someone can highlight which option is correct or maybe someone can suggest a different option entirely?
SimplerProblemLinearCoupled_problem_adaptive.edp (2.9 KB)

Many thanks for the help.

I remember originally I tried to do all of that manually with interpolate but it looks
like assigning on fespace to the other does all of that. if you add the brackets I think
that refers to the underlying array which does not know about the mesh.
The semantics take a little getting used to. I’ve tried to write some small specialized “almost” languages
and I understand from the author’s POV how these things happen but it can be confusing
to learn as a user lol.

Hi Mike, yes it’s certainly a bit confusing! So, which of the methods (options A-D) do you think is most appropriate? Are you suggesting that I use option A?

Hi Mike, did you have any thoughts on my question?