Ok, it seems to work with that switch, however, when I try to use DM sequence,
I get errors like:
0 SNES Function norm 5.067205249874e-03
1 SNES Function norm 7.983917252341e-08
2 SNES Function norm 7.291012540201e-11
0 SNES Function norm 2.228951406196e+02
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Argument out of range
[0]PETSC ERROR: New nonzero at (0,3) caused a malloc
Use MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to turn off
this check
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for
trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.5.3, unknown
[0]PETSC ERROR: ./blowup_batch2 on a arch-macports named gs_air by gideon Thu
Aug 27 23:53:19 2015
[0]PETSC ERROR: Configure options --prefix=/opt/local
--prefix=/opt/local/lib/petsc --with-valgrind=0 --with-shared-libraries
--with-debugging=0 --with-c2html-dir=/opt/local --with-x=0
--with-blas-lapack-lib=/System/Library/Frameworks/Accelerate.framework/Versions/Current/Accelerate
--with-hwloc-dir=/opt/local --with-suitesparse-dir=/opt/local
--with-superlu-dir=/opt/local --with-metis-dir=/opt/local
--with-parmetis-dir=/opt/local --with-scalapack-dir=/opt/local
--with-mumps-dir=/opt/local --with-superlu_dist-dir=/opt/local
CC=/opt/local/bin/mpicc-mpich-mp CXX=/opt/local/bin/mpicxx-mpich-mp
FC=/opt/local/bin/mpif90-mpich-mp F77=/opt/local/bin/mpif90-mpich-mp
F90=/opt/local/bin/mpif90-mpich-mp COPTFLAGS=-Os CXXOPTFLAGS=-Os FOPTFLAGS=-Os
LDFLAGS="-L/opt/local/lib -Wl,-headerpad_max_install_names"
CPPFLAGS=-I/opt/local/include CFLAGS="-Os -arch x86_64" CXXFLAGS=-Os FFLAGS=-Os
FCFLAGS=-Os F90FLAGS=-Os PETSC_ARCH=arch-macports
--with-mpiexec=mpiexec-mpich-mp
[0]PETSC ERROR: #1 MatSetValues_SeqAIJ() line 487 in
/opt/local/var/macports/build/_opt_local_var_macports_sources_rsync.macports.org_release_tarballs_ports_math_petsc/petsc/work/v3.5.3/src/mat/impls/aij/seq/aij.c
[0]PETSC ERROR: #2 MatSetValues() line 1135 in
/opt/local/var/macports/build/_opt_local_var_macports_sources_rsync.macports.org_release_tarballs_ports_math_petsc/petsc/work/v3.5.3/src/mat/interface/matrix.c
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Argument out of range
[0]PETSC ERROR: New nonzero at (1,4) caused a malloc
Use MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to turn off
this check
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for
trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.5.3, unknown
[0]PETSC ERROR: ./blowup_batch2 on a arch-macports named gs_air by gideon Thu
Aug 27 23:53:19 2015
[0]PETSC ERROR: Configure options --prefix=/opt/local
--prefix=/opt/local/lib/petsc --with-valgrind=0 --with-shared-libraries
--with-debugging=0 --with-c2html-dir=/opt/local --with-x=0
--with-blas-lapack-lib=/System/Library/Frameworks/Accelerate.framework/Versions/Current/Accelerate
--with-hwloc-dir=/opt/local --with-suitesparse-dir=/opt/local
--with-superlu-dir=/opt/local --with-metis-dir=/opt/local
--with-parmetis-dir=/opt/local --with-scalapack-dir=/opt/local
--with-mumps-dir=/opt/local --with-superlu_dist-dir=/opt/local
CC=/opt/local/bin/mpicc-mpich-mp CXX=/opt/local/bin/mpicxx-mpich-mp
FC=/opt/local/bin/mpif90-mpich-mp F77=/opt/local/bin/mpif90-mpich-mp
F90=/opt/local/bin/mpif90-mpich-mp COPTFLAGS=-Os CXXOPTFLAGS=-Os FOPTFLAGS=-Os
LDFLAGS="-L/opt/local/lib -Wl,-headerpad_max_install_names"
CPPFLAGS=-I/opt/local/include CFLAGS="-Os -arch x86_64" CXXFLAGS=-Os FFLAGS=-Os
FCFLAGS=-Os F90FLAGS=-Os PETSC_ARCH=arch-macports
--with-mpiexec=mpiexec-mpich-mp
[0]PETSC ERROR: #3 MatSetValues_SeqAIJ() line 487 in
/opt/local/var/macports/build/_opt_local_var_macports_sources_rsync.macports.org_release_tarballs_ports_math_petsc/petsc/work/v3.5.3/src/mat/impls/aij/seq/aij.c
[0]PETSC ERROR: #4 MatSetValues() line 1135 in
/opt/local/var/macports/build/_opt_local_var_macports_sources_rsync.macports.org_release_tarballs_ports_math_petsc/petsc/work/v3.5.3/src/mat/interface/matrix.c
And in my main program, I did set
MatSetOption(J, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
for the Jacobian.
-gideon
> On Aug 27, 2015, at 11:23 PM, Barry Smith <[email protected]> wrote:
>
>
>> On Aug 27, 2015, at 10:15 PM, Gideon Simpson <[email protected]>
>> wrote:
>>
>> That’s correct, I am not using the SNESGetDM. I suppose I could. Keep in
>> mind that I’m trying to solve, simultaneously, for the scalar parameters and
>> the vector field. I guess what I am unclear about is how DMRefine is to
>> know that the unknown associated with the scalar parameters can never be
>> coarsened out, but must be retained at all iterations.
>
> Nothing ever gets coarsened in grid sequencing, it only gets refined.
>
> The reason it knows not to "refine" the scalars is because the scalars are
> created with DMRedundant and the DMRedundant object knows that refinement
> means "leave as is, since there is no grid" while the DMDA knows it is a grid
> and knows how to refine itself. So when it "interpolates" the DMRedundant
> variables it just copies them (or multiples them by the matrix 1 which is
> just a copy).
>>
>> Here is my form function. I can send more code if needed.
>>
>
> Just change the user->packer that you use to be the DM obtained with
> SNESGetDM()
>
>
>> /* Form the system of equations for computing a blowup solution*/
>> PetscErrorCode form_function(SNES snes, Vec U, Vec F, void *ctx){
>>
>> blowup_ctx *user = (blowup_ctx *) ctx;
>> PetscInt i;
>> PetscScalar dx, dx2, xmax,x;
>> PetscScalar u, v, f,g, ux, vx, uxx, vxx, fx,gx, fxx, gxx;
>> DMDALocalInfo info;
>> Vec p_vec, Q_vec, Fp_vec, FQ_vec;
>> PetscScalar *p_array, *Fp_array;
>> Q *Qvals, *FQvals;
>> PetscScalar Q2sig, W2sig;
>> PetscScalar a,a2, b, u0, sigma;
>>
>> dx = user->dx; dx2 = dx *dx;
>> xmax = user->xmax;
>> sigma = user->sigma;
>>
>> /* PetscPrintf(PETSC_COMM_SELF, " dx = %g, sigma = %g\n", dx, sigma); */
>>
>> /* Extract raw arrays */
>> DMCompositeGetLocalVectors(user->packer, &p_vec, &Q_vec);
>> DMCompositeGetLocalVectors(user->packer, &Fp_vec, &FQ_vec);
>>
>> DMCompositeScatter(user->packer, U, p_vec, Q_vec);
>> /* VecView(Q_vec, PETSC_VIEWER_STDOUT_SELF); */
>>
>> VecGetArray(p_vec,&p_array);
>> VecGetArray(Fp_vec,&Fp_array);
>>
>> DMDAVecGetArray(user->Q_dm, Q_vec, &Qvals);
>> DMDAVecGetArray(user->Q_dm, FQ_vec, &FQvals);
>>
>> DMDAGetLocalInfo(user->Q_dm, &info);
>>
>> a = p_array[0]; a2 = a*a;
>> b = p_array[1];
>> u0 = p_array[2];
>>
>> /* Set boundary conditions at the origin*/
>> if(info.xs ==0){
>> set_origin_bcs(u0, Qvals);
>> }
>> /* Set boundray conditions in the far field */
>> if(info.xs+ info.xm == info.mx){
>> set_farfield_bcs(xmax,sigma, a, b, dx, Qvals,info.mx);
>> }
>>
>> /* Solve auxiliary equations */
>> if(info.xs ==0){
>> uxx = (2 * Qvals[0].u-2 * u0)/dx2;
>> vxx = (Qvals[0].v + Qvals[0].g)/dx2;
>> vx = (Qvals[0].v - Qvals[0].g)/(2*dx);
>> Fp_array[0] = Qvals[0].u - Qvals[0].f;
>> Fp_array[1] = -vxx - (1/a) * (.5/sigma) * u0;
>> Fp_array[2] = -uxx + (1/a2) * u0
>> + (1/a) * (-b * vx + PetscPowScalar(u0 * u0, sigma) * vx);
>> }
>>
>> /* Solve equations in the bulk */
>> for(i=info.xs; i < info.xs + info.xm;i++){
>>
>> u = Qvals[i].u;
>> v = Qvals[i].v;
>> f = Qvals[i].f;
>> g = Qvals[i].g;
>>
>> x = (i+1) * dx;
>>
>> Q2sig = PetscPowScalar(u*u + v*v,sigma);
>> W2sig= PetscPowScalar(f*f + g*g, sigma);
>>
>> ux = (Qvals[i+1].u-Qvals[i-1].u)/(2*dx);
>> vx = (Qvals[i+1].v-Qvals[i-1].v)/(2*dx);
>> fx = (Qvals[i+1].f-Qvals[i-1].f)/(2*dx);
>> gx = (Qvals[i+1].g-Qvals[i-1].g)/(2*dx);
>>
>> uxx = (Qvals[i+1].u+Qvals[i-1].u - 2 *u)/(dx2);
>> vxx = (Qvals[i+1].v+Qvals[i-1].v- 2 *v)/(dx2);
>> fxx = (Qvals[i+1].f+Qvals[i-1].f -2*f)/(dx2);
>> gxx = (Qvals[i+1].g+Qvals[i-1].g -2*g)/(dx2);
>>
>> FQvals[i].u = -uxx +1/a2 * u
>> + 1/a *(.5/sigma* v +x * vx- b* vx +Q2sig* vx);
>>
>> FQvals[i].v = -vxx +1/a2 * v
>> - 1/a *(.5/sigma * u +x * ux- b* ux +Q2sig* ux);
>>
>> FQvals[i].f = -fxx +1/a2 * f
>> + 1/a *(.5/sigma * g +x * gx+ b* gx -W2sig* gx);
>>
>> FQvals[i].g =-gxx +1/a2 * g
>> - 1/a *(.5/sigma * f +x * fx+ b* fx -W2sig* fx);
>> }
>>
>> /* Restore raw arrays */
>> VecRestoreArray(p_vec, &p_array);
>> VecRestoreArray(Fp_vec, &Fp_array);
>>
>> DMDAVecRestoreArray(user->Q_dm, Q_vec, &Qvals);
>> DMDAVecRestoreArray(user->Q_dm, FQ_vec, &FQvals);
>>
>> DMCompositeGather(user->packer,F,INSERT_VALUES, Fp_vec, FQ_vec);
>> DMCompositeRestoreLocalVectors(user->packer, &p_vec, &Q_vec);
>> DMCompositeRestoreLocalVectors(user->packer, &Fp_vec, &FQ_vec);
>>
>> return 0;
>> }
>>
>>
>> Here is the form function:
>>
>>
>>
>> -gideon
>>
>>> On Aug 27, 2015, at 11:09 PM, Barry Smith <[email protected]> wrote:
>>>
>>>
>>> Can you send the code, that will be the easiest way to find the problem.
>>>
>>> My guess is that you have hardwired in your function/Jacobian computation
>>> the use of the original DM for computations instead of using the current DM
>>> (with refinement there will be a new DM on the second level different than
>>> your original DM). So what you need to do in writing your FormFunction and
>>> FormJacobian is to call SNESGetDM() to get the current DM and then use
>>> DMComputeGet... to access the individual DMDA and DMRedundent for the
>>> parts. I notice you have this user.Q_dm I bet inside your form functions
>>> you use this DM? You have to remove this logic.
>>>
>>> Barry
>>>
>>>> On Aug 27, 2015, at 9:42 PM, Gideon Simpson <[email protected]>
>>>> wrote:
>>>>
>>>> I have it set up as:
>>>>
>>>> DMCompositeCreate(PETSC_COMM_WORLD, &user.packer);
>>>> DMRedundantCreate(PETSC_COMM_WORLD, 0, 3, &user.p_dm);
>>>> DMCompositeAddDM(user.packer,user.p_dm);
>>>> DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_GHOSTED,
>>>> nx, 4, 1, NULL, &user.Q_dm);
>>>> DMCompositeAddDM(user.packer,user.Q_dm);
>>>> DMCreateGlobalVector(user.packer,&U);
>>>>
>>>> where the user.packer structure has
>>>>
>>>> DM packer;
>>>> DM p_dm, Q_dm;
>>>>
>>>> Q_dm holds the field variables and p_dm holds the scalar values (the
>>>> nonlinear eigenvalues).
>>>>
>>>> Here are some of the errors that are generated:
>>>>
>>>> [0]PETSC ERROR: --------------------- Error Message
>>>> --------------------------------------------------------------
>>>> [0]PETSC ERROR: Argument out of range
>>>> [0]PETSC ERROR: New nonzero at (0,3) caused a malloc
>>>> Use MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to turn
>>>> off this check
>>>> [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
>>>> for trouble shooting.
>>>> [0]PETSC ERROR: Petsc Release Version 3.5.3, unknown
>>>> [0]PETSC ERROR: ./blowup_batch2 on a arch-macports named gs_air by gideon
>>>> Thu Aug 27 22:40:54 2015
>>>> [0]PETSC ERROR: Configure options --prefix=/opt/local
>>>> --prefix=/opt/local/lib/petsc --with-valgrind=0 --with-shared-libraries
>>>> --with-debugging=0 --with-c2html-dir=/opt/local --with-x=0
>>>> --with-blas-lapack-lib=/System/Library/Frameworks/Accelerate.framework/Versions/Current/Accelerate
>>>> --with-hwloc-dir=/opt/local --with-suitesparse-dir=/opt/local
>>>> --with-superlu-dir=/opt/local --with-metis-dir=/opt/local
>>>> --with-parmetis-dir=/opt/local --with-scalapack-dir=/opt/local
>>>> --with-mumps-dir=/opt/local --with-superlu_dist-dir=/opt/local
>>>> CC=/opt/local/bin/mpicc-mpich-mp CXX=/opt/local/bin/mpicxx-mpich-mp
>>>> FC=/opt/local/bin/mpif90-mpich-mp F77=/opt/local/bin/mpif90-mpich-mp
>>>> F90=/opt/local/bin/mpif90-mpich-mp COPTFLAGS=-Os CXXOPTFLAGS=-Os
>>>> FOPTFLAGS=-Os LDFLAGS="-L/opt/local/lib -Wl,-headerpad_max_install_names"
>>>> CPPFLAGS=-I/opt/local/include CFLAGS="-Os -arch x86_64" CXXFLAGS=-Os
>>>> FFLAGS=-Os FCFLAGS=-Os F90FLAGS=-Os PETSC_ARCH=arch-macports
>>>> --with-mpiexec=mpiexec-mpich-mp
>>>> [0]PETSC ERROR: #1 MatSetValues_MPIAIJ() line 530 in
>>>> /opt/local/var/macports/build/_opt_local_var_macports_sources_rsync.macports.org_release_tarballs_ports_math_petsc/petsc/work/v3.5.3/src/mat/impls/aij/mpi/mpiaij.c
>>>> [0]PETSC ERROR: [1]PETSC ERROR: --------------------- Error Message
>>>> --------------------------------------------------------------
>>>> [1]PETSC ERROR: Argument out of range
>>>> [1]PETSC ERROR: Inserting a new nonzero (40003, 0) into matrix
>>>> [1]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html
>>>> for trouble shooting.
>>>> [1]PETSC ERROR: Petsc Release Version 3.5.3, unknown
>>>> [1]PETSC ERROR: ./blowup_batch2 on a arch-macports named gs_air by gideon
>>>> Thu Aug 27 22:40:54 2015
>>>> [1]PETSC ERROR: Configure options --prefix=/opt/local
>>>> --prefix=/opt/local/lib/petsc --with-valgrind=0 --with-shared-libraries
>>>> --with-debugging=0 --with-c2html-dir=/opt/local --with-x=0
>>>> --with-blas-lapack-lib=/System/Library/Frameworks/Accelerate.framework/Versions/Current/Accelerate
>>>> --with-hwloc-dir=/opt/local --with-suitesparse-dir=/opt/local
>>>> --with-superlu-dir=/opt/local --with-metis-dir=/opt/local
>>>> --with-parmetis-dir=/opt/local --with-scalapack-dir=/opt/local
>>>> --with-mumps-dir=/opt/local --with-superlu_dist-dir=/opt/local
>>>> CC=/opt/local/bin/mpicc-mpich-mp CXX=/opt/local/bin/mpicxx-mpich-mp
>>>> FC=/opt/local/bin/mpif90-mpich-mp F77=/opt/local/bin/mpif90-mpich-mp
>>>> F90=/opt/local/bin/mpif90-mpich-mp COPTFLAGS=-Os CXXOPTFLAGS=-Os
>>>> FOPTFLAGS=-Os LDFLAGS="-L/opt/local/lib -Wl,-headerpad_max_install_names"
>>>> CPPFLAGS=-I/opt/local/include CFLAGS="-Os -arch x86_64" CXXFLAGS=-Os
>>>> FFLAGS=-Os FCFLAGS=-Os F90FLAGS=-Os PETSC_ARCH=arch-macports
>>>> --with-mpiexec=mpiexec-mpich-mp
>>>> [1]PETSC ERROR: #1 MatSetValues_MPIAIJ() line 561 in
>>>> /opt/local/var/macports/build/_opt_local_var_macports_sources_rsync.macports.org_release_tarballs_ports_math_petsc/petsc/work/v3.5.3/src/mat/impls/aij/mpi/mpiaij.c
>>>> [1]PETSC ERROR: #2 MatSetValues() line 1135 in
>>>> /opt/local/var/macports/build/_opt_local_var_macports_sources_rsync.macports.org_release_tarballs_ports_math_petsc/petsc/work/v3.5.3/src/mat/interface/matrix.c
>>>>
>>>>
>>>>
>>>> -gideon
>>>>
>>>>> On Aug 27, 2015, at 10:37 PM, Barry Smith <[email protected]> wrote:
>>>>>
>>>>>
>>>>> We need the full error message.
>>>>>
>>>>> But are you using a DMDA for the scalars? You should not be, you should
>>>>> be using a DMRedundant for the scalars.
>>>>>
>>>>> Barry
>>>>>
>>>>> Though you should not get this error even if you are using a DMDA there.
>>>>>
>>>>>> On Aug 27, 2015, at 9:32 PM, Gideon Simpson <[email protected]>
>>>>>> wrote:
>>>>>>
>>>>>> I’m getting the following errors:
>>>>>>
>>>>>> [1]PETSC ERROR: Argument out of range
>>>>>> [1]PETSC ERROR: Inserting a new nonzero (40003, 0) into matrix
>>>>>>
>>>>>> Could this have to do with me using the DMComposite with one da holding
>>>>>> the scalar parameters and the other holding the field variables?
>>>>>>
>>>>>> -gideon
>>>>>>
>>>>>>> On Aug 27, 2015, at 10:15 PM, Matthew Knepley <[email protected]> wrote:
>>>>>>>
>>>>>>> On Thu, Aug 27, 2015 at 9:11 PM, Gideon Simpson
>>>>>>> <[email protected]> wrote:
>>>>>>> HI Barry,
>>>>>>>
>>>>>>> Nope, I’m not doing any grid sequencing. Clearly that makes a lot of
>>>>>>> sense, to solve on a spatially coarse mesh for the field variables,
>>>>>>> interpolate onto the finer mesh, and then solve again. I’m not
>>>>>>> entirely clear on the practical implementation
>>>>>>>
>>>>>>> SNES should do this automatically using -snes_grid_sequence <k>. If
>>>>>>> this does not work, complain. Loudly.
>>>>>>>
>>>>>>> Matt
>>>>>>>
>>>>>>> -gideon
>>>>>>>
>>>>>>>> On Aug 27, 2015, at 10:02 PM, Barry Smith <[email protected]> wrote:
>>>>>>>>
>>>>>>>>
>>>>>>>> Gideon,
>>>>>>>>
>>>>>>>> Are you using grid sequencing? Simply solve on a coarse grid,
>>>>>>>> interpolate u1 and u2 to a once refined version of the grid and use
>>>>>>>> that plus the mu lam as initial guess for the next level. Repeat to as
>>>>>>>> fine a grid as you want. You can use DMRefine() and
>>>>>>>> DMGetInterpolation() to get the interpolation needed to interpolate
>>>>>>>> from the coarse to finer mesh.
>>>>>>>>
>>>>>>>> Then and only then you can use multigrid (with or without fieldsplit)
>>>>>>>> to solve the linear problems for finer meshes. Once you have the grid
>>>>>>>> sequencing working we can help you with this.
>>>>>>>>
>>>>>>>> Barry
>>>>>>>>
>>>>>>>>> On Aug 27, 2015, at 7:00 PM, Gideon Simpson
>>>>>>>>> <[email protected]> wrote:
>>>>>>>>>
>>>>>>>>> I’m working on a problem which, morally, can be posed as a system of
>>>>>>>>> coupled semi linear elliptic PDEs together with unknown nonlinear
>>>>>>>>> eigenvalue parameters, loosely, of the form
>>>>>>>>>
>>>>>>>>> -\Delta u_1 + f(u_1, u_2) = lam * u1 - mu * du2/dx
>>>>>>>>> -\Delta u_2 + g(u_1, u_2) = lam * u2 + mu * du1/dx
>>>>>>>>>
>>>>>>>>> Currently, I have it set up with a DMComposite with two sub da’s, one
>>>>>>>>> for the parameters (lam, mu), and one for the vector field (u_1, u_2)
>>>>>>>>> on the mesh. I have had success in solving this as a fully coupled
>>>>>>>>> system with SNES + sparse direct solvers (MUMPS, SuperLU).
>>>>>>>>>
>>>>>>>>> Lately, I am finding that, when the mesh resolution gets fine enough
>>>>>>>>> (i.e. 10^6-10^8 lattice points), my SNES gets stuck with the
>>>>>>>>> function norm = O(10^{-4}), eventually returning reason -6 (failed
>>>>>>>>> line search).
>>>>>>>>>
>>>>>>>>> Perhaps there is another way around the above problem, but one thing
>>>>>>>>> I was thinking of trying would be to get away from direct solvers,
>>>>>>>>> and I was hoping to use field split for this. However, it’s a bit
>>>>>>>>> beyond what I’ve seen examples for because it has 2 types of
>>>>>>>>> variables: scalar parameters which appear globally in the system and
>>>>>>>>> vector valued field variables. Any suggestions on how to get started?
>>>>>>>>>
>>>>>>>>> -gideon
>>>>>>>>>
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> --
>>>>>>> What most experimenters take for granted before they begin their
>>>>>>> experiments is infinitely more interesting than any results to which
>>>>>>> their experiments lead.
>>>>>>> -- Norbert Wiener
>>>>>>
>>>>>
>>>>
>>>
>>
>
