PETSc error when using 'adapt()' and '.setorder' is greater than one

[caseyjamesdavis]

Hi Alex,

I have been learning how to use the adaptive mesh refinement (AMR) and I ran into this while running the amr-1d2d3d example.

The problem appears to be isolated to the 2D portion of the example.

If I change the following lines:

v.setorder(s1, 5);
v.setorder(s2, 4);
v.setorder(s3, 3);
v.setorder(s4, 3);

to:

v.setorder(s1, 1);
v.setorder(s2, 1);
v.setorder(s3, 1);
v.setorder(s4, 1);

the error goes away and the example runs as expected.

I noticed this first in another simulation I'm working on where the temperature field was set to order 2 and adapt() would give the same error. If I reduced the order of the temperature field to 1 the simulation would work as expected.

Any ideas?

Casey

ps - i did a full re-install this morning. cjd.

[halbux]

Hi Casey,

Another user has had similar issues on a windows compilled petsc. I call a matrix permutation operation of petsc when the order is >1 which seems to be here again the culprit. (when order = 1 operations are straightforward and you will not get that error, that's normal).

What's your OS? Did you change the configuration settings of petsc? Did you ask petsc to install MUMPS?

[halbux]

BTW: I see you are using libopenblas from /usr/lib...: did you install openblas from the repo? I strongly recommend not, it can be extremely slow. Best is to let petsc download it and compile it

[caseyjamesdavis]

I'm using Manjaro 20.2.1 and a 5.9 kernel.

The only thing I did was run the 'install_petsc.sh' script. I did not change any settings.

[halbux]

Ok. It works on my fedora 30 though (it's part of the examples I run to validate the code every day). I'll have a closer look asap.

[halbux]

Could you please confirm that the error happens during the call

A.permute(renumtodiagblocks, renumtodiagblocks);

in src/field/rawfield.cpp (or at least another permute call there).

If yes can you also confirm the call happens during the call to

MatPermute

in src/formulation/mat.cpp

Thanks

[caseyjamesdavis]

Happy to help.

Not exactly sure how to check. Do I run the main.cpp as usual with some kind of debug flag?

[halbux]

You could edit the src code to print out a message before the call and after, recompile then run.

Alternatively you could assume that the error comes from there, in which case you could try to add before the MatPermute call

IsSetPermutation(isobject) for the 2 IS objects.

Might also be needed in src/formulation/vec.cpp in the similar permute function.

This might solve your issue.

[caseyjamesdavis]

Thanks. Got it. I'll add the print statements to your code and let you know what I find.

[halbux]

I just checked the petsc main branch and it seems 22 hours ago they changed the behavior of MatPermute to throw the erro you see, so I will have that error with a fresh enough petsc version. I will fix it now , or tomorrow latest

[caseyjamesdavis]

I did get the same error with a petsc/sparselizard install date of 2021/02/04.

[halbux]

should be fixed, could you git pull and try again? Output should be 1.33227e-15 3.26364e-11 1.35181e-12 with a 1 returned

[caseyjamesdavis]

I added the two print statements to src/field/rawfield.cpp:

And added the two print statements to src/formulation/mat.cpp

Results:

[caseyjamesdavis]

No problem. I reinstall and try again.

[halbux]

hm, just git pull, you can keep the petsc as it is :)

[caseyjamesdavis]

thanks. I was wondering if that would work.

[caseyjamesdavis]

Nice work - no error ;-) See results below.

Is it normal for it to take about 15 seconds to transition from the 2D to 3D example?

[halbux]

Yes it's not surprising it takes at least a couple of seconds, the mesh adaptation is the visible part of the work that is done, but behind the scene the fields must be automatically reinterpolated on the new mesh (here an order 4 curvature mesh with order 5 interpolation fields, that can be a little heavy since it chooses order 2*5+2 = 12 when it needs to integrate for the projection on the new mesh). You can speed it up with a decreased field.setupdateaccuracy , but make sure the projection is still good, or better field.noautomaticupdate if you set the field value again yourself on the new mesh anyways without using the value as it was on the previous mesh

[caseyjamesdavis]

That makes sense. I'll explore the two options you mention.

Thanks for digging into this so swiftly.

Cheers, Casey

[halbux]

.noautomaticupdate is generally what you want to use for linear static problems. For transient resolutions you cannot use that since the field value at the previous time step (and therefore probably previous mesh) is needed.

Alex