17 logical,
pointer :: smgnc => null()
18 logical,
pointer ::
implicit => null()
19 logical,
pointer :: i2kn => null()
20 integer(I4B),
pointer :: nexg => null()
21 integer(I4B),
pointer :: numjs => null()
24 integer(I4B),
dimension(:),
pointer,
contiguous :: nodem1 => null()
25 integer(I4B),
dimension(:),
pointer,
contiguous :: nodem2 => null()
26 integer(I4B),
dimension(:, :),
pointer,
contiguous :: nodesj => null()
27 real(dp),
dimension(:),
pointer,
contiguous :: cond => null()
28 integer(I4B),
dimension(:),
pointer,
contiguous :: idxglo => null()
29 integer(I4B),
dimension(:),
pointer,
contiguous :: idxsymglo => null()
30 real(dp),
dimension(:, :),
pointer,
contiguous :: alphasj => null()
31 integer(I4B),
dimension(:),
pointer,
contiguous :: idiagn => null()
32 integer(I4B),
dimension(:),
pointer,
contiguous :: idiagm => null()
33 integer(I4B),
dimension(:, :),
pointer,
contiguous :: jposinrown => null()
34 integer(I4B),
dimension(:, :),
pointer,
contiguous :: jposinrowm => null()
60 subroutine gnc_cr(gncobj, name_parent, inunit, iout)
63 character(len=*),
intent(in) :: name_parent
64 integer(I4B),
intent(in) :: inunit
65 integer(I4B),
intent(in) :: iout
72 call gncobj%set_names(1, name_parent,
'GNC',
'GNC')
75 call gncobj%allocate_scalars()
78 gncobj%inunit = inunit
102 if (
present(m2))
then
108 call this%parser%Initialize(this%inunit, this%iout)
111 call this%read_options()
114 call this%read_dimensions()
117 call this%allocate_arrays()
120 call this%read_data()
123 if (this%m1%idsoln /= this%m2%idsoln)
then
124 if (this%implicit)
then
125 write (
errmsg,
'(a)')
'GNC is implicit but models are in '// &
126 'different solutions.'
147 integer(I4B) :: ignc, jidx, noden, nodem, nodej
152 if (this%implicit)
then
153 do ignc = 1, this%nexg
154 noden = this%nodem1(ignc) + this%m1%moffset
155 nodem = this%nodem2(ignc) + this%m2%moffset
156 jloop:
do jidx = 1, this%numjs
157 nodej = this%nodesj(jidx, ignc)
158 if (nodej == 0) cycle
159 nodej = nodej + this%m1%moffset
160 call sparse%addconnection(nodem, nodej, 1)
161 call sparse%addconnection(nodej, nodem, 1)
162 call sparse%addconnection(noden, nodej, 1)
163 call sparse%addconnection(nodej, noden, 1)
189 integer(I4B) :: noden, nodem, ipos, ignc, jidx, nodej
191 character(len=*),
parameter :: fmterr = &
192 "('GHOST NODE ERROR. Cell ', i0, ' in model ', a, &
193 &' is not connected to cell ', i0, ' in model ', a)"
197 do ignc = 1, this%nexg
198 noden = this%nodem1(ignc) + this%m1%moffset
199 nodem = this%nodem2(ignc) + this%m2%moffset
202 this%idiagn(ignc) = matrix_sln%get_position_diag(noden)
203 this%idiagm(ignc) = matrix_sln%get_position_diag(nodem)
210 this%idxglo(ignc) = matrix_sln%get_position(noden, nodem)
211 this%idxsymglo(ignc) = matrix_sln%get_position(nodem, noden)
214 if (this%idxglo(ignc) == -1)
then
215 write (
errmsg, fmterr) this%nodem1(ignc), trim(this%m1%name), &
216 this%nodem2(ignc), trim(this%m2%name)
228 if (this%implicit)
then
229 do ignc = 1, this%nexg
230 noden = this%nodem1(ignc) + this%m1%moffset
231 nodem = this%nodem2(ignc) + this%m2%moffset
233 do jidx = 1, this%numjs
234 nodej = this%nodesj(jidx, ignc)
235 if (nodej > 0) nodej = nodej + this%m1%moffset
240 this%jposinrown(jidx, ignc) = ipos
242 this%jposinrown(jidx, ignc) = matrix_sln%get_position(noden, nodej)
248 this%jposinrowm(jidx, ignc) = ipos
250 this%jposinrowm(jidx, ignc) = matrix_sln%get_position(nodem, nodej)
268 integer(I4B),
intent(in) :: kiter
271 integer(I4B) :: ignc, ipos
276 gncloop:
do ignc = 1, this%nexg
277 ipos = this%idxglo(ignc)
279 cond = matrix%get_value_pos(ipos)
283 this%cond(ignc) = cond
300 integer(I4B),
intent(in) :: kiter
303 integer(I4B) :: ignc, j, noden, nodem, ipos, jidx, iposjn, iposjm
304 real(DP) :: cond, alpha, aterm, rterm
308 if (this%smgnc)
call this%gnc_fmsav(kiter, matrix)
312 gncloop:
do ignc = 1, this%nexg
313 noden = this%nodem1(ignc)
314 nodem = this%nodem2(ignc)
315 if (this%m1%ibound(noden) == 0 .or. &
316 this%m2%ibound(nodem) == 0) cycle gncloop
317 ipos = this%idxglo(ignc)
318 cond = this%cond(ignc)
319 jloop:
do jidx = 1, this%numjs
320 j = this%nodesj(jidx, ignc)
322 alpha = this%alphasj(jidx, ignc)
323 if (alpha ==
dzero) cycle
325 if (this%implicit)
then
326 iposjn = this%jposinrown(jidx, ignc)
327 iposjm = this%jposinrowm(jidx, ignc)
328 call matrix%add_value_pos(this%idiagn(ignc), aterm)
329 call matrix%add_value_pos(iposjn, -aterm)
330 call matrix%add_value_pos(this%idxsymglo(ignc), -aterm)
331 call matrix%add_value_pos(iposjm, aterm)
333 rterm = aterm * (this%m1%x(noden) - this%m1%x(j))
334 this%m1%rhs(noden) = this%m1%rhs(noden) - rterm
335 this%m2%rhs(nodem) = this%m2%rhs(nodem) + rterm
358 subroutine gnc_fn(this, kiter, matrix_sln, condsat, ihc_opt, &
359 ivarcv_opt, ictm1_opt, ictm2_opt)
365 integer(I4B) :: kiter
367 real(DP),
dimension(:),
intent(in) :: condsat
368 integer(I4B),
dimension(:),
optional :: ihc_opt
369 integer(I4B),
optional :: ivarcv_opt
370 integer(I4B),
dimension(:),
optional :: ictm1_opt
371 integer(I4B),
dimension(:),
optional :: ictm2_opt
373 integer(I4B) :: ignc, jidx, ipos, isympos, ihc, ivarcv
374 integer(I4B) :: nodej, noden, nodem
375 integer(I4B) :: iups, ictup
376 real(DP) :: csat, alpha, consterm, term, derv
377 real(DP) :: xup, topup, botup
382 if (
present(ivarcv_opt)) ivarcv = ivarcv_opt
384 gncloop:
do ignc = 1, this%nexg
385 noden = this%nodem1(ignc)
386 nodem = this%nodem2(ignc)
387 if (this%m1%ibound(noden) == 0 .or. &
388 this%m2%ibound(nodem) == 0) cycle gncloop
393 ipos = this%m1%dis%con%getjaindex(noden, nodem)
394 isympos = this%m1%dis%con%jas(ipos)
395 ihc = this%m1%dis%con%ihc(isympos)
396 csat = condsat(isympos)
403 if (ihc == 0 .and. ivarcv == 0) cycle
407 if (this%m2%x(nodem) > this%m1%x(noden)) iups = 1
412 topup = this%m1%dis%top(noden)
413 botup = this%m1%dis%bot(noden)
415 if (
present(ictm1_opt)) ictup = ictm1_opt(noden)
416 xup = this%m1%x(noden)
418 topup = this%m2%dis%top(nodem)
419 botup = this%m2%dis%bot(nodem)
421 if (
present(ictm2_opt)) ictup = ictm2_opt(nodem)
422 xup = this%m2%x(nodem)
426 if (ictup == 0) cycle
430 topup = min(this%m1%dis%top(noden), this%m2%dis%top(nodem))
431 botup = max(this%m1%dis%bot(noden), this%m2%dis%bot(nodem))
435 jloop:
do jidx = 1, this%numjs
436 nodej = this%nodesj(jidx, ignc)
437 if (nodej == 0) cycle
438 if (this%m1%ibound(nodej) == 0) cycle
439 alpha = this%alphasj(jidx, ignc)
440 if (alpha ==
dzero) cycle
441 consterm = csat * alpha * (this%m1%x(noden) - this%m1%x(nodej))
443 term = consterm * derv
445 call matrix_sln%add_value_pos(this%idiagn(ignc), term)
446 if (this%m2%ibound(nodem) > 0)
then
447 call matrix_sln%add_value_pos(this%idxsymglo(ignc), -term)
449 this%m1%rhs(noden) = this%m1%rhs(noden) + term * this%m1%x(noden)
450 this%m2%rhs(nodem) = this%m2%rhs(nodem) - term * this%m1%x(noden)
452 call matrix_sln%add_value_pos(this%idiagm(ignc), -term)
453 if (this%m1%ibound(noden) > 0)
then
454 call matrix_sln%add_value_pos(this%idxglo(ignc), term)
456 this%m1%rhs(noden) = this%m1%rhs(noden) + term * this%m2%x(nodem)
457 this%m2%rhs(nodem) = this%m2%rhs(nodem) - term * this%m2%x(nodem)
473 integer(I4B),
intent(in) :: ibudfl
476 real(DP) :: deltaQgnc
477 character(len=LINELENGTH) :: nodenstr, nodemstr
479 character(len=*),
parameter :: fmtgnc =
"(i10, 2a10, 2(1pg15.6))"
482 if (ibudfl /= 0 .and. this%iprflow /= 0)
then
483 write (this%iout,
'(//, a)')
'GHOST NODE CORRECTION RESULTS'
484 write (this%iout,
'(3a10, 2a15)')
'GNC NUM',
'NODEN',
'NODEM', &
485 'DELTAQGNC',
'CONDNM'
486 do ignc = 1, this%nexg
487 deltaqgnc = this%deltaQgnc(ignc)
488 call this%m1%dis%noder_to_string(this%nodem1(ignc), nodenstr)
489 call this%m2%dis%noder_to_string(this%nodem2(ignc), nodemstr)
490 write (this%iout, fmtgnc) ignc, trim(adjustl(nodenstr)), &
491 trim(adjustl(nodemstr)), &
492 deltaqgnc, this%cond(ignc)
505 real(DP),
dimension(:),
intent(inout) :: flowja
507 integer(I4B) :: ignc, n1, n2, ipos, isympos
508 real(DP) :: deltaQgnc
511 do ignc = 1, this%nexg
514 n1 = this%nodem1(ignc)
515 n2 = this%nodem2(ignc)
516 deltaqgnc = this%deltaQgnc(ignc)
519 ipos = this%m1%dis%con%getjaindex(n1, n2)
520 isympos = this%m1%dis%con%isym(ipos)
523 flowja(ipos) = flowja(ipos) + deltaqgnc
524 flowja(isympos) = flowja(isympos) - deltaqgnc
543 integer(I4B),
intent(in) :: ignc
545 integer(I4B) :: noden, nodem, nodej, jidx
546 real(dp) :: sigalj, alpha, hd, aterm, cond
552 noden = this%nodem1(ignc)
553 nodem = this%nodem2(ignc)
556 if (this%m1%ibound(noden) /= 0 .and. this%m2%ibound(nodem) /= 0)
then
557 jloop:
do jidx = 1, this%numjs
558 nodej = this%nodesj(jidx, ignc)
559 if (nodej == 0) cycle jloop
560 if (this%m1%ibound(nodej) == 0) cycle jloop
561 alpha = this%alphasj(jidx, ignc)
562 sigalj = sigalj + alpha
563 hd = hd + alpha * this%m1%x(nodej)
565 aterm = sigalj * this%m1%x(noden) - hd
566 cond = this%cond(ignc)
583 call this%NumericalPackageType%allocate_scalars()
586 call mem_allocate(this%implicit,
'IMPLICIT', this%memoryPath)
593 this%implicit = .true.
611 call mem_allocate(this%nodem1, this%nexg,
'NODEM1', this%memoryPath)
612 call mem_allocate(this%nodem2, this%nexg,
'NODEM2', this%memoryPath)
613 call mem_allocate(this%nodesj, this%numjs, this%nexg,
'NODESJ', &
615 call mem_allocate(this%alphasj, this%numjs, this%nexg,
'ALPHASJ', &
617 call mem_allocate(this%cond, this%nexg,
'COND', this%memoryPath)
618 call mem_allocate(this%idxglo, this%nexg,
'IDXGLO', this%memoryPath)
619 call mem_allocate(this%idiagn, this%nexg,
'IDIAGN', this%memoryPath)
620 call mem_allocate(this%idiagm, this%nexg,
'IDIAGM', this%memoryPath)
621 call mem_allocate(this%idxsymglo, this%nexg,
'IDXSYMGLO', this%memoryPath)
622 if (this%implicit)
then
623 call mem_allocate(this%jposinrown, this%numjs, this%nexg,
'JPOSINROWN', &
625 call mem_allocate(this%jposinrowm, this%numjs, this%nexg,
'JPOSINROWM', &
628 call mem_allocate(this%jposinrown, 0, 0,
'JPOSINROWN', this%memoryPath)
629 call mem_allocate(this%jposinrowm, 0, 0,
'JPOSINROWM', this%memoryPath)
651 if (this%inunit > 0)
then
666 call this%NumericalPackageType%da()
683 character(len=LINELENGTH) :: keyword
685 logical :: isfound, endOfBlock
688 call this%parser%GetBlock(
'OPTIONS', isfound, ierr, &
689 supportopenclose=.true., blockrequired=.false.)
693 write (this%iout,
'(1x,a)')
'PROCESSING GNC OPTIONS'
695 call this%parser%GetNextLine(endofblock)
697 call this%parser%GetStringCaps(keyword)
698 select case (keyword)
701 write (this%iout,
'(4x,a)') &
702 'THE LIST OF GHOST-NODE CORRECTIONS WILL BE PRINTED.'
705 write (this%iout,
'(4x,a)') &
706 'DELTAQGNC VALUES WILL BE PRINTED TO THE LIST FILE.'
709 write (this%iout,
'(4x,a)') &
710 'SECOND ORDER CORRECTION WILL BE APPLIED.'
712 this%implicit = .false.
713 write (this%iout,
'(4x,a)')
'GHOST NODE CORRECTION IS EXPLICIT.'
715 write (
errmsg,
'(a,a)')
'Unknown GNC option: ', &
718 call this%parser%StoreErrorUnit()
721 write (this%iout,
'(1x,a)')
'END OF GNC OPTIONS'
725 if (this%implicit) this%iasym = 1
742 character(len=LINELENGTH) :: keyword
744 logical :: isfound, endOfBlock
747 call this%parser%GetBlock(
'DIMENSIONS', isfound, ierr, &
748 supportopenclose=.true.)
752 write (this%iout,
'(1x,a)')
'PROCESSING GNC DIMENSIONS'
754 call this%parser%GetNextLine(endofblock)
756 call this%parser%GetStringCaps(keyword)
757 select case (keyword)
759 this%nexg = this%parser%GetInteger()
760 write (this%iout,
'(4x,a,i7)')
'NUMGNC = ', this%nexg
762 this%numjs = this%parser%GetInteger()
763 write (this%iout,
'(4x,a,i7)')
'NUMAPHAJ = ', this%numjs
765 write (
errmsg,
'(a,a)')
'Unknown GNC dimension: ', &
768 call this%parser%StoreErrorUnit()
771 write (this%iout,
'(1x,a)')
'END OF GNC DIMENSIONS'
773 call store_error(
'Required DIMENSIONS block not found.', terminate=.true.)
791 character(len=LINELENGTH) :: line, nodestr, fmtgnc, cellid, &
793 integer(I4B) :: lloc, ierr, ival
794 integer(I4B) :: ignc, jidx, nodeun, nodeum, nerr
795 integer(I4B),
dimension(:),
allocatable :: nodesuj
796 logical :: isfound, endOfBlock
799 write (fmtgnc,
'("(2i10,",i0,"i10,",i0, "(1pg15.6))")') this%numjs, &
804 allocate (nodesuj(this%numjs))
807 call this%parser%GetBlock(
'GNCDATA', isfound, ierr, supportopenclose=.true.)
811 write (this%iout,
'(1x,a)')
'PROCESSING GNCDATA'
812 do ignc = 1, this%nexg
813 call this%parser%GetNextLine(endofblock)
815 call this%parser%GetCurrentLine(line)
819 call this%parser%GetCellid(this%m1%dis%ndim, cellidn)
820 nodeun = this%m1%dis%nodeu_from_cellid(cellidn, this%parser%iuactive, &
824 call this%nodeu_to_noder(nodeun, this%nodem1(ignc), this%m1)
827 call this%parser%GetCellid(this%m2%dis%ndim, cellidm)
828 nodeum = this%m2%dis%nodeu_from_cellid(cellidm, this%parser%iuactive, &
832 call this%nodeu_to_noder(nodeum, this%nodem2(ignc), this%m2)
835 do jidx = 1, this%numjs
837 call this%parser%GetCellid(this%m1%dis%ndim, cellid)
838 ival = this%m1%dis%nodeu_from_cellid(cellid, this%parser%iuactive, &
839 this%iout, allow_zero=.true.)
842 call this%nodeu_to_noder(ival, this%nodesj(jidx, ignc), this%m1)
844 this%nodesj(jidx, ignc) = 0
849 do jidx = 1, this%numjs
850 this%alphasj(jidx, ignc) = this%parser%GetDouble()
854 if (this%iprpak /= 0) &
855 write (this%iout, fmtgnc) nodeun, nodeum, &
856 (nodesuj(jidx), jidx=1, this%numjs), &
857 (this%alphasj(jidx, ignc), jidx=1, this%numjs)
860 if (this%nodem1(ignc) <= 0)
then
861 call this%m1%dis%nodeu_to_string(nodeun, nodestr)
863 trim(adjustl(this%m1%name))// &
864 ' Cell is outside active grid domain: '// &
865 trim(adjustl(nodestr))
870 if (this%nodem2(ignc) <= 0)
then
871 call this%m2%dis%nodeu_to_string(nodeum, nodestr)
873 trim(adjustl(this%m2%name))// &
874 ' Cell is outside active grid domain: '// &
875 trim(adjustl(nodestr))
880 do jidx = 1, this%numjs
881 if (this%nodesj(jidx, ignc) < 0)
then
882 call this%m1%dis%nodeu_to_string(nodesuj(jidx), nodestr)
884 trim(adjustl(this%m1%name))// &
885 ' Cell is outside active grid domain: '// &
886 trim(adjustl(nodestr))
896 call store_error(
'Errors encountered in GNC input file.')
897 call this%parser%StoreErrorUnit()
900 write (this%iout,
'(1x,a)')
'END OF GNCDATA'
902 write (
errmsg,
'(a)')
'Required GNCDATA block not found.'
904 call this%parser%StoreErrorUnit()
923 integer(I4B),
intent(in) :: nodeu
924 integer(I4B),
intent(inout) :: noder
927 if (nodeu < 1 .or. nodeu > model%dis%nodesuser)
then
929 trim(adjustl(model%name))// &
930 ' node number < 0 or > model nodes: ', nodeu
933 noder = model%dis%get_nodenumber(nodeu, 0)
This module contains block parser methods.
This module contains simulation constants.
integer(i4b), parameter linelength
maximum length of a standard line
real(dp), parameter dzero
real constant zero
subroutine read_dimensions(this)
Single Model GNC Read Dimensions.
subroutine nodeu_to_noder(this, nodeu, noder, model)
Convert the user-based node number into a reduced number.
subroutine gnc_da(this)
Deallocate memory.
subroutine gnc_fn(this, kiter, matrix_sln, condsat, ihc_opt, ivarcv_opt, ictm1_opt, ictm2_opt)
Fill GNC Newton terms.
subroutine, public gnc_cr(gncobj, name_parent, inunit, iout)
Create new GNC exchange object.
subroutine read_data(this)
Read a GNCDATA block.
subroutine allocate_scalars(this)
Allocate gnc scalar variables.
subroutine gnc_fc(this, kiter, matrix)
Fill matrix terms.
subroutine gnc_df(this, m1, m2)
Initialize a gnc object.
subroutine read_options(this)
Read a gnc options block.
subroutine allocate_arrays(this)
Allocate gnc scalar variables.
subroutine gnc_ot(this, ibudfl)
Single Model GNC Output.
subroutine gnc_cq(this, flowja)
Add GNC to flowja.
subroutine gnc_mc(this, matrix_sln)
Single or Two-Model GNC Map Connections.
subroutine gnc_fmsav(this, kiter, matrix)
Store the n-m Picard conductance in cond prior to the Newton terms terms being added.
subroutine gnc_ac(this, sparse)
Single or Two-Model GNC Add Connections.
real(dp) function deltaqgnc(this, ignc)
Single Model deltaQgnc (ghost node correction flux)
This module defines variable data types.
This module contains the base numerical package type.
This module contains simulation methods.
subroutine, public store_error(msg, terminate)
Store an error message.
integer(i4b) function, public count_errors()
Return number of errors.
subroutine, public store_error_unit(iunit, terminate)
Store the file unit number.
This module contains simulation variables.
character(len=maxcharlen) errmsg
error message string
real(dp) function squadraticsaturationderivative(top, bot, x, eps)
@ brief Derivative of the quadratic saturation function