########## 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 #include #include 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 #include #include #include #include #include #include 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