Re: [petsc-users] VTS writer
Dear All
I forgot to put the 3 or more lines of code I am using. Here they are.
I am using the following 3 lines to write vts files. It is a structured FD
mesh. I now have difference space steps in the 3 directions: it is 1 and 1 in
X and Y directions, and 2 in Z direction. Can anyone say how to write this to
the vts header?
// FD spacing in each direction
#define DX 0.250/* X internode spacing
*/
#define DY 0.250/* Y internode spacing
*/
#define DZ 0.500/* Z internode spacing: See suppl of
paper. */
.
.
main(){
.
.
sprintf(str,my_3d%d.vts,file_Counter++);
PetscViewer viewer;
PetscViewerVTKOpen(PETSC_COMM_WORLD,str,FILE_MODE_WRITE, viewer);
VecView(u, viewer);
PetscViewerDestroy(viewer);
Re: [petsc-users] VTS writer
Kharche, Sanjay s.r.khar...@exeter.ac.uk writes:
Dear All
I forgot to put the 3 or more lines of code I am using. Here they are.
I am using the following 3 lines to write vts files. It is a structured FD
mesh. I now have difference space steps in the 3 directions: it is 1 and 1
in X and Y directions, and 2 in Z direction. Can anyone say how to write
this to the vts header?
// FD spacing in each direction
#define DX 0.250/* X internode spacing
*/
#define DY 0.250/* Y internode spacing
*/
#define DZ 0.500/* Z internode spacing: See suppl
of paper. */
DMDASetUniformCoordinates
.
.
main(){
.
.
sprintf(str,my_3d%d.vts,file_Counter++);
PetscViewer viewer;
PetscViewerVTKOpen(PETSC_COMM_WORLD,str,FILE_MODE_WRITE, viewer);
VecView(u, viewer);
PetscViewerDestroy(viewer);
signature.asc
Description: PGP signature
Re: [petsc-users] Tao iterations
Hi Justin, This expected behavior due to the accumulation of numerical round-offs. If this is a problem or if you just want to confirm that this is the cause, you can try configuring PETSc for quad precision (--with-precision=__float128, works with GNU compilers) and the results should match better. Jason On Tue, Apr 28, 2015 at 10:19 PM, Justin Chang jchan...@uh.edu wrote: Jason (or anyone), I am noticing that the iteration numbers reported by TaoGetSolutionStatus() for blmvm differ whenever I change the number of processes. The solution seems to remain the same though. Is there a reason why this could be happening? Thanks, On Tue, Apr 21, 2015 at 10:40 AM, Jason Sarich jason.sar...@gmail.com wrote: Justin, 1) The big difference between TRON and BLMVM is that TRON requires hessian information, BLMVM only uses gradient information. Thus TRON will usually converge faster, but requires more information, memory, and a KSP solver. GPCG (gradient projected conjugate gradient) is another gradient-only option, but usually performs worse than BLMVM. 2) TaoGetLinearSolveIterations() will get the total number of KSP iterations per solve Jason On Tue, Apr 21, 2015 at 10:33 AM, Justin Chang jychan...@gmail.com wrote: Jason, Tightening the tolerances did the trick. Thanks. Though I do have a couple more related questions: 1) Is there a general guideline for choosing tron over blmvm or vice versa? Also is there another tao type that is also suitable given only bounded constraints? 2) Is it possible to obtain the total number of KSP and/or PG iterations from tron? Thanks, Justin On Tue, Apr 21, 2015 at 9:52 AM, Jason Sarich jason.sar...@gmail.com wrote: Hi Justin, blmvm believes that it is already sufficiently close to a minimum, so it doesn't do anything. You may need to tighten some of the tolerance to force an iteration. Jason On Tue, Apr 21, 2015 at 9:48 AM, Justin Chang jychan...@gmail.com wrote: Time step 1: Tao Object: 1 MPI processes type: blmvm Gradient steps: 0 TaoLineSearch Object: 1 MPI processes type: more-thuente Active Set subset type: subvec convergence tolerances: fatol=0.0001, frtol=0.0001 convergence tolerances: gatol=0, steptol=0, gttol=0 Residual in Function/Gradient:=0.0663148 Objective value=-55.5945 total number of iterations=35, (max: 2000) total number of function/gradient evaluations=37, (max: 4000) Solution converged: estimated |f(x)-f(X*)|/|f(X*)| = frtol Time step 2: Tao Object: 1 MPI processes type: blmvm Gradient steps: 0 TaoLineSearch Object: 1 MPI processes type: more-thuente Active Set subset type: subvec convergence tolerances: fatol=0.0001, frtol=0.0001 convergence tolerances: gatol=0, steptol=0, gttol=0 Residual in Function/Gradient:=0.0682307 Objective value=-66.9675 total number of iterations=23, (max: 2000) total number of function/gradient evaluations=25, (max: 4000) Solution converged: estimated |f(x)-f(X*)|/|f(X*)| = frtol Time step 3: Tao Object: 1 MPI processes type: blmvm Gradient steps: 0 TaoLineSearch Object: 1 MPI processes type: more-thuente Active Set subset type: subvec convergence tolerances: fatol=0.0001, frtol=0.0001 convergence tolerances: gatol=0, steptol=0, gttol=0 Residual in Function/Gradient:=0.0680522 Objective value=-71.8211 total number of iterations=19, (max: 2000) total number of function/gradient evaluations=22, (max: 4000) Solution converged: estimated |f(x)-f(X*)|/|f(X*)| = frtol Time step 4: Tao Object: 1 MPI processes type: blmvm Gradient steps: 0 TaoLineSearch Object: 1 MPI processes type: more-thuente Active Set subset type: subvec convergence tolerances: fatol=0.0001, frtol=0.0001 convergence tolerances: gatol=0, steptol=0, gttol=0 Residual in Function/Gradient:=0.0551556 Objective value=-75.1252 total number of iterations=18, (max: 2000) total number of function/gradient evaluations=20, (max: 4000) Solution converged: estimated |f(x)-f(X*)|/|f(X*)| = frtol Time step 5: Tao Object: 1 MPI processes type: blmvm Gradient steps: 0 TaoLineSearch Object: 1 MPI processes type: more-thuente Active Set subset type: subvec convergence tolerances: fatol=0.0001, frtol=0.0001 convergence tolerances: gatol=0, steptol=0, gttol=0 Residual in Function/Gradient:=0.0675667 Objective value=-77.4414 total number of iterations=6, (max: 2000) total number of function/gradient evaluations=8, (max: 4000) Solution converged: estimated |f(x)-f(X*)|/|f(X*)| = frtol Time step 6: Tao Object: 1 MPI processes type: blmvm Gradient steps: 0 TaoLineSearch
Re: [petsc-users] Floating point exception in hypre BoomerAMG
Ok, run on one process in the debugger and put a break point in VecLockPush and VecLockPopthen run and each time it breaks on VecLockPush or VecLockPop type bt then type cont send all the output, this will tell us where the vector got locked readonly and did not get unlocked. Barry On Apr 29, 2015, at 2:35 PM, Danyang Su danyang...@gmail.com wrote: On 15-04-29 12:19 PM, Barry Smith wrote: Ok, your code seems fine in terms of logic. But the vector x_vec_gbl i.e. b_react should not be in a read only state. Are you somewhere calling a VecGetArray() and not calling the VecRestoreArray()? Barry I just made a triple check on VecGetArrayF90(). All the VecGetArrayF90() come with VecRestoreArrayF90(). No VecGetArray() is used in my codes. Danyang On Apr 29, 2015, at 1:52 PM, Danyang Su danyang...@gmail.com wrote: On 15-04-29 11:30 AM, Barry Smith wrote: On Apr 29, 2015, at 12:15 PM, Danyang Su danyang...@gmail.com wrote: On 15-04-28 06:50 PM, Barry Smith wrote: We started enforcing more checks on writing to vectors that you should not write to. Where are you calling DMLocalToGlobalBegin() ? It looks like you are trying to copy values into a global array that you should not because it is a read only input to a function, for example it is the right hand side of a linear system or the input into a SNESFormFunction(). So check the output to that call. Barry This is used when local RHS vec is calculated and copied to global RHS. Huhh, why are you copying anything into the RHS? Is this before you call the linear system solve? Send the code that calls the DMLocalToGlobalBegin() Yes. This is called before the linear system solver. At present, we use PETSc KSP solver. The code is first developed in totally sequential version and then ported to the parallel version. So the subdomain has its own space for jacobi matrix and rhs. We do not use global RHS vector directly to set the RHS at this moment. All the RHS values of subdomain are copied to the global RHS. ! This is the function where RHS values of subdomain are copied to the global RHS values. subroutine compute_function(rank,da,x_array_loc,x_vec_loc, x_vec_gbl,nngl,row_idx_l2pg,col_idx_l2pg, b_non_interlaced) implicit none #include finclude/petscsys.h #include finclude/petscis.h #include finclude/petscvec.h #include finclude/petscvec.h90 #include finclude/petscmat.h #include finclude/petscpc.h #include finclude/petscksp.h #include finclude/petscdmda.h #include finclude/petscdmda.h90 PetscInt :: rank DM :: da PetscReal, allocatable :: x_array_loc(:) Vec :: x_vec_loc Vec :: x_vec_gbl PetscInt :: nngl PetscInt, allocatable :: row_idx_l2pg(:) PetscInt, allocatable :: col_idx_l2pg(:) PetscBool :: b_non_interlaced PetscInt :: info_debug PetscInt :: i, j PetscErrorCode :: ierr PetscScalar, pointer :: vecpointer(:) !Zero entries call VecZeroEntries(x_vec_loc, ierr) !Get a pointer to vector data when you need access to the array call VecGetArrayF90(x_vec_loc, vecpointer, ierr) !Compute the function over the locally owned part of the grid if(b_non_interlaced) then j = nngl/2 do i = 1, j vecpointer(2*i-1) = x_array_loc(i) vecpointer(2*i) = x_array_loc(i+j) end do else do i = 1, nngl vecpointer(i) = x_array_loc(i) end do end if !Restore the vector when you no longer need access to the array call VecRestoreArrayF90(x_vec_loc,vecpointer,ierr) !Insert values into global vector call DMLocalToGlobalBegin(da,x_vec_loc,INSERT_VALUES, x_vec_gbl,ierr) !By placing code between these two statements, computations can be !done while messages are in transition. call DMLocalToGlobalEnd(da,x_vec_loc,INSERT_VALUES, x_vec_gbl,ierr) return end subroutine ! This is the function where reactive transport equations are solved. subroutine solver_dd_snes_solve_react(ilog,idetail,a_in,b_in, x_inout,ia_in,ja_in,nngl_in,itsolv, over_flow,rnorm,row_idx_l2pg,col_idx_l2pg, b_non_interlaced) use gen, only : rank, node_idx_l2lg, ittot_rt, b_output_matrix, b_enable_output use solver_snes_function, only :
Re: [petsc-users] Discrete adjoint and adaptive time stepping
On Apr 29, 2015, at 2:48 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: If we calculate the gradient using the discrete adjoint without differentiating the controller, and then calculate the same gradient using finite difference (allowing the time steps to freely change), how different these results are? Clearly I don't know but I image they could be very different. I would do as Jed suggests and use one run to get the adaptive dt and then use that same collect of dt for all the runs. Barry Miguel On Wed, Apr 29, 2015 at 2:12 PM, Barry Smith bsm...@mcs.anl.gov wrote: On Apr 29, 2015, at 1:51 PM, Jed Brown j...@jedbrown.org wrote: Barry Smith bsm...@mcs.anl.gov writes: If so, how is it done? We just keep the history of the time-step sizes and then use those time-step sizes when doing the backward integration. Seems simple to me, am I missing something? Barry, if you're using this for optimization, you might want the gradient to be exactly consistent with the objective functional. But for that, you would need to differentiate the controller, which is non-smooth in practice because the number of time steps can change and stages could be rejected (solver failure). Ahh, yes for multiple forward runs yup. One approach would be to save the timestep sequence and have the controller use that in subsequent *forward* runs. If the dynamical system behaves similarly for those steps, it would be okay to use the same timestep sequence. Presumably if that single set of dt (from the first run) is not sufficient for some later runs one could possibly use the union of the dt of several runs for all the runs. (that is run adaptively and inconsistently several runs to determine where dt needs to be controlled and then use the various smaller of the dt at the different time regions for a full set of consistent runs). Of course if the various smaller of the dt requires a tiny dt for all time steps then you are not getting an advantage of adaptive time-stepping, but ok. The idea of actually propagating the gradients through the time-step controller seems IMHO to be absurd; I won't even put it on our game plan until we have many more things done and much more practical experience. Barry -- Miguel Angel Salazar de Troya Graduate Research Assistant Department of Mechanical Science and Engineering University of Illinois at Urbana-Champaign (217) 550-2360 salaz...@illinois.edu
Re: [petsc-users] Floating point exception in hypre BoomerAMG
Ok, your code seems fine in terms of logic. But the vector x_vec_gbl i.e. b_react should not be in a read only state. Are you somewhere calling a VecGetArray() and not calling the VecRestoreArray()? Barry On Apr 29, 2015, at 1:52 PM, Danyang Su danyang...@gmail.com wrote: On 15-04-29 11:30 AM, Barry Smith wrote: On Apr 29, 2015, at 12:15 PM, Danyang Su danyang...@gmail.com wrote: On 15-04-28 06:50 PM, Barry Smith wrote: We started enforcing more checks on writing to vectors that you should not write to. Where are you calling DMLocalToGlobalBegin() ? It looks like you are trying to copy values into a global array that you should not because it is a read only input to a function, for example it is the right hand side of a linear system or the input into a SNESFormFunction(). So check the output to that call. Barry This is used when local RHS vec is calculated and copied to global RHS. Huhh, why are you copying anything into the RHS? Is this before you call the linear system solve? Send the code that calls the DMLocalToGlobalBegin() Yes. This is called before the linear system solver. At present, we use PETSc KSP solver. The code is first developed in totally sequential version and then ported to the parallel version. So the subdomain has its own space for jacobi matrix and rhs. We do not use global RHS vector directly to set the RHS at this moment. All the RHS values of subdomain are copied to the global RHS. ! This is the function where RHS values of subdomain are copied to the global RHS values. subroutine compute_function(rank,da,x_array_loc,x_vec_loc, x_vec_gbl,nngl,row_idx_l2pg,col_idx_l2pg, b_non_interlaced) implicit none #include finclude/petscsys.h #include finclude/petscis.h #include finclude/petscvec.h #include finclude/petscvec.h90 #include finclude/petscmat.h #include finclude/petscpc.h #include finclude/petscksp.h #include finclude/petscdmda.h #include finclude/petscdmda.h90 PetscInt :: rank DM :: da PetscReal, allocatable :: x_array_loc(:) Vec :: x_vec_loc Vec :: x_vec_gbl PetscInt :: nngl PetscInt, allocatable :: row_idx_l2pg(:) PetscInt, allocatable :: col_idx_l2pg(:) PetscBool :: b_non_interlaced PetscInt :: info_debug PetscInt :: i, j PetscErrorCode :: ierr PetscScalar, pointer :: vecpointer(:) !Zero entries call VecZeroEntries(x_vec_loc, ierr) !Get a pointer to vector data when you need access to the array call VecGetArrayF90(x_vec_loc, vecpointer, ierr) !Compute the function over the locally owned part of the grid if(b_non_interlaced) then j = nngl/2 do i = 1, j vecpointer(2*i-1) = x_array_loc(i) vecpointer(2*i) = x_array_loc(i+j) end do else do i = 1, nngl vecpointer(i) = x_array_loc(i) end do end if !Restore the vector when you no longer need access to the array call VecRestoreArrayF90(x_vec_loc,vecpointer,ierr) !Insert values into global vector call DMLocalToGlobalBegin(da,x_vec_loc,INSERT_VALUES, x_vec_gbl,ierr) !By placing code between these two statements, computations can be !done while messages are in transition. call DMLocalToGlobalEnd(da,x_vec_loc,INSERT_VALUES, x_vec_gbl,ierr) return end subroutine ! This is the function where reactive transport equations are solved. subroutine solver_dd_snes_solve_react(ilog,idetail,a_in,b_in, x_inout,ia_in,ja_in,nngl_in,itsolv, over_flow,rnorm,row_idx_l2pg,col_idx_l2pg, b_non_interlaced) use gen, only : rank, node_idx_l2lg, ittot_rt, b_output_matrix, b_enable_output use solver_snes_function, only : form_initial_guess, compute_function, compute_jacobian use petsc_mpi_common, only : petsc_mpi_finalize implicit none #include finclude/petscsys.h #include finclude/petscvec.h #include finclude/petscvec.h90 #include finclude/petscdmda.h PetscInt :: ilog PetscInt :: idetail PetscInt :: nngl_in PetscReal, allocatable :: a_in(:) PetscReal, allocatable :: b_in(:)
Re: [petsc-users] Discrete adjoint and adaptive time stepping
If we calculate the gradient using the discrete adjoint without differentiating the controller, and then calculate the same gradient using finite difference (allowing the time steps to freely change), how different these results are? Miguel On Wed, Apr 29, 2015 at 2:12 PM, Barry Smith bsm...@mcs.anl.gov wrote: On Apr 29, 2015, at 1:51 PM, Jed Brown j...@jedbrown.org wrote: Barry Smith bsm...@mcs.anl.gov writes: If so, how is it done? We just keep the history of the time-step sizes and then use those time-step sizes when doing the backward integration. Seems simple to me, am I missing something? Barry, if you're using this for optimization, you might want the gradient to be exactly consistent with the objective functional. But for that, you would need to differentiate the controller, which is non-smooth in practice because the number of time steps can change and stages could be rejected (solver failure). Ahh, yes for multiple forward runs yup. One approach would be to save the timestep sequence and have the controller use that in subsequent *forward* runs. If the dynamical system behaves similarly for those steps, it would be okay to use the same timestep sequence. Presumably if that single set of dt (from the first run) is not sufficient for some later runs one could possibly use the union of the dt of several runs for all the runs. (that is run adaptively and inconsistently several runs to determine where dt needs to be controlled and then use the various smaller of the dt at the different time regions for a full set of consistent runs). Of course if the various smaller of the dt requires a tiny dt for all time steps then you are not getting an advantage of adaptive time-stepping, but ok. The idea of actually propagating the gradients through the time-step controller seems IMHO to be absurd; I won't even put it on our game plan until we have many more things done and much more practical experience. Barry -- *Miguel Angel Salazar de Troya* Graduate Research Assistant Department of Mechanical Science and Engineering University of Illinois at Urbana-Champaign (217) 550-2360 salaz...@illinois.edu
Re: [petsc-users] Floating point exception in hypre BoomerAMG
On Thu, Apr 30, 2015 at 5:35 AM, Danyang Su danyang...@gmail.com wrote: On 15-04-29 12:19 PM, Barry Smith wrote: Ok, your code seems fine in terms of logic. But the vector x_vec_gbl i.e. b_react should not be in a read only state. Are you somewhere calling a VecGetArray() and not calling the VecRestoreArray()? Barry I just made a triple check on VecGetArrayF90(). All the VecGetArrayF90() come with VecRestoreArrayF90(). No VecGetArray() is used in my codes. So you can check the lock status of the Vec: http://www.mcs.anl.gov/petsc/petsc-master/docs/manualpages/Vec/VecLockGet.html Put these in your code until you can see who is locking the Vec before the call. Thanks, Matt Danyang On Apr 29, 2015, at 1:52 PM, Danyang Su danyang...@gmail.com wrote: On 15-04-29 11:30 AM, Barry Smith wrote: On Apr 29, 2015, at 12:15 PM, Danyang Su danyang...@gmail.com wrote: On 15-04-28 06:50 PM, Barry Smith wrote: We started enforcing more checks on writing to vectors that you should not write to. Where are you calling DMLocalToGlobalBegin() ? It looks like you are trying to copy values into a global array that you should not because it is a read only input to a function, for example it is the right hand side of a linear system or the input into a SNESFormFunction(). So check the output to that call. Barry This is used when local RHS vec is calculated and copied to global RHS. Huhh, why are you copying anything into the RHS? Is this before you call the linear system solve? Send the code that calls the DMLocalToGlobalBegin() Yes. This is called before the linear system solver. At present, we use PETSc KSP solver. The code is first developed in totally sequential version and then ported to the parallel version. So the subdomain has its own space for jacobi matrix and rhs. We do not use global RHS vector directly to set the RHS at this moment. All the RHS values of subdomain are copied to the global RHS. ! This is the function where RHS values of subdomain are copied to the global RHS values. subroutine compute_function(rank,da,x_array_loc,x_vec_loc, x_vec_gbl,nngl,row_idx_l2pg,col_idx_l2pg, b_non_interlaced) implicit none #include finclude/petscsys.h #include finclude/petscis.h #include finclude/petscvec.h #include finclude/petscvec.h90 #include finclude/petscmat.h #include finclude/petscpc.h #include finclude/petscksp.h #include finclude/petscdmda.h #include finclude/petscdmda.h90 PetscInt :: rank DM :: da PetscReal, allocatable :: x_array_loc(:) Vec :: x_vec_loc Vec :: x_vec_gbl PetscInt :: nngl PetscInt, allocatable :: row_idx_l2pg(:) PetscInt, allocatable :: col_idx_l2pg(:) PetscBool :: b_non_interlaced PetscInt :: info_debug PetscInt :: i, j PetscErrorCode :: ierr PetscScalar, pointer :: vecpointer(:) !Zero entries call VecZeroEntries(x_vec_loc, ierr) !Get a pointer to vector data when you need access to the array call VecGetArrayF90(x_vec_loc, vecpointer, ierr) !Compute the function over the locally owned part of the grid if(b_non_interlaced) then j = nngl/2 do i = 1, j vecpointer(2*i-1) = x_array_loc(i) vecpointer(2*i) = x_array_loc(i+j) end do else do i = 1, nngl vecpointer(i) = x_array_loc(i) end do end if !Restore the vector when you no longer need access to the array call VecRestoreArrayF90(x_vec_loc,vecpointer,ierr) !Insert values into global vector call DMLocalToGlobalBegin(da,x_vec_loc,INSERT_VALUES, x_vec_gbl,ierr) !By placing code between these two statements, computations can be !done while messages are in transition. call DMLocalToGlobalEnd(da,x_vec_loc,INSERT_VALUES, x_vec_gbl,ierr) return end subroutine ! This is the function where reactive transport equations are solved. subroutine solver_dd_snes_solve_react(ilog,idetail,a_in,b_in, x_inout,ia_in,ja_in,nngl_in,itsolv, over_flow,rnorm,row_idx_l2pg,col_idx_l2pg, b_non_interlaced) use gen, only : rank, node_idx_l2lg, ittot_rt, b_output_matrix, b_enable_output use solver_snes_function, only : form_initial_guess,
Re: [petsc-users] Discrete adjoint and adaptive time stepping
On Apr 29, 2015, at 1:51 PM, Jed Brown j...@jedbrown.org wrote: Barry Smith bsm...@mcs.anl.gov writes: If so, how is it done? We just keep the history of the time-step sizes and then use those time-step sizes when doing the backward integration. Seems simple to me, am I missing something? Barry, if you're using this for optimization, you might want the gradient to be exactly consistent with the objective functional. But for that, you would need to differentiate the controller, which is non-smooth in practice because the number of time steps can change and stages could be rejected (solver failure). Ahh, yes for multiple forward runs yup. One approach would be to save the timestep sequence and have the controller use that in subsequent *forward* runs. If the dynamical system behaves similarly for those steps, it would be okay to use the same timestep sequence. Presumably if that single set of dt (from the first run) is not sufficient for some later runs one could possibly use the union of the dt of several runs for all the runs. (that is run adaptively and inconsistently several runs to determine where dt needs to be controlled and then use the various smaller of the dt at the different time regions for a full set of consistent runs). Of course if the various smaller of the dt requires a tiny dt for all time steps then you are not getting an advantage of adaptive time-stepping, but ok. The idea of actually propagating the gradients through the time-step controller seems IMHO to be absurd; I won't even put it on our game plan until we have many more things done and much more practical experience. Barry
Re: [petsc-users] Floating point exception in hypre BoomerAMG
On 15-04-29 12:19 PM, Barry Smith wrote: Ok, your code seems fine in terms of logic. But the vector x_vec_gbl i.e. b_react should not be in a read only state. Are you somewhere calling a VecGetArray() and not calling the VecRestoreArray()? Barry I just made a triple check on VecGetArrayF90(). All the VecGetArrayF90() come with VecRestoreArrayF90(). No VecGetArray() is used in my codes. Danyang On Apr 29, 2015, at 1:52 PM, Danyang Su danyang...@gmail.com wrote: On 15-04-29 11:30 AM, Barry Smith wrote: On Apr 29, 2015, at 12:15 PM, Danyang Su danyang...@gmail.com wrote: On 15-04-28 06:50 PM, Barry Smith wrote: We started enforcing more checks on writing to vectors that you should not write to. Where are you calling DMLocalToGlobalBegin() ? It looks like you are trying to copy values into a global array that you should not because it is a read only input to a function, for example it is the right hand side of a linear system or the input into a SNESFormFunction(). So check the output to that call. Barry This is used when local RHS vec is calculated and copied to global RHS. Huhh, why are you copying anything into the RHS? Is this before you call the linear system solve? Send the code that calls the DMLocalToGlobalBegin() Yes. This is called before the linear system solver. At present, we use PETSc KSP solver. The code is first developed in totally sequential version and then ported to the parallel version. So the subdomain has its own space for jacobi matrix and rhs. We do not use global RHS vector directly to set the RHS at this moment. All the RHS values of subdomain are copied to the global RHS. ! This is the function where RHS values of subdomain are copied to the global RHS values. subroutine compute_function(rank,da,x_array_loc,x_vec_loc, x_vec_gbl,nngl,row_idx_l2pg,col_idx_l2pg, b_non_interlaced) implicit none #include finclude/petscsys.h #include finclude/petscis.h #include finclude/petscvec.h #include finclude/petscvec.h90 #include finclude/petscmat.h #include finclude/petscpc.h #include finclude/petscksp.h #include finclude/petscdmda.h #include finclude/petscdmda.h90 PetscInt :: rank DM :: da PetscReal, allocatable :: x_array_loc(:) Vec :: x_vec_loc Vec :: x_vec_gbl PetscInt :: nngl PetscInt, allocatable :: row_idx_l2pg(:) PetscInt, allocatable :: col_idx_l2pg(:) PetscBool :: b_non_interlaced PetscInt :: info_debug PetscInt :: i, j PetscErrorCode :: ierr PetscScalar, pointer :: vecpointer(:) !Zero entries call VecZeroEntries(x_vec_loc, ierr) !Get a pointer to vector data when you need access to the array call VecGetArrayF90(x_vec_loc, vecpointer, ierr) !Compute the function over the locally owned part of the grid if(b_non_interlaced) then j = nngl/2 do i = 1, j vecpointer(2*i-1) = x_array_loc(i) vecpointer(2*i) = x_array_loc(i+j) end do else do i = 1, nngl vecpointer(i) = x_array_loc(i) end do end if !Restore the vector when you no longer need access to the array call VecRestoreArrayF90(x_vec_loc,vecpointer,ierr) !Insert values into global vector call DMLocalToGlobalBegin(da,x_vec_loc,INSERT_VALUES, x_vec_gbl,ierr) !By placing code between these two statements, computations can be !done while messages are in transition. call DMLocalToGlobalEnd(da,x_vec_loc,INSERT_VALUES, x_vec_gbl,ierr) return end subroutine ! This is the function where reactive transport equations are solved. subroutine solver_dd_snes_solve_react(ilog,idetail,a_in,b_in, x_inout,ia_in,ja_in,nngl_in,itsolv, over_flow,rnorm,row_idx_l2pg,col_idx_l2pg, b_non_interlaced) use gen, only : rank, node_idx_l2lg, ittot_rt, b_output_matrix, b_enable_output use solver_snes_function, only : form_initial_guess, compute_function, compute_jacobian use petsc_mpi_common, only : petsc_mpi_finalize implicit none #include finclude/petscsys.h #include finclude/petscvec.h #include finclude/petscvec.h90 #include finclude/petscdmda.h PetscInt ::
Re: [petsc-users] Turning off TSADAPT still adapts time step
To continue beating this horse, I wrote this:
PetscErrorCode CheckStageTrue(TSAdapt adapt,TS ts,PetscBool *accept){
*accept = PETSC_TRUE;
}
TSAdaptSetCheckStage(adapt,CheckStageTrue);
But TSAdapt then still rejects the stage, but quits altogether instead
of continuing:
Nonlinear solve did not converge due to DIVERGED_LINEAR_SOLVE
iterations 4
TSAdapt 'basic': step 0 stage rejected t=0 + 1.000e-01
retrying with dt=2.500e-02
Am I not using TSAdaptSetCheckStage correctly?
On 04/27/2015 12:55 PM, Jed Brown wrote:
Mark Lohry mlo...@princeton.edu writes:
TSAdaptChoose_None is being called, yes. Output with -ts_adapt monitor
looks like this:
TSAdapt 'none': step 17 accepted t=680+ 4.000e+01
family='beuler' scheme=0:'(null)' dt=4.000e+01
TSAdapt 'none': step 18 accepted t=720+ 4.000e+01
family='beuler' scheme=0:'(null)' dt=4.000e+01
TSAdapt 'none': step 19 accepted t=760+ 4.000e+01
family='beuler' scheme=0:'(null)' dt=4.000e+01
TSAdapt 'none': step 20 stage rejected t=800+ 4.000e+01
retrying with dt=1.000e+01
Add -snes_converged_reason. What do you want to do when the solver
fails? Pretend like it succeeded and get the wrong answer? Usually
people shorten the time step and retry, which is what you see happening
here.
You can call TSAdaptSetCheckStage and have your function unconditionally
set accept to PETSC_TRUE (and you'll likely get the wrong answer).
Re: [petsc-users] Discrete adjoint and adaptive time stepping
On Apr 29, 2015, at 12:57 PM, Miguel Angel Salazar de Troya salazardetr...@gmail.com wrote: Hi I was wondering if your implementation to calculate sensitivities with discrete adjoints will account for the adaptive time stepping. The time step history depends on the parameters, so it should be included in the derivation of the discrete adjoint. Is this calculated in your implementation? Sure, If so, how is it done? We just keep the history of the time-step sizes and then use those time-step sizes when doing the backward integration. Seems simple to me, am I missing something? Barry Thanks in advance Miguel -- Miguel Angel Salazar de Troya Graduate Research Assistant Department of Mechanical Science and Engineering University of Illinois at Urbana-Champaign (217) 550-2360 salaz...@illinois.edu
Re: [petsc-users] Discrete adjoint and adaptive time stepping
Barry Smith bsm...@mcs.anl.gov writes: If so, how is it done? We just keep the history of the time-step sizes and then use those time-step sizes when doing the backward integration. Seems simple to me, am I missing something? Barry, if you're using this for optimization, you might want the gradient to be exactly consistent with the objective functional. But for that, you would need to differentiate the controller, which is non-smooth in practice because the number of time steps can change and stages could be rejected (solver failure). One approach would be to save the timestep sequence and have the controller use that in subsequent *forward* runs. If the dynamical system behaves similarly for those steps, it would be okay to use the same timestep sequence. signature.asc Description: PGP signature
[petsc-users] Discrete adjoint and adaptive time stepping
Hi I was wondering if your implementation to calculate sensitivities with discrete adjoints will account for the adaptive time stepping. The time step history depends on the parameters, so it should be included in the derivation of the discrete adjoint. Is this calculated in your implementation? If so, how is it done? Thanks in advance Miguel -- *Miguel Angel Salazar de Troya* Graduate Research Assistant Department of Mechanical Science and Engineering University of Illinois at Urbana-Champaign (217) 550-2360 salaz...@illinois.edu
Re: [petsc-users] Floating point exception in hypre BoomerAMG
On Apr 29, 2015, at 12:15 PM, Danyang Su danyang...@gmail.com wrote:
On 15-04-28 06:50 PM, Barry Smith wrote:
We started enforcing more checks on writing to vectors that you should not
write to. Where are you calling DMLocalToGlobalBegin() ? It looks like you
are trying to copy values into a global array that you should not because it
is a read only input to a function, for example it is the right hand side of
a linear system or the input into a SNESFormFunction(). So check the output
to that call.
Barry
This is used when local RHS vec is calculated and copied to global RHS.
Huhh, why are you copying anything into the RHS? Is this before you call the
linear system solve? Send the code that calls the DMLocalToGlobalBegin()
Barry
It is strange why this error does not occur at the first timestep.
Thanks,
Danyang
On Apr 28, 2015, at 7:30 PM, Danyang Su danyang...@gmail.com wrote:
Hi Barry,
There seems another bug (not pretty sure) in PETSc-dev, as shown below. The
case I used is similar with the one I mentioned recently. I have no problem
running this case using PETSc 3.5.2. But it give out the following error
using PETSc-dev, even with only one processor.
Thanks,
Danyang
timestep:2048 time: 4.615E+00 years delt: 1.460E-20 years iter: 1
max.sia: 0.000E+00 tol.sia: 0.000E+00
[0]PETSC ERROR: - Error Message
--
[0]PETSC ERROR: Object is in wrong state
[0]PETSC ERROR: Vec is locked read only, argument # 1
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for
trouble shooting.
[0]PETSC ERROR: Petsc Development GIT revision: v3.5.3-2796-g23b6c49 GIT
Date: 2015-04-27 11:32:44 -0500
[0]PETSC ERROR: ../min3p_thcm on a linux-gnu-dbg named nwmop by dsu Tue Apr
28 15:56:41 2015
[0]PETSC ERROR: Configure options PETSC_ARCH=linux-gnu-dbg --with-cc=gcc
--with-cxx=g++ --with-fc=gfortran --download-fblaslapack --download-mpich
--download-mumps --download-hypre --download-superlu_dist --download-metis
--download-parmetis --download-scalapack
[0]PETSC ERROR: #349 VecGetArray() line 1646 in
/home/dsu/Soft/PETSc/petsc-dev/src/vec/vec/interface/rvector.c
[0]PETSC ERROR: #350 VecGetArrayPair() line 434 in
/home/dsu/Soft/PETSc/petsc-dev/include/petscvec.h
[0]PETSC ERROR: #351 VecScatterBegin_SSToSS() line 649 in
/home/dsu/Soft/PETSc/petsc-dev/src/vec/vec/utils/vscat.c
[0]PETSC ERROR: #352 VecScatterBegin() line 1694 in
/home/dsu/Soft/PETSc/petsc-dev/src/vec/vec/utils/vscat.c
[0]PETSC ERROR: #353 DMLocalToGlobalBegin_DA() line 56 in
/home/dsu/Soft/PETSc/petsc-dev/src/dm/impls/da/dagtol.c
[0]PETSC ERROR: #354 DMLocalToGlobalBegin() line 1944 in
/home/dsu/Soft/PETSc/petsc-dev/src/dm/interface/dm.c
Reduce time step for reactive transport
no further time step reduction possible
Optimal time increment computed by MIN3P 2.6644769214899893E-018
Minimum time increment specified by user 3.6498E-018
Please, reduce the MINIMUM TIME INCREMENT
stop signal in time step reduction failed
On 15-04-28 11:15 AM, Barry Smith wrote:
I am forwarding this on to the hypre support email list; hopefully they
can have some advice on how to proceed.
hypre folks, this is with version hypre-2.10.0b (the previous version
had the same behavior). We are getting a divide by zero in gselim() line
4363 (this happens in a time dependent problem after many successful
solves with time dependent matrix). I looked at the code, below, and note
that there are some checks for the diagonal of A not being zero but not
for the line that causes the divide by zero; is this perhaps an oversight
in the hypre code? Any advice appreciated.
Thanks
Barry
HYPRE_Int gselim(A,x,n)
HYPRE_Real *A;
HYPRE_Real *x;
HYPRE_Int n;
{
HYPRE_Interr_flag = 0;
HYPRE_Intj,k,m;
HYPRE_Real factor;
if (n==1) /* A is 1x1 */
{
if (A[0] != 0.0)
{
x[0] = x[0]/A[0];
return(err_flag);
}
else
{
err_flag = 1;
return(err_flag);
}
}
else /* A is nxn. Forward elimination */
{
for (k = 0; k n-1; k++)
{
if (A[k*n+k] != 0.0)
{
for (j = k+1; j n; j++)
{
if (A[j*n+k] != 0.0)
{
factor = A[j*n+k]/A[k*n+k];
for (m = k+1; m n; m++)
{
A[j*n+m] -= factor * A[k*n+m];
}
/* Elimination step for rhs */
x[j] -= factor * x[k];
}
}
}
}
/* Back Substitution */
for (k = n-1; k 0; --k)
{
Re: [petsc-users] Tao iterations
Okay that's what I figured, thanks you very much On Wed, Apr 29, 2015 at 10:39 AM, Jason Sarich jason.sar...@gmail.com wrote: Hi Justin, This expected behavior due to the accumulation of numerical round-offs. If this is a problem or if you just want to confirm that this is the cause, you can try configuring PETSc for quad precision (--with-precision=__float128, works with GNU compilers) and the results should match better. Jason On Tue, Apr 28, 2015 at 10:19 PM, Justin Chang jchan...@uh.edu wrote: Jason (or anyone), I am noticing that the iteration numbers reported by TaoGetSolutionStatus() for blmvm differ whenever I change the number of processes. The solution seems to remain the same though. Is there a reason why this could be happening? Thanks, On Tue, Apr 21, 2015 at 10:40 AM, Jason Sarich jason.sar...@gmail.com wrote: Justin, 1) The big difference between TRON and BLMVM is that TRON requires hessian information, BLMVM only uses gradient information. Thus TRON will usually converge faster, but requires more information, memory, and a KSP solver. GPCG (gradient projected conjugate gradient) is another gradient-only option, but usually performs worse than BLMVM. 2) TaoGetLinearSolveIterations() will get the total number of KSP iterations per solve Jason On Tue, Apr 21, 2015 at 10:33 AM, Justin Chang jychan...@gmail.com wrote: Jason, Tightening the tolerances did the trick. Thanks. Though I do have a couple more related questions: 1) Is there a general guideline for choosing tron over blmvm or vice versa? Also is there another tao type that is also suitable given only bounded constraints? 2) Is it possible to obtain the total number of KSP and/or PG iterations from tron? Thanks, Justin On Tue, Apr 21, 2015 at 9:52 AM, Jason Sarich jason.sar...@gmail.com wrote: Hi Justin, blmvm believes that it is already sufficiently close to a minimum, so it doesn't do anything. You may need to tighten some of the tolerance to force an iteration. Jason On Tue, Apr 21, 2015 at 9:48 AM, Justin Chang jychan...@gmail.com wrote: Time step 1: Tao Object: 1 MPI processes type: blmvm Gradient steps: 0 TaoLineSearch Object: 1 MPI processes type: more-thuente Active Set subset type: subvec convergence tolerances: fatol=0.0001, frtol=0.0001 convergence tolerances: gatol=0, steptol=0, gttol=0 Residual in Function/Gradient:=0.0663148 Objective value=-55.5945 total number of iterations=35, (max: 2000) total number of function/gradient evaluations=37, (max: 4000) Solution converged: estimated |f(x)-f(X*)|/|f(X*)| = frtol Time step 2: Tao Object: 1 MPI processes type: blmvm Gradient steps: 0 TaoLineSearch Object: 1 MPI processes type: more-thuente Active Set subset type: subvec convergence tolerances: fatol=0.0001, frtol=0.0001 convergence tolerances: gatol=0, steptol=0, gttol=0 Residual in Function/Gradient:=0.0682307 Objective value=-66.9675 total number of iterations=23, (max: 2000) total number of function/gradient evaluations=25, (max: 4000) Solution converged: estimated |f(x)-f(X*)|/|f(X*)| = frtol Time step 3: Tao Object: 1 MPI processes type: blmvm Gradient steps: 0 TaoLineSearch Object: 1 MPI processes type: more-thuente Active Set subset type: subvec convergence tolerances: fatol=0.0001, frtol=0.0001 convergence tolerances: gatol=0, steptol=0, gttol=0 Residual in Function/Gradient:=0.0680522 Objective value=-71.8211 total number of iterations=19, (max: 2000) total number of function/gradient evaluations=22, (max: 4000) Solution converged: estimated |f(x)-f(X*)|/|f(X*)| = frtol Time step 4: Tao Object: 1 MPI processes type: blmvm Gradient steps: 0 TaoLineSearch Object: 1 MPI processes type: more-thuente Active Set subset type: subvec convergence tolerances: fatol=0.0001, frtol=0.0001 convergence tolerances: gatol=0, steptol=0, gttol=0 Residual in Function/Gradient:=0.0551556 Objective value=-75.1252 total number of iterations=18, (max: 2000) total number of function/gradient evaluations=20, (max: 4000) Solution converged: estimated |f(x)-f(X*)|/|f(X*)| = frtol Time step 5: Tao Object: 1 MPI processes type: blmvm Gradient steps: 0 TaoLineSearch Object: 1 MPI processes type: more-thuente Active Set subset type: subvec convergence tolerances: fatol=0.0001, frtol=0.0001 convergence tolerances: gatol=0, steptol=0, gttol=0 Residual in Function/Gradient:=0.0675667 Objective value=-77.4414 total number of iterations=6, (max: 2000) total number of function/gradient evaluations=8, (max: 4000) Solution converged:
Re: [petsc-users] Floating point exception in hypre BoomerAMG
Ok, it sounds to me like you are not checking the error codes and so when hypre fails your code continues along as if nothing happened. Thus the lock set in KSPSolve() remains in place and hits your code later when you try to write to those locations. You really should be checking error codes and doing something different if an error occurred and not just keep on running. You Barry On Apr 29, 2015, at 5:19 PM, Danyang Su danyang...@gmail.com wrote: Hi Barry, I haven't tried your suggestion but it seems the bug from HYPRE. Here are my tests: Test A: If I run with start_in_debugger, I got PETSC ERROR: PCApply_HYPRE(), as shown below. If I run without start_in_debugger, then I got Vec lock problem as I mentioned before, exactly at the same timestep. timestep:2048 time: 4.619E+00 years delt: 4.208E-03 years iter: 1 max.sia: 0.000E+00 tol.sia: 0.000E+00 Newton Iteration Convergence Summary: Newton maximumfor maximumsolver iteration update component residual iterations maxvol nexvol 1 3.D+00 o2(aq) 2.4957D-03 1 33 39 2 3.D+00 o2(aq) 6.4858D-03 1 33 39 3 3.D+00 o2(aq) 5.7149D-03 1 33 39 4 3.D+00 o2(aq) 4.9772D-03 1 33 39 5 -3.D+00 co3-25.9683D-06 180 33 39 6 3.D+00 h+1 2.5193D+04 1 1 39 7 3.D+00 co3-22.4783D+16 3 1 39 8 3.D+00 h+1 3.6929D+29 14 1 39 [0]PETSC ERROR: PCApply_HYPRE() line 174 in /home/dsu/Soft/PETSc/petsc-3.5.2/src/ksp/pc/impls/hypre/hypre.c Error in HYPRE solver, error code 1 Test B: I changed the preconditioner to bjacobi then all the errors are gone, whatever with or without start_in_debugger option. Usually we use HYPRE as preconditioner method as it is generally faster than other methods for our cases. Do I still need to test on this case to locate the error or just wait until the bug of HYPRE is fixed. Thanks, Danyang On 15-04-29 12:44 PM, Barry Smith wrote: Ok, run on one process in the debugger and put a break point in VecLockPush and VecLockPopthen run and each time it breaks on VecLockPush or VecLockPop type bt then type cont send all the output, this will tell us where the vector got locked readonly and did not get unlocked. Barry On Apr 29, 2015, at 2:35 PM, Danyang Su danyang...@gmail.com wrote: On 15-04-29 12:19 PM, Barry Smith wrote: Ok, your code seems fine in terms of logic. But the vector x_vec_gbl i.e. b_react should not be in a read only state. Are you somewhere calling a VecGetArray() and not calling the VecRestoreArray()? Barry I just made a triple check on VecGetArrayF90(). All the VecGetArrayF90() come with VecRestoreArrayF90(). No VecGetArray() is used in my codes. Danyang On Apr 29, 2015, at 1:52 PM, Danyang Su danyang...@gmail.com wrote: On 15-04-29 11:30 AM, Barry Smith wrote: On Apr 29, 2015, at 12:15 PM, Danyang Su danyang...@gmail.com wrote: On 15-04-28 06:50 PM, Barry Smith wrote: We started enforcing more checks on writing to vectors that you should not write to. Where are you calling DMLocalToGlobalBegin() ? It looks like you are trying to copy values into a global array that you should not because it is a read only input to a function, for example it is the right hand side of a linear system or the input into a SNESFormFunction(). So check the output to that call. Barry This is used when local RHS vec is calculated and copied to global RHS. Huhh, why are you copying anything into the RHS? Is this before you call the linear system solve? Send the code that calls the DMLocalToGlobalBegin() Yes. This is called before the linear system solver. At present, we use PETSc KSP solver. The code is first developed in totally sequential version and then ported to the parallel version. So the subdomain has its own space for jacobi matrix and rhs. We do not use global RHS vector directly to set the RHS at this moment. All the RHS values of subdomain are copied to the global RHS. ! This is the function where RHS values of subdomain are copied to the global RHS values. subroutine compute_function(rank,da,x_array_loc,x_vec_loc, x_vec_gbl,nngl,row_idx_l2pg,col_idx_l2pg, b_non_interlaced) implicit none #include finclude/petscsys.h #include finclude/petscis.h #include finclude/petscvec.h #include finclude/petscvec.h90 #include finclude/petscmat.h #include finclude/petscpc.h #include finclude/petscksp.h #include finclude/petscdmda.h #include finclude/petscdmda.h90 PetscInt :: rank DM :: da PetscReal, allocatable :: x_array_loc(:) Vec :: x_vec_loc
Re: [petsc-users] Floating point exception in hypre BoomerAMG
Hi Barry, I haven't tried your suggestion but it seems the bug from HYPRE. Here are my tests: Test A: If I run with start_in_debugger, I got PETSC ERROR: PCApply_HYPRE(), as shown below. If I run without start_in_debugger, then I got Vec lock problem as I mentioned before, exactly at the same timestep. timestep:2048 time: 4.619E+00 years delt: 4.208E-03 years iter: 1 max.sia: 0.000E+00 tol.sia: 0.000E+00 Newton Iteration Convergence Summary: Newton maximumfor maximumsolver iteration update component residual iterations maxvol nexvol 1 3.D+00 o2(aq) 2.4957D-03 1 33 39 2 3.D+00 o2(aq) 6.4858D-03 1 33 39 3 3.D+00 o2(aq) 5.7149D-03 1 33 39 4 3.D+00 o2(aq) 4.9772D-03 1 33 39 5 -3.D+00 co3-25.9683D-06 180 33 39 6 3.D+00 h+1 2.5193D+04 1 1 39 7 3.D+00 co3-22.4783D+16 3 1 39 8 3.D+00 h+1 3.6929D+29 14 1 39 [0]PETSC ERROR: PCApply_HYPRE() line 174 in /home/dsu/Soft/PETSc/petsc-3.5.2/src/ksp/pc/impls/hypre/hypre.c Error in HYPRE solver, error code 1 Test B: I changed the preconditioner to bjacobi then all the errors are gone, whatever with or without start_in_debugger option. Usually we use HYPRE as preconditioner method as it is generally faster than other methods for our cases. Do I still need to test on this case to locate the error or just wait until the bug of HYPRE is fixed. Thanks, Danyang On 15-04-29 12:44 PM, Barry Smith wrote: Ok, run on one process in the debugger and put a break point in VecLockPush and VecLockPopthen run and each time it breaks on VecLockPush or VecLockPop type bt then type cont send all the output, this will tell us where the vector got locked readonly and did not get unlocked. Barry On Apr 29, 2015, at 2:35 PM, Danyang Su danyang...@gmail.com wrote: On 15-04-29 12:19 PM, Barry Smith wrote: Ok, your code seems fine in terms of logic. But the vector x_vec_gbl i.e. b_react should not be in a read only state. Are you somewhere calling a VecGetArray() and not calling the VecRestoreArray()? Barry I just made a triple check on VecGetArrayF90(). All the VecGetArrayF90() come with VecRestoreArrayF90(). No VecGetArray() is used in my codes. Danyang On Apr 29, 2015, at 1:52 PM, Danyang Su danyang...@gmail.com wrote: On 15-04-29 11:30 AM, Barry Smith wrote: On Apr 29, 2015, at 12:15 PM, Danyang Su danyang...@gmail.com wrote: On 15-04-28 06:50 PM, Barry Smith wrote: We started enforcing more checks on writing to vectors that you should not write to. Where are you calling DMLocalToGlobalBegin() ? It looks like you are trying to copy values into a global array that you should not because it is a read only input to a function, for example it is the right hand side of a linear system or the input into a SNESFormFunction(). So check the output to that call. Barry This is used when local RHS vec is calculated and copied to global RHS. Huhh, why are you copying anything into the RHS? Is this before you call the linear system solve? Send the code that calls the DMLocalToGlobalBegin() Yes. This is called before the linear system solver. At present, we use PETSc KSP solver. The code is first developed in totally sequential version and then ported to the parallel version. So the subdomain has its own space for jacobi matrix and rhs. We do not use global RHS vector directly to set the RHS at this moment. All the RHS values of subdomain are copied to the global RHS. ! This is the function where RHS values of subdomain are copied to the global RHS values. subroutine compute_function(rank,da,x_array_loc,x_vec_loc, x_vec_gbl,nngl,row_idx_l2pg,col_idx_l2pg, b_non_interlaced) implicit none #include finclude/petscsys.h #include finclude/petscis.h #include finclude/petscvec.h #include finclude/petscvec.h90 #include finclude/petscmat.h #include finclude/petscpc.h #include finclude/petscksp.h #include finclude/petscdmda.h #include finclude/petscdmda.h90 PetscInt :: rank DM :: da PetscReal, allocatable :: x_array_loc(:) Vec :: x_vec_loc Vec :: x_vec_gbl PetscInt :: nngl PetscInt, allocatable :: row_idx_l2pg(:) PetscInt, allocatable :: col_idx_l2pg(:) PetscBool :: b_non_interlaced PetscInt :: info_debug PetscInt :: i, j PetscErrorCode :: ierr PetscScalar, pointer :: vecpointer(:) !Zero entries call VecZeroEntries(x_vec_loc, ierr) !Get a pointer to vector data when you need access to the array call
