Example code not working

General Discussion
1

Dear @prj,

I would like to use the ‘Augmented Lagrangian Preconditioner for Large-Scale Hydrodynamic Stability Analysis’ method based on your codes. However, when I try to run the codes, the base-flow calculation runs fine, but the eigenvalue computation does not work. It stops when the method ‘EPSSolve’ is called, with the following message:

[0]PETSC ERROR: #1 PetscCommDuplicate() line 118 in /builds/workspace/FreeFEM-sources-deployDEB-withPETSc/3rdparty/ff-petsc/petsc-3.13.0/src/sys/objects/tagm.c
[0]PETSC ERROR: #2 PetscHeaderCreate_Private() line 63 in /builds/workspace/FreeFEM-sources-deployDEB-withPETSc/3rdparty/ff-petsc/petsc-3.13.0/src/sys/objects/inherit.c
[0]PETSC ERROR: #3 MatCreate() line 83 in /builds/workspace/FreeFEM-sources-deployDEB-withPETSc/3rdparty/ff-petsc/petsc-3.13.0/src/mat/utils/gcreate.c
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: [1]PETSC ERROR: #1 PetscCommDuplicate() line 118 in /builds/workspace/FreeFEM-sources-deployDEB-withPETSc/3rdparty/ff-petsc/petsc-3.13.0/src/sys/objects/tagm.c
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://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind
[0]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS X to find memory corruption errors
[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: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Signal received
[0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.13.0, Mar 29, 2020
[0]PETSC ERROR: FreeFem+±mpi on a named DESKTOP-8VC8OSF by root Sat Sep 12 11:17:02 2020
[0]PETSC ERROR: Configure options MAKEFLAGS= --prefix=/usr/local/ff-petsc//c --with-debugging=0 COPTFLAGS="-O3 -mtune=generic" CXXOPTFLAGS="-O3 -mtune=generic" FOPTFLAGS="-O3 -mtune=generic" --with-cxx-dialect=C++11 --with-ssl=0 --with-x=0 --with-fortran-bindings=0 --with-cc=/usr/bin/mpicc --with-cxx=/usr/bin/mpic++ --with-fc=/usr/bin/mpif90 --with-scalar-type=complex --with-blaslapack-include= --with-blaslapack-lib="-llapack -lblas" --with-scalapack-dir=/usr/local/ff-petsc//r --with-metis-dir=/usr/local/ff-petsc//r --with-ptscotch-dir=/usr/local/ff-petsc//r --with-suitesparse-include=/usr/local/ff-petsc/r/include --with-suitesparse-lib="-Wl,-rpath,/usr/local/ff-petsc/r/lib -L/usr/local/ff-petsc/r/lib -lumfpack -lklu -lcholmod -lbtf -lccolamd -lcolamd -lcamd -lamd -lsuitesparseconfig" --with-mumps-dir=/usr/local/ff-petsc//r --with-parmetis-dir=/usr/local/ff-petsc//r --with-tetgen-dir=/usr/local/ff-petsc//r --download-slepc --download-hpddm --download-hpddm-commit=e8639ff PETSC_ARCH=fc
[0]PETSC ERROR: #4 User provided function() line 0 in unknown file
[DESKTOP-8VC8OSF:00232] 3 more processes have sent help message help-mpi-api.txt / mpi-abort

Do you have any idea what causes this issue? I am running the code in the Windows Ubuntu app version Ubuntu 18.04.4 LTS, FreeFem version 4.6.

Thank you for your help in advance
András

2

Could you please try with the latest FreeFEM version (4.7)?
Also, in Eigensolver.edp, could you try to replace -st_pc_type fieldsplit by -st_pc_type lu, and tell me if it is still not working?

3

Dear @prj,

thanks for the fast reply. I managed to resolve the issue: the macro_ddm.idp file was outdated (a previous reinstallation of FreeFem++ did not set the path properly).

However, I have and other question related to the variational formulation. In this example, in the variational formulation of both the residual and the jacobian the boundary conditions are the same. However, in the augmented Lagrangian preconditioner example, the boundary condition in the jacobian and in the residual are different: in the residual, the nonzero velocities are prescribed, but in the jacobian every boundary condition is zero. Why is that?

Thank you for your help
András

4

Good to know the code is now working properly. The value you set with the on keyword on a bilinear form does not matter, it’s only taken into account for the assembly of the linear form. So, for the Jacobian, on(1, u = 1) or on(1, u = 0) does not matter. However, for the residual, it does matter. Here is an example so that you can see it for yourself.

mesh Th = square(4, 4);
varf vPb1(u, v) = on(1,2,3,4, u = 1);
varf vPb2(u, v) = on(1,2,3,4, u = 1234);
fespace Vh(Th, P1);
matrix A1 = vPb1(Vh, Vh, tgv = -1);
matrix A2 = vPb2(Vh, Vh, tgv = -1);
matrix diff = A1 + (-1.0)*A2;
cout << diff << endl;

Effectively, this means that the 2D and 3D code are computing the same quantities, sorry for the confusion.

5

Thank you for the quick answer, now the boundary condition specification is clear!

6

I have two more questions about the method.

  1. It is customary in flow stability that the base-flow calculation is carried out in multiple steps: first, at a low Reynolds number, the initial guess in the solution of the nonlinear equation is ad-hoc, and then the Reynolds number is increased in multiple steps. However, then I tried to use this procedure, the nonlinear equation solution failed to converge, but the solution converged with a worse initial guess, e.g. [u1, u2, p] = [1, 0, 0]. I also noticed that in your paper, you only took one previous step in the base flow calculation (Re=50 first, then Re=100). My question is basically the following: is this normal, and when using iterative solvers, a too good initial guess can cause convergence problems?
  2. In the eigenvalue solution file , at line 74 the operator s^-1 Lp is assembled, and a boundary condition is specified. This is unclear to me: how should I specify this boundary condition for a general problem?

Thanks for the help

András

7
  1. Yes, this is customary to all nonlinear problems. You can do that, but it’s not implemented (I think? I.e., doing -Re 50, then -Re 75, and finally -Re 100, will use the baseflow for Re 50 for the Re 100 solve, but I’m not sure)
  2. There is no boundary condition specified in Lp