<html><head><meta http-equiv="Content-Type" content="text/html charset=utf-8"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space;" class="">Yes, to clarify, the problem with two scalar fields and two scalar “eigenvalues” is<div class=""><br class=""></div><div class=""><br class="">-\Delta u_1 + f(u_1, u_2) = lam * u1 - mu * du2/dx <br class="">-\Delta u_2 + g(u_1, u_2) = lam * u2 + mu * du1/dx <br class=""><br class=""><div class="">(u_1, u_2) are defined on the mesh, and (lam, mu) are unknown scalars. My actual problem has a 4 degrees of freedom at each mesh point and 3 unknown scalars, but this makes the point.</div><div class=""><br class="webkit-block-placeholder"></div><div apple-content-edited="true" class="">
<span class="Apple-style-span" style="border-collapse: separate; border-spacing: 0px;">-gideon</span>
</div>
<br class=""><div><blockquote type="cite" class=""><div class="">On Aug 28, 2015, at 6:20 AM, Matthew Knepley <<a href="mailto:knepley@gmail.com" class="">knepley@gmail.com</a>> wrote:</div><br class="Apple-interchange-newline"><div class=""><div dir="ltr" style="font-family: Helvetica; font-size: 14px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px; -webkit-text-stroke-width: 0px;" class=""><div class="gmail_extra"><div class="gmail_quote">On Thu, Aug 27, 2015 at 10:23 PM, Barry Smith<span class="Apple-converted-space"> </span><span dir="ltr" class=""><<a href="mailto:bsmith@mcs.anl.gov" target="_blank" class="">bsmith@mcs.anl.gov</a>></span><span class="Apple-converted-space"> </span>wrote:<br class=""><blockquote class="gmail_quote" style="margin: 0px 0px 0px 0.8ex; border-left-width: 1px; border-left-color: rgb(204, 204, 204); border-left-style: solid; padding-left: 1ex;"><span class=""><br class="">> On Aug 27, 2015, at 10:15 PM, Gideon Simpson <<a href="mailto:gideon.simpson@gmail.com" class="">gideon.simpson@gmail.com</a>> wrote:<br class="">><br class="">> 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.<br class=""><br class=""></span> <span class="Apple-converted-space"> </span>Nothing ever gets coarsened in grid sequencing, it only gets refined.<br class=""><br class=""> <span class="Apple-converted-space"> </span>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).<br class=""><span class=""></span></blockquote><div class=""><br class=""></div><div class="">I think you might be misunderstanding the "scalars" part. He is solving a nonlinear eigenproblem (which he did not write down) for some variables. Then he</div><div class="">uses those variable in the coupled diffusion equations he did write down. He has wrapped the whole problem in a SNES with 2 parts: the nonlinear eigenproblem</div><div class="">and the diffusion equations. He uses DMComposite to deal with all the unknowns.</div><div class=""><br class=""></div><div class="">I think Nonlinear Block Gauss-Siedel on the different problems would be a useful starting point, but we do not have that.</div><div class=""><br class=""></div><div class=""> Thanks,</div><div class=""><br class=""></div><div class=""> Matt</div><div class=""> </div><blockquote class="gmail_quote" style="margin: 0px 0px 0px 0.8ex; border-left-width: 1px; border-left-color: rgb(204, 204, 204); border-left-style: solid; padding-left: 1ex;"><span class="">><br class="">> Here is my form function. I can send more code if needed.<br class="">><br class=""><br class=""></span> <span class="Apple-converted-space"> </span>Just change the user->packer that you use to be the DM obtained with SNESGetDM()<br class=""><div class="HOEnZb"><div class="h5"><br class=""><br class="">> /* Form the system of equations for computing a blowup solution*/<br class="">> PetscErrorCode form_function(SNES snes, Vec U, Vec F, void *ctx){<br class="">><br class="">> blowup_ctx *user = (blowup_ctx *) ctx;<br class="">> PetscInt i;<br class="">> PetscScalar dx, dx2, xmax,x;<br class="">> PetscScalar u, v, f,g, ux, vx, uxx, vxx, fx,gx, fxx, gxx;<br class="">> DMDALocalInfo info;<br class="">> Vec p_vec, Q_vec, Fp_vec, FQ_vec;<br class="">> PetscScalar *p_array, *Fp_array;<br class="">> Q *Qvals, *FQvals;<br class="">> PetscScalar Q2sig, W2sig;<br class="">> PetscScalar a,a2, b, u0, sigma;<br class="">><br class="">> dx = user->dx; dx2 = dx *dx;<br class="">> xmax = user->xmax;<br class="">> sigma = user->sigma;<br class="">><br class="">> /* PetscPrintf(PETSC_COMM_SELF, " dx = %g, sigma = %g\n", dx, sigma); */<br class="">><br class="">> /* Extract raw arrays */<br class="">> DMCompositeGetLocalVectors(user->packer, &p_vec, &Q_vec);<br class="">> DMCompositeGetLocalVectors(user->packer, &Fp_vec, &FQ_vec);<br class="">><br class="">> DMCompositeScatter(user->packer, U, p_vec, Q_vec);<br class="">> /* VecView(Q_vec, PETSC_VIEWER_STDOUT_SELF); */<br class="">><br class="">> VecGetArray(p_vec,&p_array);<br class="">> VecGetArray(Fp_vec,&Fp_array);<br class="">><br class="">> DMDAVecGetArray(user->Q_dm, Q_vec, &Qvals);<br class="">> DMDAVecGetArray(user->Q_dm, FQ_vec, &FQvals);<br class="">><br class="">> DMDAGetLocalInfo(user->Q_dm, &info);<br class="">><br class="">> a = p_array[0]; a2 = a*a;<br class="">> b = p_array[1];<br class="">> u0 = p_array[2];<br class="">><br class="">> /* Set boundary conditions at the origin*/<br class="">> if(info.xs ==0){<br class="">> set_origin_bcs(u0, Qvals);<br class="">> }<br class="">> /* Set boundray conditions in the far field */<br class="">> if(info.xs+ info.xm ==<span class="Apple-converted-space"> </span><a href="http://info.mx/" rel="noreferrer" target="_blank" class="">info.mx</a>){<br class="">> set_farfield_bcs(xmax,sigma, a, b, dx, Qvals,<a href="http://info.mx/" rel="noreferrer" target="_blank" class="">info.mx</a>);<br class="">> }<br class="">><br class="">> /* Solve auxiliary equations */<br class="">> if(info.xs ==0){<br class="">> uxx = (2 * Qvals[0].u-2 * u0)/dx2;<br class="">> vxx = (Qvals[0].v + Qvals[0].g)/dx2;<br class="">> vx = (Qvals[0].v - Qvals[0].g)/(2*dx);<br class="">> Fp_array[0] = Qvals[0].u - Qvals[0].f;<br class="">> Fp_array[1] = -vxx - (1/a) * (.5/sigma) * u0;<br class="">> Fp_array[2] = -uxx + (1/a2) * u0<br class="">> + (1/a) * (-b * vx + PetscPowScalar(u0 * u0, sigma) * vx);<br class="">> }<br class="">><br class="">> /* Solve equations in the bulk */<br class="">> for(i=info.xs; i < info.xs + info.xm;i++){<br class="">><br class="">> u = Qvals[i].u;<br class="">> v = Qvals[i].v;<br class="">> f = Qvals[i].f;<br class="">> g = Qvals[i].g;<br class="">><br class="">> x = (i+1) * dx;<br class="">><br class="">> Q2sig = PetscPowScalar(u*u + v*v,sigma);<br class="">> W2sig= PetscPowScalar(f*f + g*g, sigma);<br class="">><br class="">> ux = (Qvals[i+1].u-Qvals[i-1].u)/(2*dx);<br class="">> vx = (Qvals[i+1].v-Qvals[i-1].v)/(2*dx);<br class="">> fx = (Qvals[i+1].f-Qvals[i-1].f)/(2*dx);<br class="">> gx = (Qvals[i+1].g-Qvals[i-1].g)/(2*dx);<br class="">><br class="">> uxx = (Qvals[i+1].u+Qvals[i-1].u - 2 *u)/(dx2);<br class="">> vxx = (Qvals[i+1].v+Qvals[i-1].v- 2 *v)/(dx2);<br class="">> fxx = (Qvals[i+1].f+Qvals[i-1].f -2*f)/(dx2);<br class="">> gxx = (Qvals[i+1].g+Qvals[i-1].g -2*g)/(dx2);<br class="">><br class="">> FQvals[i].u = -uxx +1/a2 * u<br class="">> + 1/a *(.5/sigma* v +x * vx- b* vx +Q2sig* vx);<br class="">><br class="">> FQvals[i].v = -vxx +1/a2 * v<br class="">> - 1/a *(.5/sigma * u +x * ux- b* ux +Q2sig* ux);<br class="">><br class="">> FQvals[i].f = -fxx +1/a2 * f<br class="">> + 1/a *(.5/sigma * g +x * gx+ b* gx -W2sig* gx);<br class="">><br class="">> FQvals[i].g =-gxx +1/a2 * g<br class="">> - 1/a *(.5/sigma * f +x * fx+ b* fx -W2sig* fx);<br class="">> }<br class="">><br class="">> /* Restore raw arrays */<br class="">> VecRestoreArray(p_vec, &p_array);<br class="">> VecRestoreArray(Fp_vec, &Fp_array);<br class="">><br class="">> DMDAVecRestoreArray(user->Q_dm, Q_vec, &Qvals);<br class="">> DMDAVecRestoreArray(user->Q_dm, FQ_vec, &FQvals);<br class="">><br class="">> DMCompositeGather(user->packer,F,INSERT_VALUES, Fp_vec, FQ_vec);<br class="">> DMCompositeRestoreLocalVectors(user->packer, &p_vec, &Q_vec);<br class="">> DMCompositeRestoreLocalVectors(user->packer, &Fp_vec, &FQ_vec);<br class="">><br class="">> return 0;<br class="">> }<br class="">><br class="">><br class="">> Here is the form function:<br class="">><br class="">><br class="">><br class="">> -gideon<br class="">><br class="">>> On Aug 27, 2015, at 11:09 PM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov" class="">bsmith@mcs.anl.gov</a>> wrote:<br class="">>><br class="">>><br class="">>> Can you send the code, that will be the easiest way to find the problem.<br class="">>><br class="">>> 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.<br class="">>><br class="">>> Barry<br class="">>><br class="">>>> On Aug 27, 2015, at 9:42 PM, Gideon Simpson <<a href="mailto:gideon.simpson@gmail.com" class="">gideon.simpson@gmail.com</a>> wrote:<br class="">>>><br class="">>>> I have it set up as:<br class="">>>><br class="">>>> DMCompositeCreate(PETSC_COMM_WORLD, &user.packer);<br class="">>>> DMRedundantCreate(PETSC_COMM_WORLD, 0, 3, &user.p_dm);<br class="">>>> DMCompositeAddDM(user.packer,user.p_dm);<br class="">>>> DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_GHOSTED,<br class="">>>> nx, 4, 1, NULL, &user.Q_dm);<br class="">>>> DMCompositeAddDM(user.packer,user.Q_dm);<br class="">>>> DMCreateGlobalVector(user.packer,&U);<br class="">>>><br class="">>>> where the user.packer structure has<br class="">>>><br class="">>>> DM packer;<br class="">>>> DM p_dm, Q_dm;<br class="">>>><br class="">>>> Q_dm holds the field variables and p_dm holds the scalar values (the nonlinear eigenvalues).<br class="">>>><br class="">>>> Here are some of the errors that are generated:<br class="">>>><br class="">>>> [0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------<br class="">>>> [0]PETSC ERROR: Argument out of range<br class="">>>> [0]PETSC ERROR: New nonzero at (0,3) caused a malloc<br class="">>>> Use MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to turn off this check<br class="">>>> [0]PETSC ERROR: See<span class="Apple-converted-space"> </span><a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" rel="noreferrer" target="_blank" class="">http://www.mcs.anl.gov/petsc/documentation/faq.html</a><span class="Apple-converted-space"> </span>for trouble shooting.<br class="">>>> [0]PETSC ERROR: Petsc Release Version 3.5.3, unknown<br class="">>>> [0]PETSC ERROR: ./blowup_batch2 on a arch-macports named gs_air by gideon Thu Aug 27 22:40:54 2015<br class="">>>> [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<br class="">>>> [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<br class="">>>> [0]PETSC ERROR: [1]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------<br class="">>>> [1]PETSC ERROR: Argument out of range<br class="">>>> [1]PETSC ERROR: Inserting a new nonzero (40003, 0) into matrix<br class="">>>> [1]PETSC ERROR: See<span class="Apple-converted-space"> </span><a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" rel="noreferrer" target="_blank" class="">http://www.mcs.anl.gov/petsc/documentation/faq.html</a><span class="Apple-converted-space"> </span>for trouble shooting.<br class="">>>> [1]PETSC ERROR: Petsc Release Version 3.5.3, unknown<br class="">>>> [1]PETSC ERROR: ./blowup_batch2 on a arch-macports named gs_air by gideon Thu Aug 27 22:40:54 2015<br class="">>>> [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<br class="">>>> [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<br class="">>>> [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<br class="">>>><br class="">>>><br class="">>>><br class="">>>> -gideon<br class="">>>><br class="">>>>> On Aug 27, 2015, at 10:37 PM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov" class="">bsmith@mcs.anl.gov</a>> wrote:<br class="">>>>><br class="">>>>><br class="">>>>> We need the full error message.<br class="">>>>><br class="">>>>> But are you using a DMDA for the scalars? You should not be, you should be using a DMRedundant for the scalars.<br class="">>>>><br class="">>>>> Barry<br class="">>>>><br class="">>>>> Though you should not get this error even if you are using a DMDA there.<br class="">>>>><br class="">>>>>> On Aug 27, 2015, at 9:32 PM, Gideon Simpson <<a href="mailto:gideon.simpson@gmail.com" class="">gideon.simpson@gmail.com</a>> wrote:<br class="">>>>>><br class="">>>>>> I’m getting the following errors:<br class="">>>>>><br class="">>>>>> [1]PETSC ERROR: Argument out of range<br class="">>>>>> [1]PETSC ERROR: Inserting a new nonzero (40003, 0) into matrix<br class="">>>>>><br class="">>>>>> Could this have to do with me using the DMComposite with one da holding the scalar parameters and the other holding the field variables?<br class="">>>>>><br class="">>>>>> -gideon<br class="">>>>>><br class="">>>>>>> On Aug 27, 2015, at 10:15 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com" class="">knepley@gmail.com</a>> wrote:<br class="">>>>>>><br class="">>>>>>> On Thu, Aug 27, 2015 at 9:11 PM, Gideon Simpson <<a href="mailto:gideon.simpson@gmail.com" class="">gideon.simpson@gmail.com</a>> wrote:<br class="">>>>>>> HI Barry,<br class="">>>>>>><br class="">>>>>>> 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<br class="">>>>>>><br class="">>>>>>> SNES should do this automatically using -snes_grid_sequence <k>. If this does not work, complain. Loudly.<br class="">>>>>>><br class="">>>>>>> Matt<br class="">>>>>>><br class="">>>>>>> -gideon<br class="">>>>>>><br class="">>>>>>>> On Aug 27, 2015, at 10:02 PM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov" class="">bsmith@mcs.anl.gov</a>> wrote:<br class="">>>>>>>><br class="">>>>>>>><br class="">>>>>>>> Gideon,<br class="">>>>>>>><br class="">>>>>>>> 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.<br class="">>>>>>>><br class="">>>>>>>> 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.<br class="">>>>>>>><br class="">>>>>>>> Barry<br class="">>>>>>>><br class="">>>>>>>>> On Aug 27, 2015, at 7:00 PM, Gideon Simpson <<a href="mailto:gideon.simpson@gmail.com" class="">gideon.simpson@gmail.com</a>> wrote:<br class="">>>>>>>>><br class="">>>>>>>>> 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<br class="">>>>>>>>><br class="">>>>>>>>> -\Delta u_1 + f(u_1, u_2) = lam * u1 - mu * du2/dx<br class="">>>>>>>>> -\Delta u_2 + g(u_1, u_2) = lam * u2 + mu * du1/dx<br class="">>>>>>>>><br class="">>>>>>>>> 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).<br class="">>>>>>>>><br class="">>>>>>>>> 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).<br class="">>>>>>>>><br class="">>>>>>>>> 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?<br class="">>>>>>>>><br class="">>>>>>>>> -gideon<br class="">>>>>>>>><br class="">>>>>>>><br class="">>>>>>><br class="">>>>>>><br class="">>>>>>><br class="">>>>>>><br class="">>>>>>> --<br class="">>>>>>> What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br class="">>>>>>> -- Norbert Wiener<br class="">>>>>><br class="">>>>><br class="">>>><br class="">>><br class="">><br class=""><br class=""></div></div></blockquote></div><br class=""><br clear="all" class=""><div class=""><br class=""></div>--<span class="Apple-converted-space"> </span><br class=""><div class="gmail_signature">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br class="">-- Norbert Wiener</div></div></div></div></blockquote></div><br class=""></div></body></html>