Hello everyone,
I hope this email finds you well.
My Name is Sawsan Shatanawi, and I was developing a Fortran code for
simulating groundwater flow in a 3D system with nonlinear behavior. I solved
the nonlinear system using the PCG solver and Picard iteration, but I did not
get good results although I checked my matrix and RHS and everything, I decided
to change my solver to Newton Rapson method.
I checked PETSc documents but I have a few questions:
1) My groundwater system is time-dependent, so should I use TS only instead of
SNES?
2) My system has its deltaT, would using deltaT as dt affect my solver, or is
it better to use TS-PETSc dt? Also, would using PETSc dt affect the simulation
of the groundwater system
3) I want my Jacobian matrix to be calculated by PETSc automatically
4) Do I need to define and calculate the residual vector?
My A-Matrix contains coefficients and external sources and my RHS vector
includes the boundary conditions
Please find the attached file contains a draft of my code
Thank you in advance for your time and help.
Best regards,
Sawsan
________________________________
From: Shatanawi, Sawsan Muhammad <[email protected]>
Sent: Tuesday, January 16, 2024 10:43 AM
To: Junchao Zhang <[email protected]>
Cc: Barry Smith <[email protected]>; Matthew Knepley <[email protected]>; Mark
Adams <[email protected]>; [email protected] <[email protected]>
Subject: Re: [petsc-users] Help with Integrating PETSc into Fortran Groundwater
Flow Simulation Code
Hello all,
Thank you for your valuable help. I will do your recommendations and hope it
will run without any issues.
Bests,
Sawsan
________________________________
From: Junchao Zhang <[email protected]>
Sent: Friday, January 12, 2024 8:46 AM
To: Shatanawi, Sawsan Muhammad <[email protected]>
Cc: Barry Smith <[email protected]>; Matthew Knepley <[email protected]>; Mark
Adams <[email protected]>; [email protected] <[email protected]>
Subject: Re: [petsc-users] Help with Integrating PETSc into Fortran Groundwater
Flow Simulation Code
[EXTERNAL EMAIL]
Hi, Sawsan,
First in test_main.F90, you need to call VecGetArrayF90(temp_solution,
H_vector, ierr) and VecRestoreArrayF90 (temp_solution, H_vector, ierr) as
Barry mentioned.
Secondly, in the loop of test_main.F90, it calls GW_solver(). Within it, it
calls PetscInitialize()/PetscFinalize(). But without MPI being initialized,
PetscInitialize()/PetscFinalize() can only be called once.
do timestep =2 , NTSP
call GW_boundary_conditions(timestep-1)
!print *,HNEW(1,1,1)
call GW_elevation()
! print *, GWTOP(2,2,2)
call GW_conductance()
! print *, CC(2,2,2)
call GW_recharge()
! print *, B_Rech(5,4)
call GW_pumping(timestep-1)
! print *, B_pump(2,2,2)
call GW_SW(timestep-1)
print *,B_RIVER (2,2,2)
call GW_solver(timestep-1,N)
call GW_deallocate_loop()
end do
A solution is to delete PetscInitialize()/PetscFinalize() in GW_solver_try.F90
and move it to test_main.F90, outside the do loop.
diff --git a/test_main.F90 b/test_main.F90
index b5997c55..107bd3ee 100644
--- a/test_main.F90
+++ b/test_main.F90
@@ -1,5 +1,6 @@
program test_GW
+#include <petsc/finclude/petsc.h>
use petsc
use GW_constants
use GW_param_by_user
@@ -8,6 +9,9 @@ program test_GW
implicit none
integer :: N
integer :: timestep
+ PetscErrorCode ierr
+
+ call PetscInitialize(ierr)
call GW_domain(N)
!print *, "N=",N
!print *, DELTAT
@@ -37,4 +41,5 @@ program test_GW
end do
print *, HNEW(NCOL,3,2)
call GW_deallocate ()
+ call PetscFinalize(ierr)
end program test_GW
With that, the MPI error will be fixed. The code could run to gw_deallocate ()
before abort. There are other memory errors. You can install/use valgrind to
fix them. Run it with valgrind ./GW.exe and look through the output
Thanks.
--Junchao Zhang
On Thu, Jan 11, 2024 at 10:49 PM Shatanawi, Sawsan Muhammad
<[email protected]<mailto:[email protected]>> wrote:
Hello,
Thank you all for your help.
I have changed VecGetArray to VecGetArrayF90, and the location of destory call.
but I want to make sure that VecGet ArrayF90 is to make a new array( vector)
that I can use in the rest of my Fortran code?
when I run it and debugged it, I got
5.2000000E-03
50.00000
10.00000
0.0000000E+00
PETSC: Attaching gdb to
/weka/data/lab/richey/sawsan/GW_CODE/code2024/SS_GWM/./GW.exe of pid 33065 on
display :0.0 on machine sn16
Unable to start debugger in xterm: No such file or directory
0.0000000E+00
Attempting to use an MPI routine after finalizing MPICH
srun: error: sn16: task 0: Exited with exit code 1
[sawsan.shatanawi@login-p2n02 SS_GWM]$ gdb ./GW/exe
GNU gdb (GDB) Red Hat Enterprise Linux 7.6.1-100.el7
Copyright (C) 2013 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later
<https://urldefense.us/v3/__http://gnu.org/licenses/gpl.html__;!!G_uCfscf7eWS!aDPVIaXxxRNvtadvU_F4xbshLyPmphwU9eHojKr-5XWWpGF7rqbuu-EKWjZxPfCL5B090Aolufvvd329KjQF3rYmFLVEpUkn6Q$
<https://urldefense.com/v3/__http://gnu.org/licenses/gpl.html__;!!JmPEgBY0HMszNaDT!o4qS1zaFLg2L8PlawWJVyYsJpnwYHuL6SIIZsmzTO98RCzbP26HwTi_0-ipCS_D2SBG0X4gtEnM13-nbzFEFKvtraKaM$>>
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law. Type "show copying"
and "show warranty" for details.
This GDB was configured as "x86_64-redhat-linux-gnu".
For bug reporting instructions, please see:
<https://urldefense.us/v3/__http://www.gnu.org/software/gdb/bugs/__;!!G_uCfscf7eWS!aDPVIaXxxRNvtadvU_F4xbshLyPmphwU9eHojKr-5XWWpGF7rqbuu-EKWjZxPfCL5B090Aolufvvd329KjQF3rYmFLWBCWQ8eA$
<https://urldefense.com/v3/__http://www.gnu.org/software/gdb/bugs/__;!!JmPEgBY0HMszNaDT!o4qS1zaFLg2L8PlawWJVyYsJpnwYHuL6SIIZsmzTO98RCzbP26HwTi_0-ipCS_D2SBG0X4gtEnM13-nbzFEFKjO8jNfs$>>...
./GW/exe: No such file or directory.
(gdb) run
Starting program:
No executable file specified.
Use the "file" or "exec-file" command.
(gdb) bt
No stack.
(gdb)
If the highlighted line is the error, I don't know why when I write gdb , it
does not show me the location of error
The code : sshatanawi/SS_GWM
(github.com)<https://urldefense.com/v3/__https://github.com/sshatanawi/SS_GWM__;!!JmPEgBY0HMszNaDT!o4qS1zaFLg2L8PlawWJVyYsJpnwYHuL6SIIZsmzTO98RCzbP26HwTi_0-ipCS_D2SBG0X4gtEnM13-nbzFEFKp8i33ur$>
I really appreciate your helps
Sawsan
________________________________
From: Barry Smith <[email protected]<mailto:[email protected]>>
Sent: Wednesday, January 10, 2024 5:35 PM
To: Junchao Zhang <[email protected]<mailto:[email protected]>>
Cc: Shatanawi, Sawsan Muhammad
<[email protected]<mailto:[email protected]>>; Mark Adams
<[email protected]<mailto:[email protected]>>;
[email protected]<mailto:[email protected]>
<[email protected]<mailto:[email protected]>>
Subject: Re: [petsc-users] Help with Integrating PETSc into Fortran Groundwater
Flow Simulation Code
[EXTERNAL EMAIL]
On Jan 10, 2024, at 6:49 PM, Junchao Zhang
<[email protected]<mailto:[email protected]>> wrote:
Hi, Sawsan,
I could build your code and I also could gdb it.
$ gdb ./GW.exe
...
$ Thread 1 "GW.exe" received signal SIGSEGV, Segmentation fault.
0x00007ffff1e6d44f in vecgetarray_ (x=0x7fffffffa718, fa=0x0,
ia=0x7fffffffa75c, ierr=0x0) at
/scratch/jczhang/petsc/src/vec/vec/interface/ftn-custom/zvectorf.c:257
257 *ierr = VecGetArray(*x, &lx);
(gdb) bt
#0 0x00007ffff1e6d44f in vecgetarray_ (x=0x7fffffffa718, fa=0x0,
ia=0x7fffffffa75c, ierr=0x0) at
/scratch/jczhang/petsc/src/vec/vec/interface/ftn-custom/zvectorf.c:257
#1 0x000000000040b6e3 in gw_solver (t_s=1.40129846e-45, n=300) at
GW_solver_try.F90:169
#2 0x000000000040c6a8 in test_gw () at test_main.F90:35
ierr=0x0 caused the segfault. See
https://urldefense.us/v3/__https://petsc.org/release/manualpages/Vec/VecGetArray/*vecgetarray__;Iw!!G_uCfscf7eWS!aDPVIaXxxRNvtadvU_F4xbshLyPmphwU9eHojKr-5XWWpGF7rqbuu-EKWjZxPfCL5B090Aolufvvd329KjQF3rYmFLWZhBG_dA$
<https://urldefense.com/v3/__https://petsc.org/release/manualpages/Vec/VecGetArray/*vecgetarray__;Iw!!JmPEgBY0HMszNaDT!tqBApprMfYxwNz4Zvnk8coNE5AeWjA9wSdAM7QJcIIVP1z0VDsVIalo4Sew2b0fW3bZtTAbPh-h0MUsZ9Km12jA$>,
you should use VecGetArrayF90 instead.
BTW, Barry, the code
https://urldefense.us/v3/__https://github.com/sshatanawi/SS_GWM/blob/main/GW_solver_try.F90*L169__;Iw!!G_uCfscf7eWS!aDPVIaXxxRNvtadvU_F4xbshLyPmphwU9eHojKr-5XWWpGF7rqbuu-EKWjZxPfCL5B090Aolufvvd329KjQF3rYmFLWGbsJMyg$
<https://urldefense.com/v3/__https://github.com/sshatanawi/SS_GWM/blob/main/GW_solver_try.F90*L169__;Iw!!JmPEgBY0HMszNaDT!tqBApprMfYxwNz4Zvnk8coNE5AeWjA9wSdAM7QJcIIVP1z0VDsVIalo4Sew2b0fW3bZtTAbPh-h0MUsZh2eAi4o$>
has "call VecGetArray(temp_solution, H_vector, ierr)". I don't find petsc
Fortran examples doing VecGetArray. Do we still support it?
This is not the correct calling sequence for VecGetArray() from Fortran.
Regardless, definitely should not be writing any new code that uses
VecGetArray() from Fortran. Should use VecGetArrayF90().
--Junchao Zhang
On Wed, Jan 10, 2024 at 2:38 PM Shatanawi, Sawsan Muhammad via petsc-users
<[email protected]<mailto:[email protected]>> wrote:
Hello all,
I hope you are doing well.
Generally, I use gdb <the name of my exe.file> to debug the code.
I got the attached error message.
I have tried to add the flag -start_in_debugger in the make file, but it didn't
work, so it seems I was doing it in the wrong way
This is the link for the whole code: sshatanawi/SS_GWM
(github.com)<https://urldefense.com/v3/__https://github.com/sshatanawi/SS_GWM__;!!JmPEgBY0HMszNaDT!tqBApprMfYxwNz4Zvnk8coNE5AeWjA9wSdAM7QJcIIVP1z0VDsVIalo4Sew2b0fW3bZtTAbPh-h0MUsZ1Veab3M$>
[https://urldefense.us/v3/__https://opengraph.githubassets.com/9eb6cd14baf12f04848ed209b6f502415eb531bdd7b3a5f9696af68663b870c0/sshatanawi/SS_GWM__;!!G_uCfscf7eWS!aDPVIaXxxRNvtadvU_F4xbshLyPmphwU9eHojKr-5XWWpGF7rqbuu-EKWjZxPfCL5B090Aolufvvd329KjQF3rYmFLVukOhM2Q$
]<https://urldefense.com/v3/__https://github.com/sshatanawi/SS_GWM__;!!JmPEgBY0HMszNaDT!tqBApprMfYxwNz4Zvnk8coNE5AeWjA9wSdAM7QJcIIVP1z0VDsVIalo4Sew2b0fW3bZtTAbPh-h0MUsZ1Veab3M$>
GitHub -
sshatanawi/SS_GWM<https://urldefense.com/v3/__https://github.com/sshatanawi/SS_GWM__;!!JmPEgBY0HMszNaDT!tqBApprMfYxwNz4Zvnk8coNE5AeWjA9wSdAM7QJcIIVP1z0VDsVIalo4Sew2b0fW3bZtTAbPh-h0MUsZ1Veab3M$>
Contribute to sshatanawi/SS_GWM development by creating an account on GitHub.
github.com<https://urldefense.com/v3/__http://github.com/__;!!JmPEgBY0HMszNaDT!tqBApprMfYxwNz4Zvnk8coNE5AeWjA9wSdAM7QJcIIVP1z0VDsVIalo4Sew2b0fW3bZtTAbPh-h0MUsZ8rcrPiA$>
You can read the description of the code in " Model Desprciption.pdf"
the compiling file is makefile_f90 where you can find the linked code files
I really appreciate your help
Bests,
Sawsan
________________________________
From: Mark Adams <[email protected]<mailto:[email protected]>>
Sent: Friday, January 5, 2024 4:53 AM
To: Shatanawi, Sawsan Muhammad
<[email protected]<mailto:[email protected]>>
Cc: Matthew Knepley <[email protected]<mailto:[email protected]>>;
[email protected]<mailto:[email protected]>
<[email protected]<mailto:[email protected]>>
Subject: Re: [petsc-users] Help with Integrating PETSc into Fortran Groundwater
Flow Simulation Code
[EXTERNAL EMAIL]
This is a segv. As Matt said, you need to use a debugger for this or add print
statements to narrow down the place where this happens.
You will need to learn how to use debuggers to do your project so you might as
well start now.
If you have a machine with a GUI debugger that is easier but command line
debuggers are good to learn anyway.
I tend to run debuggers directly (eg, lldb ./a.out -- program-args ...) and use
a GUI debugger (eg, Totalview or DDT) if available.
Mark
On Wed, Dec 20, 2023 at 10:02 PM Shatanawi, Sawsan Muhammad via petsc-users
<[email protected]<mailto:[email protected]>> wrote:
Hello Matthew,
Thank you for your help. I am sorry that I keep coming back with my error
messages, but I reached a point that I don't know how to fix them, and I don't
understand them easily.
The list of errors is getting shorter, now I am getting the attached error
messages
Thank you again,
Sawsan
________________________________
From: Matthew Knepley <[email protected]<mailto:[email protected]>>
Sent: Wednesday, December 20, 2023 6:54 PM
To: Shatanawi, Sawsan Muhammad
<[email protected]<mailto:[email protected]>>
Cc: Barry Smith <[email protected]<mailto:[email protected]>>;
[email protected]<mailto:[email protected]>
<[email protected]<mailto:[email protected]>>
Subject: Re: [petsc-users] Help with Integrating PETSc into Fortran Groundwater
Flow Simulation Code
[EXTERNAL EMAIL]
On Wed, Dec 20, 2023 at 9:49 PM Shatanawi, Sawsan Muhammad via petsc-users
<[email protected]<mailto:[email protected]>> wrote:
Hello Barry,
Thank you a lot for your help, Now I am getting the attached error message.
Do not destroy the PC from KSPGetPC()
THanks,
Matt
Bests,
Sawsan
________________________________
From: Barry Smith <[email protected]<mailto:[email protected]>>
Sent: Wednesday, December 20, 2023 6:32 PM
To: Shatanawi, Sawsan Muhammad
<[email protected]<mailto:[email protected]>>
Cc: Mark Adams <[email protected]<mailto:[email protected]>>;
[email protected]<mailto:[email protected]>
<[email protected]<mailto:[email protected]>>
Subject: Re: [petsc-users] Help with Integrating PETSc into Fortran Groundwater
Flow Simulation Code
[EXTERNAL EMAIL]
Instead of
call PCCreate(PETSC_COMM_WORLD, pc, ierr)
call PCSetType(pc, PCILU,ierr) ! Choose a preconditioner type (ILU)
call KSPSetPC(ksp, pc,ierr) ! Associate the preconditioner with the KSP
solver
do
call KSPGetPC(ksp,pc,ierr)
call PCSetType(pc, PCILU,ierr)
Do not call KSPSetUp(). It will be taken care of automatically during the solve
On Dec 20, 2023, at 8:52 PM, Shatanawi, Sawsan Muhammad via petsc-users
<[email protected]<mailto:[email protected]>> wrote:
Hello,
I don't think that I set preallocation values when I created the matrix, would
you please have look at my code. It is just the petsc related part from my code.
I was able to fix some of the error messages. Now I have a new set of error
messages related to the KSP solver (attached)
I appreciate your help
Sawsan
________________________________
From: Mark Adams <[email protected]<mailto:[email protected]>>
Sent: Wednesday, December 20, 2023 6:44 AM
To: Shatanawi, Sawsan Muhammad
<[email protected]<mailto:[email protected]>>
Cc: [email protected]<mailto:[email protected]>
<[email protected]<mailto:[email protected]>>
Subject: Re: [petsc-users] Help with Integrating PETSc into Fortran Groundwater
Flow Simulation Code
[EXTERNAL EMAIL]
Did you set preallocation values when you created the matrix?
Don't do that.
On Wed, Dec 20, 2023 at 9:36 AM Shatanawi, Sawsan Muhammad
<[email protected]<mailto:[email protected]>> wrote:
Hello,
I am trying to create a sparse matrix( which is as I believe a zero matrix)
then adding some nonzero elements to it over a loop, then assembling it
Get Outlook for
iOS<https://urldefense.com/v3/__https://aka.ms/o0ukef__;!!JmPEgBY0HMszNaDT!uUJ_jeYf45gcXDGR_PeMjhU7hbd_fKcXJPn0pM9eb-YQihKNYuXMYM9x-hglsbXsCFIwNBWgHXdetHODupsOloE$>
________________________________
From: Mark Adams <[email protected]<mailto:[email protected]>>
Sent: Wednesday, December 20, 2023 2:48 AM
To: Shatanawi, Sawsan Muhammad
<[email protected]<mailto:[email protected]>>
Cc: [email protected]<mailto:[email protected]>
<[email protected]<mailto:[email protected]>>
Subject: Re: [petsc-users] Help with Integrating PETSc into Fortran Groundwater
Flow Simulation Code
[EXTERNAL EMAIL]
I am guessing that you are creating a matrix, adding to it, finalizing it
("assembly"), and then adding to it again, which is fine, but you are adding
new non-zeros to the sparsity pattern.
If this is what you want then you can tell the matrix to let you do that.
Otherwise you have a bug.
Mark
On Tue, Dec 19, 2023 at 9:50 PM Shatanawi, Sawsan Muhammad via petsc-users
<[email protected]<mailto:[email protected]>> wrote:
Hello everyone,
I hope this email finds you well.
My Name is Sawsan Shatanawi, and I am currently working on developing a
Fortran code for simulating groundwater flow in a 3D system. The code involves
solving a nonlinear system, and I have created the matrix to be solved using
the PCG solver and Picard iteration. However, when I tried to assign it as a
PETSc matrix I started getting a lot of error messages.
I am kindly asking if someone can help me, I would be happy to share my code
with him/her.
Please find the attached file contains a list of errors I have gotten
Thank you in advance for your time and assistance.
Best regards,
Sawsan
<Matrix_RHS.F90><out.txt><solver.F90>
--
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
https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!aDPVIaXxxRNvtadvU_F4xbshLyPmphwU9eHojKr-5XWWpGF7rqbuu-EKWjZxPfCL5B090Aolufvvd329KjQF3rYmFLX0Bmayww$
<https://urldefense.com/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!JmPEgBY0HMszNaDT!uskvAyF0pMMWDbMIexr9g4qN46V7Rea17GQdNIVG2vH_HMaX7mXgie4ZYgusmPpss_DS7H1_8vn8arGQNSkC$>
########## AMatrix subroutine ###########################
!This subroutine to set up the coefficient matrix, RHS vector, and the boundary
conditions.
! It contains the vertical flow correction when the confined aquifer acts as
! unconfined aquifer ( The hydraulic head becomes lower than the top of confined
! aquifer)
subroutine GW_AMatrix_RHS_BC(t_s,N,RHS_vector, &
temp_solution,A_Matrix)
!subroutine GW_AMatrix_RHS_BC(t_s,N,RECH,DEEPRECH,RHS_vector, &
! temp_solution,A_Matrix)
#include <petsc/finclude/petsc.h>
#include <petsc/finclude/petscvec.h>
#include <petsc/finclude/petscmat.h>
use petsc
use petscvec
use petscmat
use GW_constants
use GW_param_by_user
use GW_variables
implicit none
integer,intent(in) :: N
integer,intent(in) :: t_s
! real,allocatable,dimension(:,:),intent(in) :: RECH
! real,allocatable,dimension(:,:),intent(in) :: DEEPRECH
Vec, intent(inout) ::RHS_vector
Vec, intent(inout) ::temp_solution
Mat, intent(inout) ::A_Matrix
real,allocatable,dimension(:,:,:) :: HCOF ! Head Coefficent
real,allocatable,dimension(:,:,:) :: B !The constant cofficient from all
exteral sources
real,allocatable,dimension(:,:,:) :: P !Head dependent cofficient from all
external sources
! Declare PETSc objects
PetscErrorCode ierr
PetscInt i, j, k, indexGW
PetscScalar A_value,h,RHS
if (.not. allocated (HCOF))then
allocate(HCOF(NCOL,NROW,NLAY))
HCOF = 0.0
end if
if (.not. allocated(B))then
allocate (B(NCOL,NROW,NLAY))
B = 0.0
end if
if (.not. allocated(P))then
allocate(P(NCOL,NROW,NLAY))
P = 0.0
end if
! Initialize PETSc
call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
! Update the matrix A and vector RHS
do k = 1, NLAY
do i = 1, NROW
do j = 1, NCOL
! assign the old HNEW to HOLD
! HOLD(j,i,k)=HNEW(j,i,k) ! it is already in elevation
! Calculate the global index of the current cell
indexGW= (k-1)*NCOL*NROW+(i-1)*NCOL+j
! Calculate P(j,i,k) and B(j,i,k)
! No need for P_Recharge because it does not depend on
! head
P(j, i, k) = -1.0*P_RIVER(j, i, k) !- P_BOUND(j,i,k)
B(j, i, k) = B_RIVER(j, i, k) + B_Rech(j, i) + B_pump(j,i,k)!
+B_BOUND(j,i,k)
! Calculate HCOF(j,i,k)
if (k == 1) then
HCOF(j, i, k) = P(j, i, k) - (SY(j, i) * DELTAC *DELTAR *
DELTAV(j, i, k) / DELTAT)
else
HCOF(j, i, k) = P(j, i, k) - (SS(j, i) * DELTAC *DELTAR *
DELTAV(j, i, k) / DELTAT)
end if
! update diagonal element of A_Matrix
if(indexGW>=1 .and. indexGW<=N)then
A_value =
-(CR_left(j,i,k)+CR_right(j,i,k)+CC_left(j,i,k)+CC_right
(j,i,k)+CV_left(j,i,k))+HCOF(j,i,k)
call
MatSetValue(A_Matrix,indexGW-1,indexGW-1,A_value,INSERT_VALUES,ierr)
! end if
! Update the neighboring elements of A_Matrix
if (j < NCOL) then
A_value = CR_right(j,i,k)
call
MatSetValue(A_Matrix,indexGW-1,indexGW+1-1,A_value,INSERT_VALUES,ierr)
end if
if (j > 1) then
A_value = CR_left (j,i,k)
call
MatSetValue(A_Matrix,indexGW-1,indexGW-1-1,A_value,INSERT_VALUES,ierr)
end if
if (i < NROW) then
A_value = CC_right(j,i,k)
call
MatSetValue(A_Matrix,indexGW-1,indexGW+NCOL-1,A_value,INSERT_VALUES,ierr)
end if
if (i > 1) then
A_value =CC_right(j,i,k)
call
MatSetValue(A_Matrix,indexGW-1,indexGW-NCOL-1,A_value,INSERT_VALUES,ierr)
end if
if (k < NLAY ) then
A_value = CV_left(j,i,k)
call
MatSetValue(A_Matrix,indexGW-1,indexGW+NCOL*NROW-1,A_value,INSERT_VALUES,ierr)
end if
! if (k > 1) then
! A_value = CV(j, i, k - 1)
! call
MatSetValue(A_Matrix,indexGW-1,indexGW-NCOL*NROW-1,A_value,INSERT_VALUES,ierr)
end if
! Calculate RHS based on dry/wet conditions
! If the confined aquifer is fully saturated (the
! piezometric head higher than the top of the confined
! aquifer)
if (HOLD(j, i, 2) > (GWBOT(j, i, 1) - BED_THK))then
if (k == 1) then
if (GW_CELL(j,i,k) == 0 .or. GW_CELL(j,i,k) ==2 ) then
RHS = - H_constant (j,i,k) !-HNOFLOW (j,i,k)
!elseif (GW_CELL (j,i,k) == 2) then
! RHS = -HBOUND(j,i,k)
else
RHS = -B(j, i, k) - ((SY(j, i) *DELTAC * DELTAR
* DELTAV(j, i, k) * HOLD(j, i, k) / DELTAT))
end if
call VecSetValue(RHS_vector, indexGW-1, RHS, INSERT_VALUES,
ierr)
else
if (GW_CELL(j,i,k) == 0 .or. GW_CELL (j,i,k) == 2 ) then
RHS = - H_constant( j,i,k) !-HNOFLOW (j,i,k)
! elseif (GW_CELL (j,i,k) == 2) then
! RHS = -HBOUND(j,i,k)
else
RHS = -B(j,i,k)
-((SS(j,i)*DELTAC*DELTAR*DELTAV(j,i,k)*HOLD(j,i,k)/DELTAT))
endif
call VecSetValue(RHS_vector, indexGW-1, RHS,
INSERT_VALUES, ierr)
end if
end if
!! if the confined aquifer has dewatered (the peizomteric
!head is below the top of the confined aquifer, start
!behave as unconfined)
if (HOLD (j,i,2) .LT. (GWBOT (j,i,1)-BED_THK)) then
if ( k .EQ. 1 ) then
if (GW_CELL(j,i,k) == 0 .or. GW_CELL(j,i,k) == 2) then
RHS = - H_constant (j,i,k) !-HNOFLOW (j,i,k)
!elseif (GW_CELL (j,i,k) == 2) then
! RHS = -HBOUND(j,i,k)
else
RHS = -B(j,i,k)
-((SY(j,i)*DELTAC*DELTAR*DELTAV(j,i,k)*HOLD(j,i,k)/DELTAT))+
CV_left(j,i,k)*(HOLD (j,i,k+1)-GWTOP(j,i,k+1))
call VecSetValue(RHS_vector, indexGW-1,
RHS,INSERT_VALUES, ierr)
endif
else
if (GW_CELL(j,i,k) == 0 .or. GW_CELL(j,i,k) == 2) then
RHS = -H_constant (j,i,k) !-HNOFLOW (j,i,k)
!elseif (GW_CELL (j,i,k) == 2) then
! RHS = -HBOUND(j,i,k)
else
RHS = -B(j,i,k)
-((SS(j,i)*DELTAC*DELTAR*DELTAV(j,i,k)*HOLD(j,i,k)/DELTAT))+
CV_left(j,i,k)*(GWTOP(j,i,k)-HOLD(j,i,k))
endif
call VecSetValue(RHS_vector, indexGW-1, RHS,INSERT_VALUES,
ierr)
end if
end if
! this the temporary h vector that we want to find
h =HOLD(j,i,k)
call VecSetValue(temp_solution, indexGW-1, h,INSERT_VALUES,
ierr)
!end if
end do
end do
end do
! Set values in A_Martrix (A_Matrix)
call MatAssemblyBegin(A_Matrix,MAT_FINAL_ASSEMBLY, ierr)
call MatAssemblyEnd(A_Matrix,MAT_FINAL_ASSEMBLY ,ierr)
! Set values in the solution vector (temp_solution)
call VecAssemblyBegin(temp_solution, ierr)
call VecAssemblyEnd(temp_solution, ierr)
! Set values in the right-hand side vector (RHS_vector)
call VecAssemblyBegin(RHS_vector, ierr)
call VecAssemblyEnd(RHS_vector, ierr)
! deallocate
deallocate (HCOF)
deallocate (B)
deallocate (P)
end subroutine GW_AMatrix_RHS_BC
################################### SOLVER Subroutine ##########################
! This subroutine to solve the system of equation numerically by finite
difference method
! We used PETSc toolkit for scientfic compution
! PCG (to solve the system of linear equations) with Picard iteration to deal
with non-linearity in the equations
! Now I want to try Newton Raphson solver
subroutine GW_solver(t_s,N)
!subroutine GW_solver(t_s,N,RECH,DEEPRECH)
#include <petsc/finclude/petsc.h>
#include <petsc/finclude/petscksp.h>
#include <petsc/finclude/petscvec.h>
#include <petsc/finclude/petscmat.h>
#include <petsc/finclude/petscts.h>
#include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h>
use GW_constants
use GW_param_by_user
use GW_variables
use petsc
use petscvec
use petscmat
use petscts
use petscsnes
implicit none
! Declare arguments for the solver subroutine
integer, intent(in) ::t_s
integer, intent(in) ::N
!real, allocatable,dimension(:,:),intent(inout) ::RECH
!real, allocatable,dimension(:,:),intent(inout) ::DEEPRECH
PetscScalar,pointer :: H_vector(:)
Mat :: A_Matrix
Vec :: RHS_vector
Vec :: old_solution
Vec :: temp_solution
Vec :: r_vector
Vec :: diff_vector
Vec :: norm_vector
Vec :: Residual_vector
PetscInt iter,i,j,k, indexGW
KSP ksp
PC pc
PetscReal tol,norm
PetscErrorCode ierr
PetscReal norm_picard
PetscReal diff_norm
PetscScalar neg_one
PetscScalar zero
!Set tolerance
tol = 1.0e-6
ierr = 0.0
zero = 0.0
neg_one = -1.0
! Initialize PETSc
!call PetscInitialize(ierr)
! create AMatrix and RHS_vector and temp_solution vectors
! print *, N
call MatCreate(PETSC_COMM_WORLD, A_Matrix, ierr)
! print *, "matrix created"
call MatSetSizes(A_Matrix, PETSC_DECIDE, PETSC_DECIDE, N, N, ierr)
call MatSetType(A_Matrix, MATAIJ, ierr) ! as sparse matrix
call MatSetFromOptions(A_Matrix,ierr)
call MatSetUp(A_Matrix,ierr)
call VecCreate(PETSC_COMM_WORLD, RHS_vector, ierr)
call VecSetSizes(RHS_vector, PETSC_DECIDE, N, ierr)
call VecSetFromOptions(RHS_vector,ierr)
call VecSet(RHS_vector, zero, ierr) ! initialize RHS_vector with zero
call VecCreate(PETSC_COMM_WORLD, temp_solution, ierr)
call VecSetSizes(temp_solution, PETSC_DECIDE, N, ierr)
call VecSetFromOptions(temp_solution,ierr)
call VecSet(temp_solution, zero, ierr) ! initialize temp_solution with zero
call VecCreate(PETSC_COMM_WORLD, Residual_vector, ierr)
call VecSetSizes(Residual_vector, PETSC_DECIDE, N, ierr)
call VecSetFromOptions(Residual_vector,ierr)
call VecSet(Residual_vector, zero, ierr) ! initialize Residual_vector with
zero
! Create vectors for old solution and residual
call VecDuplicate(temp_solution, old_solution, ierr)
call VecDuplicate(RHS_vector, r_vector, ierr)
!!!!_________________NEWTON
__________________________________________________
call TSCreate (PETSC_COMM_WORLD, ts, ierr)
call TSSetProblemType (ts, TS_NONLINEAR, ierr)
call TSSetType (ts, TS_BEULER, ierr)
call TSSetFromOption (ts)
call TSSetTimeStep (ts,DELTAT) ! specify time step to DELTAT
call TSSetIJacobian (ts, A_Matrix, A_Matrix, ComputeJacobian, user_data,
ierr)
call TSSetSolution (ts,temp_solution, ierr)
call TSSolve(ts, temp_solution, ierr)
!!!____________________________________________________________________
! ! Create the KSP solver cotext
! call KSPCreate(PETSC_COMM_WORLD, ksp, ierr)
! call KSPSetOperators(ksp, A_Matrix, A_Matrix, ierr)
!
! ! Set the solver type (CG for Conjugate Gradient)
! call KSPSetType(ksp, KSPCG, ierr)
! ! assign the preconditioner type
! call KSPGetPC(ksp,pc,ierr)
! !call PCSetType(pc,PCILU,ierr)
! call PCSetType(pc,PCICC,ierr)
!
! ! Set solver options (e.g tolerance, max_iteration for the solver (default))
! call
KSPSetTolerances(ksp,tol,PETSC_DEFAULT_REAL,PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr)
!
! call MatGetSize(A_matrix,N,N,ierr)
!
! ! Nonlinear Picard Iteration
! do iter = 1, MAX_ITER
!
! call VecCopy(temp_solution, old_solution, ierr)
!
! ! Update the RHS_vector and A-Matrix based on the current
! ! temp solution_vector
!
! ! Call GW_AMatrix_RHS to updqate A_matrix and RHS_vector
!
! call GW_AMatrix_RHS_BC(t_s,N,RHS_vector, &
! temp_solution,A_Matrix)
! !call GW_AMatrix_RHS(t_s,N,RECH,DEEPRECH,RHS_vector, &
! ! temp_solution,A_Matrix)
!
!
!
! ! Solve the linear system using the PCG solver
! call KSPSolve(ksp, RHS_vector, temp_solution, ierr)
! !_______________________ Check for
convergence_________________________!
!
! ! Compute the difference between the residual and the RHS
vector
! ! Calculate the residual vector
! call MatMult(A_Matrix, temp_solution, r_vector,ierr)
!
! ! Compute the difference between the residual and the RHS
vector
! call VecAXPY(r_vector, neg_one, RHS_vector,ierr)
!
! ! Compute the L2 norm of the residual vector
!
! call VecNorm(r_vector, NORM_2,norm_picard, ierr)
!
! ! Calculate the absolute difference between temp_solution and
! ! old_solution
!
!
! ! Now, norm_picard contains the maximum absolute difference
! if (norm_picard .LT. 1.0E-4) then
!
! exit
! else
! continue
! endif
! end do
!
!
! View the matrix
! call MatView(A_Matrix, PETSC_VIEWER_STDOUT_WORLD, ierr)
call VecView(RHS_vector, PETSC_VIEWER_STDOUT_WORLD, ierr)
! store the values from temp_solution into H_vector as array
call VecGetArrayF90(temp_solution, H_vector, ierr)
call VecRestoreArray (temp_solution, H_vector, ierr)
! Now, you can use the 'H_vector' array in your main program or another
! subroutine
do k=1, NLAY
do i=1, NROW
do j= 1, NCOL
indexGW= (k-1)*NCOL*NROW+(i-1)*NCOL+j
HNEW (j,i,k) = H_vector (indexGW)
if (k==1) then
UNCONF_H(j,i) = HNEW(j,i,k)
WTD(j,i) = SurElev(j,i)-UNCONF_H(j,i)
end if
end do
end do
end do
GW_HEAD(:,:,:,t_s+1)=HNEW(:,:,:)
! store the values from temp_solution into H_vector as array
!call VecRestoreArray (temp_solution, H_vector, ierr)
call KSPDestroy(ksp,ierr)
!call PCDestroy(pc, ierr) ! Destroy the preconditioner object
call VecDestroy(old_solution, ierr)
call VecDestroy(r_vector, ierr)
call VecDestroy(Residual_vector, ierr)
call VecDestroy(temp_solution, ierr)
call VecDestroy(RHS_vector, ierr)
call MatDestroy(A_Matrix,ierr)
! call PetscFinalize(ierr)
! print *, " Done first iteration"
end subroutine GW_solver
