numericalsolutionmodule Module Reference

Data Types
type	numericalsolutiontype
interface	synchronize_iface

Functions/Subroutines
subroutine, public	create_numerical_solution (num_sol, filename, id)
	@ brief Create a new solution More...
subroutine	allocate_scalars (this)
	@ brief Allocate scalars More...
subroutine	allocate_arrays (this)
	@ brief Allocate arrays More...
subroutine	sln_df (this)
	@ brief Define the solution More...
subroutine	sln_ar (this)
	@ brief Allocate and read data More...
subroutine	sln_dt (this)
	@ brief Calculate delt More...
subroutine	sln_ad (this)
	@ brief Advance solution More...
subroutine	sln_ot (this)
	@ brief Output solution More...
subroutine	sln_fp (this)
	@ brief Finalize solution More...
subroutine	sln_da (this)
	@ brief Deallocate solution More...
subroutine	sln_ca (this, isgcnvg, isuppress_output)
	@ brief Solve solution More...
subroutine	writecsvheader (this)
	@ brief CSV header More...
subroutine	writeptcinfotofile (this, kper)
	@ brief PTC header More...
subroutine	preparesolve (this)
	@ brief prepare to solve More...
subroutine	solve (this, kiter)
	@ brief Build and solve the simulation More...
subroutine	finalizesolve (this, kiter, isgcnvg, isuppress_output)
	@ brief finalize a solution More...
subroutine	sln_buildsystem (this, kiter, inewton)
subroutine	convergence_summary (this, iu, im, itertot_timestep)
	@ brief Solution convergence summary More...
subroutine	csv_convergence_summary (this, iu, totim, kper, kstp, kouter, niter, istart, kstart)
	@ brief Solution convergence CSV summary More...
subroutine	save (this, filename)
	@ brief Save solution data to a file More...
subroutine	add_model (this, mp)
	@ brief Add a model More...
type(listtype) function, pointer	get_models (this)
	Get a list of models. More...
subroutine	add_exchange (this, exchange)
	Add exchange. More...
type(listtype) function, pointer	get_exchanges (this)
	Returns a pointer to the list of exchanges in this solution. More...
subroutine	sln_connect (this)
	@ brief Assign solution connections More...
subroutine	sln_reset (this)
	@ brief Reset the solution More...
subroutine	sln_ls (this, kiter, kstp, kper, in_iter, iptc, ptcf)
	@ brief Solve the linear system of equations More...
subroutine	sln_setouter (this, ifdparam)
	@ brief Set default Picard iteration variables More...
subroutine	sln_backtracking (this, mp, cp, kiter)
	@ brief Perform backtracking More...
subroutine	sln_backtracking_xupdate (this, bt_flag)
	@ brief Backtracking update of the dependent variable More...
integer(i4b) function	get_backtracking_flag (this)
	Check if backtracking should be applied for this solution,. More...
integer(i4b) function	sln_get_idvscale (this)
	Check if dependent variable scalining should be applied for this solution,. More...
subroutine	apply_backtracking (this)
	Update x with backtracking. More...
subroutine	sln_l2norm (this, l2norm)
	@ brief Calculate the solution L-2 norm for all active cells using More...
subroutine	sln_maxval (this, nsize, v, vmax)
	@ brief Get the maximum value from a vector More...
subroutine	sln_calcdx (this, neq, active, x, xtemp, dx)
	@ brief Calculate dependent-variable change More...
subroutine	sln_calc_ptc (this, iptc, ptcf)
	Calculate pseudo-transient continuation factor. More...
subroutine	sln_calc_residual (this, vec_resid)
	Calculate the current residual vector r = A*x - b,. More...
subroutine	sln_underrelax (this, kiter, bigch, neq, active, x, xtemp)
	@ brief Under-relaxation More...
subroutine	sln_get_dxmax (this, hncg, lrch)
	@ brief Determine maximum dependent-variable change More...
logical(lgp) function	sln_has_converged (this, max_dvc)
integer(i4b) function	sln_package_convergence (this, dpak, cpakout, iend)
	Check package convergence. More...
integer(i4b) function	sln_sync_newtonur_flag (this, inewtonur)
	Synchronize Newton Under-relaxation flag. More...
logical(lgp) function	sln_nur_has_converged (this, dxold_max, hncg)
	Custom convergence check for when Newton UR has been applied. More...
subroutine	sln_get_loc (this, nodesln, str)
	@ brief Get cell location string More...
subroutine	sln_get_nodeu (this, nodesln, im, nodeu)
	@ brief Get user node number More...
class(numericalsolutiontype) function, pointer, public	castasnumericalsolutionclass (obj)
	@ brief Cast a object as a Numerical Solution More...
class(numericalsolutiontype) function, pointer, public	getnumericalsolutionfromlist (list, idx)
	@ brief Get a numerical solution More...

Variables
integer(i4b), parameter	ims_solver = 1
integer(i4b), parameter	petsc_solver = 2

Function/Subroutine Documentation

add_exchange()

subroutine numericalsolutionmodule::add_exchange ( class(numericalsolutiontype)  this, class(baseexchangetype), intent(in), pointer  exchange  )

private

Add and exchange to thisexchangelist.

Parameters

	this	NumericalSolutionType instance
[in]	exchange	model exchange instance

Definition at line 2340 of file NumericalSolution.f90.

    ! -- dummy variables
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    class(BaseExchangeType), pointer, intent(in) :: exchange !< model exchange instance
    ! -- local variables
    class(NumericalExchangeType), pointer :: num_ex => null()
    !
    ! -- add exchange
    select type (exchange)
    class is (numericalexchangetype)
      num_ex => exchange
      call addnumericalexchangetolist(this%exchangelist, num_ex)
    end select

Here is the call graph for this function:

add_model()

subroutine numericalsolutionmodule::add_model	(	class(numericalsolutiontype)	this,
		class(basemodeltype), intent(in), pointer	mp
	)

Add a model to solution.

Parameters

	this	NumericalSolutionType instance
[in]	mp	model instance

Definition at line 2305 of file NumericalSolution.f90.

    ! -- dummy variables
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    class(BaseModelType), pointer, intent(in) :: mp !< model instance
    ! -- local variables
    class(NumericalModelType), pointer :: m => null()
    !
    ! -- add a model
    select type (mp)
    class is (numericalmodeltype)
      m => mp
      call addnumericalmodeltolist(this%modellist, m)
    end select

Here is the call graph for this function:

allocate_arrays()

subroutine numericalsolutionmodule::allocate_arrays

(

class(numericalsolutiontype)

this

)

Allocate arrays for a new solution.

Parameters

this	NumericalSolutionType instance

Definition at line 387 of file NumericalSolution.f90.

    ! -- modules
    use memorymanagermodule, only: mem_allocate
    ! -- dummy variables
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    ! -- local variables
    class(NumericalModelType), pointer :: mp => null()
    integer(I4B) :: i
    integer(I4B) :: ieq
    !
    ! -- initialize the number of models in the solution
    this%convnmod = this%modellist%Count()
    !
    ! -- allocate arrays
    call mem_allocate(this%active, this%neq, 'IACTIVE', this%memory_path)
    call mem_allocate(this%xtemp, this%neq, 'XTEMP', this%memory_path)
    call mem_allocate(this%dxold, this%neq, 'DXOLD', this%memory_path)
    call mem_allocate(this%hncg, 0, 'HNCG', this%memory_path)
    call mem_allocate(this%lrch, 3, 0, 'LRCH', this%memory_path)
    call mem_allocate(this%wsave, 0, 'WSAVE', this%memory_path)
    call mem_allocate(this%hchold, 0, 'HCHOLD', this%memory_path)
    call mem_allocate(this%deold, 0, 'DEOLD', this%memory_path)
    call mem_allocate(this%convmodstart, this%convnmod + 1, 'CONVMODSTART', &
                      this%memory_path)
    !
    ! -- initialize allocated arrays
    do i = 1, this%neq
      this%xtemp(i) = dzero
      this%dxold(i) = dzero
      this%active(i) = 1 !default is active
    end do
    !
    ! -- initialize convmodstart
    ieq = 1
    this%convmodstart(1) = ieq
    do i = 1, this%modellist%Count()
      mp => getnumericalmodelfromlist(this%modellist, i)
      ieq = ieq + mp%neq
      this%convmodstart(i + 1) = ieq
    end do

Here is the call graph for this function:

allocate_scalars()

subroutine numericalsolutionmodule::allocate_scalars

(

class(numericalsolutiontype)

this

)

Allocate scalars for a new solution.

Definition at line 285 of file NumericalSolution.f90.

    ! -- modules
    use memorymanagermodule, only: mem_allocate
    ! -- dummy variables
    class(NumericalSolutionType) :: this
    !
    ! -- allocate scalars
    call mem_allocate(this%id, 'ID', this%memory_path)
    call mem_allocate(this%iu, 'IU', this%memory_path)
    call mem_allocate(this%ttform, 'TTFORM', this%memory_path)
    call mem_allocate(this%ttsoln, 'TTSOLN', this%memory_path)
    call mem_allocate(this%isymmetric, 'ISYMMETRIC', this%memory_path)
    call mem_allocate(this%neq, 'NEQ', this%memory_path)
    call mem_allocate(this%matrix_offset, 'MATRIX_OFFSET', this%memory_path)
    call mem_allocate(this%dvclose, 'DVCLOSE', this%memory_path)
    call mem_allocate(this%bigchold, 'BIGCHOLD', this%memory_path)
    call mem_allocate(this%bigch, 'BIGCH', this%memory_path)
    call mem_allocate(this%relaxold, 'RELAXOLD', this%memory_path)
    call mem_allocate(this%res_prev, 'RES_PREV', this%memory_path)
    call mem_allocate(this%res_new, 'RES_NEW', this%memory_path)
    call mem_allocate(this%icnvg, 'ICNVG', this%memory_path)
    call mem_allocate(this%itertot_timestep, 'ITERTOT_TIMESTEP', this%memory_path)
    call mem_allocate(this%iouttot_timestep, 'IOUTTOT_TIMESTEP', this%memory_path)
    call mem_allocate(this%itertot_sim, 'INNERTOT_SIM', this%memory_path)
    call mem_allocate(this%mxiter, 'MXITER', this%memory_path)
    call mem_allocate(this%linsolver, 'LINSOLVER', this%memory_path)
    call mem_allocate(this%nonmeth, 'NONMETH', this%memory_path)
    call mem_allocate(this%iprims, 'IPRIMS', this%memory_path)
    call mem_allocate(this%theta, 'THETA', this%memory_path)
    call mem_allocate(this%akappa, 'AKAPPA', this%memory_path)
    call mem_allocate(this%gamma, 'GAMMA', this%memory_path)
    call mem_allocate(this%amomentum, 'AMOMENTUM', this%memory_path)
    call mem_allocate(this%breduc, 'BREDUC', this%memory_path)
    call mem_allocate(this%btol, 'BTOL', this%memory_path)
    call mem_allocate(this%res_lim, 'RES_LIM', this%memory_path)
    call mem_allocate(this%numtrack, 'NUMTRACK', this%memory_path)
    call mem_allocate(this%ibflag, 'IBFLAG', this%memory_path)
    call mem_allocate(this%icsvouterout, 'ICSVOUTEROUT', this%memory_path)
    call mem_allocate(this%icsvinnerout, 'ICSVINNEROUT', this%memory_path)
    call mem_allocate(this%nitermax, 'NITERMAX', this%memory_path)
    call mem_allocate(this%convnmod, 'CONVNMOD', this%memory_path)
    call mem_allocate(this%iallowptc, 'IALLOWPTC', this%memory_path)
    call mem_allocate(this%iptcopt, 'IPTCOPT', this%memory_path)
    call mem_allocate(this%iptcout, 'IPTCOUT', this%memory_path)
    call mem_allocate(this%l2norm0, 'L2NORM0', this%memory_path)
    call mem_allocate(this%ptcdel, 'PTCDEL', this%memory_path)
    call mem_allocate(this%ptcdel0, 'PTCDEL0', this%memory_path)
    call mem_allocate(this%ptcexp, 'PTCEXP', this%memory_path)
    call mem_allocate(this%atsfrac, 'ATSFRAC', this%memory_path)
    call mem_allocate(this%idv_scale, 'IDV_SCALE', this%memory_path)
    call mem_allocate(this%dscale, 'DSCALE', this%memory_path)
    !
    ! -- initialize scalars
    this%isymmetric = 0
    this%id = 0
    this%iu = 0
    this%ttform = dzero
    this%ttsoln = dzero
    this%neq = 0
    this%dvclose = dzero
    this%bigchold = dzero
    this%bigch = dzero
    this%relaxold = dzero
    this%res_prev = dzero
    this%icnvg = 0
    this%itertot_timestep = 0
    this%iouttot_timestep = 0
    this%itertot_sim = 0
    this%mxiter = 0
    this%linsolver = ims_solver
    this%nonmeth = 0
    this%iprims = 0
    this%theta = done
    this%akappa = dzero
    this%gamma = done
    this%amomentum = dzero
    this%breduc = dzero
    this%btol = 0
    this%res_lim = dzero
    this%numtrack = 0
    this%ibflag = 0
    this%icsvouterout = 0
    this%icsvinnerout = 0
    this%nitermax = 0
    this%convnmod = 0
    this%iallowptc = 1
    this%iptcopt = 0
    this%iptcout = 0
    this%l2norm0 = dzero
    this%ptcdel = dzero
    this%ptcdel0 = dzero
    this%ptcexp = done
    this%atsfrac = donethird
    this%idv_scale = 0
    this%dscale = done

apply_backtracking()

subroutine numericalsolutionmodule::apply_backtracking ( class(numericalsolutiontype)  this )

private

Parameters

this	NumericalSolutionType instance

Definition at line 2891 of file NumericalSolution.f90.

    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    ! local
    integer(I4B) :: n
    real(DP) :: delx
 
    do n = 1, this%neq
      if (this%active(n) < 1) cycle
      delx = this%breduc * (this%x(n) - this%xtemp(n))
      this%x(n) = this%xtemp(n) + delx
    end do
 

castasnumericalsolutionclass()

class(numericalsolutiontype) function, pointer, public numericalsolutionmodule::castasnumericalsolutionclass

(

class(*), intent(inout), pointer

obj

)

Get a numerical solution from a list.

Parameters

[in,out]

obj

generic object

Returns

output NumericalSolutionType

Definition at line 3365 of file NumericalSolution.f90.

    ! -- dummy variables
    class(*), pointer, intent(inout) :: obj !< generic object
    ! -- return variable
    class(NumericalSolutionType), pointer :: res !< output NumericalSolutionType
    !
    ! -- initialize return variable
    res => null()
    !
    ! -- determine if obj is associated
    if (.not. associated(obj)) return
    !
    ! -- set res
    select type (obj)
    class is (numericalsolutiontype)
      res => obj
    end select

Here is the caller graph for this function:

convergence_summary()

subroutine numericalsolutionmodule::convergence_summary ( class(numericalsolutiontype)  this, integer(i4b), intent(in)  iu, integer(i4b), intent(in)  im, integer(i4b), intent(in)  itertot_timestep  )

private

Save convergence summary to a File.

Parameters

	this	NumericalSolutionType instance
[in]	iu	file unit number for summary file
[in]	im	model number
[in]	itertot_timestep	total iteration for the time step

Definition at line 2036 of file NumericalSolution.f90.

    ! -- modules
    use inputoutputmodule, only: getunit
    ! -- dummy variables
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    integer(I4B), intent(in) :: iu !< file unit number for summary file
    integer(I4B), intent(in) :: im !< model number
    integer(I4B), intent(in) :: itertot_timestep !< total iteration for the time step
    ! -- local variables
    character(len=LINELENGTH) :: title
    character(len=LINELENGTH) :: tag
    character(len=LENPAKLOC) :: loc_dvmax_str
    character(len=LENPAKLOC) :: loc_rmax_str
    integer(I4B) :: ntabrows
    integer(I4B) :: ntabcols
    integer(I4B) :: iinner
    integer(I4B) :: i0
    integer(I4B) :: iouter
    integer(I4B) :: j
    integer(I4B) :: k
    integer(I4B) :: locdv
    integer(I4B) :: locdr
    real(DP) :: dv !< the maximum change in the dependent variable
    real(DP) :: res !< the maximum value of the residual vector
    !
    ! -- initialize local variables
    loc_dvmax_str = ''
    loc_rmax_str = ''
    iouter = 1
    locdv = 0
    locdr = 0
    !
    ! -- initialize inner iteration summary table
    if (.not. associated(this%innertab)) then
      !
      ! -- create outer iteration table
      ! -- table dimensions
      ntabrows = itertot_timestep
      ntabcols = 7
      !
      ! -- initialize table and define columns
      title = trim(this%memory_path)//' INNER ITERATION SUMMARY'
      call table_cr(this%innertab, this%name, title)
      call this%innertab%table_df(ntabrows, ntabcols, iu)
      tag = 'TOTAL ITERATION'
      call this%innertab%initialize_column(tag, 10, alignment=tabright)
      tag = 'OUTER ITERATION'
      call this%innertab%initialize_column(tag, 10, alignment=tabright)
      tag = 'INNER ITERATION'
      call this%innertab%initialize_column(tag, 10, alignment=tabright)
      tag = 'MAXIMUM CHANGE'
      call this%innertab%initialize_column(tag, 15, alignment=tabright)
      tag = 'MAXIMUM CHANGE MODEL-(CELLID)'
      call this%innertab%initialize_column(tag, lenpakloc, alignment=tabright)
      tag = 'MAXIMUM RESIDUAL'
      call this%innertab%initialize_column(tag, 15, alignment=tabright)
      tag = 'MAXIMUM RESIDUAL MODEL-(CELLID)'
      call this%innertab%initialize_column(tag, lenpakloc, alignment=tabright)
      !
      ! -- reset the output unit and the number of rows (maxbound)
    else
      call this%innertab%set_maxbound(itertot_timestep)
      call this%innertab%set_iout(iu)
    end if
    !
    ! -- write the inner iteration summary to unit iu
    i0 = 0
    do k = 1, itertot_timestep
      iinner = this%cnvg_summary%itinner(k)
      if (iinner <= i0) then
        iouter = iouter + 1
      end if
      if (im > this%convnmod) then
        dv = dzero
        res = dzero
        do j = 1, this%convnmod
          if (abs(this%cnvg_summary%convdvmax(j, k)) > abs(dv)) then
            locdv = this%cnvg_summary%convlocdv(j, k)
            dv = this%cnvg_summary%convdvmax(j, k)
          end if
          if (abs(this%cnvg_summary%convrmax(j, k)) > abs(res)) then
            locdr = this%cnvg_summary%convlocr(j, k)
            res = this%cnvg_summary%convrmax(j, k)
          end if
        end do
      else
        locdv = this%cnvg_summary%convlocdv(im, k)
        locdr = this%cnvg_summary%convlocr(im, k)
        dv = this%cnvg_summary%convdvmax(im, k)
        res = this%cnvg_summary%convrmax(im, k)
      end if
      call this%sln_get_loc(locdv, loc_dvmax_str)
      call this%sln_get_loc(locdr, loc_rmax_str)
      !
      ! -- add data to innertab
      call this%innertab%add_term(k)
      call this%innertab%add_term(iouter)
      call this%innertab%add_term(iinner)
      call this%innertab%add_term(dv)
      call this%innertab%add_term(adjustr(trim(loc_dvmax_str)))
      call this%innertab%add_term(res)
      call this%innertab%add_term(adjustr(trim(loc_rmax_str)))
      !
      ! -- update i0
      i0 = iinner
    end do

Here is the call graph for this function:

create_numerical_solution()

subroutine, public numericalsolutionmodule::create_numerical_solution	(	class(numericalsolutiontype), pointer	num_sol,
		character(len=*), intent(in)	filename,
		integer(i4b), intent(in)	id
	)

Create a new solution using the data in filename, assign this new solution an id number and store the solution in the basesolutionlist. Also open the filename for later reading.

Parameters

	num_sol	the create solution
[in]	filename	solution input file name
[in]	id	solution id

Definition at line 238 of file NumericalSolution.f90.

    ! -- modules
    use simvariablesmodule, only: iout
    use inputoutputmodule, only: getunit, openfile
    ! -- dummy variables
    class(NumericalSolutionType), pointer :: num_sol !< the create solution
    character(len=*), intent(in) :: filename !< solution input file name
    integer(I4B), intent(in) :: id !< solution id
    ! -- local variables
    integer(I4B) :: inunit
    class(BaseSolutionType), pointer :: solbase => null()
    character(len=LENSOLUTIONNAME) :: solutionname
    !
    ! -- Create a new solution and add it to the basesolutionlist container
    solbase => num_sol
    write (solutionname, '(a, i0)') 'SLN_', id
    !
    num_sol%name = solutionname
    num_sol%memory_path = create_mem_path(solutionname)
    allocate (num_sol%modellist)
    allocate (num_sol%exchangelist)
    !
    call num_sol%allocate_scalars()
    !
    call addbasesolutiontolist(basesolutionlist, solbase)
    !
    num_sol%id = id
    !
    ! -- Open solution input file for reading later after problem size is known
    !    Check to see if the file is already opened, which can happen when
    !    running in single model mode
    inquire (file=filename, number=inunit)
 
    if (inunit < 0) inunit = getunit()
    num_sol%iu = inunit
    write (iout, '(/a,a)') ' Creating solution: ', num_sol%name
    call openfile(num_sol%iu, iout, filename, 'IMS')
    !
    ! -- Initialize block parser
    call num_sol%parser%Initialize(num_sol%iu, iout)

Here is the call graph for this function:

Here is the caller graph for this function:

csv_convergence_summary()

subroutine numericalsolutionmodule::csv_convergence_summary	(	class(numericalsolutiontype)	this,
		integer(i4b), intent(in)	iu,
		real(dp), intent(in)	totim,
		integer(i4b), intent(in)	kper,
		integer(i4b), intent(in)	kstp,
		integer(i4b), intent(in)	kouter,
		integer(i4b), intent(in)	niter,
		integer(i4b), intent(in)	istart,
		integer(i4b), intent(in)	kstart
	)

Save convergence summary to a comma-separated value file.

Parameters

	this	NumericalSolutionType instance
[in]	iu	file unit number
[in]	totim	total simulation time
[in]	kper	stress period number
[in]	kstp	time step number
[in]	kouter	number of outer (Picard) iterations
[in]	niter	number of inner iteration in this time step
[in]	istart	starting iteration number for this time step
[in]	kstart	starting position in the conv* arrays

Definition at line 2149 of file NumericalSolution.f90.

    ! -- modules
    use inputoutputmodule, only: getunit
    ! -- dummy variables
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    integer(I4B), intent(in) :: iu !< file unit number
    real(DP), intent(in) :: totim !< total simulation time
    integer(I4B), intent(in) :: kper !< stress period number
    integer(I4B), intent(in) :: kstp !< time step number
    integer(I4B), intent(in) :: kouter !< number of outer (Picard) iterations
    integer(I4B), intent(in) :: niter !< number of inner iteration in this time step
    integer(I4B), intent(in) :: istart !< starting iteration number for this time step
    integer(I4B), intent(in) :: kstart !< starting position in the conv* arrays
    ! -- local
    integer(I4B) :: itot
    integer(I4B) :: m_idx, j, k
    integer(I4B) :: kpos
    integer(I4B) :: loc_dvmax !< solution node number (row) of max. dep. var. change
    integer(I4B) :: loc_rmax !< solution node number (row) of max. residual
    integer(I4B) :: model_id, node_user
    real(DP) :: dvmax !< maximum dependent variable change
    real(DP) :: rmax !< maximum residual
    class(NumericalModelType), pointer :: num_mod => null()
    !
    ! -- initialize local variables
    itot = istart
    !
    ! -- write inner iteration results to the inner csv output file
    do k = 1, niter
      kpos = kstart + k - 1
      write (iu, '(*(G0,:,","))', advance='NO') &
        itot, totim, kper, kstp, kouter, k
      !
      ! -- solution summary
      dvmax = dzero
      rmax = dzero
      do j = 1, this%convnmod
        if (abs(this%cnvg_summary%convdvmax(j, kpos)) > abs(dvmax)) then
          loc_dvmax = this%cnvg_summary%convlocdv(j, kpos)
          dvmax = this%cnvg_summary%convdvmax(j, kpos)
        end if
        if (abs(this%cnvg_summary%convrmax(j, kpos)) > abs(rmax)) then
          loc_rmax = this%cnvg_summary%convlocr(j, kpos)
          rmax = this%cnvg_summary%convrmax(j, kpos)
        end if
      end do
      !
      ! -- no change, could be anywhere
      if (dvmax == dzero) loc_dvmax = 0
      if (rmax == dzero) loc_rmax = 0
      !
      ! -- get model number and user node number for max. dep. var. change
      if (loc_dvmax > 0) then
        call this%sln_get_nodeu(loc_dvmax, m_idx, node_user)
        num_mod => getnumericalmodelfromlist(this%modellist, m_idx)
        model_id = num_mod%id
      else
        model_id = 0
        node_user = 0
      end if
      write (iu, '(*(G0,:,","))', advance='NO') '', dvmax, model_id, node_user
      !
      ! -- get model number and user node number for max. residual
      if (loc_rmax > 0) then
        call this%sln_get_nodeu(loc_rmax, m_idx, node_user)
        num_mod => getnumericalmodelfromlist(this%modellist, m_idx)
        model_id = num_mod%id
      else
        model_id = 0
        node_user = 0
      end if
      write (iu, '(*(G0,:,","))', advance='NO') '', rmax, model_id, node_user
      !
      ! -- write ims acceleration parameters
      if (this%linsolver == ims_solver) then
        write (iu, '(*(G0,:,","))', advance='NO') &
          '', trim(adjustl(this%caccel(kpos)))
      end if
      !
      ! -- write information for each model
      if (this%convnmod > 1 .or. simulation_mode == "PARALLEL") then
        do j = 1, this%cnvg_summary%convnmod
          loc_dvmax = this%cnvg_summary%convlocdv(j, kpos)
          dvmax = this%cnvg_summary%convdvmax(j, kpos)
          loc_rmax = this%cnvg_summary%convlocr(j, kpos)
          rmax = this%cnvg_summary%convrmax(j, kpos)
          !
          ! -- get model number and user node number for max. dep. var. change
          if (loc_dvmax > 0) then
            call this%sln_get_nodeu(loc_dvmax, m_idx, node_user)
          else
            node_user = 0
          end if
          write (iu, '(*(G0,:,","))', advance='NO') '', dvmax, node_user
          !
          ! -- get model number and user node number for max. residual
          if (loc_rmax > 0) then
            call this%sln_get_nodeu(loc_rmax, m_idx, node_user)
          else
            node_user = 0
          end if
          write (iu, '(*(G0,:,","))', advance='NO') '', rmax, node_user
        end do
      end if
      !
      ! -- write line
      write (iu, '(a)') ''
      !
      ! -- update itot
      itot = itot + 1
    end do
    !
    ! -- flush file
    flush (iu)

Here is the call graph for this function:

finalizesolve()

subroutine numericalsolutionmodule::finalizesolve	(	class(numericalsolutiontype)	this,
		integer(i4b), intent(in)	kiter,
		integer(i4b), intent(inout)	isgcnvg,
		integer(i4b), intent(in)	isuppress_output
	)

Finalize the solution. Called after the outer iteration loop.

Parameters

	this	NumericalSolutionType instance
[in]	kiter	Picard iteration number after convergence or failure
[in,out]	isgcnvg	solution group convergence flag
[in]	isuppress_output	flag for suppressing output

Definition at line 1882 of file NumericalSolution.f90.

    ! -- modules
    use tdismodule, only: kper, kstp
    ! -- dummy variables
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    integer(I4B), intent(in) :: kiter !< Picard iteration number after convergence or failure
    integer(I4B), intent(inout) :: isgcnvg !< solution group convergence flag
    integer(I4B), intent(in) :: isuppress_output !< flag for suppressing output
    ! -- local variables
    integer(I4B) :: ic, im
    class(NumericalModelType), pointer :: mp => null()
    class(NumericalExchangeType), pointer :: cp => null()
    ! -- formats for convergence info
    character(len=*), parameter :: fmtnocnvg = &
      "(1X,'Solution ', i0, ' did not converge for stress period ', i0, &
      &' and time step ', i0)"
    character(len=*), parameter :: fmtcnvg = &
      "(1X, I0, ' CALLS TO NUMERICAL SOLUTION ', 'IN TIME STEP ', I0, &
      &' STRESS PERIOD ',I0,/1X,I0,' TOTAL ITERATIONS')"
 
    ! start timer
    call g_prof%start("Finalize solve"//this%id_postfix, this%tmr_final_solve)
 
    !
    ! -- finalize the outer iteration table
    if (this%iprims > 0) then
      call this%outertab%finalize_table()
    end if
    !
    ! -- write convergence info
    !
    ! -- convergence was achieved
    if (this%icnvg /= 0) then
      if (this%iprims > 0) then
        write (iout, fmtcnvg) kiter, kstp, kper, this%itertot_timestep
      end if
      !
      ! -- convergence was not achieved
    else
      write (iout, fmtnocnvg) this%id, kper, kstp
    end if
    !
    ! -- write inner iteration convergence summary
    if (this%iprims == 2) then
      !
      ! -- write summary for each model
      do im = 1, this%modellist%Count()
        mp => getnumericalmodelfromlist(this%modellist, im)
        call this%convergence_summary(mp%iout, im, this%itertot_timestep)
      end do
      !
      ! -- write summary for entire solution
      call this%convergence_summary(iout, this%convnmod + 1, &
                                    this%itertot_timestep)
    end if
    !
    ! -- set solution group convergence flag
    if (this%icnvg == 0) isgcnvg = 0
 
    call g_prof%start("Calculate flows", this%tmr_flows)
 
    !
    ! -- Calculate flow for each model
    do im = 1, this%modellist%Count()
      mp => getnumericalmodelfromlist(this%modellist, im)
      call mp%model_cq(this%icnvg, isuppress_output)
    end do
    !
    ! -- Calculate flow for each exchange
    do ic = 1, this%exchangelist%Count()
      cp => getnumericalexchangefromlist(this%exchangelist, ic)
      call cp%exg_cq(isgcnvg, isuppress_output, this%id)
    end do
 
    call g_prof%stop(this%tmr_flows)
    call g_prof%start("Calculate budgets", this%tmr_budgets)
 
    !
    ! -- Budget terms for each model
    do im = 1, this%modellist%Count()
      mp => getnumericalmodelfromlist(this%modellist, im)
      call mp%model_bd(this%icnvg, isuppress_output)
    end do
    !
    ! -- Budget terms for each exchange
    do ic = 1, this%exchangelist%Count()
      cp => getnumericalexchangefromlist(this%exchangelist, ic)
      call cp%exg_bd(isgcnvg, isuppress_output, this%id)
    end do
 
    ! stop timer
    call g_prof%stop(this%tmr_budgets)
    call g_prof%stop(this%tmr_final_solve)
 

Here is the call graph for this function:

get_backtracking_flag()

integer(i4b) function numericalsolutionmodule::get_backtracking_flag ( class(numericalsolutiontype)  this )

private

Parameters

this	NumericalSolutionType instance

Returns

backtracking flag (1) backtracking performed (0) backtracking not performed

Definition at line 2838 of file NumericalSolution.f90.

    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    integer(I4B) :: bt_flag !< backtracking flag (1) backtracking performed (0) backtracking not performed
    ! local
    integer(I4B) :: n
    real(DP) :: dx
    real(DP) :: dx_abs
    real(DP) :: dx_abs_max
 
    ! default is off
    bt_flag = 0
 
    ! find max. change
    dx_abs_max = 0.0
    do n = 1, this%neq
      if (this%active(n) < 1) cycle
      dx = this%x(n) - this%xtemp(n)
      dx_abs = abs(dx)
      if (dx_abs > dx_abs_max) dx_abs_max = dx_abs
    end do
 
    ! if backtracking, set flag
    if (this%breduc * dx_abs_max >= this%dvclose) then
      bt_flag = 1
    end if
 

get_exchanges()

type(listtype) function, pointer numericalsolutionmodule::get_exchanges ( class(numericalsolutiontype)  this )

private

Parameters

this	instance of the numerical solution

Returns

pointer to the exchange list

Definition at line 2357 of file NumericalSolution.f90.

    class(NumericalSolutionType) :: this !< instance of the numerical solution
    type(ListType), pointer :: exchanges !< pointer to the exchange list
 
    exchanges => this%exchangelist
 

get_models()

type(listtype) function, pointer numericalsolutionmodule::get_models ( class(numericalsolutiontype)  this )

private

Returns a pointer to the list of models in this solution.

Returns

pointer to the model list

Parameters

this	NumericalSolutionType instance

Definition at line 2325 of file NumericalSolution.f90.

    ! -- return variable
    type(ListType), pointer :: models !< pointer to the model list
    ! -- dummy variables
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
 
    models => this%modellist
 

getnumericalsolutionfromlist()

class(numericalsolutiontype) function, pointer, public numericalsolutionmodule::getnumericalsolutionfromlist	(	type(listtype), intent(inout)	list,
		integer(i4b), intent(in)	idx
	)

Get a numerical solution from a list.

Parameters

[in,out]	list	list of numerical solutions
[in]	idx	value to retrieve from the list

Returns

numerical solution

Definition at line 3389 of file NumericalSolution.f90.

    ! -- dummy variables
    type(ListType), intent(inout) :: list !< list of numerical solutions
    integer(I4B), intent(in) :: idx !< value to retrieve from the list
    ! -- return variables
    class(NumericalSolutionType), pointer :: res !< numerical solution
    ! -- local variables
    class(*), pointer :: obj
    !
    obj => list%GetItem(idx)
    res => castasnumericalsolutionclass(obj)

Here is the call graph for this function:

preparesolve()

subroutine numericalsolutionmodule::preparesolve ( class(numericalsolutiontype)  this )

private

Prepare for the system solve by advancing the simulation.

Parameters

this	NumericalSolutionType instance

Definition at line 1467 of file NumericalSolution.f90.

    ! -- dummy variables
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    ! -- local variables
    integer(I4B) :: ic
    integer(I4B) :: im
    class(NumericalExchangeType), pointer :: cp => null()
    class(NumericalModelType), pointer :: mp => null()
 
    ! start timer
    call g_prof%start("Prepare solve"//this%id_postfix, this%tmr_prep_solve)
 
    ! synchronize for AD
    call this%synchronize(stg_bfr_exg_ad, this%synchronize_ctx)
 
    ! -- Exchange advance
    do ic = 1, this%exchangelist%Count()
      cp => getnumericalexchangefromlist(this%exchangelist, ic)
      call cp%exg_ad()
    end do
 
    ! -- Model advance
    do im = 1, this%modellist%Count()
      mp => getnumericalmodelfromlist(this%modellist, im)
      call mp%model_ad()
    end do
 
    ! advance solution
    call this%sln_ad()
 
    ! stop timer
    call g_prof%stop(this%tmr_prep_solve)
 

Here is the call graph for this function:

save()

subroutine numericalsolutionmodule::save	(	class(numericalsolutiontype)	this,
		character(len=*), intent(in)	filename
	)

Save solution ia vector, ja vector , coefficient matrix, right-hand side vector, and the dependent-variable vector to a file.

Parameters

	this	NumericalSolutionType instance
[in]	filename	filename to save solution data

Definition at line 2272 of file NumericalSolution.f90.

    use sparsematrixmodule
    ! -- modules
    use inputoutputmodule, only: getunit
    ! -- dummy variables
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    character(len=*), intent(in) :: filename !< filename to save solution data
    ! -- local variables
    integer(I4B) :: inunit
    !
    select type (spm => this%system_matrix)
    class is (sparsematrixtype)
      inunit = getunit()
      open (unit=inunit, file=filename, status='unknown')
      write (inunit, *) 'ia'
      write (inunit, *) spm%ia
      write (inunit, *) 'ja'
      write (inunit, *) spm%ja
      write (inunit, *) 'amat'
      write (inunit, *) spm%amat
      write (inunit, *) 'rhs'
      write (inunit, *) this%rhs
      write (inunit, *) 'x'
      write (inunit, *) this%x
      close (inunit)
    end select

Here is the call graph for this function:

sln_ad()

subroutine numericalsolutionmodule::sln_ad

(

class(numericalsolutiontype)

this

)

Advance solution.

Parameters

this	NumericalSolutionType instance

Definition at line 1140 of file NumericalSolution.f90.

    ! -- modules
    use tdismodule, only: kstp, kper
    ! -- dummy variables
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    !
    ! -- write headers to CSV file
    if (kper == 1 .and. kstp == 1) then
      call this%writeCSVHeader()
    end if
 
    ! write PTC info on models to iout
    call this%writePTCInfoToFile(kper)
 
    ! reset convergence flag and inner solve counter
    this%icnvg = 0
    this%itertot_timestep = 0
    this%iouttot_timestep = 0

sln_ar()

subroutine numericalsolutionmodule::sln_ar

(

class(numericalsolutiontype)

this

)

Allocate and read data for a solution.

Parameters

this	NumericalSolutionType instance

Definition at line 535 of file NumericalSolution.f90.

    ! -- modules
    use memorymanagermodule, only: mem_reallocate
    use simvariablesmodule, only: iout
    use simmodule, only: store_error, count_errors, deprecation_warning
    use inputoutputmodule, only: getunit, openfile
    ! -- dummy variables
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    ! -- local variables
    class(NumericalModelType), pointer :: mp => null()
    class(NumericalExchangeType), pointer :: cp => null()
    character(len=linelength) :: warnmsg
    character(len=linelength) :: keyword
    character(len=linelength) :: fname
    character(len=linelength) :: msg
    integer(I4B) :: i
    integer(I4B) :: ifdparam, mxvl, npp
    integer(I4B) :: ierr
    logical(LGP) :: isfound, endOfBlock
    integer(I4B) :: ival
    real(DP) :: rval
    character(len=*), parameter :: fmtcsvout = &
      "(4x, 'CSV OUTPUT WILL BE SAVED TO FILE: ', a, &
      &/4x, 'OPENED ON UNIT: ', I7)"
    character(len=*), parameter :: fmtptcout = &
      "(4x, 'PTC OUTPUT WILL BE SAVED TO FILE: ', a, &
      &/4x, 'OPENED ON UNIT: ', I7)"
    character(len=*), parameter :: fmterrasym = &
      "(a,' **',a,'** PRODUCES AN ASYMMETRIC COEFFICIENT MATRIX, BUT THE       &
      &CONJUGATE GRADIENT METHOD WAS SELECTED. USE BICGSTAB INSTEAD. ')"
    !
    ! identify package and initialize.
    WRITE (iout, 1) this%iu
00001 FORMAT(1x, /1x, 'IMS -- ITERATIVE MODEL SOLUTION PACKAGE, VERSION 6', &
           ', 4/28/2017', /, 9x, 'INPUT READ FROM UNIT', i5)
    !
    ! -- initialize
    i = 1
    ifdparam = 1
    npp = 0
    mxvl = 0
    !
    ! -- get options block
    call this%parser%GetBlock('OPTIONS', isfound, ierr, &
                              supportopenclose=.true., blockrequired=.false.)
    !
    ! -- parse options block if detected
    if (isfound) then
      write (iout, '(/1x,a)') 'PROCESSING IMS OPTIONS'
      do
        call this%parser%GetNextLine(endofblock)
        if (endofblock) exit
        call this%parser%GetStringCaps(keyword)
        select case (keyword)
        case ('PRINT_OPTION')
          call this%parser%GetStringCaps(keyword)
          if (keyword .eq. 'NONE') then
            this%iprims = 0
          else if (keyword .eq. 'SUMMARY') then
            this%iprims = 1
          else if (keyword .eq. 'ALL') then
            this%iprims = 2
          else
            write (errmsg, '(3a)') &
              'Unknown IMS print option (', trim(keyword), ').'
            call store_error(errmsg)
          end if
        case ('COMPLEXITY')
          call this%parser%GetStringCaps(keyword)
          if (keyword .eq. 'SIMPLE') then
            ifdparam = 1
            WRITE (iout, 21)
          else if (keyword .eq. 'MODERATE') then
            ifdparam = 2
            WRITE (iout, 23)
          else if (keyword .eq. 'COMPLEX') then
            ifdparam = 3
            WRITE (iout, 25)
          else
            write (errmsg, '(3a)') &
              'Unknown IMS COMPLEXITY option (', trim(keyword), ').'
            call store_error(errmsg)
          end if
        case ('CSV_OUTER_OUTPUT')
          call this%parser%GetStringCaps(keyword)
          if (keyword == 'FILEOUT') then
            call this%parser%GetString(fname)
            if (nr_procs > 1) then
              call append_processor_id(fname, proc_id)
            end if
            this%icsvouterout = getunit()
            call openfile(this%icsvouterout, iout, fname, 'CSV_OUTER_OUTPUT', &
                          filstat_opt='REPLACE')
            write (iout, fmtcsvout) trim(fname), this%icsvouterout
          else
            write (errmsg, '(a)') 'Optional CSV_OUTER_OUTPUT '// &
              'keyword must be followed by FILEOUT'
            call store_error(errmsg)
          end if
        case ('CSV_INNER_OUTPUT')
          call this%parser%GetStringCaps(keyword)
          if (keyword == 'FILEOUT') then
            call this%parser%GetString(fname)
            if (nr_procs > 1) then
              call append_processor_id(fname, proc_id)
            end if
            this%icsvinnerout = getunit()
            call openfile(this%icsvinnerout, iout, fname, 'CSV_INNER_OUTPUT', &
                          filstat_opt='REPLACE')
            write (iout, fmtcsvout) trim(fname), this%icsvinnerout
          else
            write (errmsg, '(a)') 'Optional CSV_INNER_OUTPUT '// &
              'keyword must be followed by FILEOUT'
            call store_error(errmsg)
          end if
        case ('NO_PTC')
          call this%parser%GetStringCaps(keyword)
          select case (keyword)
          case ('ALL')
            ival = 0
            msg = 'ALL'
          case ('FIRST')
            ival = -1
            msg = 'THE FIRST'
          case default
            ival = 0
            msg = 'ALL'
          end select
          this%iallowptc = ival
          write (iout, '(3x,A)') 'PSEUDO-TRANSIENT CONTINUATION DISABLED FOR'// &
            ' '//trim(adjustl(msg))//' STRESS-PERIOD(S)'
        case ('ATS_OUTER_MAXIMUM_FRACTION')
          rval = this%parser%GetDouble()
          if (rval < dzero .or. rval > dhalf) then
            write (errmsg, '(a,G0)') 'Value for ATS_OUTER_MAXIMUM_FRAC must be &
              &between 0 and 0.5.  Found ', rval
            call store_error(errmsg)
          end if
          this%atsfrac = rval
          write (iout, '(3x,A,G0)') 'ADAPTIVE TIME STEP SETTING FOUND.  FRACTION &
            &OF OUTER MAXIMUM USED TO INCREASE OR DECREASE TIME STEP SIZE IS ',&
            &this%atsfrac
          !
          ! -- DEPRECATED OPTIONS
        case ('CSV_OUTPUT')
          call this%parser%GetStringCaps(keyword)
          if (keyword == 'FILEOUT') then
            call this%parser%GetString(fname)
            this%icsvouterout = getunit()
            call openfile(this%icsvouterout, iout, fname, 'CSV_OUTPUT', &
                          filstat_opt='REPLACE')
            write (iout, fmtcsvout) trim(fname), this%icsvouterout
            !
            ! -- create warning message
            write (warnmsg, '(a)') &
              'OUTER ITERATION INFORMATION WILL BE SAVED TO '//trim(fname)
            !
            ! -- create deprecation warning
            call deprecation_warning('OPTIONS', 'CSV_OUTPUT', '6.1.1', &
                                     warnmsg, this%parser%GetUnit())
          else
            write (errmsg, '(a)') 'Optional CSV_OUTPUT '// &
              'keyword must be followed by FILEOUT'
            call store_error(errmsg)
          end if
          !
          ! -- right now these are options that are only available in the
          !    development version and are not included in the documentation.
          !    These options are only available when IDEVELOPMODE in
          !    constants module is set to 1
        case ('DEV_PTC')
          call this%parser%DevOpt()
          this%iallowptc = 1
          write (iout, '(1x,A)') 'PSEUDO-TRANSIENT CONTINUATION ENABLED'
        case ('DEV_PTC_OUTPUT')
          call this%parser%DevOpt()
          this%iallowptc = 1
          call this%parser%GetStringCaps(keyword)
          if (keyword == 'FILEOUT') then
            call this%parser%GetString(fname)
            if (nr_procs > 1) then
              call append_processor_id(fname, proc_id)
            end if
            this%iptcout = getunit()
            call openfile(this%iptcout, iout, fname, 'PTC-OUT', &
                          filstat_opt='REPLACE')
            write (iout, fmtptcout) trim(fname), this%iptcout
          else
            write (errmsg, '(a)') &
              'Optional PTC_OUTPUT keyword must be followed by FILEOUT'
            call store_error(errmsg)
          end if
        case ('DEV_PTC_OPTION')
          call this%parser%DevOpt()
          this%iallowptc = 1
          this%iptcopt = 1
          write (iout, '(1x,A)') &
            'PSEUDO-TRANSIENT CONTINUATION USES BNORM AND L2NORM TO '// &
            'SET INITIAL VALUE'
        case ('DEV_PTC_EXPONENT')
          call this%parser%DevOpt()
          rval = this%parser%GetDouble()
          if (rval < dzero) then
            write (errmsg, '(a)') 'PTC_EXPONENT must be > 0.'
            call store_error(errmsg)
          else
            this%iallowptc = 1
            this%ptcexp = rval
            write (iout, '(1x,A,1x,g15.7)') &
              'PSEUDO-TRANSIENT CONTINUATION EXPONENT', this%ptcexp
          end if
        case ('DEV_PTC_DEL0')
          call this%parser%DevOpt()
          rval = this%parser%GetDouble()
          if (rval < dzero) then
            write (errmsg, '(a)') 'IMS sln_ar: PTC_DEL0 must be > 0.'
            call store_error(errmsg)
          else
            this%iallowptc = 1
            this%ptcdel0 = rval
            write (iout, '(1x,A,1x,g15.7)') &
              'PSEUDO-TRANSIENT CONTINUATION INITIAL TIMESTEP', this%ptcdel0
          end if
        case default
          write (errmsg, '(a,2(1x,a))') &
            'Unknown IMS option  (', trim(keyword), ').'
          call store_error(errmsg)
        end select
      end do
      write (iout, '(1x,a)') 'END OF IMS OPTIONS'
    else
      write (iout, '(1x,a)') 'NO IMS OPTION BLOCK DETECTED.'
    end if
 
00021 FORMAT(1x, 'SIMPLE OPTION:', /, &
           1x, 'DEFAULT SOLVER INPUT VALUES FOR FAST SOLUTIONS')
00023 FORMAT(1x, 'MODERATE OPTION:', /, 1x, 'DEFAULT SOLVER', &
           ' INPUT VALUES REFLECT MODERATELY NONLINEAR MODEL')
00025 FORMAT(1x, 'COMPLEX OPTION:', /, 1x, 'DEFAULT SOLVER', &
           ' INPUT VALUES REFLECT STRONGLY NONLINEAR MODEL')
 
    !-------READ NONLINEAR ITERATION PARAMETERS AND LINEAR SOLVER SELECTION INDEX
    ! -- set default nonlinear parameters
    call this%sln_setouter(ifdparam)
    !
    ! -- get NONLINEAR block
    call this%parser%GetBlock('NONLINEAR', isfound, ierr, &
                              supportopenclose=.true., blockrequired=.false.)
    !
    ! -- parse NONLINEAR block if detected
    if (isfound) then
      write (iout, '(/1x,a)') 'PROCESSING IMS NONLINEAR'
      do
        call this%parser%GetNextLine(endofblock)
        if (endofblock) exit
        call this%parser%GetStringCaps(keyword)
        ! -- parse keyword
        select case (keyword)
        case ('OUTER_DVCLOSE')
          this%dvclose = this%parser%GetDouble()
        case ('OUTER_MAXIMUM')
          this%mxiter = this%parser%GetInteger()
        case ('UNDER_RELAXATION')
          call this%parser%GetStringCaps(keyword)
          ival = 0
          if (keyword == 'NONE') then
            ival = 0
          else if (keyword == 'SIMPLE') then
            ival = 1
          else if (keyword == 'COOLEY') then
            ival = 2
          else if (keyword == 'DBD') then
            ival = 3
          else
            write (errmsg, '(3a)') &
              'Unknown UNDER_RELAXATION specified (', trim(keyword), ').'
            call store_error(errmsg)
          end if
          this%nonmeth = ival
        case ('LINEAR_SOLVER')
          call this%parser%GetStringCaps(keyword)
          ival = ims_solver
          if (keyword .eq. 'DEFAULT' .or. &
              keyword .eq. 'LINEAR') then
            ival = ims_solver
          else
            write (errmsg, '(3a)') &
              'Unknown LINEAR_SOLVER specified (', trim(keyword), ').'
            call store_error(errmsg)
          end if
          this%linsolver = ival
        case ('UNDER_RELAXATION_THETA')
          this%theta = this%parser%GetDouble()
        case ('UNDER_RELAXATION_KAPPA')
          this%akappa = this%parser%GetDouble()
        case ('UNDER_RELAXATION_GAMMA')
          this%gamma = this%parser%GetDouble()
        case ('UNDER_RELAXATION_MOMENTUM')
          this%amomentum = this%parser%GetDouble()
        case ('BACKTRACKING_NUMBER')
          this%numtrack = this%parser%GetInteger()
          IF (this%numtrack > 0) this%ibflag = 1
        case ('BACKTRACKING_TOLERANCE')
          this%btol = this%parser%GetDouble()
        case ('BACKTRACKING_REDUCTION_FACTOR')
          this%breduc = this%parser%GetDouble()
        case ('BACKTRACKING_RESIDUAL_LIMIT')
          this%res_lim = this%parser%GetDouble()
          !
          ! -- deprecated variables
        case ('OUTER_HCLOSE')
          this%dvclose = this%parser%GetDouble()
          !
          ! -- create warning message
          write (warnmsg, '(a)') &
            'SETTING OUTER_DVCLOSE TO OUTER_HCLOSE VALUE'
          !
          ! -- create deprecation warning
          call deprecation_warning('NONLINEAR', 'OUTER_HCLOSE', '6.1.1', &
                                   warnmsg, this%parser%GetUnit())
        case ('OUTER_RCLOSEBND')
          !
          ! -- create warning message
          write (warnmsg, '(a)') &
            'OUTER_DVCLOSE IS USED TO EVALUATE PACKAGE CONVERGENCE'
          !
          ! -- create deprecation warning
          call deprecation_warning('NONLINEAR', 'OUTER_RCLOSEBND', '6.1.1', &
                                   warnmsg, this%parser%GetUnit())
        case default
          write (errmsg, '(3a)') &
            'Unknown IMS NONLINEAR keyword (', trim(keyword), ').'
          call store_error(errmsg)
        end select
      end do
      write (iout, '(1x,a)') 'END OF IMS NONLINEAR DATA'
    else
      if (ifdparam .EQ. 0) then
        write (errmsg, '(a)') 'NO IMS NONLINEAR block detected.'
        call store_error(errmsg)
      end if
    end if
    !
    if (this%theta < dem3) then
      this%theta = dem3
    end if
    !
    ! -- backtracking should only be used if this%nonmeth > 0
    if (this%nonmeth < 1) then
      this%ibflag = 0
    end if
    !
    ! -- check that MXITER is greater than zero
    if (this%mxiter <= 0) then
      write (errmsg, '(a)') 'Outer iteration number must be > 0.'
      call store_error(errmsg)
    END IF
    !
    ! -- write under-relaxation option
    if (this%nonmeth > 0) then
      WRITE (iout, *) '**UNDER-RELAXATION WILL BE USED***'
      WRITE (iout, *)
    elseif (this%nonmeth == 0) then
      WRITE (iout, *) '***UNDER-RELAXATION WILL NOT BE USED***'
      WRITE (iout, *)
    else
      WRITE (errmsg, '(a)') &
        'Incorrect value for variable NONMETH was specified.'
      call store_error(errmsg)
    end if
    !
    ! -- ensure gamma is > 0 for simple
    if (this%nonmeth == 1) then
      if (this%gamma == 0) then
        WRITE (errmsg, '(a)') &
          'GAMMA must be greater than zero if SIMPLE under-relaxation is used.'
        call store_error(errmsg)
      end if
    end if
 
    if (this%solver_mode == 'PETSC') then
      this%linsolver = petsc_solver
    end if
 
    ! configure linear settings
    call this%linear_settings%init(this%memory_path)
    call this%linear_settings%preset_config(ifdparam)
    call this%linear_settings%read_from_file(this%parser, iout)
    !
    if (this%linear_settings%ilinmeth == cg_method) then
      this%isymmetric = 1
    end if
    !
    ! -- call secondary subroutine to initialize and read linear
    !    solver parameters IMSLINEAR solver
    if (this%solver_mode == "IMS") then
      allocate (this%imslinear)
      WRITE (iout, *) '***IMS LINEAR SOLVER WILL BE USED***'
      call this%imslinear%imslinear_allocate(this%name, iout, this%iprims, &
                                             this%mxiter, this%neq, &
                                             this%system_matrix, this%rhs, &
                                             this%x, this%linear_settings)
      !
      ! -- petsc linear solver flag
    else if (this%solver_mode == "PETSC") then
      call this%linear_solver%initialize(this%system_matrix, &
                                         this%linear_settings, &
                                         this%cnvg_summary)
      !
      ! -- incorrect linear solver flag
    else
      write (errmsg, '(a)') &
        'Incorrect value for linear solution method specified.'
      call store_error(errmsg)
    end if
    !
    ! -- write message about matrix symmetry
    if (this%isymmetric == 1) then
      write (iout, '(1x,a,/)') 'A symmetric matrix will be solved'
    else
      write (iout, '(1x,a,/)') 'An asymmetric matrix will be solved'
    end if
    !
    ! -- If CG, then go through each model and each exchange and check
    !    for asymmetry
    if (this%isymmetric == 1) then
      !
      ! -- Models
      do i = 1, this%modellist%Count()
        mp => getnumericalmodelfromlist(this%modellist, i)
        if (mp%get_iasym() /= 0) then
          write (errmsg, fmterrasym) 'MODEL', trim(adjustl(mp%name))
          call store_error(errmsg)
        end if
      end do
      !
      ! -- Exchanges
      do i = 1, this%exchangelist%Count()
        cp => getnumericalexchangefromlist(this%exchangelist, i)
        if (cp%get_iasym() /= 0) then
          write (errmsg, fmterrasym) 'EXCHANGE', trim(adjustl(cp%name))
          call store_error(errmsg)
        end if
      end do
      !
    end if
 
    !
    ! determine if the x and rhs should be scaled
    this%idv_scale = this%sln_get_idvscale()
 
    if (this%idv_scale > 0) then
      write (iout, '(2(1x,a,/),1x,a,/,6x,a,/)') &
        'X and RHS will be scaled to avoid very large positive or negative', &
        'dependent variable values in the model IMS package.', &
        'NOTE: Specified outer and inner DVCLOSE values in the model IMS &
        &package', 'will be relative closure criteria.'
    else if (this%idv_scale < 0) then
      write (errmsg, '(2(a,1x))') &
        'dependent_variable_scaling must be specified for all models in', &
        'the solution and can only be used with GWT and GWE models. '
      call store_error(errmsg)
    end if
    !
    !
    ! -- write solver data to output file
    !
    ! -- non-linear solver data
    WRITE (iout, 9002) this%dvclose, this%mxiter, &
      this%iprims, this%nonmeth, this%linsolver
    !
    ! -- standard outer iteration formats
9002 FORMAT(1x, 'OUTER ITERATION CONVERGENCE CRITERION    (DVCLOSE) = ', e15.6, &
           /1x, 'MAXIMUM NUMBER OF OUTER ITERATIONS        (MXITER) = ', i0, &
           /1x, 'SOLVER PRINTOUT INDEX                     (IPRIMS) = ', i0, &
           /1x, 'NONLINEAR ITERATION METHOD            (NONLINMETH) = ', i0, &
           /1x, 'LINEAR SOLUTION METHOD                   (LINMETH) = ', i0)
    !
    if (this%nonmeth == 1) then ! simple
      write (iout, 9003) this%gamma
    else if (this%nonmeth == 2) then ! cooley
      write (iout, 9004) this%gamma
    else if (this%nonmeth == 3) then ! delta bar delta
      write (iout, 9005) this%theta, this%akappa, this%gamma, this%amomentum
    end if
    !
    ! -- write backtracking information
    if (this%numtrack /= 0) write (iout, 9006) this%numtrack, this%btol, &
      this%breduc, this%res_lim
    !
    ! -- under-relaxation formats (simple, cooley, dbd)
9003 FORMAT(1x, 'UNDER-RELAXATION FACTOR                    (GAMMA) = ', e15.6)
9004 FORMAT(1x, 'UNDER-RELAXATION PREVIOUS HISTORY FACTOR   (GAMMA) = ', e15.6)
9005 FORMAT(1x, 'UNDER-RELAXATION WEIGHT REDUCTION FACTOR   (THETA) = ', e15.6, &
           /1x, 'UNDER-RELAXATION WEIGHT INCREASE INCREMENT (KAPPA) = ', e15.6, &
           /1x, 'UNDER-RELAXATION PREVIOUS HISTORY FACTOR   (GAMMA) = ', e15.6, &
           /1x, 'UNDER-RELAXATION MOMENTUM TERM         (AMOMENTUM) = ', e15.6)
    !
    ! -- backtracking formats
9006 FORMAT(1x, 'MAXIMUM NUMBER OF BACKTRACKS            (NUMTRACK) = ', i0, &
           /1x, 'BACKTRACKING TOLERANCE FACTOR               (BTOL) = ', e15.6, &
           /1x, 'BACKTRACKING REDUCTION FACTOR             (BREDUC) = ', e15.6, &
           /1x, 'BACKTRACKING RESIDUAL LIMIT              (RES_LIM) = ', e15.6)
    !
    ! -- linear solver data
    if (this%linsolver == ims_solver) then
      call this%imslinear%imslinear_summary(this%mxiter)
    else
      call this%linear_solver%print_summary()
    end if
 
    ! -- write summary of solver error messages
    ierr = count_errors()
    if (ierr > 0) then
      call this%parser%StoreErrorUnit()
    end if
    !
    ! reallocate space for nonlinear arrays and initialize
    call mem_reallocate(this%hncg, this%mxiter, 'HNCG', this%name)
    call mem_reallocate(this%lrch, 3, this%mxiter, 'LRCH', this%name)
 
    ! delta-bar-delta under-relaxation
    if (this%nonmeth == 3) then
      call mem_reallocate(this%wsave, this%neq, 'WSAVE', this%name)
      call mem_reallocate(this%hchold, this%neq, 'HCHOLD', this%name)
      call mem_reallocate(this%deold, this%neq, 'DEOLD', this%name)
      do i = 1, this%neq
        this%wsave(i) = dzero
        this%hchold(i) = dzero
        this%deold(i) = dzero
      end do
    end if
    this%hncg = dzero
    this%lrch = 0
 
    ! allocate space for saving solver convergence history
    if (this%iprims == 2 .or. this%icsvinnerout > 0) then
      this%nitermax = this%linear_settings%iter1 * this%mxiter
    else
      this%nitermax = 1
    end if
 
    allocate (this%caccel(this%nitermax))
 
    !
    ! -- resize convergence report
    call this%cnvg_summary%reinit(this%nitermax)
    !
    ! -- check for numerical solution errors
    ierr = count_errors()
    if (ierr > 0) then
      call this%parser%StoreErrorUnit()
    end if
    !
    ! -- close ims input file
    call this%parser%Clear()

Here is the call graph for this function:

sln_backtracking()

subroutine numericalsolutionmodule::sln_backtracking ( class(numericalsolutiontype), intent(inout)  this, class(numericalmodeltype), pointer  mp, class(numericalexchangetype), pointer  cp, integer(i4b), intent(in)  kiter  )

private

Perform backtracking on the solution in the maximum number of backtrack iterations (nbtrack) is greater than 0 and the backtracking criteria are exceeded.

Parameters

[in,out]	this	NumericalSolutionType instance
	mp	model pointer (currently null())
	cp	exchange pointer (currently null())
[in]	kiter	Picard iteration number

Definition at line 2718 of file NumericalSolution.f90.

    ! -- dummy variables
    class(NumericalSolutionType), intent(inout) :: this !< NumericalSolutionType instance
    class(NumericalModelType), pointer :: mp !< model pointer (currently null())
    class(NumericalExchangeType), pointer :: cp !< exchange pointer (currently null())
    integer(I4B), intent(in) :: kiter !< Picard iteration number
    ! -- local variables
    character(len=7) :: cmsg
    integer(I4B) :: nb
    integer(I4B) :: btflag
    integer(I4B) :: ibflag
    integer(I4B) :: ibtcnt
    real(DP) :: resin
    !
    ! -- initialize local variables
    ibflag = 0
 
    !
    ! -- refill amat and rhs with standard conductance
    call this%sln_buildsystem(kiter, inewton=0)
 
    !
    ! -- calculate initial l2 norm
    if (kiter == 1) then
      call this%sln_l2norm(this%res_prev)
      resin = this%res_prev
      ibflag = 0
    else
      call this%sln_l2norm(this%res_new)
      resin = this%res_new
    end if
    ibtcnt = 0
    if (kiter > 1) then
      if (this%res_new > this%res_prev * this%btol) then
        !
        ! -- iterate until backtracking complete
        btloop: do nb = 1, this%numtrack
          !
          ! -- backtrack the dependent variable
          call this%sln_backtracking_xupdate(btflag)
          !
          ! -- dependent-variable change less than dvclose
          if (btflag == 0) then
            ibflag = 4
            exit btloop
          end if
          !
          ibtcnt = nb
 
          ! recalculate linear system (amat and rhs)
          call this%sln_buildsystem(kiter, inewton=0)
 
          !
          ! -- calculate updated l2norm
          call this%sln_l2norm(this%res_new)
          !
          ! -- evaluate if back tracking can be terminated
          if (nb == this%numtrack) then
            ibflag = 2
            exit btloop
          end if
          if (this%res_new < this%res_prev * this%btol) then
            ibflag = 1
            exit btloop
          end if
          if (this%res_new < this%res_lim) then
            exit btloop
          end if
        end do btloop
      end if
      ! -- save new residual
      this%res_prev = this%res_new
    end if
    !
    ! -- write back backtracking results
    if (this%iprims > 0) then
      if (ibtcnt > 0) then
        cmsg = ' '
      else
        cmsg = '*'
      end if
      !
      ! -- add data to outertab
      call this%outertab%add_term('Backtracking')
      call this%outertab%add_term(kiter)
      call this%outertab%add_term(' ')
      if (this%numtrack > 0) then
        call this%outertab%add_term(ibflag)
        call this%outertab%add_term(ibtcnt)
        call this%outertab%add_term(resin)
        call this%outertab%add_term(this%res_prev)
      end if
      call this%outertab%add_term(' ')
      call this%outertab%add_term(cmsg)
      call this%outertab%add_term(' ')
    end if

sln_backtracking_xupdate()

subroutine numericalsolutionmodule::sln_backtracking_xupdate ( class(numericalsolutiontype), intent(inout)  this, integer(i4b), intent(inout)  bt_flag  )

private

Backtracking update of the dependent variable if the calculated backtracking update exceeds the dependent variable closure criteria.

Parameters

[in,out]	this	NumericalSolutionType instance
[in,out]	bt_flag	backtracking flag (1) backtracking performed (0) backtracking not performed

Definition at line 2822 of file NumericalSolution.f90.

    ! -- dummy variables
    class(NumericalSolutionType), intent(inout) :: this !< NumericalSolutionType instance
    integer(I4B), intent(inout) :: bt_flag !< backtracking flag (1) backtracking performed (0) backtracking not performed
 
    bt_flag = this%get_backtracking_flag()
 
    ! perform backtracking if ...
    if (bt_flag > 0) then
      call this%apply_backtracking()
    end if
 

sln_buildsystem()

subroutine numericalsolutionmodule::sln_buildsystem	(	class(numericalsolutiontype)	this,
		integer(i4b), intent(in)	kiter,
		integer(i4b), intent(in)	inewton
	)

Definition at line 1979 of file NumericalSolution.f90.

    class(NumericalSolutionType) :: this
    integer(I4B), intent(in) :: kiter
    integer(I4B), intent(in) :: inewton
    ! local
    integer(I4B) :: im, ic
    class(NumericalModelType), pointer :: mp
    class(NumericalExchangeType), pointer :: cp
    !
    ! -- Set amat and rhs to zero
    call this%sln_reset()
 
    ! reset models
    do im = 1, this%modellist%Count()
      mp => getnumericalmodelfromlist(this%modellist, im)
      call mp%model_reset()
    end do
 
    ! synchronize for CF
    call this%synchronize(stg_bfr_exg_cf, this%synchronize_ctx)
 
    !
    ! -- Calculate the matrix terms for each exchange
    do ic = 1, this%exchangelist%Count()
      cp => getnumericalexchangefromlist(this%exchangelist, ic)
      call cp%exg_cf(kiter)
    end do
    !
    ! -- Calculate the matrix terms for each model
    do im = 1, this%modellist%Count()
      mp => getnumericalmodelfromlist(this%modellist, im)
      call mp%model_cf(kiter)
    end do
 
    ! synchronize for FC
    call this%synchronize(stg_bfr_exg_fc, this%synchronize_ctx)
 
    !
    ! -- Add exchange coefficients to the solution
    do ic = 1, this%exchangelist%Count()
      cp => getnumericalexchangefromlist(this%exchangelist, ic)
      call cp%exg_fc(kiter, this%system_matrix, this%rhs, inewton)
    end do
    !
    ! -- Add model coefficients to the solution
    do im = 1, this%modellist%Count()
      mp => getnumericalmodelfromlist(this%modellist, im)
      call mp%model_fc(kiter, this%system_matrix, inewton)
    end do
 

Here is the call graph for this function:

sln_ca()

subroutine numericalsolutionmodule::sln_ca	(	class(numericalsolutiontype)	this,
		integer(i4b), intent(inout)	isgcnvg,
		integer(i4b), intent(in)	isuppress_output
	)

Solve the models in this solution for kper and kstp.

Parameters

	this	NumericalSolutionType instance
[in,out]	isgcnvg	solution group convergence flag
[in]	isuppress_output	flag for suppressing output

Definition at line 1317 of file NumericalSolution.f90.

    ! -- dummy variables
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    integer(I4B), intent(inout) :: isgcnvg !< solution group convergence flag
    integer(I4B), intent(in) :: isuppress_output !< flag for suppressing output
    ! -- local variables
    class(NumericalModelType), pointer :: mp => null()
    character(len=LINELENGTH) :: line
    character(len=LINELENGTH) :: fmt
    integer(I4B) :: im
    integer(I4B) :: kiter ! non-linear iteration counter
 
    ! advance the models, exchanges, and solution
    call this%prepareSolve()
 
    select case (isim_mode)
    case (mvalidate)
      line = 'mode="validation" -- Skipping matrix assembly and solution.'
      fmt = "(/,1x,a,/)"
      do im = 1, this%modellist%Count()
        mp => getnumericalmodelfromlist(this%modellist, im)
        call mp%model_message(line, fmt=fmt)
      end do
    case (mnormal)
      ! nonlinear iteration loop for this solution
      outerloop: do kiter = 1, this%mxiter
 
        ! perform a single iteration
        call this%solve(kiter)
 
        ! exit if converged
        if (this%icnvg == 1) then
          exit outerloop
        end if
 
      end do outerloop
 
      ! finish up, write convergence info, CSV file, budgets and flows, ...
      call this%finalizeSolve(kiter, isgcnvg, isuppress_output)
    end select

Here is the call graph for this function:

sln_calc_ptc()

subroutine numericalsolutionmodule::sln_calc_ptc ( class(numericalsolutiontype)  this, integer(i4b)  iptc, real(dp)  ptcf  )

private

Parameters

this	NumericalSolutionType instance
iptc	PTC (1) or not (0)
ptcf	the PTC factor calculated

Definition at line 2997 of file NumericalSolution.f90.

    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    integer(I4B) :: iptc !< PTC (1) or not (0)
    real(DP) :: ptcf !< the PTC factor calculated
    ! local
    integer(I4B) :: im
    class(NumericalModelType), pointer :: mp
    class(VectorBaseType), pointer :: vec_resid
 
    iptc = 0
    ptcf = dzero
 
    ! calc. residual vector
    vec_resid => this%system_matrix%create_vec(this%neq)
    call this%sln_calc_residual(vec_resid)
 
    ! determine ptc
    do im = 1, this%modellist%Count()
      mp => getnumericalmodelfromlist(this%modellist, im)
      call mp%model_ptc(vec_resid, iptc, ptcf)
    end do
 
    ! clean up temp. vector
    call vec_resid%destroy()
    deallocate (vec_resid)
 

Here is the call graph for this function:

sln_calc_residual()

subroutine numericalsolutionmodule::sln_calc_residual ( class(numericalsolutiontype)  this, class(vectorbasetype), pointer  vec_resid  )

private

Parameters

this	NumericalSolutionType instance
vec_resid	the residual vector

Definition at line 3027 of file NumericalSolution.f90.

    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    class(VectorBaseType), pointer :: vec_resid !< the residual vector
    ! local
    integer(I4B) :: n
 
    call this%system_matrix%multiply(this%vec_x, vec_resid) ! r = A*x
 
    call vec_resid%axpy(-1.0_dp, this%vec_rhs) ! r = r - b
 
    do n = 1, this%neq
      if (this%active(n) < 1) then
        call vec_resid%set_value_local(n, 0.0_dp) ! r_i = 0 if inactive
      end if
    end do
 

sln_calcdx()

subroutine numericalsolutionmodule::sln_calcdx ( class(numericalsolutiontype), intent(inout)  this, integer(i4b), intent(in)  neq, integer(i4b), dimension(neq), intent(in)  active, real(dp), dimension(neq), intent(in)  x, real(dp), dimension(neq), intent(in)  xtemp, real(dp), dimension(neq), intent(inout)  dx  )

private

Calculate the dependent-variable change for every cell.

Parameters

[in,out]	this	NumericalSolutionType instance
[in]	neq	number of equations
[in]	active	active cell flag (1)
[in]	x	current dependent-variable
[in]	xtemp	previous dependent-variable
[in,out]	dx	dependent-variable change

Definition at line 2973 of file NumericalSolution.f90.

    ! -- dummy variables
    class(NumericalSolutionType), intent(inout) :: this !< NumericalSolutionType instance
    integer(I4B), intent(in) :: neq !< number of equations
    integer(I4B), dimension(neq), intent(in) :: active !< active cell flag (1)
    real(DP), dimension(neq), intent(in) :: x !< current dependent-variable
    real(DP), dimension(neq), intent(in) :: xtemp !< previous dependent-variable
    real(DP), dimension(neq), intent(inout) :: dx !< dependent-variable change
    ! -- local
    integer(I4B) :: n
    !
    ! -- calculate dependent-variable change
    do n = 1, neq
      ! -- skip inactive nodes
      if (active(n) < 1) then
        dx(n) = dzero
      else
        dx(n) = x(n) - xtemp(n)
      end if
    end do

sln_connect()

subroutine numericalsolutionmodule::sln_connect ( class(numericalsolutiontype)  this )

private

Assign solution connections. This is the main workhorse method for a solution. The method goes through all the models and all the connections and builds up the sparse matrix. Steps are (1) add internal model connections, (2) add cross terms, (3) allocate solution arrays, (4) create mapping arrays, and (5) fill cross term values if necessary.

Parameters

this	NumericalSolutionType instance

Definition at line 2374 of file NumericalSolution.f90.

    ! -- modules
    use memorymanagermodule, only: mem_allocate
    ! -- dummy variables
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    ! -- local variables
    class(NumericalModelType), pointer :: mp => null()
    class(NumericalExchangeType), pointer :: cp => null()
    integer(I4B) :: im
    integer(I4B) :: ic
    !
    ! -- Add internal model connections to sparse
    do im = 1, this%modellist%Count()
      mp => getnumericalmodelfromlist(this%modellist, im)
      call mp%model_ac(this%sparse)
    end do
    !
    ! -- synchronize before AC
    call this%synchronize(stg_bfr_exg_ac, this%synchronize_ctx)
    !
    ! -- Add the cross terms to sparse
    do ic = 1, this%exchangelist%Count()
      cp => getnumericalexchangefromlist(this%exchangelist, ic)
      call cp%exg_ac(this%sparse)
    end do
    !
    ! -- The number of non-zero array values are now known so
    ! -- ia and ja can be created from sparse. then destroy sparse
    call this%sparse%sort()
    call this%system_matrix%init(this%sparse, this%name)
    call this%sparse%destroy()
    !
    ! -- Create mapping arrays for each model.  Mapping assumes
    ! -- that each row has the diagonal in the first position,
    ! -- however, rows do not need to be sorted.
    do im = 1, this%modellist%Count()
      mp => getnumericalmodelfromlist(this%modellist, im)
      call mp%model_mc(this%system_matrix)
    end do
    !
    ! -- Create arrays for mapping exchange connections to global solution
    do ic = 1, this%exchangelist%Count()
      cp => getnumericalexchangefromlist(this%exchangelist, ic)
      call cp%exg_mc(this%system_matrix)
    end do

Here is the call graph for this function:

sln_da()

subroutine numericalsolutionmodule::sln_da

(

class(numericalsolutiontype)

this

)

Deallocate a solution.

Parameters

this	NumericalSolutionType instance

Definition at line 1199 of file NumericalSolution.f90.

    ! -- modules
    use memorymanagermodule, only: mem_deallocate
    ! -- dummy variables
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    !
    ! -- IMSLinearModule
    if (this%linsolver == ims_solver) then
      call this%imslinear%imslinear_da()
      deallocate (this%imslinear)
    end if
    !
    ! -- lists
    call this%modellist%Clear()
    call this%exchangelist%Clear()
    deallocate (this%modellist)
    deallocate (this%exchangelist)
 
    call this%system_matrix%destroy()
    deallocate (this%system_matrix)
    call this%vec_x%destroy()
    deallocate (this%vec_x)
    call this%vec_rhs%destroy()
    deallocate (this%vec_rhs)
 
    !
    ! -- character arrays
    deallocate (this%caccel)
    !
    ! -- inner iteration table object
    if (associated(this%innertab)) then
      call this%innertab%table_da()
      deallocate (this%innertab)
      nullify (this%innertab)
    end if
    !
    ! -- outer iteration table object
    if (associated(this%outertab)) then
      call this%outertab%table_da()
      deallocate (this%outertab)
      nullify (this%outertab)
    end if
    !
    ! -- arrays
    call mem_deallocate(this%active)
    call mem_deallocate(this%xtemp)
    call mem_deallocate(this%dxold)
    call mem_deallocate(this%hncg)
    call mem_deallocate(this%lrch)
    call mem_deallocate(this%wsave)
    call mem_deallocate(this%hchold)
    call mem_deallocate(this%deold)
    call mem_deallocate(this%convmodstart)
    !
    ! -- convergence report
    call this%cnvg_summary%destroy()
    deallocate (this%cnvg_summary)
    !
    ! -- linear solver
    call this%linear_solver%destroy()
    deallocate (this%linear_solver)
    !
    ! -- linear solver settings
    call this%linear_settings%destroy()
    deallocate (this%linear_settings)
    !
    ! -- Scalars
    call mem_deallocate(this%id)
    call mem_deallocate(this%iu)
    call mem_deallocate(this%ttform)
    call mem_deallocate(this%ttsoln)
    call mem_deallocate(this%isymmetric)
    call mem_deallocate(this%neq)
    call mem_deallocate(this%matrix_offset)
    call mem_deallocate(this%dvclose)
    call mem_deallocate(this%bigchold)
    call mem_deallocate(this%bigch)
    call mem_deallocate(this%relaxold)
    call mem_deallocate(this%res_prev)
    call mem_deallocate(this%res_new)
    call mem_deallocate(this%icnvg)
    call mem_deallocate(this%itertot_timestep)
    call mem_deallocate(this%iouttot_timestep)
    call mem_deallocate(this%itertot_sim)
    call mem_deallocate(this%mxiter)
    call mem_deallocate(this%linsolver)
    call mem_deallocate(this%nonmeth)
    call mem_deallocate(this%iprims)
    call mem_deallocate(this%theta)
    call mem_deallocate(this%akappa)
    call mem_deallocate(this%gamma)
    call mem_deallocate(this%amomentum)
    call mem_deallocate(this%breduc)
    call mem_deallocate(this%btol)
    call mem_deallocate(this%res_lim)
    call mem_deallocate(this%numtrack)
    call mem_deallocate(this%ibflag)
    call mem_deallocate(this%icsvouterout)
    call mem_deallocate(this%icsvinnerout)
    call mem_deallocate(this%nitermax)
    call mem_deallocate(this%convnmod)
    call mem_deallocate(this%iallowptc)
    call mem_deallocate(this%iptcopt)
    call mem_deallocate(this%iptcout)
    call mem_deallocate(this%l2norm0)
    call mem_deallocate(this%ptcdel)
    call mem_deallocate(this%ptcdel0)
    call mem_deallocate(this%ptcexp)
    call mem_deallocate(this%atsfrac)
    call mem_deallocate(this%idv_scale)
    call mem_deallocate(this%dscale)

sln_df()

subroutine numericalsolutionmodule::sln_df

(

class(numericalsolutiontype)

this

)

Define a new solution. Must be called after the models and exchanges have been added to solution. The order of the steps is (1) Allocate neq and nja, (2) Assign model offsets and solution ids, (3) Allocate and initialize the solution arrays, (4) Point each model's x and rhs arrays, and (5) Initialize the sparsematrix instance

Parameters

this	NumericalSolutionType instance

Definition at line 438 of file NumericalSolution.f90.

    ! modules
    use memorymanagermodule, only: mem_allocate
    use simvariablesmodule, only: simulation_mode
    ! -- dummy variables
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    ! -- local variables
    class(NumericalModelType), pointer :: mp => null()
    integer(I4B) :: i
    integer(I4B), allocatable, dimension(:) :: rowmaxnnz
    integer(I4B) :: ncol, irow_start, irow_end
    integer(I4B) :: mod_offset
    !
    ! -- set sol id and determine nr. of equation in this solution
    do i = 1, this%modellist%Count()
      mp => getnumericalmodelfromlist(this%modellist, i)
      call mp%set_idsoln(this%id)
      this%neq = this%neq + mp%neq
    end do
    !
    ! -- set up the (possibly parallel) linear system
    if (simulation_mode == 'PARALLEL') then
      this%solver_mode = 'PETSC'
    else
      this%solver_mode = 'IMS'
    end if
    !
    ! -- allocate settings structure
    allocate (this%linear_settings)
    !
    ! -- create linear system matrix and compatible vectors
    this%linear_solver => create_linear_solver(this%solver_mode, this%name)
    this%system_matrix => this%linear_solver%create_matrix()
    this%vec_x => this%system_matrix%create_vec_mm(this%neq, 'X', &
                                                   this%memory_path)
    this%x => this%vec_x%get_array()
    this%vec_rhs => this%system_matrix%create_vec_mm(this%neq, 'RHS', &
                                                     this%memory_path)
    this%rhs => this%vec_rhs%get_array()
    !
    call this%vec_rhs%get_ownership_range(irow_start, irow_end)
    ncol = this%vec_rhs%get_size()
    !
    ! -- calculate and set offsets
    mod_offset = irow_start - 1
    this%matrix_offset = irow_start - 1
    do i = 1, this%modellist%Count()
      mp => getnumericalmodelfromlist(this%modellist, i)
      call mp%set_moffset(mod_offset)
      mod_offset = mod_offset + mp%neq
    end do
    !
    ! -- Allocate and initialize solution arrays
    call this%allocate_arrays()
    !
    ! -- Create convergence summary report
    allocate (this%cnvg_summary)
    call this%cnvg_summary%init(this%modellist%Count(), this%convmodstart, &
                                this%memory_path)
    !
    ! -- Go through each model and point x, ibound, and rhs to solution
    do i = 1, this%modellist%Count()
      mp => getnumericalmodelfromlist(this%modellist, i)
      call mp%set_xptr(this%x, this%matrix_offset, 'X', this%name)
      call mp%set_rhsptr(this%rhs, this%matrix_offset, 'RHS', this%name)
      call mp%set_iboundptr(this%active, this%matrix_offset, 'IBOUND', this%name)
    end do
    !
    ! -- Create the sparsematrix instance
    allocate (rowmaxnnz(this%neq))
    do i = 1, this%neq
      rowmaxnnz(i) = 4
    end do
    call this%sparse%init(this%neq, ncol, rowmaxnnz)
    this%sparse%offset = this%matrix_offset
    deallocate (rowmaxnnz)
    !
    ! -- Assign connections, fill ia/ja, map connections
    call this%sln_connect()
 
    ! add timers
    write (this%id_postfix, '(a,i0,a)') " (", this%id, ")"
    this%tmr_prep_solve = -1
    this%tmr_solve = -1
    this%tmr_final_solve = -1
    this%tmr_formulate = -1
    this%tmr_linsolve = -1
    this%tmr_flows = -1
    this%tmr_budgets = -1
 

Here is the call graph for this function:

sln_dt()

subroutine numericalsolutionmodule::sln_dt

(

class(numericalsolutiontype)

this

)

Calculate time step length.

Parameters

this	NumericalSolutionType instance

Definition at line 1098 of file NumericalSolution.f90.

    ! -- modules
    use tdismodule, only: kstp, kper, delt
    use adaptivetimestepmodule, only: ats_submit_delt
    use constantsmodule, only: dtwo, dthree
    ! -- dummy variables
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    ! -- local variables
    integer(I4B) :: idir
    real(DP) :: delt_temp
    real(DP) :: fact_lower
    real(DP) :: fact_upper
    !
    ! -- increase or decrease delt based on kiter fraction.  atsfrac should be
    !    a value of about 1/3.  If the number of outer iterations is less than
    !    1/3 of mxiter, then increase step size.  If the number of outer
    !    iterations is greater than 2/3 of mxiter, then decrease step size.
    if (this%atsfrac > dzero) then
      delt_temp = delt
      fact_lower = this%mxiter * this%atsfrac
      fact_upper = this%mxiter - fact_lower
      if (this%iouttot_timestep < int(fact_lower)) then
        ! -- increase delt according to tsfactats
        idir = 1
      else if (this%iouttot_timestep > int(fact_upper)) then
        ! -- decrease delt according to tsfactats
        idir = -1
      else
        ! -- do not change delt
        idir = 0
      end if
      !
      ! -- submit stable dt for upcoming step
      call ats_submit_delt(kstp, kper, delt_temp, this%memory_path, idir=idir)
    end if

Here is the call graph for this function:

sln_fp()

subroutine numericalsolutionmodule::sln_fp ( class(numericalsolutiontype)  this )

private

Finalize a solution.

Parameters

this	NumericalSolutionType instance

Definition at line 1177 of file NumericalSolution.f90.

    use simvariablesmodule, only: iout
    ! -- dummy variables
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    !
    ! -- write timer output
    if (idevelopmode == 1) then
      write (iout, '(//1x,a,1x,a,1x,a)') &
        'Solution', trim(adjustl(this%name)), 'summary'
      write (iout, "(1x,70('-'))")
      write (iout, '(1x,a,1x,g0,1x,a)') &
        'Total formulate time: ', this%ttform, 'seconds'
      write (iout, '(1x,a,1x,g0,1x,a,/)') &
        'Total solution time:  ', this%ttsoln, 'seconds'
    end if

sln_get_dxmax()

subroutine numericalsolutionmodule::sln_get_dxmax ( class(numericalsolutiontype), intent(inout)  this, real(dp), intent(inout)  hncg, integer(i4b), intent(inout)  lrch  )

private

Determine the maximum dependent-variable change at the end of a Picard iteration.

Parameters

[in,out]	this	NumericalSolutionType instance
[in,out]	hncg	maximum dependent-variable change
[in,out]	lrch	location of the maximum dependent-variable change

Definition at line 3183 of file NumericalSolution.f90.

    ! -- dummy variables
    class(NumericalSolutionType), intent(inout) :: this !< NumericalSolutionType instance
    real(DP), intent(inout) :: hncg !< maximum dependent-variable change
    integer(I4B), intent(inout) :: lrch !< location of the maximum dependent-variable change
    ! -- local variables
    integer(I4B) :: nb
    real(DP) :: bigch
    real(DP) :: abigch
    integer(I4B) :: n
    real(DP) :: hdif
    real(DP) :: ahdif
    !
    ! -- determine the maximum change
    nb = 0
    bigch = dzero
    abigch = dzero
    do n = 1, this%neq
      if (this%active(n) < 1) cycle
      hdif = this%x(n) - this%xtemp(n)
      ahdif = abs(hdif)
      if (ahdif > abigch) then
        bigch = hdif
        abigch = ahdif
        nb = n
      end if
    end do
    !
    !-----store maximum change value and location
    hncg = bigch
    lrch = nb

sln_get_idvscale()

integer(i4b) function numericalsolutionmodule::sln_get_idvscale ( class(numericalsolutiontype)  this )

private

Parameters

this	NumericalSolutionType instance

Returns

backtracking flag (1) backtracking performed (0) backtracking not performed

Definition at line 2868 of file NumericalSolution.f90.

    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    integer(I4B) :: idv_scale !< backtracking flag (1) backtracking performed (0) backtracking not performed
    ! local
    class(NumericalModelType), pointer :: mp => null()
    integer(I4B) :: i
 
    idv_scale = 0
    do i = 1, this%modellist%Count()
      mp => getnumericalmodelfromlist(this%modellist, i)
      if (mp%get_idv_scale() /= 0) then
        idv_scale = 1
      else
        if (idv_scale == 1) then
          idv_scale = -1
        end if
      end if
    end do
 

Here is the call graph for this function:

sln_get_loc()

subroutine numericalsolutionmodule::sln_get_loc ( class(numericalsolutiontype), intent(inout)  this, integer(i4b), intent(in)  nodesln, character(len=*), intent(inout)  str  )

private

Get the cell location string for the provided solution node number.

Parameters

[in,out]	this	NumericalSolutionType instance
[in]	nodesln	solution node number
[in,out]	str	string with user node number

Definition at line 3285 of file NumericalSolution.f90.

    ! -- dummy variables
    class(NumericalSolutionType), intent(inout) :: this !< NumericalSolutionType instance
    integer(I4B), intent(in) :: nodesln !< solution node number
    character(len=*), intent(inout) :: str !< string with user node number
    ! -- local variables
    class(NumericalModelType), pointer :: mp => null()
    integer(I4B) :: i
    integer(I4B) :: istart
    integer(I4B) :: iend
    integer(I4B) :: noder
    integer(I4B) :: nglo
    !
    ! -- initialize dummy variables
    str = ''
    !
    ! -- initialize local variables
    noder = 0
    !
    ! -- when parallel, account for offset
    nglo = nodesln + this%matrix_offset
    !
    ! -- calculate and set offsets
    do i = 1, this%modellist%Count()
      mp => getnumericalmodelfromlist(this%modellist, i)
      istart = 0
      iend = 0
      call mp%get_mrange(istart, iend)
      if (nglo >= istart .and. nglo <= iend) then
        noder = nglo - istart + 1
        call mp%get_mcellid(noder, str)
        exit
      end if
    end do

Here is the call graph for this function:

sln_get_nodeu()

subroutine numericalsolutionmodule::sln_get_nodeu ( class(numericalsolutiontype), intent(inout)  this, integer(i4b), intent(in)  nodesln, integer(i4b), intent(inout)  im, integer(i4b), intent(inout)  nodeu  )

private

Get the user node number from a model for the provided solution node number.

Parameters

[in,out]	this	NumericalSolutionType instance
[in]	nodesln	solution node number
[in,out]	im	solution model index (index in model list for this solution)
[in,out]	nodeu	user node number

Definition at line 3326 of file NumericalSolution.f90.

    ! -- dummy variables
    class(NumericalSolutionType), intent(inout) :: this !< NumericalSolutionType instance
    integer(I4B), intent(in) :: nodesln !< solution node number
    integer(I4B), intent(inout) :: im !< solution model index (index in model list for this solution)
    integer(I4B), intent(inout) :: nodeu !< user node number
    ! -- local variables
    class(NumericalModelType), pointer :: mp => null()
    integer(I4B) :: i
    integer(I4B) :: istart
    integer(I4B) :: iend
    integer(I4B) :: noder, nglo
    !
    ! -- initialize local variables
    noder = 0
    !
    ! -- when parallel, account for offset
    nglo = nodesln + this%matrix_offset
    !
    ! -- calculate and set offsets
    do i = 1, this%modellist%Count()
      mp => getnumericalmodelfromlist(this%modellist, i)
      istart = 0
      iend = 0
      call mp%get_mrange(istart, iend)
      if (nglo >= istart .and. nglo <= iend) then
        noder = nglo - istart + 1
        call mp%get_mnodeu(noder, nodeu)
        im = i
        exit
      end if
    end do

Here is the call graph for this function:

sln_has_converged()

logical(lgp) function numericalsolutionmodule::sln_has_converged ( class(numericalsolutiontype)  this, real(dp)  max_dvc  )

private

Parameters

this	NumericalSolutionType instance
max_dvc	the maximum dependent variable change

Returns

True, when converged

Definition at line 3216 of file NumericalSolution.f90.

    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    real(DP) :: max_dvc !< the maximum dependent variable change
    logical(LGP) :: has_converged !< True, when converged
 
    has_converged = .false.
    if (abs(max_dvc) <= this%dvclose) then
      has_converged = .true.
    end if
 

sln_l2norm()

subroutine numericalsolutionmodule::sln_l2norm ( class(numericalsolutiontype)  this, real(dp)  l2norm  )

private

A = the linear system matrix x = the dependent variable vector b = the right-hand side vector

 r = A * x - b

r_i = 0 if cell i is inactive L2norm = || r ||_2

Parameters

this	NumericalSolutionType instance
l2norm	calculated L-2 norm

Definition at line 2916 of file NumericalSolution.f90.

    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    real(DP) :: l2norm !< calculated L-2 norm
    ! local
    class(VectorBaseType), pointer :: vec_resid
 
    ! calc. residual vector
    vec_resid => this%system_matrix%create_vec(this%neq)
    call this%sln_calc_residual(vec_resid)
 
    ! 2-norm
    l2norm = vec_resid%norm2()
 
    ! clean up temp. vector
    call vec_resid%destroy()
    deallocate (vec_resid)

sln_ls()

subroutine numericalsolutionmodule::sln_ls ( class(numericalsolutiontype), intent(inout)  this, integer(i4b), intent(in)  kiter, integer(i4b), intent(in)  kstp, integer(i4b), intent(in)  kper, integer(i4b), intent(inout)  in_iter, integer(i4b), intent(inout)  iptc, real(dp), intent(in)  ptcf  )

private

Solve the linear system of equations. Steps include (1) matrix cleanup, (2) add pseudo-transient continuation terms, and (3) residual reduction.

Parameters

[in,out]

this

NumericalSolutionType instance

Definition at line 2442 of file NumericalSolution.f90.

    ! -- dummy variables
    class(NumericalSolutionType), intent(inout) :: this !< NumericalSolutionType instance
    integer(I4B), intent(in) :: kiter
    integer(I4B), intent(in) :: kstp
    integer(I4B), intent(in) :: kper
    integer(I4B), intent(inout) :: in_iter
    integer(I4B), intent(inout) :: iptc
    real(DP), intent(in) :: ptcf
    ! -- local variables
    logical(LGP) :: lsame
    integer(I4B) :: ieq
    integer(I4B) :: irow_glo
    integer(I4B) :: itestmat
    integer(I4B) :: ipos
    integer(I4B) :: icol_s
    integer(I4B) :: icol_e
    integer(I4B) :: jcol
    integer(I4B) :: iptct
    integer(I4B) :: iallowptc
    real(DP) :: adiag
    real(DP) :: diagval
    real(DP) :: l2norm
    real(DP) :: ptcval
    real(DP) :: bnorm
    character(len=50) :: fname
    character(len=*), parameter :: fmtfname = "('mf6mat_', i0, '_', i0, &
      &'_', i0, '_', i0, '.txt')"
    !
    ! -- take care of loose ends for all nodes before call to solver
    do ieq = 1, this%neq
      !
      ! -- get (global) cell id
      irow_glo = ieq + this%matrix_offset
      !
      ! -- store x in temporary location
      this%xtemp(ieq) = this%x(ieq)
      !
      ! -- make adjustments to the continuity equation for the node
      ! -- adjust small diagonal coefficient in an active cell
      if (this%active(ieq) > 0) then
        diagval = -done
        adiag = abs(this%system_matrix%get_diag_value(irow_glo))
        if (adiag < dem15) then
          call this%system_matrix%set_diag_value(irow_glo, diagval)
          this%rhs(ieq) = this%rhs(ieq) + diagval * this%x(ieq)
        end if
        ! -- Dirichlet boundary or no-flow cell
      else
        call this%system_matrix%set_diag_value(irow_glo, done)
        call this%system_matrix%zero_row_offdiag(irow_glo)
        this%rhs(ieq) = this%x(ieq)
      end if
    end do
    !
    ! -- complete adjustments for Dirichlet boundaries for a symmetric matrix
    if (this%isymmetric == 1 .and. simulation_mode == "SEQUENTIAL") then
      do ieq = 1, this%neq
        if (this%active(ieq) > 0) then
          icol_s = this%system_matrix%get_first_col_pos(ieq)
          icol_e = this%system_matrix%get_last_col_pos(ieq)
          do ipos = icol_s, icol_e
            jcol = this%system_matrix%get_column(ipos)
            if (jcol == ieq) cycle
            if (this%active(jcol) < 0) then
              this%rhs(ieq) = this%rhs(ieq) - &
                              (this%system_matrix%get_value_pos(ipos) * &
                               this%x(jcol))
              call this%system_matrix%set_value_pos(ipos, dzero)
            end if
 
          end do
        end if
      end do
    end if
    !
    ! -- pseudo transient continuation
    !
    ! -- set iallowptc
    ! -- no_ptc_option is FIRST
    if (this%iallowptc < 0) then
      if (kper > 1) then
        iallowptc = 1
      else
        iallowptc = 0
      end if
      !
      ! -- no_ptc_option is ALL (0) or using PTC (1)
    else
      iallowptc = this%iallowptc
    end if
    !
    ! -- set iptct
    iptct = iptc * iallowptc
    !
    ! -- calculate or modify pseudo transient continuation terms and add
    !    to amat diagonals
    if (iptct /= 0) then
      call this%sln_l2norm(l2norm)
      ! -- confirm that the l2norm exceeds previous l2norm
      !    if not, there is no need to add ptc terms
      if (kiter == 1) then
        if (kper > 1 .or. kstp > 1) then
          if (l2norm <= this%l2norm0) then
            iptc = 0
          end if
        end if
      else
        lsame = is_close(l2norm, this%l2norm0)
        if (lsame) then
          iptc = 0
        end if
      end if
    end if
    iptct = iptc * iallowptc
    if (iptct /= 0) then
      if (kiter == 1) then
        if (this%iptcout > 0) then
          write (this%iptcout, '(A10,6(1x,A15))') 'OUTER ITER', &
            '         PTCDEL', '        L2NORM0', '         L2NORM', &
            '        RHSNORM', '       1/PTCDEL', ' RHSNORM/L2NORM'
        end if
        if (this%ptcdel0 > dzero) then
          this%ptcdel = this%ptcdel0
        else
          if (this%iptcopt == 0) then
            !
            ! -- ptcf is the reciprocal of the pseudo-time step
            this%ptcdel = done / ptcf
          else
            bnorm = dzero
            do ieq = 1, this%neq
              if (this%active(ieq) .gt. 0) then
                bnorm = bnorm + this%rhs(ieq) * this%rhs(ieq)
              end if
            end do
            bnorm = sqrt(bnorm)
            this%ptcdel = bnorm / l2norm
          end if
        end if
      else
        if (l2norm > dzero) then
          this%ptcdel = this%ptcdel * (this%l2norm0 / l2norm)**this%ptcexp
        else
          this%ptcdel = dzero
        end if
      end if
      if (this%ptcdel > dzero) then
        ptcval = done / this%ptcdel
      else
        ptcval = done
      end if
      bnorm = dzero
      do ieq = 1, this%neq
        irow_glo = ieq + this%matrix_offset
        if (this%active(ieq) > 0) then
          diagval = abs(this%system_matrix%get_diag_value(irow_glo))
          bnorm = bnorm + this%rhs(ieq) * this%rhs(ieq)
          call this%system_matrix%add_diag_value(irow_glo, -ptcval)
          this%rhs(ieq) = this%rhs(ieq) - ptcval * this%x(ieq)
        end if
      end do
      bnorm = sqrt(bnorm)
      if (this%iptcout > 0) then
        write (this%iptcout, '(i10,5(1x,e15.7),1(1x,f15.6))') &
          kiter, this%ptcdel, this%l2norm0, l2norm, bnorm, &
          ptcval, bnorm / l2norm
      end if
      this%l2norm0 = l2norm
    end if
    !
    ! -- save rhs, amat to a file
    !    to enable set itestmat to 1 and recompile
    !-------------------------------------------------------
    itestmat = 0
    if (itestmat == 1) then
      write (fname, fmtfname) this%id, kper, kstp, kiter
      print *, 'Saving amat to: ', trim(adjustl(fname))
 
      itestmat = getunit()
      open (itestmat, file=trim(adjustl(fname)))
      write (itestmat, *) 'NODE, RHS, AMAT FOLLOW'
      do ieq = 1, this%neq
        irow_glo = ieq + this%matrix_offset
        icol_s = this%system_matrix%get_first_col_pos(irow_glo)
        icol_e = this%system_matrix%get_last_col_pos(irow_glo)
        write (itestmat, '(*(G0,:,","))') &
          irow_glo, &
          this%rhs(ieq), &
          (this%system_matrix%get_column(ipos), ipos=icol_s, icol_e), &
          (this%system_matrix%get_value_pos(ipos), ipos=icol_s, icol_e)
      end do
      close (itestmat)
      !stop
    end if
    !-------------------------------------------------------
    !
    ! -- call appropriate linear solver
    !
    ! -- ims linear solver - linmeth option 1
    if (this%linsolver == ims_solver) then
      call this%imslinear%imslinear_apply(this%icnvg, kstp, kiter, in_iter, &
                                          this%nitermax, this%convnmod, &
                                          this%convmodstart, this%caccel, &
                                          this%cnvg_summary)
    else if (this%linsolver == petsc_solver) then
      call this%linear_solver%solve(kiter, this%vec_rhs, &
                                    this%vec_x, this%cnvg_summary)
      in_iter = this%linear_solver%iteration_number
      this%icnvg = this%linear_solver%is_converged
    end if

Here is the call graph for this function:

sln_maxval()

subroutine numericalsolutionmodule::sln_maxval ( class(numericalsolutiontype)  this, integer(i4b), intent(in)  nsize, real(dp), dimension(nsize), intent(in)  v, real(dp), intent(inout)  vmax  )

private

Return the maximum value in a vector using a normalized form.

Parameters

	this	NumericalSolutionType instance
[in]	nsize	length of vector
[in]	v	input vector
[in,out]	vmax	maximum value

Definition at line 2939 of file NumericalSolution.f90.

    ! -- dummy variables
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    integer(I4B), intent(in) :: nsize !< length of vector
    real(DP), dimension(nsize), intent(in) :: v !< input vector
    real(DP), intent(inout) :: vmax !< maximum value
    ! -- local variables
    integer(I4B) :: n
    real(DP) :: d
    real(DP) :: denom
    real(DP) :: dnorm
    !
    ! -- determine maximum value
    vmax = v(1)
    do n = 2, nsize
      d = v(n)
      denom = abs(vmax)
      if (denom == dzero) then
        denom = dprec
      end if
      !
      ! -- calculate normalized value
      dnorm = abs(d) / denom
      if (dnorm > done) then
        vmax = d
      end if
    end do

sln_nur_has_converged()

logical(lgp) function numericalsolutionmodule::sln_nur_has_converged ( class(numericalsolutiontype)  this, real(dp), intent(in)  dxold_max, real(dp), intent(in)  hncg  )

private

Parameters

	this	NumericalSolutionType instance
[in]	dxold_max	the maximum dependent variable change for unrelaxed cells
[in]	hncg	largest dep. var. change at end of Picard iteration

Returns

True, when converged

Definition at line 3266 of file NumericalSolution.f90.

    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    real(DP), intent(in) :: dxold_max !< the maximum dependent variable change for unrelaxed cells
    real(DP), intent(in) :: hncg !< largest dep. var. change at end of Picard iteration
    logical(LGP) :: has_converged !< True, when converged
 
    has_converged = .false.
    if (abs(dxold_max) <= this%dvclose .and. &
        abs(hncg) <= this%dvclose) then
      has_converged = .true.
    end if
 

sln_ot()

subroutine numericalsolutionmodule::sln_ot

(

class(numericalsolutiontype)

this

)

Output solution data. Currently does nothing.

Parameters

this	NumericalSolutionType instance

Definition at line 1165 of file NumericalSolution.f90.

    ! -- dummy variables
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    !
    ! -- Nothing to do here

sln_package_convergence()

integer(i4b) function numericalsolutionmodule::sln_package_convergence ( class(numericalsolutiontype)  this, real(dp), intent(in)  dpak, character(len=lenpakloc), intent(in)  cpakout, integer(i4b), intent(in)  iend  )

private

Parameters

	this	NumericalSolutionType instance
[in]	dpak	Newton Under-relaxation flag
[in]	cpakout	string with package that caused failure
[in]	iend	flag indicating if last inner iteration (iend=1)

Definition at line 3230 of file NumericalSolution.f90.

    ! dummy
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    real(DP), intent(in) :: dpak !< Newton Under-relaxation flag
    character(len=LENPAKLOC), intent(in) :: cpakout !< string with package that caused failure
    integer(I4B), intent(in) :: iend !< flag indicating if last inner iteration (iend=1)
    ! local
    integer(I4B) :: ivalue
    ivalue = 1
    if (abs(dpak) > this%dvclose) then
      ivalue = 0
      ! -- write message to stdout
      if (iend /= 0) then
        write (errmsg, '(3a)') &
          'PACKAGE (', trim(cpakout), ') CAUSED CONVERGENCE FAILURE'
        call write_message(errmsg)
      end if
    end if
 

Here is the call graph for this function:

sln_reset()

subroutine numericalsolutionmodule::sln_reset

(

class(numericalsolutiontype)

this

)

Reset the solution by setting the coefficient matrix and right-hand side vectors to zero.

Parameters

this	NumericalSolutionType instance

Definition at line 2427 of file NumericalSolution.f90.

    ! -- dummy variables
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    !
    ! -- reset the solution
    call this%system_matrix%zero_entries()
    call this%vec_rhs%zero_entries()

sln_setouter()

subroutine numericalsolutionmodule::sln_setouter ( class(numericalsolutiontype), intent(inout)  this, integer(i4b), intent(in)  ifdparam  )

private

Set default Picard iteration variables based on passed complexity option.

Parameters

[in,out]	this	NumericalSolutionType instance
[in]	ifdparam	complexity option (1) simple (2) moderate (3) complex

Definition at line 2661 of file NumericalSolution.f90.

    ! -- dummy variables
    class(NumericalSolutionType), intent(inout) :: this !< NumericalSolutionType instance
    integer(I4B), intent(in) :: ifdparam !< complexity option (1) simple (2) moderate (3) complex
    !
    ! -- simple option
    select case (ifdparam)
    case (1)
      this%dvclose = dem3
      this%mxiter = 25
      this%nonmeth = 0
      this%theta = done
      this%akappa = dzero
      this%gamma = done
      this%amomentum = dzero
      this%numtrack = 0
      this%btol = dzero
      this%breduc = dzero
      this%res_lim = dzero
      !
      ! -- moderate
    case (2)
      this%dvclose = dem2
      this%mxiter = 50
      this%nonmeth = 3
      this%theta = 0.9d0
      this%akappa = 0.0001d0
      this%gamma = dzero
      this%amomentum = dzero
      this%numtrack = 0
      this%btol = dzero
      this%breduc = dzero
      this%res_lim = dzero
      !
      ! -- complex
    case (3)
      this%dvclose = dem1
      this%mxiter = 100
      this%nonmeth = 3
      this%theta = 0.8d0
      this%akappa = 0.0001d0
      this%gamma = dzero
      this%amomentum = dzero
      this%numtrack = 20
      this%btol = 1.05d0
      this%breduc = 0.1d0
      this%res_lim = 0.002d0
    end select

sln_sync_newtonur_flag()

integer(i4b) function numericalsolutionmodule::sln_sync_newtonur_flag ( class(numericalsolutiontype)  this, integer(i4b), intent(in)  inewtonur  )

private

Parameters

	this	NumericalSolutionType instance
[in]	inewtonur	Newton Under-relaxation flag

Returns

Default is set to current value (1 = under-relaxation applied)

Definition at line 3253 of file NumericalSolution.f90.

    ! dummy
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    integer(I4B), intent(in) :: inewtonur !< Newton Under-relaxation flag
    ! local
    integer(I4B) :: ivalue !< Default is set to current value (1 = under-relaxation applied)
 
    ivalue = inewtonur
 

sln_underrelax()

subroutine numericalsolutionmodule::sln_underrelax ( class(numericalsolutiontype)  this, integer(i4b), intent(in)  kiter, real(dp), intent(in)  bigch, integer(i4b), intent(in)  neq, integer(i4b), dimension(neq), intent(in)  active, real(dp), dimension(neq), intent(inout)  x, real(dp), dimension(neq), intent(in)  xtemp  )

private

Under relax using the simple, cooley, or delta-bar-delta methods.

Parameters

	this	NumericalSolutionType instance
[in]	kiter	Picard iteration number
[in]	bigch	maximum dependent-variable change
[in]	neq	number of equations
[in]	active	active cell flag (1)
[in,out]	x	current dependent-variable
[in]	xtemp	previous dependent-variable

Definition at line 3050 of file NumericalSolution.f90.

    ! -- dummy variables
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    integer(I4B), intent(in) :: kiter !< Picard iteration number
    real(DP), intent(in) :: bigch !< maximum dependent-variable change
    integer(I4B), intent(in) :: neq !< number of equations
    integer(I4B), dimension(neq), intent(in) :: active !< active cell flag (1)
    real(DP), dimension(neq), intent(inout) :: x !< current dependent-variable
    real(DP), dimension(neq), intent(in) :: xtemp !< previous dependent-variable
    ! -- local variables
    integer(I4B) :: n
    real(DP) :: ww
    real(DP) :: delx
    real(DP) :: relax
    real(DP) :: es
    real(DP) :: aes
    real(DP) :: amom
    !
    ! -- option for using simple dampening (as done by MODFLOW-2005 PCG)
    if (this%nonmeth == 1) then
      do n = 1, neq
        !
        ! -- skip inactive nodes
        if (active(n) < 1) cycle
        !
        ! -- compute step-size (delta x)
        delx = x(n) - xtemp(n)
        this%dxold(n) = delx
 
        ! -- dampen dependent variable solution
        x(n) = xtemp(n) + this%gamma * delx
      end do
      !
      ! -- option for using cooley underrelaxation
    else if (this%nonmeth == 2) then
      !
      ! -- set bigch
      this%bigch = bigch
      !
      ! -- initialize values for first iteration
      if (kiter == 1) then
        relax = done
        this%relaxold = done
        this%bigchold = bigch
      else
        !
        ! -- compute relaxation factor
        es = this%bigch / (this%bigchold * this%relaxold)
        aes = abs(es)
        if (es < -done) then
          relax = dhalf / aes
        else
          relax = (dthree + es) / (dthree + aes)
        end if
      end if
      this%relaxold = relax
      !
      ! -- modify cooley to use weighted average of past changes
      this%bigchold = (done - this%gamma) * this%bigch + this%gamma * &
                      this%bigchold
      !
      ! -- compute new dependent variable after under-relaxation
      if (relax < done) then
        do n = 1, neq
          !
          ! -- skip inactive nodes
          if (active(n) < 1) cycle
          !
          ! -- update dependent variable
          delx = x(n) - xtemp(n)
          this%dxold(n) = delx
          x(n) = xtemp(n) + relax * delx
        end do
      end if
      !
      ! -- option for using delta-bar-delta scheme to under-relax for all equations
    else if (this%nonmeth == 3) then
      do n = 1, neq
        !
        ! -- skip inactive nodes
        if (active(n) < 1) cycle
        !
        ! -- compute step-size (delta x) and initialize d-b-d parameters
        delx = x(n) - xtemp(n)
        !
        ! -- initialize values for first iteration
        if (kiter == 1) then
          this%wsave(n) = done
          this%hchold(n) = dem20
          this%deold(n) = dzero
        end if
        !
        ! -- compute new relaxation term as per delta-bar-delta
        ww = this%wsave(n)
        !
        ! for flip-flop condition, decrease factor
        if (this%deold(n) * delx < dzero) then
          ww = this%theta * this%wsave(n)
          ! -- when change is of same sign, increase factor
        else
          ww = this%wsave(n) + this%akappa
        end if
        if (ww > done) ww = done
        this%wsave(n) = ww
        !
        ! -- compute weighted average of past changes in hchold
        if (kiter == 1) then
          this%hchold(n) = delx
        else
          this%hchold(n) = (done - this%gamma) * delx + &
                           this%gamma * this%hchold(n)
        end if
        !
        ! -- store slope (change) term for next iteration
        this%deold(n) = delx
        this%dxold(n) = delx
        !
        ! -- compute accepted step-size and new dependent variable
        amom = dzero
        if (kiter > 4) amom = this%amomentum
        delx = delx * ww + amom * this%hchold(n)
        x(n) = xtemp(n) + delx
      end do
      !
    end if

solve()

subroutine numericalsolutionmodule::solve ( class(numericalsolutiontype)  this, integer(i4b), intent(in)  kiter  )

private

Builds and solve the system for this numerical solution. It roughly consists of the following steps (1) backtracking, (2) reset amat and rhs (3) calculate matrix terms (*_cf), (4) add coefficients to matrix (*_fc), (6) newton-raphson, (6) PTC, (7) linear solve, (8) convergence checks, (9) write output, and (10) underrelaxation

Parameters

	this	NumericalSolutionType instance
[in]	kiter	Picard iteration number

Definition at line 1512 of file NumericalSolution.f90.

    ! -- modules
    use tdismodule, only: kstp, kper, totim
    ! -- dummy variables
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    integer(I4B), intent(in) :: kiter !< Picard iteration number
    ! -- local variables
    class(NumericalModelType), pointer :: mp => null()
    class(NumericalExchangeType), pointer :: cp => null()
    character(len=LINELENGTH) :: title
    character(len=LINELENGTH) :: tag
    character(len=LENPAKLOC) :: cmod
    character(len=LENPAKLOC) :: cpak
    character(len=LENPAKLOC) :: cpakout
    character(len=LENPAKLOC) :: strh
    character(len=25) :: cval
    character(len=7) :: cmsg
    integer(I4B) :: ic
    integer(I4B) :: im, m_idx, model_id
    integer(I4B) :: icsv0
    integer(I4B) :: kcsv0
    integer(I4B) :: ntabrows
    integer(I4B) :: ntabcols
    integer(I4B) :: i0, i1
    integer(I4B) :: itestmat, n
    integer(I4B) :: iter
    integer(I4B) :: inewtonur
    integer(I4B) :: locmax_nur
    integer(I4B) :: iend
    integer(I4B) :: icnvgmod
    integer(I4B) :: iptc
    integer(I4B) :: node_user
    integer(I4B) :: ipak
    integer(I4B) :: ipos0
    integer(I4B) :: ipos1
    real(DP) :: dxmax_nur
    real(DP) :: dxold_max
    real(DP) :: ptcf
    real(DP) :: ttform
    real(DP) :: ttsoln
    real(DP) :: dpak
    real(DP) :: outer_hncg
 
    ! start timer
    call g_prof%start("Solve"//this%id_postfix, this%tmr_solve)
 
    !
    ! -- initialize local variables
    icsv0 = max(1, this%itertot_sim + 1)
    kcsv0 = max(1, this%itertot_timestep + 1)
    !
    ! -- create header for outer iteration table
    if (this%iprims > 0) then
      if (.not. associated(this%outertab)) then
        !
        ! -- create outer iteration table
        ! -- table dimensions
        ntabrows = 1
        ntabcols = 6
        if (this%numtrack > 0) then
          ntabcols = ntabcols + 4
        end if
        !
        ! -- initialize table and define columns
        title = trim(this%memory_path)//' OUTER ITERATION SUMMARY'
        call table_cr(this%outertab, this%name, title)
        call this%outertab%table_df(ntabrows, ntabcols, iout, &
                                    finalize=.false.)
        tag = 'OUTER ITERATION STEP'
        call this%outertab%initialize_column(tag, 25, alignment=tableft)
        tag = 'OUTER ITERATION'
        call this%outertab%initialize_column(tag, 10, alignment=tabright)
        tag = 'INNER ITERATION'
        call this%outertab%initialize_column(tag, 10, alignment=tabright)
        if (this%numtrack > 0) then
          tag = 'BACKTRACK FLAG'
          call this%outertab%initialize_column(tag, 10, alignment=tabright)
          tag = 'BACKTRACK ITERATIONS'
          call this%outertab%initialize_column(tag, 10, alignment=tabright)
          tag = 'INCOMING RESIDUAL'
          call this%outertab%initialize_column(tag, 15, alignment=tabright)
          tag = 'OUTGOING RESIDUAL'
          call this%outertab%initialize_column(tag, 15, alignment=tabright)
        end if
        tag = 'MAXIMUM CHANGE'
        call this%outertab%initialize_column(tag, 15, alignment=tabright)
        tag = 'STEP SUCCESS'
        call this%outertab%initialize_column(tag, 7, alignment=tabright)
        tag = 'MAXIMUM CHANGE MODEL-(CELLID) OR MODEL-PACKAGE-(NUMBER)'
        call this%outertab%initialize_column(tag, 34, alignment=tabright)
      end if
    end if
    !
    ! -- backtracking
    if (this%numtrack > 0) then
      call this%sln_backtracking(mp, cp, kiter)
    end if
    !
    call code_timer(0, ttform, this%ttform)
    call g_prof%start("Formulate", this%tmr_formulate)
    !
    ! -- (re)build the solution matrix
    call this%sln_buildsystem(kiter, inewton=1)
    !
    ! -- Calculate pseudo-transient continuation factor for each model
    call this%sln_calc_ptc(iptc, ptcf)
    !
    ! -- Add model Newton-Raphson terms to solution
    do im = 1, this%modellist%Count()
      mp => getnumericalmodelfromlist(this%modellist, im)
      call mp%model_nr(kiter, this%system_matrix, 1)
    end do
    call code_timer(1, ttform, this%ttform)
    call g_prof%stop(this%tmr_formulate)
 
    ! x and rhs scaling
    if (this%idv_scale /= 0) then
      call this%sln_maxval(this%neq, this%x, this%dscale)
      call ims_misc_dvscale(0, this%neq, this%dscale, this%x, this%rhs)
    end if
    !
    ! -- linear solve
    call code_timer(0, ttsoln, this%ttsoln)
    call g_prof%start("Linear solve", this%tmr_linsolve)
    call this%sln_ls(kiter, kstp, kper, iter, iptc, ptcf)
    call g_prof%stop(this%tmr_linsolve)
    call code_timer(1, ttsoln, this%ttsoln)
    !
    ! -- increment counters storing the total number of linear and
    !    non-linear iterations for this timestep and the total
    !    number of linear iterations for all timesteps
    this%itertot_timestep = this%itertot_timestep + iter
    this%iouttot_timestep = this%iouttot_timestep + 1
    this%itertot_sim = this%itertot_sim + iter
    !
    ! -- save matrix to a file
    !    to enable set itestmat to 1 and recompile
    !-------------------------------------------------------
    itestmat = 0
    if (itestmat /= 0) then
      open (99, file='sol_MF6.TXT')
      WRITE (99, *) 'MATRIX SOLUTION FOLLOWS'
      WRITE (99, '(10(I8,G15.4))') (n, this%x(n), n=1, this%NEQ)
      close (99)
      call pstop()
    end if
    !-------------------------------------------------------
    !
    ! -- check convergence of solution
    call this%sln_get_dxmax(this%hncg(kiter), this%lrch(1, kiter))
    if (this%icnvg /= 0) then
      this%icnvg = 0
      if (this%sln_has_converged(this%hncg(kiter))) then
        this%icnvg = 1
      end if
    end if
    !
    ! -- set failure flag
    if (this%icnvg == 0) then
      cmsg = ' '
    else
      cmsg = '*'
    end if
    !
    ! -- set flag if this is the last outer iteration
    iend = 0
    if (kiter == this%mxiter) then
      iend = 1
    end if
    !
    ! -- write maximum dependent-variable change from linear solver to list file
    if (this%iprims > 0) then
      cval = 'Model'
      call this%sln_get_loc(this%lrch(1, kiter), strh)
      !
      ! -- add data to outertab
      call this%outertab%add_term(cval)
      call this%outertab%add_term(kiter)
      call this%outertab%add_term(iter)
      if (this%numtrack > 0) then
        call this%outertab%add_term(' ')
        call this%outertab%add_term(' ')
        call this%outertab%add_term(' ')
        call this%outertab%add_term(' ')
      end if
      call this%outertab%add_term(this%hncg(kiter))
      call this%outertab%add_term(cmsg)
      call this%outertab%add_term(trim(strh))
    end if
    !
    ! -- Additional convergence check for exchanges
    do ic = 1, this%exchangelist%Count()
      cp => getnumericalexchangefromlist(this%exchangelist, ic)
      call cp%exg_cc(this%icnvg)
    end do
    !
    ! -- additional convergence check for model packages
    icnvgmod = this%icnvg
    cpak = ' '
    ipak = 0
    dpak = dzero
    do im = 1, this%modellist%Count()
      mp => getnumericalmodelfromlist(this%modellist, im)
      call mp%get_mcellid(0, cmod)
      call mp%model_cc(this%itertot_sim, kiter, iend, icnvgmod, &
                       cpak, ipak, dpak)
      if (ipak /= 0) then
        ipos0 = index(cpak, '-', back=.true.)
        ipos1 = len_trim(cpak)
        write (cpakout, '(a,a,"-(",i0,")",a)') &
          trim(cmod), cpak(1:ipos0 - 1), ipak, cpak(ipos0:ipos1)
      else
        cpakout = ' '
      end if
    end do
    !
    ! -- evaluate package convergence - only done if convergence is achieved
    if (this%icnvg == 1) then
      this%icnvg = this%sln_package_convergence(dpak, cpakout, iend)
      !
      ! -- write maximum change in package convergence check
      if (this%iprims > 0) then
        cval = 'Package'
        if (this%icnvg /= 1) then
          cmsg = ' '
        else
          cmsg = '*'
        end if
        if (len_trim(cpakout) > 0) then
          !
          ! -- add data to outertab
          call this%outertab%add_term(cval)
          call this%outertab%add_term(kiter)
          call this%outertab%add_term(' ')
          if (this%numtrack > 0) then
            call this%outertab%add_term(' ')
            call this%outertab%add_term(' ')
            call this%outertab%add_term(' ')
            call this%outertab%add_term(' ')
          end if
          call this%outertab%add_term(dpak)
          call this%outertab%add_term(cmsg)
          call this%outertab%add_term(cpakout)
        end if
      end if
    end if
    !
    ! -- under-relaxation - only done if convergence not achieved
    if (this%icnvg /= 1) then
      if (this%nonmeth > 0) then
        call this%sln_underrelax(kiter, this%hncg(kiter), this%neq, &
                                 this%active, this%x, this%xtemp)
      else
        call this%sln_calcdx(this%neq, this%active, &
                             this%x, this%xtemp, this%dxold)
      end if
      !
      ! -- adjust heads by newton under-relaxation, if necessary
      inewtonur = 0
      dxmax_nur = dzero
      locmax_nur = 0
      do im = 1, this%modellist%Count()
        mp => getnumericalmodelfromlist(this%modellist, im)
        i0 = mp%moffset + 1 - this%matrix_offset
        i1 = i0 + mp%neq - 1
        call mp%model_nur(mp%neq, this%x(i0:i1), this%xtemp(i0:i1), &
                          this%dxold(i0:i1), inewtonur, dxmax_nur, locmax_nur)
      end do
      !
      ! -- synchronize Newton Under-relaxation flag
      inewtonur = this%sln_sync_newtonur_flag(inewtonur)
      !
      ! -- check for convergence if newton under-relaxation applied
      if (inewtonur /= 0) then
        !
        ! -- calculate maximum change in heads in cells that have
        !    not been adjusted by newton under-relxation
        call this%sln_maxval(this%neq, this%dxold, dxold_max)
        !
        ! -- evaluate convergence
        if (this%sln_nur_has_converged(dxold_max, this%hncg(kiter))) then
          !
          ! -- converged
          this%icnvg = 1
          !
          ! -- reset outer dependent-variable change and location for output
          call this%sln_get_dxmax(this%hncg(kiter), this%lrch(1, kiter))
          !
          ! -- write revised dependent-variable change data after
          !    newton under-relaxation
          if (this%iprims > 0) then
            cval = 'Newton under-relaxation'
            cmsg = '*'
            call this%sln_get_loc(this%lrch(1, kiter), strh)
            !
            ! -- add data to outertab
            call this%outertab%add_term(cval)
            call this%outertab%add_term(kiter)
            call this%outertab%add_term(iter)
            if (this%numtrack > 0) then
              call this%outertab%add_term(' ')
              call this%outertab%add_term(' ')
              call this%outertab%add_term(' ')
              call this%outertab%add_term(' ')
            end if
            call this%outertab%add_term(this%hncg(kiter))
            call this%outertab%add_term(cmsg)
            call this%outertab%add_term(trim(strh))
          end if
        end if
      end if
    end if
    !
    ! -- write to outer iteration csv file
    if (this%icsvouterout > 0) then
      !
      ! -- set outer dependent-variable change variable
      outer_hncg = this%hncg(kiter)
      !
      ! -- model convergence error
      if (abs(outer_hncg) > abs(dpak)) then
        !
        ! -- get model number and user node number
        call this%sln_get_nodeu(this%lrch(1, kiter), m_idx, node_user)
        mp => getnumericalmodelfromlist(this%modellist, m_idx)
        model_id = mp%id
        cpakout = ''
      else if (outer_hncg == dzero .and. dpak == dzero) then ! zero change, location could be any
        model_id = 0
        node_user = 0
        !
        ! -- then it's a package convergence error
      else
        !
        ! -- set convergence error, model number, user node number,
        !    and package name
        outer_hncg = dpak
        ipos0 = index(cmod, '_')
        read (cmod(1:ipos0 - 1), *) model_id
        node_user = ipak
        ipos0 = index(cpak, '-', back=.true.)
        cpakout = cpak(1:ipos0 - 1)
      end if
 
      write (this%icsvouterout, '(*(G0,:,","))') &
        this%itertot_sim, totim, kper, kstp, kiter, iter, &
        outer_hncg, model_id, trim(cpakout), node_user
    end if
    !
    ! -- write to inner iteration csv file
    if (this%icsvinnerout > 0) then
      call this%csv_convergence_summary(this%icsvinnerout, totim, kper, kstp, &
                                        kiter, iter, icsv0, kcsv0)
    end if
 
    ! undo x and rhs scaling
    if (this%idv_scale /= 0) then
      call ims_misc_dvscale(1, this%neq, this%dscale, this%x, this%rhs)
    end if
 
    ! stop timer
    call g_prof%stop(this%tmr_solve)
 

Here is the call graph for this function:

writecsvheader()

subroutine numericalsolutionmodule::writecsvheader ( class(numericalsolutiontype)  this )

private

Write header for solver output to comma-separated value files.

Parameters

this	NumericalSolutionType instance

Definition at line 1364 of file NumericalSolution.f90.

    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    ! local variables
    integer(I4B) :: im
    class(NumericalModelType), pointer :: mp => null()
    !
    ! -- outer iteration csv header
    if (this%icsvouterout > 0) then
      write (this%icsvouterout, '(*(G0,:,","))') &
        'total_inner_iterations', 'totim', 'kper', 'kstp', 'nouter', &
        'inner_iterations', 'solution_outer_dvmax', &
        'solution_outer_dvmax_model', 'solution_outer_dvmax_package', &
        'solution_outer_dvmax_node'
    end if
    !
    ! -- inner iteration csv header
    if (this%icsvinnerout > 0) then
      write (this%icsvinnerout, '(*(G0,:,","))', advance='NO') &
        'total_inner_iterations', 'totim', 'kper', 'kstp', 'nouter', &
        'ninner', 'solution_inner_dvmax', 'solution_inner_dvmax_model', &
        'solution_inner_dvmax_node'
      write (this%icsvinnerout, '(*(G0,:,","))', advance='NO') &
        '', 'solution_inner_rmax', 'solution_inner_rmax_model', &
        'solution_inner_rmax_node'
      ! solver items specific to ims solver
      if (this%linsolver == ims_solver) then
        write (this%icsvinnerout, '(*(G0,:,","))', advance='NO') &
          '', 'solution_inner_alpha'
        if (this%imslinear%ilinmeth == 2) then
          write (this%icsvinnerout, '(*(G0,:,","))', advance='NO') &
            '', 'solution_inner_omega'
        end if
      end if
      ! -- check for more than one model - ims only
      if (this%convnmod > 1 .or. simulation_mode == "PARALLEL") then
        do im = 1, this%modellist%Count()
          mp => getnumericalmodelfromlist(this%modellist, im)
          write (this%icsvinnerout, '(*(G0,:,","))', advance='NO') &
            '', trim(adjustl(mp%name))//'_inner_dvmax', &
            trim(adjustl(mp%name))//'_inner_dvmax_node', &
            trim(adjustl(mp%name))//'_inner_rmax', &
            trim(adjustl(mp%name))//'_inner_rmax_node'
        end do
      end if
      write (this%icsvinnerout, '(a)') ''
    end if

Here is the call graph for this function:

writeptcinfotofile()

subroutine numericalsolutionmodule::writeptcinfotofile ( class(numericalsolutiontype)  this, integer(i4b), intent(in)  kper  )

private

Write header for pseudo-transient continuation information to a file.

Parameters

	this	NumericalSolutionType instance
[in]	kper	current stress period number

Definition at line 1417 of file NumericalSolution.f90.

    ! -- dummy variables
    class(NumericalSolutionType) :: this !< NumericalSolutionType instance
    integer(I4B), intent(in) :: kper !< current stress period number
    ! -- local variable
    integer(I4B) :: n, im, iallowptc, iptc
    class(NumericalModelType), pointer :: mp => null()
 
    ! -- determine if PTC will be used in any model
    n = 1
    do im = 1, this%modellist%Count()
      !
      ! -- set iallowptc
      ! -- no_ptc_option is FIRST
      if (this%iallowptc < 0) then
        if (kper > 1) then
          iallowptc = 1
        else
          iallowptc = 0
        end if
        ! -- no_ptc_option is ALL (0) or using PTC (1)
      else
        iallowptc = this%iallowptc
      end if
 
      if (iallowptc > 0) then
        mp => getnumericalmodelfromlist(this%modellist, im)
        call mp%model_ptcchk(iptc)
      else
        iptc = 0
      end if
 
      if (iptc /= 0) then
        if (n == 1) then
          write (iout, '(//)')
          n = 0
        end if
        write (iout, '(1x,a,1x,i0,1x,3a)') &
          'PSEUDO-TRANSIENT CONTINUATION WILL BE APPLIED TO MODEL', im, '("', &
          trim(adjustl(mp%name)), '") DURING THIS TIME STEP'
      end if
    end do
 

Here is the call graph for this function:

Variable Documentation

ims_solver

integer(i4b), parameter numericalsolutionmodule::ims_solver = 1

private

Definition at line 55 of file NumericalSolution.f90.

  integer(I4B), parameter :: IMS_SOLVER = 1

petsc_solver

integer(i4b), parameter numericalsolutionmodule::petsc_solver = 2

private

Definition at line 56 of file NumericalSolution.f90.

  integer(I4B), parameter :: PETSC_SOLVER = 2