From 69d3b6c0046d0f5ac9f547dde7642a59150e00d6 Mon Sep 17 00:00:00 2001 From: "Neil N. Carlson" <nnc@lanl.gov> Date: Thu, 18 Feb 2021 12:04:24 -0700 Subject: [PATCH 1/6] WIP: Prototype HT solver using vector-class-based DAE integrator (#450) This adds a new HT solver that uses a redesigned DAE integrator that uses abstract vectors rather than intrinsic rank-1 real arrays. This pure HT solver is used in place of the HTSD solver whenever species diffusion is not enabled. This is working for all HT and HT+flow test problems that had used the HTSD solver, *except* for those including view factor radiation. Next step is to incorporate VF radiation. --- src/truchas/ode/CMakeLists.txt | 4 + src/truchas/ode/new_idaesol_type.F90 | 742 ++++++++++++++++++ src/truchas/ode/new_nka_type.F90 | 523 ++++++++++++ src/truchas/ode/new_state_history_type.F90 | 433 ++++++++++ src/truchas/ode/vector_class.F90 | 277 +++++++ .../heat_species_transport/CMakeLists.txt | 9 + .../diffusion_solver.F90 | 59 +- .../ht_ic_solver_type.F90 | 473 +++++++++++ .../ht_idaesol_model_type.F90 | 112 +++ .../ht_model_factory.F90 | 357 +++++++++ .../heat_species_transport/ht_model_type.F90 | 295 +++++++ .../heat_species_transport/ht_norm_type.F90 | 168 ++++ .../heat_species_transport/ht_precon_type.F90 | 337 ++++++++ .../ht_solver_factory.F90 | 78 ++ .../heat_species_transport/ht_solver_type.F90 | 319 ++++++++ .../heat_species_transport/ht_vector_type.F90 | 229 ++++++ test/enthalpy/test2.py | 4 +- 17 files changed, 4415 insertions(+), 4 deletions(-) create mode 100644 src/truchas/ode/new_idaesol_type.F90 create mode 100644 src/truchas/ode/new_nka_type.F90 create mode 100644 src/truchas/ode/new_state_history_type.F90 create mode 100644 src/truchas/ode/vector_class.F90 create mode 100644 src/truchas/physics/heat_species_transport/ht_ic_solver_type.F90 create mode 100644 src/truchas/physics/heat_species_transport/ht_idaesol_model_type.F90 create mode 100644 src/truchas/physics/heat_species_transport/ht_model_factory.F90 create mode 100644 src/truchas/physics/heat_species_transport/ht_model_type.F90 create mode 100644 src/truchas/physics/heat_species_transport/ht_norm_type.F90 create mode 100644 src/truchas/physics/heat_species_transport/ht_precon_type.F90 create mode 100644 src/truchas/physics/heat_species_transport/ht_solver_factory.F90 create mode 100644 src/truchas/physics/heat_species_transport/ht_solver_type.F90 create mode 100644 src/truchas/physics/heat_species_transport/ht_vector_type.F90 diff --git a/src/truchas/ode/CMakeLists.txt b/src/truchas/ode/CMakeLists.txt index e3af740c6..50f6a8122 100644 --- a/src/truchas/ode/CMakeLists.txt +++ b/src/truchas/ode/CMakeLists.txt @@ -6,4 +6,8 @@ target_sources(truchas PRIVATE bdf2/bdf2_profiling.F90 bdf2/solution_history.F90 idaesol_type.F90 + vector_class.F90 + new_nka_type.F90 + new_state_history_type.F90 + new_idaesol_type.F90 ) diff --git a/src/truchas/ode/new_idaesol_type.F90 b/src/truchas/ode/new_idaesol_type.F90 new file mode 100644 index 000000000..224786dd9 --- /dev/null +++ b/src/truchas/ode/new_idaesol_type.F90 @@ -0,0 +1,742 @@ +!! +!! IDAESOL_TYPE +!! +!! A solver for index-1 DAE in implicit form using the BDF2 method. +!! +!! This F2008 version is adapted from much earlier F95 and F77 implementations +!! and was obtained from http://sourceforge.net/projects/calliope.mfeproject.p/ +!! +!! Neil N. Carlson <nnc@lanl.gov> +!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! +!! Copyright (c) 2013 Neil N. Carlson <neil.n.carlson@gmail.com> +!! +!! Permission is hereby granted, free of charge, to any person obtaining a +!! copy of this software and associated documentation files (the "Software"), +!! to deal in the Software without restriction, including without limitation +!! the rights to use, copy, modify, merge, publish, distribute, sublicense, +!! and/or sell copies of the Software, and to permit persons to whom the +!! Software is furnished to do so, subject to the following conditions: +!! +!! The above copyright notice and this permission notice shall be included +!! in all copies or substantial portions of the Software. +!! +!! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +!! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +!! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +!! THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +!! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +!! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +!! DEALINGS IN THE SOFTWARE. +!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +#include "f90_assert.fpp" + +module new_idaesol_type + + use,intrinsic :: iso_fortran_env, only: r8 => real64 + use new_state_history_type + use new_nka_type + use vector_class + implicit none + private + + type, public :: idaesol + private + class(idaesol_model), pointer :: model => null() + integer :: seq = -1 ! number of steps taken + real(r8) :: hlast ! last step size + real(r8) :: hpc ! step size built into the current preconditioner + logical :: usable_pc = .false. ! whether the current preconditioner is usable + integer :: pc_age, pc_freq + integer :: freeze_count = 0 ! don't increase step size for this number of steps + integer :: mitr = 5 ! maximum number of nonlinear iterations + real(r8) :: ntol = 0.1_r8 ! nonlinear solver error tolerance (relative to 1) + type(nka), allocatable :: nka ! nonlinear Krylov accelerator + type(state_history) :: uhist ! solution history structure + + !! Persistent temporary workspace + class(vector), allocatable :: u, u0, up + class(vector), allocatable :: du, udot ! local to bce_step + + !! Perfomance counters + integer :: pcfun_calls = 0 ! number of calls to PCFUN + integer :: updpc_calls = 0 ! number of calls to UPDPC + integer :: updpc_failed = 0 ! number of UPDPC calls returning an error + integer :: retried_bce = 0 ! number of retried BCE steps + integer :: failed_bce = 0 ! number of completely failed BCE steps + integer :: rejected_steps = 0 ! number of steps rejected on error tolerance + real(r8) :: hmin = huge(1.0_r8) ! minimum step size used on a successful step + real(r8) :: hmax = tiny(1.0_r8) ! maximum step size used on a successful step + + !! Diagnostics + integer :: unit = 0 + logical :: verbose = .false. + contains + procedure :: init + procedure :: set_initial_state + procedure :: integrate => bdf2_step_driver + procedure :: step + procedure :: commit_state + procedure :: get_interpolated_state + procedure :: get_last_state_copy + procedure :: get_last_state_view + procedure :: last_time + procedure :: last_step_size + procedure :: set_verbose_stepping + procedure :: set_quiet_stepping + generic :: write_metrics => write_metrics_unit, write_metrics_string + procedure, private :: write_metrics_unit + procedure, private :: write_metrics_string + procedure :: get_stepping_statistics + end type idaesol + + type, abstract, public :: idaesol_model + contains + procedure(alloc_vector), deferred :: alloc_vector + procedure(compute_f), deferred :: compute_f + procedure(apply_precon), deferred :: apply_precon + procedure(compute_precon), deferred :: compute_precon + procedure(du_norm), deferred :: du_norm + end type + + abstract interface + subroutine alloc_vector(this, vec) + import idaesol_model, vector + class(idaesol_model), intent(in) :: this + class(vector), allocatable, intent(out) :: vec + end subroutine + subroutine compute_f(this, t, u, udot, f) + import idaesol_model, vector, r8 + class(idaesol_model) :: this + real(r8), intent(in) :: t + class(vector), intent(inout) :: u, udot + class(vector), intent(inout) :: f + end subroutine + subroutine apply_precon(this, t, u, f) + import idaesol_model, vector, r8 + class(idaesol_model) :: this + real(r8), intent(in) :: t + class(vector), intent(inout) :: u + class(vector), intent(inout) :: f + end subroutine + subroutine compute_precon(this, t, u, dt) + import idaesol_model, vector, r8 + class(idaesol_model) :: this + real(r8), intent(in) :: t, dt + class(vector), intent(inout) :: u + end subroutine + subroutine du_norm(this, t, u, du, error) + import :: idaesol_model, vector, r8 + class(idaesol_model) :: this + real(r8), intent(in) :: t + class(vector), intent(in) :: u, du + real(r8), intent(out) :: error + end subroutine + end interface + + real(r8), parameter, private :: RMIN = 0.25_r8 + real(r8), parameter, private :: RMAX = 4.0_r8 + real(r8), parameter, private :: MARGIN = 3.0_r8 + + !! Successful STATUS return codes: + integer, parameter, public :: SOLVED_TO_TOUT = 1 + integer, parameter, public :: SOLVED_TO_NSTEP = 2 + + !! Unsuccessful STATUS return codes: + integer, parameter, public :: BAD_INPUT = -1 + integer, parameter, public :: STEP_FAILED = -2 + integer, parameter, public :: STEP_SIZE_TOO_SMALL = -3 + +contains + + subroutine init(this, model, params, stat, errmsg) + + use parameter_list_type + + class(idaesol), intent(out) :: this + class(idaesol_model), intent(in), target :: model + type(parameter_list) :: params + integer, intent(out) :: stat + character(:), allocatable, intent(out) :: errmsg + + character(:), allocatable :: context + integer :: maxv + real(r8) :: vtol + + this%model => model + + context = 'processing ' // params%name() // ': ' + call params%get('nlk-max-iter', this%mitr, default=5, stat=stat, errmsg=errmsg) + if (stat /= 0) then + errmsg = context//errmsg + return + end if + if (this%mitr < 2) then + stat = 1 + errmsg = context//'"nlk-max-iter" must be > 1' + return + end if + + call params%get('nlk-tol', this%ntol, default=0.1_r8, stat=stat, errmsg=errmsg) + if (stat /= 0) then + errmsg = context//errmsg + return + end if + if (this%ntol <= 0.0_r8 .or. this%ntol > 1.0_r8) then + stat = 1 + errmsg = context//'"nlk-tol" must be > 0.0 and <= 1.0' + return + end if + + call params%get('nlk-max-vec', maxv, default=this%mitr-1, stat=stat, errmsg=errmsg) + if (stat /= 0) then + errmsg = context//errmsg + return + end if + if (maxv < 0) then + stat = 1 + errmsg = context//'"nlk-max-vec" must be >= 0' + return + end if + maxv = min(maxv, this%mitr-1) + + call params%get('nlk-vec-tol', vtol, default=0.01_r8, stat=stat, errmsg=errmsg) + if (stat /= 0) then + errmsg = context//errmsg + return + end if + if (vtol <= 0.0_r8) then + stat = 1 + errmsg = context//'"nlk-vec-tol" must be > 0.0' + return + end if + + call params%get('pc-freq', this%pc_freq, default=huge(this%pc_freq), stat=stat, errmsg=errmsg) + if (stat /= 0) then + errmsg = context//errmsg + return + end if + if (this%pc_freq < 1) then + stat = 1 + errmsg = context//'"pc-freq" must be > 0' + return + end if + + call model%alloc_vector(this%u) + call model%alloc_vector(this%u0) + call model%alloc_vector(this%up) + call model%alloc_vector(this%du) + call model%alloc_vector(this%udot) + + !! Initialize the NKA structure. + if (maxv > 0) then + allocate(this%nka) + call this%nka%init(this%u, maxv) + call this%nka%set_vec_tol(vtol) + end if + + !! We need to maintain 3 solution vectors for quadratic extrapolation. + call this%uhist%init(3, this%u) + + end subroutine init + + + subroutine set_initial_state(this, t, u, udot) + class(idaesol), intent(inout) :: this + real(r8), intent(in) :: t + class(vector), intent(in) :: u, udot + call this%uhist%flush(t, u, udot) + this%seq = 0 + this%usable_pc = .false. + end subroutine set_initial_state + + + subroutine bdf2_step_driver(this, hnext, status, nstep, tout, hmin, hmax, mtry) + + class(idaesol), intent(inout) :: this + real(r8), intent(inout) :: hnext + integer, intent(out) :: status + integer, intent(in) :: nstep, mtry + real(r8), intent(in) :: tout, hmin, hmax + optional :: nstep, tout, hmin, hmax, mtry + + integer :: max_step, max_try, step, errc + real(r8) :: tlast, h, t_out, h_min, h_max + + ASSERT(this%seq >= 0) + + !!! + !!! PROCESS THE INPUT AND CHECK IT FOR CORRECTNESS + + status = 0 + + !! Set the maximum number of time steps; default is unlimited. + max_step = huge(1) + if (present(nstep)) max_step = nstep + if (max_step < 1) status = BAD_INPUT + + !! Set the target integration time; default is +infinity. + t_out = huge(1.0_r8) + if (present(tout)) t_out = tout + if (t_out <= this%uhist%last_time()) status = BAD_INPUT + + !! Verify that at least one of NSTEP and TOUT were specified. + if (.not.present(nstep) .and. .not.present(tout)) status = BAD_INPUT + + !! Set the bounds on the step size; default is none. + h_min = tiny(1.0_r8) + h_max = huge(1.0_r8) + if (present(hmin)) h_min = hmin + if (present(hmax)) h_max = hmax + if (h_min < 0.0_r8 .or. hnext < h_min .or. hnext > h_max) status = BAD_INPUT + + !! Set the maximum number of attempts allowed for a step. + max_try = 10 + if (present(mtry)) max_try = mtry + if (max_try < 1) status = BAD_INPUT + + if (status == BAD_INPUT) return + + !!! + !!! BEGIN TIME STEPPING + + step = 0 + STEP_LOOP: do + + h = hnext + tlast = this%uhist%last_time() + + !! Check for a normal return before proceeding. + if (t_out <= tlast) then + status = SOLVED_TO_TOUT + exit STEP_LOOP + end if + if (step >= max_step) then + status = SOLVED_TO_NSTEP + exit STEP_LOOP + end if + + step = step + 1 + + call bdf2_step(this, h, h_min, max_try, this%u, hnext, errc) + if (errc /= 0) then + status = STEP_FAILED + exit STEP_LOOP ! step failed + end if + + !! BDF2 step was successful; commit the solution. + call commit_state(this, tlast+h, this%u) + + !! Set the next step size. + hnext = min(h_max, hnext) + if (this%verbose) write(this%unit,fmt=1) hnext/h + + end do STEP_LOOP + + 1 format(2x,'Changing H by a factor of ',f6.3) + + end subroutine bdf2_step_driver + + + subroutine commit_state(this, t, u) + + class(idaesol), intent(inout) :: this + real(r8), intent(in) :: t + class(vector), intent(in) :: u + + real(r8) :: h + + h = t - this%uhist%last_time() + ASSERT(h > 0.0_r8) + call this%uhist%record_state(t, u) + + this%hlast = h + this%seq = this%seq + 1 + this%freeze_count = max(0, this%freeze_count - 1) + + this%hmin = min(h, this%hmin) + this%hmax = max(h, this%hmax) + + end subroutine commit_state + + + subroutine bdf2_step(this, h, hmin, mtry, u, hnext, errc) + + type(idaesol), intent(inout) :: this + real(r8), intent(inout) :: h + real(r8), intent(in) :: hmin + integer, intent(in) :: mtry + class(vector), intent(inout) :: u + real(r8), intent(out) :: hnext + integer, intent(out) :: errc + + integer :: try + + errc = 0 + if (hmin < 0.0_r8 .or. h < hmin) errc = BAD_INPUT + if (mtry < 1) errc = BAD_INPUT + if (errc == BAD_INPUT) return + + try = 0 + do + try = try + 1 + + !! Check for too many attempts at a single step. + if (try > mtry) then + errc = STEP_FAILED + exit + end if + + !! Check for a too-small step size. + if (h < hmin) then + errc = STEP_SIZE_TOO_SMALL + exit + end if + + !! Attempt a BDF2 step. + call step(this, this%uhist%last_time()+h, u, hnext, errc) + if (errc == 0) exit + + !! Step failed; try again with the suggested step size. + if (this%verbose) write(this%unit,fmt=1) hnext/h + h = hnext + + end do + + 1 format(2x,'Changing H by a factor of ',f6.3) + + end subroutine bdf2_step + + subroutine step(this, t, u, hnext, errc) + + class(idaesol), intent(inout) :: this + real(r8), intent(in) :: t + class(vector), intent(inout) :: u ! data is intent(out) + real(r8), intent(out) :: hnext + integer, intent(out) :: errc + + real(r8) :: eta, etah, h, t0, tlast, perr, dt(3) + logical :: fresh_pc, predictor_error + + ASSERT(this%seq >= 0) + + tlast = this%uhist%last_time() + h = t - tlast + INSIST(h > 0) + + !! Predicted solution and base point for BCE step. + if (this%uhist%depth() == 2) then ! trapezoid method + if (this%verbose) write(this%unit,fmt=1) this%seq+1, tlast, h + etah = 0.5_r8 * h + t0 = tlast + etah + call this%uhist%interp_state(t, this%up, order=1) + call this%uhist%interp_state(t0, this%u0, order=1) + else ! BDF2 + if (this%verbose) write(this%unit,fmt=2) this%seq+1, tlast, h + eta = (this%hlast + h) / (this%hlast + 2.0_r8 * h) + etah = eta * h + t0 = tlast + (1.0_r8 - eta)*h + call this%uhist%interp_state(t, this%up, order=2) + call this%uhist%interp_state(t0, this%u0, order=1) + end if + + fresh_pc = .false. + + if (this%usable_pc) then + this%pc_age = this%pc_age + 1 + if (this%pc_age >= this%pc_freq) this%usable_pc = .false. + end if + + !! If the PC step size is too different than the current step size we tag + !! it as unusable in order to preempt a possible nonlinear solve failure. + if (this%usable_pc) then + if (this%hpc/etah > 1.0_r8 + MARGIN) this%usable_pc = .false. + if (etah/this%hpc > 1.0_r8 + MARGIN) this%usable_pc = .false. + end if + + BCE: do + + !! Update the preconditioner if necessary. + if (.not.this%usable_pc) then + this%updpc_calls = this%updpc_calls + 1 + call this%model%compute_precon(t, this%up, etah) + if (this%verbose) write(this%unit,fmt=3) t + this%hpc = etah + this%usable_pc = .true. + this%pc_age = 0 + fresh_pc = .true. + end if + + !! Solve the nonlinear BCE system. + call u%copy(this%up) ! Initial solution guess is the predictor. + call bce_step(this, t, etah, this%u0, u, errc) + if (errc == 0) exit BCE ! the BCE step was successful. + + if (fresh_pc) then ! preconditioner was fresh; cut h and return error condition. + this%failed_bce = this%failed_bce + 1 + hnext = 0.5_r8 * h + this%freeze_count = 1 + errc = 1 + return + else ! update the preconditioner and retry the nonlinear solve. + this%retried_bce = this%retried_bce + 1 + this%usable_pc = .false. + cycle BCE + end if + + end do BCE + + predictor_error = (this%seq >= 3) + + if (predictor_error) then + + !! Predictor error control. + !du = u - up + call this%du%copy(u) + call this%du%update(-1.0_r8, this%up) + call this%model%du_norm(t, u, this%du, perr) + if (perr < 4.0_r8) then ! accept the step. + if (this%verbose) write(this%unit,fmt=4) perr + errc = 0 + else ! reject the step; cut h and return error condition. + this%rejected_steps = this%rejected_steps + 1 + if (this%verbose) write(this%unit,fmt=5) perr + hnext = 0.5_r8 * h + this%freeze_count = 1 + errc = 2 + return + end if + + !! Select the next step size based on the predictor error and past step + !! size history, but don't change the step size by too great a factor. + dt(1) = h + dt(2:) = h + this%uhist%time_deltas() + call select_step_size(dt, perr, hnext) + hnext = max(RMIN*h, min(RMAX*h, hnext)) + if (this%freeze_count /= 0) hnext = min(h, hnext) + + else + + if (this%verbose) write(this%unit,fmt=6) + hnext = h + errc = 0 + + end if + + 1 format(/,'TRAP step ',i6,': T=',es22.15,', H=',es22.15) + 2 format(/,'BDF2 step ',i6,': T=',es22.15,', H=',es22.15) + 3 format(2x,'Preconditioner updated at T=',es22.15) + 4 format(2x,'Step accepted: perr=',es22.15) + 5 format(2x,'Step REJECTED: perr=',es22.15) + 6 format(2x,'Step accepted: no local error control') + + end subroutine step + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !! + !! SELECT_STEP_SIZE -- Choose a new time step. + !! + + subroutine select_step_size(dt, perr, h) + + real(r8), intent(in) :: dt(:), perr + real(r8), intent(out) :: h + + real(r8), parameter :: tol = 0.001_r8 + real(r8) :: a, dh, phi, dphi + + ASSERT( size(dt) == 3 ) + + a = 0.5_r8*dt(1)*dt(2)*dt(3)/max(perr,0.001_r8) + h = dt(1) + + do ! until converged -- DANGEROUS! + + phi = h*(h + dt(1))*(h + dt(2)) - a + dphi = (2.0_r8*h + dt(1))*(h + dt(2)) + h*(h + dt(1)) + + dh = phi / dphi + h = h - dh + if (abs(dh) / h < tol) exit + + end do + + end subroutine select_step_size + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !! + !! SOLVE_BCE_AIN -- Solve the Backward Cauchy-Euler system using AIN. + !! + !! The backward Cauchy-Euler (BCE) method applied to the implicit DAE + !! + !! f(t,u,u') = 0 + !! + !! yields a nonlinear system of equations for advancing the solution from a + !! given state u0 at time t - h to the unknown solution u at time t, + !! + !! f(t,u,(u-u0)/h) = 0. + !! + !! This subroutine solves this nonlinear system using an accelerated fixed + !! point iteration [1] for the preconditioned system + !! g(u) = pc(f(t,u,(u-u0)/h)) = 0: + !! + !! u given + !! Do until converged: + !! du <-- g(u) + !! du <-- NKA(du) + !! u <-- u - du + !! End do + !! + !! The procedure NKA uses information about g' gleaned from the unaccelerated + !! correction du=g(u) and previous g values to compute an improved correction. + !! The preconditioning function pc() is typically an approximate solution + !! of the Newton correction equation J*du = f(t,u,(u-u0)/h) where J is an + !! approximation to the Jacobian of f(t,u,(u-u0)/h) as a function of u. Thus + !! this method can be regarded as an accelerated inexact Newton (AIN) method. + !! + !! The dummy procedure PCFUN evaluates the preconditioned function g. + !! + !! [1] N.N.Carlson and K.Miller, "Design and application of a gradient- + !! weighted moving finite element code I: in one dimension", SIAM J. + !! Sci. Comput;, 19 (1998), pp. 728-765.. + !! + + subroutine bce_step(this, t, h, u0, u, errc) + + type(idaesol), intent(inout) :: this + real(r8), intent(in) :: t, h + class(vector), intent(in) :: u0 + class(vector), intent(inout) :: u + integer, intent(out) :: errc + + integer :: itr + real(r8) :: error + + if (allocated(this%nka)) call this%nka%restart + + itr = 0 + do + + if (itr >= this%mitr) then ! too many nonlinear iterations + if (this%verbose) write(this%unit,fmt=1) itr, error + errc = 1 + exit + end if + + itr = itr + 1 + + !! Evaluate the nonlinear function and precondition it. + this%pcfun_calls = this%pcfun_calls + 1 + + !call this%model%compute_f(t, u, (u-u0)/h, du) + call this%udot%copy(u) + call this%udot%update(-1.0_r8, u0) + call this%udot%scale(1.0_r8/h) + call this%model%compute_f(t, u, this%udot, this%du) + call this%model%apply_precon(t, u, this%du) + + !! NKA accelerated correction. + if (allocated(this%nka)) call this%nka%accel_update(this%du) + + !! Next solution iterate. + !u = u - du + call u%update(-1.0_r8, this%du) + + !! Error estimate. + call this%model%du_norm(t, u, this%du, error) + if (this%verbose) write(this%unit,fmt=3) itr, error + + !! Check for convergence. + if (((error < this%ntol) .and. (itr > 1)) .or. (error < 0.01_r8 * this%ntol)) then + if (this%verbose) write(this%unit,fmt=2) itr, error + errc = 0 + exit + end if + + end do + + 1 format(2x,'NLK BCE solve FAILED: ',i3,' iterations (max), error=',es22.15) + 2 format(2x,'NLK BCE solve succeeded: ',i3,' iterations, error=',es22.15) + 3 format(2x,i3,': error=',es22.15) + + end subroutine bce_step + + subroutine set_verbose_stepping(this, unit) + class(idaesol), intent(inout) :: this + integer, intent(in) :: unit + this%unit = unit + this%verbose = .true. + end subroutine + + subroutine set_quiet_stepping(this) + class(idaesol), intent(inout) :: this + this%verbose = .false. + end subroutine + + subroutine get_last_state_view(this, view) + class(idaesol), intent(in) :: this + class(vector), pointer :: view + call this%uhist%get_last_state_view(view) + end subroutine + + subroutine get_last_state_copy(this, copy) + class(idaesol), intent(in) :: this + class(vector), intent(inout) :: copy ! data intent(out) + call this%uhist%get_last_state_copy(copy) + end subroutine + + function last_time(this) result(t) + class(idaesol), intent(in) :: this + real(r8) :: t + t = this%uhist%last_time() + end function + + function last_step_size(this) result(h) + class(idaesol), intent(in) :: this + real(r8) :: h + h = this%hlast + end function + + subroutine get_interpolated_state(this, t, u) + class(idaesol), intent(in) :: this + real(r8), intent(in) :: t + class(vector), intent(inout) :: u ! data intent(out) + call this%uhist%interp_state(t, u) + end subroutine + + subroutine write_metrics_unit(this, unit) + class(idaesol), intent(in) :: this + integer, intent(in) :: unit + character(80) :: string(2) + call write_metrics_string(this, string) + write(unit,'(/,(a))') string + end subroutine + + subroutine write_metrics_string(this, string) + class(idaesol), intent(in) :: this + character(*), intent(out) :: string(:) + ASSERT(size(string) == 2) + write(string(1),fmt='(a,i6,a,es11.5,a,es9.3)') & + 'STEP=', this%seq, ', T=', this%uhist%last_time(), ', H=', this%hlast + write(string(2),fmt='(a,i7.7,":",i5.5,a,5(i4.4,:,":"))') & + 'NFUN:NPC=', this%pcfun_calls, this%updpc_calls, & + ', NPCF:NNR:NNF:NSR=', this%updpc_failed, & + this%retried_bce, this%failed_bce, this%rejected_steps + end subroutine + + subroutine get_stepping_statistics(this, counters) + class(idaesol), intent(in) :: this + integer, intent(out) :: counters(:) + ASSERT(size(counters) == 6) + counters(1) = this%pcfun_calls + counters(2) = this%updpc_calls + counters(3) = this%updpc_failed + counters(4) = this%retried_bce + counters(5) = this%failed_bce + counters(6) = this%rejected_steps + end subroutine + +end module new_idaesol_type diff --git a/src/truchas/ode/new_nka_type.F90 b/src/truchas/ode/new_nka_type.F90 new file mode 100644 index 000000000..45e97ddd2 --- /dev/null +++ b/src/truchas/ode/new_nka_type.F90 @@ -0,0 +1,523 @@ +!! +!! NKA_TYPE +!! +!! Neil N. Carlson <neil.n.carlson@gmail.com> +!! +!! This module implements the nonlinear Krylov accelerator (NKA) introduced +!! in [1] for fixed point or Picard iterations. Placed in the iteration loop, +!! this black-box accelerator listens to the sequence of solution updates and +!! replaces them with accelerated updates. More generally, NKA can accelerate +!! typical quasi-Newton iterations, which can usually be viewed as a fixed +!! point iteration for a preconditioned function. +!! +!! This is a Fortran 2008 adaptation of the original Fortran 95 version that +!! implements the methods as type bound procedures. +!! +!! [1] N.N.Carlson and K.Miller, "Design and application of a gradient- +!! weighted moving finite element code I: in one dimension", SIAM J. +!! Sci. Comput;, 19 (1998), pp. 728-765. See section 9. +!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! +!! Copyright (c) 2010, 2011 Neil N. Carlson +!! +!! Permission is hereby granted, free of charge, to any person obtaining a +!! copy of this software and associated documentation files (the "Software"), +!! to deal in the Software without restriction, including without limitation +!! the rights to use, copy, modify, merge, publish, distribute, sublicense, +!! and/or sell copies of the Software, and to permit persons to whom the +!! Software is furnished to do so, subject to the following conditions: +!! +!! The above copyright notice and this permission notice shall be included +!! in all copies or substantial portions of the Software. +!! +!! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +!! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +!! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +!! THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +!! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +!! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +!! DEALINGS IN THE SOFTWARE. +!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! +!! PROGRAMING INTERFACE +!! +!! This module defines the derived data type NKA, which encapsulates the state +!! of the acceleration procedure, with the following methods as type bound +!! procedures. +!! +!! INIT(VLEN, MVEC) initializes the object to handle as many as MVEC vectors +!! of length VLEN. +!! +!! SET_VEC_TOL(VTOL) sets the vector drop tolerance. A vector is dropped +!! from the acceleration subspace when the sine of the angle between the +!! vector and the subspace spanned by the preceding vectors is less than +!! this value. If not set, the default value 0.01 is used. +!! +!! SET_DOT_PROD(DOT_PROD) sets the procedure used to compute vector dot +!! products. DOT_PROD is a procedure pointer with the same interface +!! as the intrinsic function DOT_PRODUCT with real rank-1 array arguments. +!! If not set, the default is to use the instrinsic DOT_PRODUCT. In a +!! parallel context where the vector components are distributed across +!! processes, a global dot product procedure needs to be supplied that +!! performs the necessary communication internally. +!! +!! ACCEL_UPDATE(F) takes the function value F, which would be the update +!! vector in a fixed point iteration, and overwrites it with the accelerated +!! update computed from the acceleration subspace stored in the object. +!! This subspace is updated prior to computing the update using F and +!! the previous function value and update that were cached on the +!! preceding call to ACCEL_UPDATE, if any. The input F and returned update +!! are cached for use by the next call to ACCEL_UPDATE. +!! +!! RESTART() flushes the acceleration subspace, returning the object to its +!! state after the call to INIT; the next call to ACCEL_UPDATE begins the +!! process of accumulating a new acceleration subspace. Typical usage is +!! to call RESTART at the start of each nonlinear solve in a sequence of +!! solves. This allows the object to be reused and eliminates the overhead +!! of repeated memory allocation and deallocation that would otherwise occur. +!! +!! RELAX() deletes the pending vectors that were cached by the preceding +!! call to ACCEL_UPDATE, if any. This modifies the behavior of the next +!! call to ACCEL_UPDATE in that the acceleration subspace will not be +!! updated prior to computing the accelerated update. This could be used, +!! for example, to carry over the subspace from one nonlinear solve to +!! another. (Whether this is an effective strategy is an open question.) +!! ACCEL_UPDATE expects that the passed function value is connected to the +!! preceding update (if it exists), but this is not normally true for +!! the first call in a subsequent nonlinear solve, and would result in the +!! subspace being updated with bogus information. A call to RELAX at the +!! end of a nonlinear solve prevents this from occuring. +!! +!! NUM_VEC() returns the number of vectors in the acceleration subspace. +!! +!! MAX_VEC() returns the max number of vectors in the acceleration subspace. +!! +!! VEC_LEN() returns the length of the vectors. +!! +!! VEC_TOL() returns the vector drop tolerance. +!! +!! REAL_KIND() returns the kind parameter value expected of all real arguments. +!! +!! DEFINED() returns the value true if the object is well-defined; otherwise +!! it returns the value false. Defined means that the data components of +!! the object are properly and consistently defined. Due to the significant +!! effort this function goes through to examine the object, it is primarily +!! intended to be used in debugging situations. Note that this function is +!! used internally when this module is compiled without the preprocessor +!! -DNDEBUG flag. +!! +!! USAGE EXAMPLE +!! +!! The following simple example shows the usage of this acceleration +!! procedure. For more details, see the associated documentation. +!! Consider a quasi-Newton iteration for solving the nonlinear system +!! F(x) = 0. Suppose PC(y) is some preconditioning procedure that applies +!! some approximation to the inverse of the Jacobian of F(x) to the vector y. +!! The original quasi-Newton iteration (equivalent to the fixed point +!! iteration for PC(F(x)) = 0) would look something like +!! +!! x = 0 +!! do <until converged> +!! dx = PC(F(x)) +!! x = x - dx +!! end do +!! +!! The accelerated iteration would look something like +!! +!! type(nka) :: accel +!! call accel%init (size(v), mvec=5) +!! x = 0 +!! do <until converged> +!! dx = PC(F(x)) +!! call accel%accel_update(dx) +!! x = x - dx +!! end do +!! +!! The INIT call can be moved outside any nonlinear solution procedure +!! containing this iteration, and a single NKA-type variable used for repeated +!! calls to the procedure. This avoids the repeated allocation and deallocation +!! of memory associated with the accelerator. In this case, one should either +!! include a call to RESTART before the loop so that each iterative solve starts +!! with clean slate, or include a call to RELAX after the loop so that first +!! call to ACCEL_UPDATE in the next iterative solve doesn't update the +!! acceleration subspace with bogus information. +!! + +#include "f90_assert.fpp" + +module new_nka_type + + use,intrinsic :: iso_fortran_env, only: r8 => real64 + use vector_class + implicit none + private + + type, public :: nka + private + logical :: subspace = .false. + logical :: pending = .false. + integer :: mvec = 0 ! maximum number of vectors + real(r8) :: vtol = 0.01_r8 ! vector drop tolerance + !! Subspace storage. + class(vector), allocatable :: v(:) ! update vectors + class(vector), allocatable :: w(:) ! function difference vectors + real(r8), allocatable :: h(:,:) ! matrix of inner products + !! Linked-list organization of the vector storage. + integer :: first, last, free + integer, allocatable :: next(:), prev(:) + contains + procedure :: init => nka_init + procedure :: set_vec_tol => nka_set_vec_tol + procedure :: num_vec => nka_num_vec + procedure :: max_vec => nka_max_vec + procedure :: vec_tol => nka_vec_tol + procedure :: real_kind => nka_real_kind + procedure :: accel_update => nka_accel_update + procedure :: relax => nka_relax + procedure :: restart => nka_restart + procedure :: defined => nka_defined + end type nka + +contains + + subroutine nka_init (this, vec, mvec) + class(nka), intent(out) :: this + class(vector), intent(in) :: vec + integer, intent(in) :: mvec + integer :: n + ASSERT(mvec > 0) + this%mvec = mvec + n = mvec + 1 + call vec%clone(this%v, n) + call vec%clone(this%w, n) + allocate(this%h(n,n), this%next(n), this%prev(n)) + call nka_restart (this) + ASSERT(nka_defined(this)) + end subroutine nka_init + + subroutine nka_set_vec_tol (this, vtol) + class(nka), intent(inout) :: this + real(r8), intent(in) :: vtol + ASSERT(vtol > 0.0_r8) + this%vtol = vtol + end subroutine nka_set_vec_tol + + integer function nka_num_vec (this) + class(nka), intent(in) :: this + integer :: k + nka_num_vec = 0 + k = this%first + do while (k /= 0) + nka_num_vec = nka_num_vec + 1 + k = this%next(k) + end do + if (this%pending) nka_num_vec = nka_num_vec - 1 + end function nka_num_vec + + integer function nka_max_vec (this) + class(nka), intent(in) :: this + nka_max_vec = this%mvec + end function nka_max_vec + + real(r8) function nka_vec_tol (this) + class(nka), intent(in) :: this + nka_vec_tol = this%vtol + end function nka_vec_tol + + integer function nka_real_kind (this) + class(nka), intent(in) :: this + nka_real_kind = kind(this%h) + end function nka_real_kind + + + subroutine nka_accel_update (this, f) + + class(nka), intent(inout) :: this + class(vector), intent(inout) :: f + + ! local variables. + integer :: i, j, k, new, nvec + real(r8) :: s, hkk, hkj, cj, c(this%mvec+1) + + ASSERT(nka_defined(this)) + + !!! + !!! UPDATE THE ACCELERATION SUBSPACE + + if (this%pending) then + + !! Next function difference w_1. + !this%w(:,this%first) = this%w(:,this%first) - f + !s = sqrt(this%dp(this%w(:,this%first), this%w(:,this%first))) + call this%w(this%first)%update(-1.0_r8, f) + s = this%w(this%first)%norm2() + + !! If the function difference is 0, we can't update the subspace with + !! this data; so we toss it out and continue. In this case it is likely + !! that the outer iterative solution procedure has gone badly awry + !! (unless the function value is itself 0), and we merely want to do + !! something reasonable here and hope that situation is detected on the + !! outside. + if (s == 0.0_r8) call nka_relax (this) + + end if + + if (this%pending) then + + !! Normalize w_1 and apply same factor to v_1. + !this%v(:,this%first) = this%v(:,this%first) / s + !this%w(:,this%first) = this%w(:,this%first) / s + call this%v(this%first)%scale(1.0_r8/s) + call this%w(this%first)%scale(1.0_r8/s) + + !! Update H. + k = this%next(this%first) + do while (k /= 0) + !this%h(this%first,k) = this%dp(this%w(:,this%first), this%w(:,k)) + this%h(this%first,k) = this%w(this%first)%dot(this%w(k)) + k = this%next(k) + end do + + !!! + !!! CHOLESKI FACTORIZATION OF H + + this%h(this%first,this%first) = 1.0_r8 + k = this%next(this%first) + nvec = 1 + + do while (k /= 0) + nvec = nvec + 1 + if (nvec > this%mvec) then ! Maintain at most MVEC vectors: + !! Drop the last vector and update the free storage list. + ASSERT(this%last == k) + this%next(this%last) = this%free + this%free = k + this%last = this%prev(k) + this%next(this%last) = 0 + exit + end if + + hkk = 1.0_r8 ! Single stage of Choleski factorization. + j = this%first ! Original matrix kept in lower triangle (unit diagonal). + do while (j /= k) ! Upper triangle holds the factorization. + hkj = this%h(j,k) + i = this%first + do while (i /= j) + hkj = hkj - this%h(k,i) * this%h(j,i) + i = this%next(i) + end do + hkj = hkj / this%h(j,j) + hkk = hkk - hkj**2 + this%h(k,j) = hkj + j = this%next(j) + end do + + if (hkk > this%vtol**2) then + this%h(k,k) = sqrt(hkk) + else ! The current w nearly lies in the span of the previous vectors. + + !! Drop this vector + ASSERT(this%prev(k) /= 0) + this%next(this%prev(k)) = this%next(k) + if (this%next(k) == 0) then + this%last = this%prev(k) + else + this%prev(this%next(k)) = this%prev(k) + end if + + this%next(k) = this%free ! update the free storage list, + this%free = k + + k = this%prev(k) ! and back-up. + nvec = nvec - 1 + + end if + k = this%next(k) + end do + + ASSERT(this%first /= 0) + this%subspace = .true. + this%pending = .false. + + end if + + !! Locate storage for the new vectors. + ASSERT(this%free /= 0) + new = this%free + this%free = this%next(this%free) + + !! Save the original f for the next call. + !this%w(:,new) = f + call this%w(new)%copy(f) + + !!! + !!! ACCELERATED UPDATE + + if (this%subspace) then + + !! Project f onto the span of the w vectors: forward substitution + j = this%first + do while (j /= 0) + !cj = this%dp(f, this%w(:,j)) + cj = f%dot(this%w(j)) + i = this%first + do while (i /= j) + cj = cj - this%h(j,i) * c(i) + i = this%next(i) + end do + c(j) = cj / this%h(j,j) + j = this%next(j) + end do + + !! Project f onto the span of the w vectors: backward substitution + j = this%last + do while (j /= 0) + cj = c(j) + i = this%last + do while (i /= j) + cj = cj - this%h(i,j) * c(i) + i = this%prev(i) + end do + c(j) = cj / this%h(j,j) + j = this%prev(j) + end do + + !! The accelerated update + k = this%first + do while (k /= 0) + !f = f - c(k) * this%w(:,k) + c(k) * this%v(:,k) + call f%update(-c(k), this%w(k), c(k), this%v(k)) + k = this%next(k) + end do + + end if + + !! Save the update for the next call. + !this%v(:,new) = f + call this%v(new)%copy(f) + + !! Prepend the new vectors to the list. + this%prev(new) = 0 + this%next(new) = this%first + if (this%first == 0) then + this%last = new + else + this%prev(this%first) = new + end if + this%first = new + + !! The original f and accelerated update are cached for the next call. + this%pending = .true. + + end subroutine nka_accel_update + + + subroutine nka_restart (this) + class(nka), intent(inout) :: this + integer :: k + this%subspace = .false. + this%pending = .false. + !! No vectors are stored. + this%first = 0 + this%last = 0 + !! Initialize the free storage linked list. + this%free = 1 + do k = 1, size(this%next)-1 + this%next(k) = k + 1 + end do + this%next(size(this%next)) = 0 + end subroutine nka_restart + + + subroutine nka_relax (this) + class(nka), intent(inout) :: this + integer :: new + if (this%pending) then + ASSERT(this%first /= 0) + !! Drop the pending vectors. + new = this%first + this%first = this%next(this%first) + if (this%first == 0) then + this%last = 0 + else + this%prev(this%first) = 0 + end if + !! Update the free storage list. + this%next(new) = this%free + this%free = new + this%pending = .false. + end if + end subroutine nka_relax + + + logical function nka_defined (this) + + class(nka), intent(in) :: this + + integer :: n + logical, allocatable :: tag(:) + + CHECKLIST: do + nka_defined = .false. + if (this%mvec < 1) exit + if (.not.allocated(this%v)) exit + if (.not.allocated(this%w)) exit + if (size(this%v) /= this%mvec+1) exit + if (size(this%w) /= this%mvec+1) exit + if (.not.allocated(this%h)) exit + if (size(this%h,dim=1) /= this%mvec+1) exit + if (size(this%h,dim=2) /= this%mvec+1) exit + if (.not.allocated(this%next)) exit + if (size(this%next) /= this%mvec+1) exit + if (.not.allocated(this%prev)) exit + if (size(this%prev) /= this%mvec+1) exit + + if (this%vtol <= 0.0_r8) exit + + n = size(this%next) + if (any(this%next < 0) .or. any(this%next > n)) exit + if (this%first < 0 .or. this%first > n) exit + if (this%free < 0 .or. this%free > n) exit + + !! Tag array: each location is either in the free list or vector list. + allocate(tag(size(this%next))) + tag = .false. + + !! Check the vector list for consistency. + if (this%first == 0) then + if (this%last /= 0) exit + else + n = this%first + if (this%prev(n) /= 0) exit + tag(n) = .true. + do while (this%next(n) /= 0) + if (this%prev(this%next(n)) /= n) exit CHECKLIST + n = this%next(n) + if (tag(n)) exit CHECKLIST + tag(n) = .true. + end do + if (this%last /= n) exit + end if + + !! Check the free list. + n = this%free + do while (n /= 0) + if (tag(n)) exit CHECKLIST + tag(n) = .true. + n = this%next(n) + end do + + !! All locations accounted for? + if (.not.all(tag)) exit + + nka_defined = .true. + exit + end do CHECKLIST + + if (allocated(tag)) deallocate(tag) + + end function nka_defined + +end module new_nka_type diff --git a/src/truchas/ode/new_state_history_type.F90 b/src/truchas/ode/new_state_history_type.F90 new file mode 100644 index 000000000..9a37266cb --- /dev/null +++ b/src/truchas/ode/new_state_history_type.F90 @@ -0,0 +1,433 @@ +!! +!! STATE_HISTORY_TYPE +!! +!! This module provides a structure for maintaining the recent history of an +!! a solution procedure that is characterized by a (time) sequence of state +!! vectors, and methods for performing polynomial interpolation based on that +!! history. +!! +!! Neil N. Carlson <neil.n.carlson@gmail.com> 7 Jul 2003 +!! Revision for Fortran 2003, 2 Apr 2011, 27 Dec 2011 +!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! +!! Copyright (c) 2003, 2011, 2013 Neil N. Carlson +!! +!! Permission is hereby granted, free of charge, to any person obtaining a +!! copy of this software and associated documentation files (the "Software"), +!! to deal in the Software without restriction, including without limitation +!! the rights to use, copy, modify, merge, publish, distribute, sublicense, +!! and/or sell copies of the Software, and to permit persons to whom the +!! Software is furnished to do so, subject to the following conditions: +!! +!! The above copyright notice and this permission notice shall be included +!! in all copies or substantial portions of the Software. +!! +!! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +!! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +!! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +!! THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +!! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +!! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +!! DEALINGS IN THE SOFTWARE. +!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! +!! PROGRAMMING INTERFACE +!! +!! This module defines the derived data type STATE_HISTORY (with private +!! components) with the following methods as type-bound procedures. +!! +!! INIT (MVEC, {T, X [, XDOT] | VLEN}) initializes the history structure +!! to maintain up to MVEC state vectors. In the first variant, the +!! vector X with time index T is recorded as the initial vector of a new +!! history. If the optional vector XDOT is also specified, it is recorded +!! as the state vector time derivative at the same time index T. In the +!! second variant, VLEN specifies the length of the vectors to be maintained +!! but no state vector is recorded. It is permissible to re-initialize a +!! history structure that had been previously initialized. +!! +!! FLUSH (T, X [, XDOT]) flushes the accumulated state vectors and records the +!! state vector X with time index T as the initial state vector of a new +!! history. If XDOT is specified, it is also recorded as the state vector +!! time derivative at time index T. This differs from re-initializing in +!! that the maximum number of vectors and their length are not changed. +!! +!! REAL_KIND() returns the kind parameter value expected of all real arguments. +!! +!! DEPTH() returns the number of state vectors currently stored. This number +!! will vary between 0 and the value of MVEC used to initialize the object. +!! +!! MAX_DEPTH() returns the maximum number state vectors that can be stored. +!! This number is the value of MVEC used to initialize the object. +!! +!! STATE_SIZE() returns the length of the state vectors stored by the object. +!! This number will equal the size of the state vector or the value of VLEN +!! that was used to initialize the object. +!! +!! LAST_TIME() returns the time index of the last recorded state vector. +!! +!! GET_LAST_STATE_COPY (COPY) copies the last recorded state vector into +!! the array COPY, whose length should equal STATE_SIZE(). +!! +!! GET_LAST_STATE_VIEW (VIEW) associates the pointer VIEW with the last +!! recorded state vector. THIS SHOULD BE USED WITH GREAT CAUTION. The +!! target of this pointer should never be modified or deallocated. The +!! pointer will cease to reference the last recorded state vector when +!! the history is subsequently modified through calls to RECORD_STATE, +!! FLUSH or INIT. +!! +!! TIME_DELTAS() returns an array containing the time index differences: the +!! first element is the difference between the last and penultimate times; +!! the second element is the difference between the last and antepenultimate +!! times, and so forth. The length of the result equals DEPTH()-1. It is +!! an error to call this method if DEPTH() is less than 2. +!! +!! RECORD_STATE (T, X [, XDOT]) records the vector X with time index T as +!! the last state vector in the history. If the vector XDOT is present, +!! it is recorded as the state vector time derivative with the same time +!! index. The oldest state vector (or the two oldest in the case XDOT +!! is present) is discarded once the history is fully populated with MVEC +!! vectors. Note that when only one of a X/XDOT pair of vectors is +!! discarded, it is effectively the derivative vector that gets discarded. +!! +!! INTERP_STATE (T, X [, FIRST] [, ORDER]) computes the interpolated state +!! vector at time T from the set of vectors maintained by the history, and +!! returns the result in the user-supplied array X. Polynomial interpolation +!! is used, and ORDER, if present, specifies the order using the ORDER+1 +!! most recent vectors; 1 for linear interpolation, 2 for quadratic, etc. +!! It is an error to request an order for which there is insufficient data. +!! If not specified, the maximal interpolation order is used given the +!! available data; once the history is fully populated, the interpolation +!! order is MVEC-1. Typically the array X would have the same size as the +!! stored state vectors, but more generally X may return any contiguous +!! segment of the interpolated state starting at index FIRST (default 1) +!! and length the size of X. +!! +!! DEFINED() returns true if the object is well-defined; otherwise it returns +!! false. Defined means that the data components of the object are properly +!! and consistently defined. DEFINED() should return true at any time after +!! the INIT method has been called. This method is primarily intended to be +!! used in debugging situations and is used internally in assertion checks. +!! +!! REMARKS +!! +!! Oldest/newest refers to the sequence in which vectors are recorded and not +!! to the time values supplied. Typically, the corresponding time sequence +!! will be strictly increasing; anything else would be perverse. +!! +!! When recording the time derivative (XDOT) of a solution vector, this vector +!! counts as an additional vector when considered by DEPTH() and MAX_DEPTH(). +!! +!! IMPLEMENTATION NOTES +!! +!! The sequence of solution vectors is actually stored as a table of +!! divided differences, which is easily updated and which makes polynomial +!! interpolation particularly simple to express. +!! +!! History shifting (induced by recording a new vector) takes place through +!! shifting pointers; vectors themselves are not copied. The natural +!! implementation would use an array of pointers to vectors, but this is not +!! possible to do directly. The indirect means is to box up the pointer in +!! an auxillary type (PTR_BOX) and declare an array of this type. +!! + +#include "f90_assert.fpp" + +module new_state_history_type + + use,intrinsic :: iso_fortran_env, only: r8 => real64 + use vector_class + implicit none + private + + type :: vec_box + class(vector), allocatable :: vec + end type vec_box + + type, public :: state_history + private + integer :: nvec = 0 ! number of vectors + integer :: mvec = 0 ! maximum number of vectors + real(r8), allocatable :: t(:) ! vector times + type(vec_box), allocatable :: d(:) ! divided differences + contains + procedure, private :: init_vector, init_state, init_copy + generic :: init => init_vector, init_state, init_copy + procedure :: defined + procedure :: flush + procedure :: real_kind + procedure :: depth + procedure :: max_depth + procedure :: last_time + procedure :: get_last_state_view + procedure :: get_last_state_copy + procedure :: time_deltas + procedure :: record_state + procedure :: interp_state + generic :: assignment(=) => copy + procedure, private :: copy + end type state_history + +contains + + subroutine init_vector(this, mvec, vec) + class(state_history), intent(out) :: this + integer, intent(in) :: mvec + class(vector), intent(in) :: vec + integer :: j + ASSERT(mvec > 0) + this%mvec = mvec + this%nvec = 0 + allocate(this%t(mvec), this%d(mvec)) + do j = 1, mvec + call vec%clone(this%d(j)%vec) + end do + end subroutine init_vector + + subroutine init_state(this, mvec, t, x, xdot) + class(state_history), intent(out) :: this + integer, intent(in) :: mvec + real(r8), intent(in) :: t + class(vector), intent(in) :: x + class(vector), intent(in), optional :: xdot + call init_vector(this, mvec, x) + call record_state(this, t, x, xdot) + end subroutine init_state + + subroutine init_copy(this, src) + class(state_history), intent(out) :: this + class(state_history), intent(in) :: src + integer :: j + this%mvec = src%mvec + this%nvec = src%nvec + allocate(this%t(this%mvec), this%d(this%mvec)) + do j = 1, this%mvec + call src%d(1)%vec%clone(this%d(j)%vec) + end do + do j = 1, this%nvec + this%t(j) = src%t(j) + call this%d(j)%vec%copy(src%d(j)%vec) + end do + end subroutine init_copy + + subroutine flush (this, t, x, xdot) + class(state_history), intent(inout) :: this + real(r8), intent(in) :: t + class(vector), intent(in) :: x + class(vector), intent(in), optional :: xdot + this%nvec = 0 + call record_state(this, t, x, xdot) + end subroutine flush + + pure integer function real_kind (this) + class(state_history), intent(in) :: this + real_kind = kind(this%t) + end function real_kind + + integer function depth (this) + class(state_history), intent(in) :: this + depth = this%nvec + end function depth + + integer function max_depth (this) + class(state_history), intent(in) :: this + max_depth = 0 + if (allocated(this%t)) max_depth = size(this%t) + end function max_depth + + subroutine get_last_state_view (this, view) + class(state_history), intent(in), target :: this + class(vector), pointer :: view + ASSERT(this%nvec > 0) + view => this%d(1)%vec + end subroutine get_last_state_view + + subroutine get_last_state_copy (this, copy) + class(state_history), intent(in) :: this + class(vector), intent(inout) :: copy + ASSERT(this%nvec > 0) + !ASSERT(size(copy) == this%vlen) + !copy = this%d(1)%vec + call copy%copy(this%d(1)%vec) + end subroutine get_last_state_copy + + function last_time (this) result (t) + class(state_history), intent(in) :: this + real(r8) :: t + ASSERT(this%nvec > 0) + t = this%t(1) + end function last_time + + function time_deltas (this) result (deltas) + class(state_history), intent(in) :: this + real(r8) :: deltas(this%nvec-1) + ASSERT(this%nvec > 1) + deltas = this%t(1) - this%t(2:this%nvec) + end function time_deltas + + subroutine interp_state (this, t, x, order) + + class(state_history), intent(in) :: this + real(r8), intent(in) :: t + class(vector), intent(inout) :: x + integer, intent(in), optional :: order + + integer :: k, interp_order + + ASSERT(defined(this)) + ASSERT(this%nvec > 0) + + !! Set the interpolation order. + if (present(order)) then + ASSERT(order >= 0 .and. order < this%nvec) + interp_order = order + else ! do maximal order of interpolation with available data. + interp_order = this%nvec - 1 + end if + +! do j = 1, size(x) + !value = this%d(interp_order+1)%vec(j) + call x%copy(this%d(interp_order+1)%vec) + do k = interp_order, 1, -1 + !value = this%d(k)%vec(j) + (t-this%t(k)) * value + call x%update(1.0_r8, this%d(k)%vec, (t-this%t(k))) + end do +! x(j) = value +! end do + + end subroutine interp_state + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !! + !! RECORD_STATE + !! + !! This subroutine records the vector X with time index T as the most recent + !! state vector in the history structure THIS. If the vector XDOT is + !! present, it is recorded as the solution vector time derivative at the same + !! time index. The oldest solution vector (or two oldest in the case XDOT is + !! present) is discarded once the history is fully populated. Note that it + !! is actually the table of vector divided differences that is stored, rather + !! than the vectors themselves. + !! + !! The treatment of the derivative data XDOT can be viewed as recording a + !! second state vector at a slightly greater time index than T such that + !! the first divided difference is XDOT, and then taking the limit as this + !! time index approaches T. This is implemented here by first recording X + !! as a new state vector (and updating all the divided differences) and + !! then introducing X again as a new state vector at the same time index + !! but with the first divided difference set to XDOT, and then computing the + !! the higher order divided differences as usual. + !! + !! Interpolation using the resulting history works as expected: evaluating the + !! interpolant at the vector's time index yields the vector, and evaluating + !! the derivative of the interpolant at that time index would yield XDOT. + !! + !! Note that when only one of a X/XDOT pair of vectors is discarded, it is the + !! derivative vector (which exists in the higher order divided difference) that + !! gets discarded. + !! + + subroutine record_state (this, t, x, xdot) + + class(state_history), intent(inout) :: this + real(r8), intent(in) :: t + class(vector), intent(in) :: x + class(vector), intent(in), optional :: xdot + + integer :: j + type(vec_box) :: tmp + + ASSERT(defined(this)) + + this%nvec = min(1+this%nvec, this%mvec) ! update the number of vectors + + !! Shift the divided differences + !tmp = this%d(this%nvec) ! storage for oldest gets recycled for newest + call move_alloc(this%d(this%nvec)%vec, tmp%vec) + do j = this%nvec, 2, -1 + this%t(j) = this%t(j-1) + !this%d(j) = this%d(j-1) + call move_alloc(this%d(j-1)%vec, this%d(j)%vec) + end do + + !! Insert the new vector + this%t(1) = t + !this%d(1) = tmp + !this%d(1)%vec = x + call move_alloc(tmp%vec, this%d(1)%vec) + call this%d(1)%vec%copy(x) + + !! Update the divided differences + do j = 2, this%nvec + !this%d(j)%vec = (this%d(j-1)%vec - this%d(j)%vec) / (this%t(1) - this%t(j)) + !a = 1.0_r8 / (this%t(1) - this%t(j)) + !call this%d(j)%vec%update(a, this%d(j-1)%vec, -a) + call this%d(j)%vec%update(-1.0_r8, this%d(j-1)%vec) + call this%d(j)%vec%scale(1.0_r8/(this%t(j) - this%t(1))) + end do + + if (.not.present(xdot)) return + if (this%mvec == 1) return ! no room to store derivative info + + this%nvec = min(1+this%nvec, this%mvec) ! update the number of vectors + + !! Shift the divided differences, except the first; the new vector and + !! time index are the same as the most recent. + !tmp = this%d(this%nvec) ! storage for oldest gets recycled for newest + call move_alloc(this%d(this%nvec)%vec, tmp%vec) ! storage for oldest gets recycled for newest + do j = this%nvec, 3, -1 + this%t(j) = this%t(j-1) + !this%d(j) = this%d(j-1) + call move_alloc(this%d(j-1)%vec, this%d(j)%vec) + end do + + !! The first divided difference (same time index) is the specified derivative. + this%t(2) = this%t(1) + !this%d(2) = tmp + !this%d(2)%vec = xdot + call move_alloc(tmp%vec, this%d(2)%vec) + call this%d(2)%vec%copy(xdot) + + !! Update the rest of the divided differences. + do j = 3, this%nvec + !this%d(j)%vec = (this%d(j-1)%vec - this%d(j)%vec) / (this%t(1) - this%t(j)) + call this%d(j)%vec%update(-1.0_r8, this%d(j-1)%vec) + call this%d(j)%vec%scale(1.0_r8/(this%t(j) - this%t(1))) + end do + + end subroutine record_state + + logical function defined (this) + class(state_history), intent(in) :: this + integer :: j + defined = .false. + if (.not.allocated(this%t)) return + if (.not.allocated(this%d)) return + if (size(this%t) /= this%mvec) return + if (size(this%d) /= this%mvec) return + if (this%nvec < 0 .or. this%nvec > this%mvec) return + do j = 1, this%mvec + if (.not.allocated(this%d(j)%vec)) return + end do + defined = .true. + end function defined + + !! Defined assignment subroutine + subroutine copy(lhs, rhs) + class(state_history), intent(inout) :: lhs + class(state_history), intent(in) :: rhs + integer :: j + if (lhs%mvec /= rhs%mvec) then + call init_copy(lhs, rhs) + else + lhs%nvec = rhs%nvec + do j = 1, lhs%nvec + lhs%t(j) = rhs%t(j) + !lhs%d(j)%vec = rhs%d(j)%vec + call lhs%d(j)%vec%copy(rhs%d(j)%vec) + end do + end if + end subroutine copy + +end module new_state_history_type diff --git a/src/truchas/ode/vector_class.F90 b/src/truchas/ode/vector_class.F90 new file mode 100644 index 000000000..f505747ea --- /dev/null +++ b/src/truchas/ode/vector_class.F90 @@ -0,0 +1,277 @@ +!! +!! VECTOR_CLASS +!! +!! Neil N. Carlson <nnc@lanl.gov> +!! March 2020 +!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! +!! This file is part of Truchas. 3-Clause BSD license; see the LICENSE file. +!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! +!! CLASS METHODS +!! +!! CLASS(VECTOR) :: V +!! +!! CALL V%CLONE(C [,N]) creates a clone C of V. The allocatable CLASS(VECTOR) +!! variable C is allocated with the same dynamic type and internal structure +!! as V, but the values of its vector elements are not defined (see COPY). +!! C may be a rank-1 array, in which case N must also be specified, and its +!! value is the size to which C is allocated. +!! +!! This is a generic interface with deferred specific procedures CLONE1 and +!! CLONE2 that subclasses must implement. +!! +!! CALL V%COPY(SRC) copies the vector element values from the CLASS(VECTOR) +!! variable SRC to V. V and SRC must have the same dynamic type and be +!! compatibly defined, as is the case if V were a clone of SRC. +!! +!! This is implemented using the non-virtual interface (NVI) pattern. +!! Subclasses implement the private deferred procedure COPY_ and are +!! assured that SRC and V will have the same dynamic type. +!! +!! CALL V%SETVAL(VAL) sets the vector element values to the real scalar VAL. +!! +!! This is a deferred procedure implemented by subclasses. +!! +!! CALL V%SCALE(A) multiplies the vector elements of V by the real scalar A. +!! +!! This is a deferred procedure implemented by subclasses. +!! +!! CALL V%UPDATE(A, X [,B [,Y [,C]]]) performs SAXPY-like vector updates: +!! V <-- A*X + V +!! V <-- A*X + B*V +!! V <-- A*X + B*Y + V +!! V <-- A*X + B*Y + C*V +!! V, X, and Y must all have the same dynamic type and be compatibly defined. +!! A, B, and C are real scalars. +!! +!! This is a generic procedure with specific procedures implemented using +!! the NVI pattern. Subclasses implement the private deferred procedures +!! UPDATE1_, UPDATE2_, UPDATE3_, UPDATE4_ corresponding to the different +!! versions of the update, and are assured that V, X, and Y will all have +!! the same dynamic type. +!! +!! V%DOT(Y) is the dot product of V with the CLASS(VECTOR) variable Y. V and +!! Y must have the same dynamic type and be compatibly defined. +!! +!! This is a deferred procedure implemented by subclasses. +!! +!! V%NORM2() is the l2 norm of V; i.e., SQRT(V%DOT(V)) +!! +!! This is a deferred procedure implemented by subclasses. +!! +!! V%CHECKSUM([FULL]) is a string that is a checksum of the vector elements +!! of V. The default behavior is to include only the significant elements in +!! the checksum. These are the elements that would be included in a reduction +!! operation like NORM2, for example. Otherwise, when the optional argument +!! FULL is present with value TRUE, the checksum may include additional data, +!! like duplicated ghost values, at the discretion of the implementation. +!! + +module vector_class + + use,intrinsic :: iso_fortran_env, only: r8 => real64 + implicit none + private + + type, abstract, public :: vector + contains + generic :: clone => clone1, clone2 + procedure(clone1), deferred :: clone1 ! private except to subclass + procedure(clone2), deferred :: clone2 ! private except to subclass + procedure :: copy + procedure(setval), deferred :: setval + procedure(scale), deferred :: scale + generic :: update => update1, update2, update3, update4 + procedure, private :: update1, update2, update3, update4 + procedure :: dot + procedure(norm2), deferred :: norm2 + procedure(checksum), deferred :: checksum + ! remaining are all private except to subclass +#if defined(INTEL_SR04329123) || defined(__GFORTRAN__) + procedure(cpy), deferred :: copy_ + procedure(upd1), deferred :: update1_ + procedure(upd2), deferred :: update2_ + procedure(upd3), deferred :: update3_ + procedure(upd4), deferred :: update4_ + procedure(dt), deferred :: dot_ +#else + procedure(copy), deferred :: copy_ + procedure(update1), deferred :: update1_ + procedure(update2), deferred :: update2_ + procedure(update3), deferred :: update3_ + procedure(update4), deferred :: update4_ + procedure(dot), deferred :: dot_ +#endif + end type + + abstract interface + subroutine clone1(this, clone) + import vector + class(vector), intent(in) :: this + class(vector), allocatable, intent(out) :: clone + end subroutine + subroutine clone2(this, clone, n) + import vector + class(vector), intent(in) :: this + class(vector), allocatable, intent(out) :: clone(:) + integer, intent(in) :: n + end subroutine + end interface + + abstract interface + subroutine setval(this, val) + import vector, r8 + class(vector), intent(inout) :: this + real(r8), intent(in) :: val + end subroutine + end interface + + abstract interface + subroutine scale(this, a) + import vector, r8 + class(vector), intent(inout) :: this + real(r8), intent(in) :: a + end subroutine + end interface + + abstract interface + function norm2(this) + import vector, r8 + class(vector), intent(in) :: this + real(r8) :: norm2 + end function + end interface + + abstract interface + function checksum(this, full) result(string) + import vector + class(vector), intent(in) :: this + logical, intent(in), optional :: full ! default is FALSE + character(:), allocatable :: string + end function + end interface + +#if defined(INTEL_SR04329123) || defined(__GFORTRAN__) + abstract interface + subroutine cpy(dest, src) + import vector + class(vector), intent(inout) :: dest + class(vector), intent(in) :: src + end subroutine cpy + subroutine upd1(this, a, x) + import vector, r8 + class(vector), intent(inout) :: this + class(vector), intent(in) :: x + real(r8), intent(in) :: a + end subroutine + subroutine upd2(this, a, x, b) + import vector, r8 + class(vector), intent(inout) :: this + class(vector), intent(in) :: x + real(r8), intent(in) :: a, b + end subroutine + subroutine upd3(this, a, x, b, y) + import vector, r8 + class(vector), intent(inout) :: this + class(vector), intent(in) :: x, y + real(r8), intent(in) :: a, b + end subroutine + subroutine upd4(this, a, x, b, y, c) + import vector, r8 + class(vector), intent(inout) :: this + class(vector), intent(in) :: x, y + real(r8), intent(in) :: a, b, c + end subroutine + function dt(x, y) + import vector, r8 + class(vector), intent(in) :: x, y + real(r8) :: dt + end function + end interface +#endif + +contains + + recursive subroutine copy(dest, src) + class(vector), intent(inout) :: dest + class(vector), intent(in) :: src + if (same_type_as(dest, src)) then + call dest%copy_(src) + else + error stop 'incompatible arguments to VECTOR%COPY' + end if + end subroutine + + recursive function dot(x, y) + class(vector), intent(in) :: x, y + real(r8) :: dot + if (same_type_as(x, y)) then + dot = x%dot_(y) + else + error stop 'incompatible arguments to VECTOR%DOT' + end if + end function + + !! Conventional SAXPY procedure: y <-- y + a*x + recursive subroutine update1(this, a, x) + class(vector), intent(inout) :: this + class(vector), intent(in) :: x + real(r8), intent(in) :: a + if (a == 0.0_r8) return + if (same_type_as(this, x)) then + call this%update1_(a, x) + else + error stop 'incompatible arguments to VECTOR%UPDATE' + end if + end subroutine + + !! SAXPY-like procedure: y <-- a*x + b*y + recursive subroutine update2(this, a, x, b) + class(vector), intent(inout) :: this + class(vector), intent(in) :: x + real(r8), intent(in) :: a, b + if (a == 0.0_r8) then + call this%scale(b) + else if (same_type_as(this, x)) then + call this%update2_(a, x, b) + else + error stop 'incompatible arguments to VECTOR%UPDATE' + end if + end subroutine + + !! SAXPY-like procedure: z <-- a*x + b*y + z + recursive subroutine update3(this, a, x, b, y) + class(vector), intent(inout) :: this + class(vector), intent(in) :: x, y + real(r8), intent(in) :: a, b + if (a == 0.0_r8) then + call update1(this, b, y) + else if (b == 0.0_r8) then + call update1(this, a, x) + else if (same_type_as(this, x) .and. same_type_as(this, y)) then + call this%update3_(a, x, b, y) + else + error stop 'incompatible arguments to VECTOR%UPDATE' + end if + end subroutine + + !! SAXPY-like procedure: z <-- a*x + b*y + c*z + recursive subroutine update4(this, a, x, b, y, c) + class(vector), intent(inout) :: this + class(vector), intent(in) :: x, y + real(r8), intent(in) :: a, b, c + if (a == 0.0_r8) then + call update2(this, b, y, c) + else if (b == 0.0_r8) then + call update2(this, a, x, c) + else if (same_type_as(this, x) .and. same_type_as(this, y)) then + call this%update4_(a, x, b, y, c) + else + error stop 'incompatible arguments to VECTOR%UPDATE' + end if + end subroutine + +end module vector_class diff --git a/src/truchas/physics/heat_species_transport/CMakeLists.txt b/src/truchas/physics/heat_species_transport/CMakeLists.txt index 696470521..6155e5bf7 100644 --- a/src/truchas/physics/heat_species_transport/CMakeLists.txt +++ b/src/truchas/physics/heat_species_transport/CMakeLists.txt @@ -20,6 +20,15 @@ target_sources(truchas PRIVATE HTSD_solver_factory.F90 HTSD_solver_type.F90 HTSD_init_cond_type.F90 + ht_vector_type.F90 + ht_model_type.F90 + ht_model_factory.F90 + ht_norm_type.F90 + ht_precon_type.F90 + ht_idaesol_model_type.F90 + ht_solver_type.F90 + ht_solver_factory.F90 + ht_ic_solver_type.F90 TofH_type.F90 data_layout_type.F90 mfd_disc_type.F90 diff --git a/src/truchas/physics/heat_species_transport/diffusion_solver.F90 b/src/truchas/physics/heat_species_transport/diffusion_solver.F90 index 2d70b35f6..d3b22fdb2 100644 --- a/src/truchas/physics/heat_species_transport/diffusion_solver.F90 +++ b/src/truchas/physics/heat_species_transport/diffusion_solver.F90 @@ -24,9 +24,11 @@ module diffusion_solver use string_utilities, only: i_to_c use source_mesh_function use mesh_interop + use ht_model_type + use ht_solver_type use FHT_model_type use FHT_solver_type - use HTSD_model_type + use HTSD_model_type, only: htsd_model, htsd_model_delete use HTSD_solver_type use truchas_timers use truchas_logging_services @@ -74,6 +76,9 @@ module diffusion_solver !! The special model and solver that works with transient void. type(FHT_model), pointer :: mod2 => null() type(FHT_solver), pointer :: sol2 => null() + !! New prototype HT-only solver + type(ht_model), pointer :: mod3 => null() + type(ht_solver), pointer :: sol3 => null() class(enthalpy_advector), allocatable :: hadv real(r8) :: cutvof type(parameter_list) :: ds_params, bc_params, species_bc_params, thermal_source_params @@ -82,6 +87,7 @@ module diffusion_solver integer, parameter :: SOLVER1 = 1 ! the standard solver integer, parameter :: SOLVER2 = 2 ! special solver that works with transient void + integer, parameter :: SOLVER3 = 3 ! new HT prototype contains @@ -138,6 +144,9 @@ contains t = h + FHT_solver_last_time(this%sol2) call FHT_solver_advance_state(this%sol2, t, errc) hnext = merge(h/2, huge(1.0_r8), (errc /= 0)) + case (SOLVER3) ! new HT prototype + t = h + this%sol3%last_time() + call this%sol3%step(t, hnext, errc) case default INSIST(.false.) end select @@ -218,6 +227,8 @@ contains call HTSD_solver_commit_pending_state(this%sol1) case (SOLVER2) ! HT solver with transient void call FHT_solver_commit_pending_state(this%sol2) + case (SOLVER3) ! new HT prototype + call this%sol3%commit_pending_state case default INSIST(.false.) end select @@ -246,6 +257,10 @@ contains hlast = FHT_solver_last_step_size(this%sol2) call FHT_solver_get_stepping_stats (this%sol2, counters) write(message,2) hlast, counters(2:4) + case (SOLVER3) + hlast = this%sol3%last_step_size() + call this%sol3%get_stepping_stats(counters) + write(message,1) hlast, counters(1:2), counters(4:6) case default INSIST(.false.) end select @@ -274,6 +289,8 @@ contains end do case (SOLVER2) call FHT_solver_get_cell_temp_copy(this%sol2, state(:,0)) + case (SOLVER3) + call this%sol3%get_cell_temp_copy(state(:,0)) case default INSIST(.false.) end select @@ -290,6 +307,8 @@ contains call HTSD_solver_get_cell_temp_copy (this%sol1, array) case (SOLVER2) call FHT_solver_get_cell_temp_copy (this%sol2, array) + case (SOLVER3) + call this%sol3%get_cell_temp_copy(array) case default INSIST(.false.) end select @@ -304,6 +323,8 @@ contains call HTSD_solver_get_cell_heat_copy (this%sol1, array) case (SOLVER2) call FHT_solver_get_cell_heat_copy (this%sol2, array) + case (SOLVER3) + call this%sol3%get_cell_heat_copy(array) case default INSIST(.false.) end select @@ -332,6 +353,8 @@ contains call HTSD_solver_get_cell_temp_grad (this%sol1, array(:,:this%mesh%ncell_onP)) case (SOLVER2) call FHT_solver_get_cell_temp_grad (this%sol2, array(:,:this%mesh%ncell_onP)) + case (SOLVER3) + call this%sol3%get_cell_temp_grad(array(:,:this%mesh%ncell_onP)) case default INSIST(.false.) end select @@ -348,6 +371,8 @@ contains call HTSD_solver_get_cell_temp_copy (this%sol1, array) case (SOLVER2) call FHT_solver_get_cell_temp_copy (this%sol2, array) + case (SOLVER3) + call this%sol3%get_cell_temp_copy(array) case default INSIST(.false.) end select @@ -363,6 +388,8 @@ contains call HTSD_solver_get_face_temp_copy (this%sol1, array) case (SOLVER2) call FHT_solver_get_face_temp_copy (this%sol2, array) + case (SOLVER3) + call this%sol3%get_face_temp_copy(array) case default INSIST(.false.) end select @@ -377,6 +404,8 @@ contains call HTSD_solver_get_face_temp_view (this%sol1, view) case (SOLVER2) call FHT_solver_get_face_temp_view (this%sol2, view) + case (SOLVER3) + call this%sol3%get_face_temp_view(view) case default INSIST(.false.) end select @@ -393,6 +422,8 @@ contains use FHT_solver_factory use HTSD_model_factory use HTSD_solver_factory + use ht_model_factory + use ht_solver_factory use physics_module, only: flow use enthalpy_advector1_type use thermal_bc_factory1_type @@ -458,7 +489,11 @@ contains end if else !! Void (if any) is fixed; use standard solver - this%solver_type = SOLVER1 + if (this%have_species_diffusion) then + this%solver_type = SOLVER1 + else + this%solver_Type = SOLVER3 + end if if (integrator /= 'adaptive-bdf2') then call TLS_fatal ('DS_INIT: diffusion system characteristics are incompatible with STEPPING_METHOD choice.') end if @@ -493,6 +528,18 @@ contains this%sol2 => create_FHT_solver(this%mmf, this%mod2, this%ds_params) call move_alloc(this%hadv, this%sol2%hadv) + case (SOLVER3) + block + type(thermal_bc_factory1) :: tbc_fac + type(thermal_source_factory) :: tsrc_fac + call tbc_fac%init(this%mesh, stefan_boltzmann, absolute_zero, this%bc_params) + call tsrc_fac%init(this%mesh, this%thermal_source_params) + this%mod3 => create_ht_model(tinit, this%disc, this%mmf, tbc_fac, tsrc_fac, stat, errmsg2) + end block + if (stat /= 0) call TLS_fatal('DS_INIT: ' // trim(errmsg2)) + this%ht_source => this%mod3%source + this%sol3 => create_ht_solver(this%mmf, this%mod3, this%ds_params, stat, errmsg2) + if (stat /= 0) call TLS_fatal('DS_INIT: ' // trim(errmsg2)) case default INSIST(.false.) end select @@ -564,6 +611,8 @@ contains call HTSD_solver_set_initial_state (this%sol1, t, dt, temp, conc) case (SOLVER2) call FHT_solver_set_initial_state (this%sol2, t, temp) + case (SOLVER3) + call this%sol3%set_initial_state(t, temp, dt) case default INSIST(.false.) end select @@ -587,6 +636,8 @@ contains call HTSD_solver_restart (this%sol1, dt) case (SOLVER2) ! nothing to do here yet + case (SOLVER3) + call this%sol3%restart(dt) case default INSIST(.false.) end select @@ -747,6 +798,8 @@ contains call this%mod1%update_moving_vf case (SOLVER2) call this%mod2%update_moving_vf + case (SOLVER3) + call this%mod3%update_moving_vf case default INSIST(.false.) end select @@ -760,6 +813,8 @@ contains call this%mod1%add_moving_vf_events(eventq) case (SOLVER2) call this%mod2%add_moving_vf_events(eventq) + case (SOLVER3) + call this%mod3%add_moving_vf_events(eventq) case default INSIST(.false.) end select diff --git a/src/truchas/physics/heat_species_transport/ht_ic_solver_type.F90 b/src/truchas/physics/heat_species_transport/ht_ic_solver_type.F90 new file mode 100644 index 000000000..8d15b9b1d --- /dev/null +++ b/src/truchas/physics/heat_species_transport/ht_ic_solver_type.F90 @@ -0,0 +1,473 @@ +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! +!! This file is part of Truchas. 3-Clause BSD license; see the LICENSE file. +!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +#include "f90_assert.fpp" + +module ht_ic_solver_type + + use,intrinsic :: iso_fortran_env, only: r8 => real64 + use unstr_mesh_type + use index_partitioning + use ht_vector_type + use ht_model_type + use truchas_logging_services + use parallel_communication + use parameter_list_type + implicit none + private + + type, public :: ht_ic_solver + private + type(unstr_mesh), pointer :: mesh => null() ! reference only -- do not own + type(ht_model), pointer :: model => null() ! reference only -- do not own + type(parameter_list), pointer :: params => null() ! reference only -- do not own + contains + procedure :: init + procedure :: compute + procedure :: compute_udot + end type ht_ic_solver + +contains + + subroutine init (this, model, params) + + class(ht_ic_solver), intent(out) :: this + type(ht_model), intent(in), target :: model + type(parameter_list), pointer :: params + + this%model => model + this%mesh => model%mesh + this%params => params + + end subroutine init + + + subroutine compute(this, t, temp, u, udot) + + class(ht_ic_solver), intent(inout) :: this + real(r8), intent(in) :: t + real(r8), intent(in) :: temp(:) + type(ht_vector), intent(inout) :: u, udot ! data is intent(out) + target :: u + + real(r8), pointer :: state(:,:) + + call TLS_info('') + call TLS_info('Computing consistent initial state for HT solver ...') + + call u%setval(0.0_r8) + + !! Set cell temperatures (U%TC) + ASSERT(size(temp) >= this%mesh%ncell_onP) + u%tc(:this%mesh%ncell_onP) = temp(:this%mesh%ncell_onP) + if (associated(this%model%void_cell)) & + where (this%model%void_cell) u%tc = this%model%void_temp + call gather_boundary(this%mesh%cell_ip, u%tc) + + !! Solve for consistent face temperatures and radiosities. + !! Averaging adjacent cell temps provides a cheap initial guess. + call average_to_faces(this%mesh, u%tc, u%tf, this%model%void_cell) + if (allocated(this%model%bc_dir)) then + call this%model%bc_dir%compute(t) + u%tf(this%model%bc_dir%index) = this%model%bc_dir%value + end if + if (associated(this%model%void_face)) & + where (this%model%void_face) u%tf = this%model%void_temp + state(1:this%mesh%ncell,1:1) => u%tc + call compute_face_temp(this%model, t, state, u, this%params) + + !! Compute the cell enthalpy density (U%HC) + call this%model%H_of_T%compute_value(state, u%hc) + if (associated(this%model%void_cell)) & + where (this%model%void_cell) u%hc = 0.0_r8 + + call compute_udot(this, t, u, udot) + + end subroutine compute + + + subroutine compute_udot (this, t, u, udot) + + class(ht_ic_solver), intent(inout) :: this + real(r8), intent(in) :: t + type(ht_vector), intent(inout) :: u ! data is intent(in) + type(ht_vector), intent(inout) :: udot ! data is intent(out) + target :: u, udot + + integer :: j + real(r8), pointer :: state(:,:) + real(r8) :: dt, Tmin, Tmax + + type(ht_vector), target :: f + + call TLS_info('') + call TLS_info('Computing consistent initial state derivative for HT solver ...') + + !TODO: The existing prop_mesh_func%compute_value function expects a rank-2 + ! state array. This is a workaround until prop_mesh_func is redesigned. + state(1:this%mesh%ncell,1:1) => u%tc + + call udot%setval(0.0_r8) + call f%init(u) + call this%model%compute_f(t, u, udot, f) ! only care about f%tc + + !! Compute Hdot + udot%hc = -f%tc/this%mesh%volume + if (associated(this%model%void_cell)) & + where (this%model%void_cell) udot%hc = 0.0_r8 + + !! Forward Euler time step (only cell enthalpy is advanced) + call this%params%get('dt', dt) + call udot%update(1.0_r8, u, dt) ! udot = u + dt*udot + + !! Compute the advanced cell temps from the advanced cell enthalpies + if (associated(this%model%void_cell)) then + do j = 1, size(u%tc) + if (this%model%void_cell(j)) cycle + Tmin = u%tc(j) - 1 + Tmax = u%tc(j) + 1 + call this%model%T_of_H%compute(j, udot%hc(j), Tmin, Tmax, udot%tc(j)) + end do + else + do j = 1, size(u%tc) + Tmin = u%tc(j) - 1 + Tmax = u%tc(j) + 1 + call this%model%T_of_H%compute(j, udot%hc(j), Tmin, Tmax, udot%tc(j)) + end do + end if + state(1:this%mesh%ncell,1:1) => udot%tc ! at advanced solution + + !! Set consistent advanced face temps and radiosities. Use their initial + !! conditions as the initial guess for the solution procedure. + call compute_face_temp(this%model, t+dt, state, udot, this%params) + + !! Forward Euler approximation to the time derivative at T. + !f = (f - u) / dt + call udot%update(-1.0_r8, u) + call udot%scale(1.0_r8/dt) + +!TODO !! Finite difference approx of face temp and radiosity derivatives. +!TODO if (associated(this%model%ht)) then +!TODO call HTSD_model_get_face_temp_view (this%model, f, var) +!TODO call HTSD_model_set_face_temp (this%model, var, udot) +!TODO if (associated(this%model%vf_rad_prob)) then +!TODO do n = 1, size(this%model%vf_rad_prob) +!TODO call HTSD_model_get_radiosity_view (this%model, n, f, var) +!TODO call HTSD_model_set_radiosity (this%model, n, var, udot) +!TODO end do +!TODO end if +!TODO end if + + end subroutine compute_udot + + !! Given a cell-based field, this auxiliary subroutine produces a face-based + !! field whose value on each face is the average of the values on adjacent + !! cells. SKIP is an optional cell-based mask array identifying cells whose + !! values are to be ignored. Faces with no adjacent cells have the value 0. + + subroutine average_to_faces(mesh, ucell, uface, skip) + + type(unstr_mesh), intent(in) :: mesh + real(r8), intent(in) :: ucell(:) + real(r8), intent(out) :: uface(:) + logical, pointer, intent(in) :: skip(:) + + integer :: j + integer :: scale(size(uface)) + + ASSERT(size(ucell) == mesh%ncell) + ASSERT(size(uface) == mesh%nface) + + uface = 0.0_r8 + scale = 0.0_r8 + do j = 1, mesh%ncell + if (associated(skip)) then + if (skip(j)) cycle + end if + associate (cface => mesh%cface(mesh%xcface(j):mesh%xcface(j+1)-1)) + uface(cface) = uface(cface) + ucell(j) + scale(cface) = scale(cface) + 1 + end associate + end do + call gather_boundary(mesh%face_ip, uface) + call gather_boundary(mesh%face_ip, scale) + + where (scale == 0) + uface = 0.0_r8 + elsewhere + uface = uface / scale + end where + + end subroutine average_to_faces + + !! This auxillary procedure solves for the consistent face temperatures given + !! cell temperatures. It also computes the consistent radiosities if enclosure + !! radiation problems are present. + + subroutine compute_face_temp(this, t, state, u, params) + + use mfd_diff_matrix_type + use pcsr_precon_class + use pcsr_precon_factory + use parameter_list_type + use nka_type + use parallel_communication, only: global_sum, global_maxval, global_all + + type(ht_model), intent(inout) :: this + real(r8), intent(in) :: t, state(:,:) ! lower bound is now 1 + type(ht_vector), intent(inout) :: u + type(parameter_list) :: params + + integer :: n, j, stat, num_itr, max_iter, iter + type(mfd_diff_matrix), target :: matrix + real(r8) :: atol, rtol, error, r0_err, r_err, dT_max + real(r8), allocatable :: z(:) + type(ht_vector) :: udot, f + type(nka) :: nka_obj + procedure(pardp), pointer :: dp + character(80) :: string + character(:), allocatable :: errmsg + class(pcsr_precon), allocatable :: precon + + call TLS_info (' Computing consistent face temperatures and radiosities ...') + +!TODO !! Solve for the radiosity components given the (approx) face temperatures. +!TODO if (associated(this%vf_rad_prob)) then +!TODO do n = 1, size(this%vf_rad_prob) +!TODO call HTSD_model_get_radiosity_view (this, n, u, var) +!TODO var = 0.0_r8 +!TODO associate (faces => this%vf_rad_prob(n)%faces) +!TODO call this%vf_rad_prob(n)%solve_radiosity (t, u%tf(faces), var, stat, num_itr, error) +!TODO end associate +!TODO if (TLS_VERBOSITY >= TLS_VERB_NOISY) then +!TODO write(string,'(4x,a,i0,a,es9.2," (",i0,")")') 'radiosity[', n, ']: |r|/|b|=', error, num_itr +!TODO call TLS_info (string) +!TODO end if +!TODO if (stat /= 0) call TLS_info (' WARNING: radiosities not converged') +!TODO end do +!TODO end if + + !! Initial residual of the face temperature equations. + call udot%init(u) + call udot%setval(0.0_r8) + call f%init(u) + call this%compute_f(t, u, udot, f) + associate (r => f%tf(1:this%mesh%nface_onP)) + r0_err = sqrt(global_sum(norm2(r)**2)) + end associate + if (r0_err == 0.0_r8) return + if (TLS_VERBOSITY >= TLS_VERB_NOISY) then + write(string,'(4x,a,es9.2)') '||Rface(0)||_2 = ', r0_err + call TLS_info (string) + end if + + !! Form the Jacobian (approx) of the face temperature system. + call matrix%init(this%disc) + call make_matrix(matrix) + + call alloc_pcsr_precon(precon, matrix%a22, params, stat, errmsg) + if (stat /= 0) call TLS_fatal('COMPUTE_FACE_TEMP: ' // errmsg) + call precon%compute + + call params%get('atol-temp', atol, default=0.0_r8) + call params%get('rtol-temp', rtol, default=0.0_r8) + call params%get('max-iter', max_iter) + + !! Simple Picard iteration. + allocate(z(this%mesh%nface_onP)) + call nka_obj%init(size(z), 10) + dp => pardp + call nka_obj%set_dot_prod(dp) + + associate (Tface => u%tf(1:this%mesh%nface_onP), & + r => f%tf(1:this%mesh%nface_onP)) + iter = 0 + do ! until converged + + iter = iter + 1 + if (iter > max_iter) exit + + z = r + call precon%apply(z) + call nka_obj%accel_update(z) ! nonlinear krylov acceleration (NKA) + + Tface = Tface - z + + dT_max = global_maxval(abs(z)) + +!TODO !! Solve the radiosity components given the new face temperatures. +!TODO if (associated(this%vf_rad_prob)) then +!TODO do n = 1, size(this%vf_rad_prob) +!TODO call HTSD_model_get_radiosity_view(this, n, u, var) +!TODO associate (faces => this%vf_rad_prob(n)%faces) +!TODO call this%vf_rad_prob(n)%solve_radiosity(t, Tface(faces), var, stat, num_itr, error) +!TODO end associate +!TODO if (TLS_VERBOSITY >= TLS_VERB_NOISY) then +!TODO write(string,'(4x,a,i0,a,es9.2," (",i0,")")') 'radiosity[', n, ']: |r|/|b|=', error, num_itr +!TODO call TLS_info (string) +!TODO end if +!TODO if (stat /= 0) call TLS_info (' WARNING: radiosities not converged') +!TODO end do +!TODO end if + + !! Recompute the face temp residual. + call this%compute_f(t, u, udot, f) + r_err = sqrt(global_dot_product(r,r)) + if (TLS_VERBOSITY >= TLS_VERB_NOISY) then + write(string,'(4x,a,i0,2(a,es9.2))') & + '||Rface(', iter, ')||=', r_err, ', ||ΔTface||_max=', dT_max + call TLS_info (string) + end if + + !! Check for convergence. + if (global_all(abs(z) <= atol + rtol * abs(Tface))) exit + end do + end associate + + write(string,'(4x,a,i0,3(a,es9.2))') & + '||Rface(', iter, ')||=', r_err, ', ||Rface(0)||=', r0_err, & + ', ||ΔTface||_max=', dT_max + call TLS_info (string) + + if (iter > max_iter) call TLS_info (' WARNING: face temperatures not converged') + + contains + + !! This auxillary routine computes an approximate Jacobian matrix for the + !! time-independent heat equation (diffusion operator). This is a cell and + !! face system matrix, but we are only interested in the (2,2) block (face- + !! face). That part is approximate in that it ignores the dependence of + !! the radiosity on face temperatures (assuming the radiosities have been + !! solved for). Moreover the Jacobian would be constant except for the T^4 + !! terms present in radiation BC. + !! + !! Note that this is nearly identical to code in HTSD_precon_type:compute + !! that computes the preconditioner matrix. CAN WE AVOID DUPLICATION? + + subroutine make_matrix (matrix) + + type(mfd_diff_matrix), intent(inout) :: matrix + + integer :: j, n, n1, n2, index + integer, allocatable :: more_dir_faces(:) + real(r8) :: D(this%mesh%ncell), term + real(r8), allocatable :: values(:), values2(:,:) + + !! Generate list of void faces; these are treated like Dirichlet BC. + if (associated(this%void_face)) then + n = count(this%void_face) + allocate(more_dir_faces(n)) + n = 0 + do j = 1, this%mesh%nface + if (this%void_face(j)) then + n = n + 1 + more_dir_faces(n) = j + end if + end do + else + allocate(more_dir_faces(0)) + end if + + !! Compute the diffusion coefficient. + call this%conductivity%compute_value(state, D) + if (associated(this%void_cell)) then + where (this%void_cell) D = 0.0_r8 + end if + + !! Compute the basic diffusion matrix + call matrix%compute (D) + + !! Dirichlet boundary condition fixups. + if (allocated(this%bc_dir)) then + call this%bc_dir%compute(t) + call matrix%set_dir_faces(this%bc_dir%index) + end if + + !! External HTC boundary condition contribution. + if (allocated(this%bc_htc)) then + call this%bc_htc%compute(t, u%tf) + associate (index => this%bc_htc%index, & + deriv => this%bc_htc%deriv) + call matrix%incr_face_diag(index, deriv) + end associate + end if + + !! Simple radiation boundary condition contribution. + if (allocated(this%bc_rad)) then + call this%bc_rad%compute(t, u%tf) + associate (index => this%bc_rad%index, & + deriv => this%bc_rad%deriv) + call matrix%incr_face_diag(index, deriv) + end associate + end if + + !! Experimental evaporation heat flux contribution. + if (allocated(this%evap_flux)) then + call this%evap_flux%compute_deriv(t, u%tf) + associate (index => this%evap_flux%index, & + deriv => this%evap_flux%deriv) + call matrix%incr_face_diag(index, this%mesh%area(index)*deriv) + end associate + endif + + !! Internal HTC interface condition contribution. + if (allocated(this%ic_htc)) then + call this%ic_htc%compute(t, u%tf) + associate (index => this%ic_htc%index, & + deriv => this%ic_htc%deriv) + if (associated(this%void_face)) then + do j = 1, size(index,2) !FIXME? Bad form to modify deriv? + if (any(this%void_face(index(:,j)))) deriv(:,j) = 0.0_r8 + end do + end if + call matrix%incr_interface_flux3(index, deriv) !TODO: rename these methods + end associate + end if + + !! Internal gap radiation condition contribution. + if (allocated(this%ic_rad)) then + call this%ic_rad%compute(t, u%tf) + associate (index => this%ic_rad%index, & + deriv => this%ic_rad%deriv) + if (associated(this%void_face)) then + do j = 1, size(index,2) !FIXME? Bad form to modify deriv? + if (any(this%void_face(index(:,j)))) deriv(:,j) = 0.0_r8 + end do + end if + call matrix%incr_interface_flux3(index, deriv) !TODO: rename these methods + end associate + end if + + !! Dirichlet fix-ups for void faces. + call matrix%set_dir_faces (more_dir_faces) + + !! Enclosure radiation contributions to the preconditioner. + !! This captures T^4 emitted heat (simple) but ignores absorbed heat + !! from other faces (complex and non-local). + if (associated(this%vf_rad_prob)) then + do index = 1, size(this%vf_rad_prob) + associate (faces => this%vf_rad_prob(index)%faces) + allocate(values(size(faces))) + call this%vf_rad_prob(index)%rhs_deriv (t, u%tf(faces), values) + where (.not.this%vf_rad_prob(index)%fmask) values = 0 + call matrix%incr_face_diag (faces, this%mesh%area(faces) * values) + deallocate(values) + end associate + end do + end if + + end subroutine make_matrix + + end subroutine compute_face_temp + + function pardp (a, b) result (dp) + use parallel_communication, only: global_dot_product + real(r8), intent(in) :: a(:), b(:) + real(r8) :: dp + dp = global_dot_product(a, b) + end function pardp + +end module ht_ic_solver_type diff --git a/src/truchas/physics/heat_species_transport/ht_idaesol_model_type.F90 b/src/truchas/physics/heat_species_transport/ht_idaesol_model_type.F90 new file mode 100644 index 000000000..0708e5356 --- /dev/null +++ b/src/truchas/physics/heat_species_transport/ht_idaesol_model_type.F90 @@ -0,0 +1,112 @@ +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! +!! This file is part of Truchas. 3-Clause BSD license; see the LICENSE file. +!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +#include "f90_assert.fpp" + +module ht_idaesol_model_type + + use,intrinsic :: iso_fortran_env, only: r8 => real64 + use new_idaesol_type, only: idaesol_model + use vector_class + use ht_vector_type + use ht_model_type + use ht_precon_type + use ht_norm_type + implicit none + private + + type, extends(idaesol_model), public :: ht_idaesol_model + type(ht_model), pointer :: model => null() ! reference only -- not owned + type(ht_precon), pointer :: precon => null() ! reference only -- not owned + type(ht_norm), pointer :: norm => null() ! reference only -- not owned + contains + procedure :: init + !! Deferred procedures from IDAESOL_MODEL + procedure :: alloc_vector + procedure :: compute_f + procedure :: apply_precon + procedure :: compute_precon + procedure :: du_norm + end type + +contains + + subroutine init(this, model, precon, norm) + class(ht_idaesol_model), intent(out) :: this + type(ht_model), intent(in), target :: model + type(ht_precon), intent(in), target :: precon + type(ht_norm), intent(in), target :: norm + this%model => model + this%precon => precon + this%norm => norm + ASSERT(associated(this%model, precon%model)) + end subroutine + + subroutine alloc_vector(this, vec) + class(ht_idaesol_model), intent(in) :: this + class(vector), allocatable, intent(out) :: vec + type(ht_vector), allocatable :: tmp + allocate(tmp) + call tmp%init(this%model%mesh) + call move_alloc(tmp, vec) + end subroutine + + subroutine compute_f(this, t, u, udot, f) + class(ht_idaesol_model) :: this + real(r8), intent(in) :: t + class(vector), intent(inout) :: u, udot ! data is intent(in) + class(vector), intent(inout) :: f ! data is intent(out) + select type (u) + class is (ht_vector) + select type (udot) + class is (ht_vector) + select type (f) + class is (ht_vector) + call this%model%compute_f(t, u, udot, f) + end select + end select + end select + end subroutine + + subroutine apply_precon(this, t, u, f) + class(ht_idaesol_model) :: this + real(r8), intent(in) :: t + class(vector), intent(inout) :: u ! data is intent(in) + class(vector), intent(inout) :: f ! data is intent(inout) + select type (u) + class is (ht_vector) + select type (f) + class is (ht_vector) + call this%precon%apply(t, u, f) + end select + end select + end subroutine + + subroutine compute_precon(this, t, u, dt) + class(ht_idaesol_model) :: this + real(r8), intent(in) :: t, dt + class(vector), intent(inout) :: u + select type (u) + class is (ht_vector) + call this%precon%compute(t, u, dt) + end select + end subroutine + + subroutine du_norm(this, t, u, du, error) + class(ht_idaesol_model) :: this + real(r8), intent(in) :: t + class(vector), intent(in) :: u, du + real(r8), intent(out) :: error + select type (u) + class is (ht_vector) + select type (du) + class is (ht_vector) + call this%norm%compute(t, u, du, error) + end select + end select + end subroutine + +end module ht_idaesol_model_type diff --git a/src/truchas/physics/heat_species_transport/ht_model_factory.F90 b/src/truchas/physics/heat_species_transport/ht_model_factory.F90 new file mode 100644 index 000000000..2827abf31 --- /dev/null +++ b/src/truchas/physics/heat_species_transport/ht_model_factory.F90 @@ -0,0 +1,357 @@ +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! +!! This file is part of Truchas. 3-Clause BSD license; see the LICENSE file. +!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +#include "f90_assert.fpp" + +module ht_model_factory + + use,intrinsic :: iso_fortran_env, only: r8 => real64 + use ht_model_type + use unstr_mesh_type + use mfd_disc_type + use matl_mesh_func_type + use thermal_bc_factory_class + use thermal_source_factory_type + implicit none + private + + public :: create_ht_model + +contains + + function create_ht_model(tinit, disc, mmf, bc_fac, src_fac, stat, errmsg) result(model) + + use diffusion_solver_data, only: void_temperature + use rad_problem_type + + real(r8), intent(in) :: tinit + type(mfd_disc), intent(in), target :: disc + type(matl_mesh_func), intent(in), target :: mmf + class(thermal_bc_factory), intent(inout) :: bc_fac + type(thermal_source_factory), intent(inout) :: src_fac + integer, intent(out) :: stat + character(:), allocatable , intent(out):: errmsg + + character(128) :: errmsg2 + type(ht_model), pointer :: model + + allocate(model) + + call model%init(disc) + + !! Initializes the VF_RAD_PROB components of MODEL. + model%vf_rad_prob => create_vf_rad_prob(tinit, disc%mesh, stat, errmsg2) + if (stat /= 0) then + errmsg = trim(errmsg2) + return + end if + + !! Defines the equation parameter components of MODEL. + call define_system_parameters(disc%mesh, mmf, model, stat, errmsg) + if (stat /= 0) return + + !! Defines the boundary condition components of MODEL. + call define_system_bc(disc%mesh, model, stat, errmsg) + if (stat /= 0) return + + model%void_temp = void_temperature + + contains + + subroutine define_system_parameters (mesh, mmf, model, stat, errmsg) + + use matl_mesh_func_type + use ds_source_input, only: define_external_source + use parallel_communication, only: global_any + use material_model_driver, only: matl_model + use material_utilities + + type(unstr_mesh), intent(in), target :: mesh + type(matl_mesh_func), intent(in), target :: mmf + type(HT_model), intent(inout) :: model + integer, intent(out) :: stat + character(:), allocatable, intent(out) :: errmsg + + !integer, allocatable :: matid(:) + character(128) :: errmsg2 + + !! Retrieve a list of all the material IDs that may be encountered. + !call mmf%get_all_matl(matid, drop_void=.true.) + + !! Enthalpy density. + call required_property_check(matl_model, 'enthalpy', stat, errmsg) + if (stat /= 0) return + call model%H_of_T%init(mmf, 'enthalpy', stat, errmsg2) + if (global_any(stat /= 0)) then + stat = -1 + errmsg = 'unexpected error defining H_of_T: ' // trim(errmsg2) + return + end if + + !! Thermal conductivity. + call required_property_check(matl_model, 'conductivity', stat, errmsg) + if (stat /= 0) return + call model%conductivity%init(mmf, 'conductivity', stat, errmsg2) + !call pmf_set_harmonic_average (model%conductivity) + if (global_any(stat /= 0)) then + stat = -1 + errmsg = 'unexpected error defining conductivity: ' // trim(errmsg2) + return + end if + + !! External heat source. + call define_external_source (mesh, 'temperature', model%source) + + !! Additional heat sources + call src_fac%alloc_source_func1(model%src, stat, errmsg) + if (stat /= 0) return + + stat = 0 + errmsg = '' + + end subroutine define_system_parameters + + subroutine define_system_bc(mesh, model, stat, errmsg) + + use evaporation_namelist, only: evap_params => params + use evap_heat_flux_type + use index_partitioning, only: gather_boundary + use bitfield_type + use parallel_communication, only: global_all, global_any, global_count + use string_utilities, only: i_to_c + use f08_intrinsics, only: findloc + + type(unstr_mesh), intent(in), target :: mesh + type(HT_model), intent(inout) :: model + integer, intent(out) :: stat + character(:), allocatable, intent(out) :: errmsg + + integer :: j, n + logical, allocatable :: mask(:), rmask(:), fmask(:) + integer, allocatable :: setids(:) + character(160) :: string1, string2 + type(bitfield) :: bitmask + type(evap_heat_flux), allocatable :: evap_flux + + allocate(mask(mesh%nface)) + + mask = .false. + + !! Define the internal HTC interface conditions. + call bc_fac%alloc_htc_ic(model%ic_htc, stat, errmsg) + if (stat /= 0) return + if (allocated(model%ic_htc)) then + mask(model%ic_htc%index(1,:)) = .true. + mask(model%ic_htc%index(2,:)) = .true. + call gather_boundary(mesh%face_ip, mask) + end if + + !! Define the gap radiation interface conditions; + !! may be superimposed with HTC conditions. + call bc_fac%alloc_rad_ic(model%ic_rad, stat, errmsg) + if (stat /= 0) return + if (allocated(model%ic_rad)) then + mask(model%ic_rad%index(1,:)) = .true. + mask(model%ic_rad%index(2,:)) = .true. + call gather_boundary(mesh%face_ip, mask) + end if + + !! Flux-type boundary conditions. These may be superimposed. + allocate(fmask(mesh%nface)) + fmask = .false. + + !! Define the simple flux boundary conditions. + call bc_fac%alloc_flux_bc(model%bc_flux, stat, errmsg) + if (stat /= 0) return + if (allocated(model%bc_flux)) then + if (global_any(mask(model%bc_flux%index))) then + stat = -1 + errmsg = 'temperature flux boundary condition overlaps with interface conditions' + return + end if + fmask(model%bc_flux%index) = .true. ! mark the simple flux faces + end if + + !! Define the external HTC boundary conditions. + call bc_fac%alloc_htc_bc(model%bc_htc, stat, errmsg) + if (stat /= 0) return + if (allocated(model%bc_htc)) then + if (global_any(mask(model%bc_htc%index))) then + stat = -1 + errmsg = 'temperature HTC boundary condition overlaps with interface conditions' + return + end if + fmask(model%bc_htc%index) = .true. ! mark the HTC faces + end if + + !! Tag all the enclosure radation (view factor) faces. They were already + !! verified to be boundary faces and non-overlapping with each other. + allocate(rmask(mesh%nface)) + rmask = .false. + if (associated(model%vf_rad_prob)) then + do j = 1, size(model%vf_rad_prob) + rmask(model%vf_rad_prob(j)%faces) = .true. + !fmask(model%vf_rad_prob(j)%faces) = .true. + end do + call gather_boundary (mesh%face_ip, rmask) + !call gather_boundary (mesh%face_ip, fmask) + end if + + !! Define the (simple) radiation boundary conditions. + call bc_fac%alloc_rad_bc(model%bc_rad, stat, errmsg) + if (stat /= 0) return + if (allocated(model%bc_rad)) then + if (global_any(rmask(model%bc_rad%index))) then + stat = -1 + errmsg = 'temperature radiation boundary condition overlaps with enclosure radiation' + return + end if + end if + if (allocated(model%bc_rad)) then + if (global_any(mask(model%bc_rad%index))) then + stat = -1 + errmsg = 'temperature radiation boundary condition overlaps with interface conditions' + return + end if + fmask(model%bc_rad%index) = .true. ! mark the radiation faces + end if + + !! Define the evaporation heat flux boundary condition + if (allocated(evap_params)) then + allocate(evap_flux) + call evap_flux%init(mesh, evap_params, stat, errmsg) + if (stat /= 0) return + if (.not.global_all(fmask(evap_flux%index))) then + stat = -1 + errmsg = 'evaporation heat flux applied to non-flux boundary' + return + end if + call move_alloc(evap_flux, model%evap_flux) + end if + + !! Merge flux mask into main mask. + mask = mask .or. fmask + deallocate(fmask) + + !! Define the Dirichlet boundary conditions. + call bc_fac%alloc_dir_bc(model%bc_dir, stat, errmsg) + if (stat /= 0) return + if (allocated(model%bc_dir)) then + if (global_any(mask(model%bc_dir%index))) then + stat = -1 + errmsg = 'temperature dirichlet boundary condition overlaps with other conditions' + return + end if + mask(model%bc_dir%index) = .true. ! mark the dirichlet faces + end if + + !! Update enclosure radiation face mask for temperature Dirichlet faces. + if (associated(model%vf_rad_prob) .and. allocated(model%bc_dir)) then + do n = 1, size(model%vf_rad_prob) + associate (faces => model%vf_rad_prob(n)%faces, fmask => model%vf_rad_prob(n)%fmask) + do j = 1, size(faces) + if (findloc(model%bc_dir%index, faces(j), dim=1) > 0) fmask(j) = .false. + end do + end associate + end do + end if + mask = mask .or. rmask + + !! Finally verify that a condition has been applied to every boundary face. + mask = mask .neqv. btest(mesh%face_set_mask,0) + if (global_any(mask)) then + call mesh%get_face_set_ids(pack([(j,j=1,mesh%nface)], mask), setids) + if (size(setids) == 0) then + string1 = '(none)' + else + write(string1,'(i0,*(:,", ",i0))') setids + end if + call mesh%get_link_set_ids(mask, setids) + if (size(setids) == 0) then + string2 = '(none)' + else + write(string2,'(i0,*(:,", ",i0))') setids + end if + errmsg = 'incomplete temperature boundary/interface specification;' // & + ' remaining boundary faces belong to face sets ' // trim(string1) // & + '; and interface link sets ' // trim(string2) + bitmask = ibset(ZERO_BITFIELD, 0) + mask = mask .and. (mesh%face_set_mask == bitmask) + mask(mesh%lface(1,:)) = .false. + mask(mesh%lface(2,:)) = .false. + n = global_count(mask(:mesh%nface_onP)) + if (n > 0) errmsg = errmsg // '; ' // i_to_c(n) // ' faces belong to neither' + stat = -1 + return + end if + + stat = 0 + errmsg = '' + + end subroutine define_system_bc + + end function create_HT_model + + function create_vf_rad_prob (tinit, mesh, stat, errmsg) result (vf_rad_prob) + + use rad_problem_type + use enclosure_radiation_namelist, only: params + use bitfield_type, only: btest + use parallel_communication, only: global_any, global_all + use parameter_list_type + + real(r8), intent(in) :: tinit + type(unstr_mesh), intent(in) :: mesh + integer, intent(out) :: stat + character(len=*), intent(out) :: errmsg + type(rad_problem), pointer :: vf_rad_prob(:) + + integer :: n, j + logical, allocatable :: mask(:) + character(len=31), allocatable :: encl_name(:) + type(parameter_list_iterator) :: piter + type(parameter_list), pointer :: plist + + stat = 0 + errmsg = '' + + !! Initialize the enclosure radiation (view factor) problems, if any. + piter = parameter_list_iterator(params, sublists_only=.true.) + n = piter%count() + if (n > 0) then + allocate(vf_rad_prob(n), encl_name(n)) + do j = 1, n + plist => piter%sublist() + call vf_rad_prob(j)%init (mesh, piter%name(), plist, tinit) + !! Verify that these enclosure faces are boundary faces. + if (.not.global_all(btest(mesh%face_set_mask(vf_rad_prob(j)%faces),0))) then + stat = -1 + errmsg = 'some enclosure faces are not boundary faces' + exit + else if (n > 1) then ! multiple enclosures + if (j == 1) then + allocate(mask(mesh%nface)) + mask = .false. + !! Verify that these faces don't overlap with other enclosure faces. + else if (global_any(mask(vf_rad_prob(j)%faces))) then + stat = -1 + errmsg = 'some enclosure faces belong to other enclosures' + exit + !! Tag the faces belonging to this enclosure. + else + mask(vf_rad_prob(j)%faces) = .true. + end if + end if + call piter%next + end do + if (allocated(mask)) deallocate(mask) + else + vf_rad_prob => null() + end if + + end function create_vf_rad_prob + +end module ht_model_factory diff --git a/src/truchas/physics/heat_species_transport/ht_model_type.F90 b/src/truchas/physics/heat_species_transport/ht_model_type.F90 new file mode 100644 index 000000000..489d0d51f --- /dev/null +++ b/src/truchas/physics/heat_species_transport/ht_model_type.F90 @@ -0,0 +1,295 @@ +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! +!! This file is part of Truchas. 3-Clause BSD license; see the LICENSE file. +!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +#include "f90_assert.fpp" + +module ht_model_type + + use,intrinsic :: iso_fortran_env, only: r8 => real64 + use unstr_mesh_type + use ht_vector_type + use mfd_disc_type + use prop_mesh_func_type + use source_mesh_function + use scalar_mesh_func_class + use bndry_func1_class + use bndry_func2_class + use intfc_func2_class + use rad_problem_type + use index_partitioning + use TofH_type + use truchas_timers + use parameter_list_type + implicit none + private + + type, public :: ht_model + type(unstr_mesh), pointer :: mesh => null() + type(mfd_disc), pointer :: disc => null() + logical, pointer :: void_cell(:) => null(), void_face(:) => null() + real(r8) :: void_temp = 0.0_r8 + !! Equation parameters + type(prop_mesh_func) :: conductivity ! thermal conductivity + type(prop_mesh_func) :: H_of_T ! enthalpy as a function of temperature + type(TofH) :: T_of_H + type(source_mf) :: source ! external heat source + class(scalar_mesh_func), allocatable :: src ! another external heat source + !! Boundary condition data + class(bndry_func1), allocatable :: bc_dir ! Dirichlet + class(bndry_func1), allocatable :: bc_flux ! simple flux + class(bndry_func2), allocatable :: bc_htc ! external HTC + class(bndry_func2), allocatable :: bc_rad ! simple radiation + class(intfc_func2), allocatable :: ic_htc ! internal HTC + class(intfc_func2), allocatable :: ic_rad ! internal gap radiation + class(bndry_func2), allocatable :: evap_flux + !! Enclosure radiation problems + type(rad_problem), pointer :: vf_rad_prob(:) => null() + contains + procedure :: init + procedure :: compute_f + procedure :: update_moving_vf + procedure :: add_moving_vf_events + final :: ht_model_delete + end type ht_model + +contains + + subroutine ht_model_delete(this) + type(ht_model), intent(inout) :: this + call smf_destroy(this%source) + if (associated(this%vf_rad_prob)) deallocate(this%vf_rad_prob) + end subroutine + + subroutine init(this, disc) + class(ht_model), intent(out), target :: this + type(mfd_disc), intent(in), target :: disc + this%disc => disc + this%mesh => disc%mesh + block + !TODO: pass parameters + !TODO: redesign the TofH type after redesign of prop_mesh_func type + real(r8) :: eps, delta + integer :: max_try + type(parameter_list) :: params + call params%get('tofh-tol', eps, default=0.0_r8) + call params%get('tofh-max-try', max_try, default=50) + call params%get('tofh-delta', delta, default=1.0_r8) + call this%T_of_H%init (this%H_of_T, eps=eps, max_try=max_try, delta=delta) + end block + end subroutine + + + subroutine compute_f(this, t, u, udot, f) + + use mfd_disc_type + + class(ht_model), intent(inout) :: this + real(r8), intent(in) :: t + type(ht_vector), intent(inout) :: u, udot ! data is intent(in) + type(ht_vector), intent(inout) :: f ! data is intent(out) + target :: u + + integer :: j, n, n1, n2 + real(r8) :: term + real(r8), pointer :: qrad(:) + real(r8), dimension(this%mesh%ncell) :: value + real(r8), allocatable :: Tdir(:), flux(:) + integer, pointer :: faces(:) + logical, allocatable :: void_link(:) + real(r8), pointer :: state(:,:) + + call start_timer('ht-function') + + call gather_boundary(this%mesh%cell_ip, u%tc) + call gather_boundary(this%mesh%face_ip, u%tf) + + !TODO: The existing prop_mesh_func%compute_value function expects a rank-2 + ! state array. This is a workaround until prop_mesh_func is redesigned. + state(1:this%mesh%ncell,1:1) => u%tc + + + !!!! RESIDUAL OF THE ALGEBRAIC ENTHALPY-TEMPERATURE RELATION !!!!!!!!!!!!!!!!! + + call this%H_of_T%compute_value(state, value) + f%hc = u%hc - value + + !! Overwrite function value on void cells with dummy equation H=0. + if (associated(this%void_cell)) where (this%void_cell) f%hc = u%hc + + !!!! RESIDUAL OF THE HEAT EQUATION !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + !! Overwrite the temperature on Dirichlet faces with the boundary + !! data, saving the original values to restore them later. + if (allocated(this%bc_dir)) then + call this%bc_dir%compute(t) + allocate(Tdir(size(this%bc_dir%index))) + do j = 1, size(this%bc_dir%index) + n = this%bc_dir%index(j) + Tdir(j) = u%tf(n) + u%tf(n) = this%bc_dir%value(j) + end do + end if + + !! Compute the generic heat equation residual. + call this%conductivity%compute_value(state, value) + if (associated(this%void_cell)) where (this%void_cell) value = 0.0_r8 + call this%disc%apply_diff(value, u%tc, u%tf, f%tc, f%tf) + call smf_eval(this%source, t, value) + call gather_boundary(this%mesh%cell_ip, udot%hc) + f%tc = f%tc + this%mesh%volume*(udot%hc - value) + + !! Additional heat source + if (allocated(this%src)) then + call this%src%compute(t) + f%tc = f%tc - this%mesh%volume*this%src%value + end if + + !! Overwrite face residuals with the Dirichlet BC residual and + !! restore the face temperatures to their original input values. + if (allocated(this%bc_dir)) then + do j = 1, size(this%bc_dir%index) + n = this%bc_dir%index(j) + f%tf(n) = Tdir(j) - this%bc_dir%value(j) + u%tf(n) = Tdir(j) + end do + deallocate(Tdir) + end if + + !! Simple flux BC contribution. + if (allocated(this%bc_flux)) then + call this%bc_flux%compute(t) + do j = 1, size(this%bc_flux%index) + n = this%bc_flux%index(j) + f%tf(n) = f%tf(n) + this%mesh%area(n) * this%bc_flux%value(j) + end do + end if + + !! External HTC flux contribution. + if (allocated(this%bc_htc)) then + call this%bc_htc%compute(t, u%tf) + do j = 1, size(this%bc_htc%index) + n = this%bc_htc%index(j) + f%tf(n) = f%tf(n) + this%bc_htc%value(j) + end do + end if + + !! Ambient radiation BC flux contribution. + if (allocated(this%bc_rad)) then + call this%bc_rad%compute(t, u%tf) + do j = 1, size(this%bc_rad%index) + n = this%bc_rad%index(j) + f%tf(n) = f%tf(n) + this%bc_rad%value(j) + end do + end if + + !! Experimental evaporation heat flux + if (allocated(this%evap_flux)) then + call this%evap_flux%compute_value(t, u%tf) + do j = 1, size(this%evap_flux%index) + n = this%evap_flux%index(j) + f%tf(n) = f%tf(n) + this%mesh%area(n)*this%evap_flux%value(j) + end do + end if + + !! Internal HTC flux contribution. + if (allocated(this%ic_htc)) then + call this%ic_htc%compute(t, u%tf) + allocate(void_link(size(this%ic_htc%index,2))) + if (associated(this%void_face)) then + do j = 1, size(void_link) + void_link(j) = any(this%void_face(this%ic_htc%index(:,j))) + end do + else + void_link = .false. + end if + do j = 1, size(this%ic_htc%index,2) + if (void_link(j)) cycle + n1 = this%ic_htc%index(1,j) + n2 = this%ic_htc%index(2,j) + f%tf(n1) = f%tf(n1) + this%ic_htc%value(j) + f%tf(n2) = f%tf(n2) - this%ic_htc%value(j) + end do + if (allocated(void_link)) deallocate(void_link) + end if + + !! Internal gap radiation contribution. + if (allocated(this%ic_rad)) then + call this%ic_rad%compute(t, u%tf) + allocate(void_link(size(this%ic_rad%index,2))) + if (associated(this%void_face)) then + do j = 1, size(void_link) + void_link(j) = any(this%void_face(this%ic_rad%index(:,j))) + end do + else + void_link = .false. + end if + do j = 1, size(this%ic_rad%index,2) + if (void_link(j)) cycle + n1 = this%ic_rad%index(1,j) + n2 = this%ic_rad%index(2,j) + f%tf(n1) = f%tf(n1) + this%ic_rad%value(j) + f%tf(n2) = f%tf(n2) - this%ic_rad%value(j) + end do + if (allocated(void_link)) deallocate(void_link) + end if + + !! Overwrite function value on void cells and faces with dummy equation T=0. + if (associated(this%void_cell)) where (this%void_cell) f%tc = u%tc - this%void_temp + if (associated(this%void_face)) where (this%void_face) f%tf = u%tf - this%void_temp + + !!!! RESIDUALS OF THE ENCLOSURE RADIATION SYSTEMS !!!!!!!!!!!!!!!!!!!!!!!!!!!! + +!TODO if (associated(this%vf_rad_prob)) then +!TODO do n = 1, size(this%vf_rad_prob) +!TODO call HTSD_model_get_radiosity_view (this, n, u, qrad) +!TODO faces => this%vf_rad_prob(n)%faces +!TODO !! Radiative heat flux contribution to the heat conduction face residual. +!TODO allocate(flux(size(faces))) +!TODO call this%vf_rad_prob(n)%heat_flux (t, qrad, Tface(faces), flux) +!TODO do j = 1, size(faces) +!TODO if (this%vf_rad_prob(n)%fmask(j)) & +!TODO Fface(faces(j)) = Fface(faces(j)) + this%mesh%area(faces(j)) * flux(j) +!TODO end do +!TODO deallocate(flux) +!TODO !! Residual of the algebraic radiosity system. +!TODO call HTSD_model_get_radiosity_view (this, n, f, fptr) +!TODO call this%vf_rad_prob(n)%residual (t, qrad, Tface(faces), fptr) +!TODO fptr = -fptr +!TODO end do +!TODO end if + + !TODO: is this necessary? Off-process values are not needed, but may be + ! used in dummy vector operations, and we don't want fp exceptions. + call gather_boundary(this%mesh%cell_ip, f%hc) + call gather_boundary(this%mesh%cell_ip, f%tc) + call gather_boundary(this%mesh%face_ip, f%tf) + + call stop_timer('ht-function') + + end subroutine compute_f + + + subroutine update_moving_vf(this) + class(ht_model), intent(inout) :: this + integer :: n + if (.not.associated(this%vf_rad_prob)) return + do n = 1, size(this%vf_rad_prob) + call this%vf_rad_prob(n)%update_moving_vf + end do + end subroutine + + subroutine add_moving_vf_events(this, eventq) + use sim_event_queue_type + class(ht_model), intent(inout) :: this + type(sim_event_queue), intent(inout) :: eventq + integer :: n + if (.not.associated(this%vf_rad_prob)) return + do n = 1, size(this%vf_rad_prob) + call this%vf_rad_prob(n)%add_moving_vf_events(eventq) + end do + end subroutine + +end module ht_model_type diff --git a/src/truchas/physics/heat_species_transport/ht_norm_type.F90 b/src/truchas/physics/heat_species_transport/ht_norm_type.F90 new file mode 100644 index 000000000..a9e902c52 --- /dev/null +++ b/src/truchas/physics/heat_species_transport/ht_norm_type.F90 @@ -0,0 +1,168 @@ +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! +!! This file is part of Truchas. 3-Clause BSD license; see the LICENSE file. +!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +#include "f90_assert.fpp" + +module ht_norm_type + + use,intrinsic :: iso_fortran_env, only: r8 => real64 + use ht_vector_type + use ht_model_type + use rad_problem_type + use parallel_communication + implicit none + private + + type, public :: ht_norm + private + type(ht_model), pointer :: model => null() + real(r8) :: temp_atol + real(r8) :: temp_rtol + real(r8) :: enth_atol + real(r8) :: enth_rtol + contains + procedure :: init + procedure :: compute + end type ht_norm + +contains + + subroutine init(this, model, params, stat, errmsg) + + use parameter_list_type + + class(ht_norm), intent(out) :: this + type(ht_model), intent(in), target :: model + type(parameter_list) :: params + integer, intent(out) :: stat + character(:), allocatable, intent(out) :: errmsg + + character(:), allocatable :: context + + this%model => model + + context = 'processing ' // params%name() // ': ' + call params%get('abs-t-tol', this%temp_atol, default=0.0_r8, stat=stat, errmsg=errmsg) + if (stat /= 0) then + errmsg = context // errmsg + return + else if (this%temp_atol < 0) then + stat = 1 + errmsg = context // '"abs-t-tol" must be >= 0.0' + return + end if + + call params%get('rel-t-tol', this%temp_rtol, default=0.0_r8, stat=stat, errmsg=errmsg) + if (stat /= 0) then + errmsg = context // errmsg + return + else if (this%temp_rtol < 0) then + stat = 1 + errmsg = context // '"rel-t-tol" must be >= 0.0' + return + end if + + if (this%temp_atol == 0 .and. this%temp_rtol == 0) then + stat = 1 + errmsg = context // '"abs-t-tol" and "rel-t-tol" cannot both be 0.0' + return + end if + + call params%get('abs-h-tol', this%enth_atol, default=0.0_r8, stat=stat, errmsg=errmsg) + if (stat /= 0) then + errmsg = context // errmsg + return + else if (this%enth_atol < 0) then + stat = 1 + errmsg = context // '"abs-h-tol" must be >= 0.0' + return + end if + + call params%get('rel-h-tol', this%enth_rtol, default=this%temp_rtol, stat=stat, errmsg=errmsg) + if (stat /= 0) then + errmsg = context // errmsg + return + else if (this%enth_rtol < 0) then + stat = 1 + errmsg = context // '"rel-h-tol" must be >= 0.0' + return + end if + + if (this%enth_atol == 0 .and. this%enth_rtol == 0) then + stat = 1 + errmsg = context // '"abs-h-tol" and "rel-h-tol" cannot both be 0.0' + return + end if + + end subroutine init + + + subroutine compute(this, t, u, du, du_norm) + + class(ht_norm), intent(in) :: this + real(r8), intent(in) :: t + type(ht_vector), intent(in) :: u, du + real(r8), intent(out) :: du_norm + + !integer :: n + !real(r8) :: qerror + !real(r8), pointer :: qrad(:) + !integer, pointer :: faces(:) + !real(r8), allocatable :: res(:), rhs(:) + integer :: ncell, nface + + ncell = this%model%mesh%ncell_onP + nface = this%model%mesh%nface_onP + + du_norm = 0.0_r8 + du_norm = max(du_norm, maxerr(u%hc(:ncell), du%hc(:ncell), this%enth_atol, this%enth_rtol, this%model%void_cell)) + du_norm = max(du_norm, maxerr(u%tc(:ncell), du%tc(:ncell), this%temp_atol, this%temp_rtol, this%model%void_cell)) + du_norm = max(du_norm, maxerr(u%tf(:nface), du%tf(:nface), this%temp_atol, this%temp_rtol, this%model%void_face)) + du_norm = global_maxval(du_norm) + !! Enclosure radiation system error norms + !TODO! This is a quick hack that needs to be fixed. The tolerance + !TODO! is hardwired, and we are (re)computing the residual. The BDF2 + !TODO! integrator thinks it's computing the norm of the correction du but + !TODO! in this case it is the actual residual norm. We need more general + !TODO! norms in the integrator. +!TODO if (associated(this%model%ht%vf_rad_prob)) then +!TODO !if (is_IOP) write(*,'(a)',advance='no')'ER error: ||res||/||rhs||=' +!TODO do n = 1, size(this%model%ht%vf_rad_prob) +!TODO faces => this%model%ht%vf_rad_prob(n)%faces +!TODO allocate(res(size(faces)), rhs(size(faces))) +!TODO call HTSD_model_get_radiosity_view (this%model, n, u, qrad) +!TODO call HTSD_model_get_face_temp_view (this%model, u, temp) +!TODO call this%model%ht%vf_rad_prob(n)%residual (t, qrad, temp(faces), res) +!TODO call this%model%ht%vf_rad_prob(n)%rhs (t, temp(faces), rhs) +!TODO qerror = sqrt(global_sum(res**2)) / sqrt(global_sum(rhs**2)) +!TODO !if (is_IOP) write(*,'(e10.3)',advance='no') qerror +!TODO ht_du_norm = max(ht_du_norm, qerror/1.0d-3) +!TODO deallocate(res, rhs) +!TODO end do +!TODO !if (is_IOP) write(*,*) +!TODO end if + + contains + + real(r8) function maxerr (u, du, atol, rtol, void) + real(r8), intent(in) :: u(:), du(:), atol, rtol + logical, pointer :: void(:) + real(r8) :: array(size(du)) + if (associated(void)) then + where (void(:size(du))) + array = 0.0_r8 + elsewhere + array = abs(du) / (atol + rtol*abs(u)) + end where + else + array = abs(du) / (atol + rtol*abs(u)) + end if + maxerr = maxval(array) + end function maxerr + + end subroutine compute + +end module ht_norm_type diff --git a/src/truchas/physics/heat_species_transport/ht_precon_type.F90 b/src/truchas/physics/heat_species_transport/ht_precon_type.F90 new file mode 100644 index 000000000..0dc4d563c --- /dev/null +++ b/src/truchas/physics/heat_species_transport/ht_precon_type.F90 @@ -0,0 +1,337 @@ +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! +!! This file is part of Truchas. 3-Clause BSD license; see the LICENSE file. +!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +#include "f90_assert.fpp" + +module ht_precon_type + + use,intrinsic :: iso_fortran_env, only: r8 => real64 + use ht_vector_type + use ht_model_type + use rad_problem_type + use unstr_mesh_type + use mfd_diff_precon_type + use mfd_diff_matrix_type + use index_partitioning + use truchas_timers + implicit none + private + + type, public :: ht_precon + type(ht_model), pointer :: model => null() + type(unstr_mesh), pointer :: mesh => null() + real(r8) :: dt + real(r8), allocatable :: dHdT(:) + integer, allocatable :: vfr_precon_coupling(:) + type(mfd_diff_precon) :: pc + contains + procedure :: init + procedure :: compute + procedure :: apply + end type ht_precon + + !! Methods of coupling heat conduction and radiosity preconditioning. + integer, parameter :: VFR_JAC = 1 ! Jacobi (radiosity and conduction decoupled) + integer, parameter :: VFR_FGS = 2 ! Forward Gauss-Seidel (radiosity, then conduction) + integer, parameter :: VFR_BGS = 3 ! Backward Gauss-Seidel (conduction, then radiosity) + integer, parameter :: VFR_FAC = 4 ! Factorization (radiosity, conduction, radiosity) + +contains + + subroutine init(this, model, params, stat, errmsg) + + use parameter_list_type + + class(ht_precon), intent(out) :: this + type(ht_model), intent(in), target :: model + type(parameter_list), intent(inout) :: params + integer, intent(out) :: stat + character(:), allocatable, intent(out) :: errmsg + + type(mfd_diff_matrix), allocatable :: matrix + + this%model => model + this%mesh => model%mesh + + allocate(this%dHdT(this%mesh%ncell)) + allocate(matrix) + call matrix%init(model%disc) + call this%pc%init(matrix, params, stat, errmsg) + + if (associated(model%vf_rad_prob)) then + block + integer :: n, j + character(:), allocatable :: array(:) + n = size(model%vf_rad_prob) + allocate(this%vfr_precon_coupling(n)) + call params%get('vfr-precon-coupling', array) + INSIST(size(array) == n) + do j = 1, n + select case (array(j)) + case ('JACOBI') + this%vfr_precon_coupling(j) = VFR_JAC + case ('FORWARD GS') + this%vfr_precon_coupling(j) = VFR_FGS + case ('BACKWARD GS') + this%vfr_precon_coupling(j) = VFR_BGS + case ('FACTORIZATION') + this%vfr_precon_coupling(j) = VFR_FAC + case default + INSIST(.false.) + end select + end do + end block + end if + + end subroutine init + + + subroutine compute(this, t, u, dt) + + class(ht_precon), intent(inout) :: this + real(r8), intent(in) :: t, dt + type(ht_vector), intent(inout) :: u + target :: u + + real(r8), pointer :: state(:,:) => null() + integer, allocatable :: more_dir_faces(:) + + integer :: index, j, n, n1, n2 + real(r8) :: D(this%mesh%ncell), A(this%mesh%ncell), term + real(r8), allocatable :: values(:), values2(:,:) + type(mfd_diff_matrix), pointer :: dm + integer, pointer :: faces(:) => null() + + ASSERT(dt > 0.0_r8) + + call start_timer('ht-precon-compute') + + !! Generate list of void faces; these are treated like Dirichlet BC. + if (associated(this%model%void_face)) then + n = count(this%model%void_face) + allocate(more_dir_faces(n)) + n = 0 + do j = 1, this%mesh%nface + if (this%model%void_face(j)) then + n = n + 1 + more_dir_faces(n) = j + end if + end do + else + allocate(more_dir_faces(0)) + end if + + !TODO: The existing prop_mesh_func%compute_value function expects a rank-2 + ! state array. This is a workaround until prop_mesh_func is redesigned. + state(1:this%mesh%ncell,1:1) => u%tc + + call gather_boundary(this%mesh%cell_ip, u%tc) + call gather_boundary(this%mesh%face_ip, u%tf) + + this%dt = dt + + !! Hardwired assumption that T is the first component of state -- FIXME! + call this%model%H_of_T%compute_deriv(state, 1, this%dHdT) + call this%model%conductivity%compute_value(state, D) + A = this%mesh%volume * this%dHdT / dt + + !! Correct data on void cells. + if (associated(this%model%void_cell)) then + where (this%model%void_cell) + this%dHdT = 0.0_r8 + D = 0.0_r8 + A = 1.0_r8 + end where + end if + + !! Jacobian of the basic heat equation that ignores nonlinearities + !! in the conductivity. This has the H/T relation eliminated. + dm => this%pc%matrix_ref() + call dm%compute(D) + call dm%incr_cell_diag(A) + + !! Dirichlet boundary condition fixups. + if (allocated(this%model%bc_dir)) then + call this%model%bc_dir%compute(t) + call dm%set_dir_faces(this%model%bc_dir%index) + end if + + !! External HTC boundary condition contribution. + if (allocated(this%model%bc_htc)) then + call this%model%bc_htc%compute_deriv(t, u%tf) + associate (index => this%model%bc_htc%index, & + deriv => this%model%bc_htc%deriv) + call dm%incr_face_diag(index, deriv) + end associate + end if + + !! Simple radiation boundary condition contribution. + if (allocated(this%model%bc_rad)) then + call this%model%bc_rad%compute(t, u%tf) + associate (index => this%model%bc_rad%index, & + deriv => this%model%bc_rad%deriv) + call dm%incr_face_diag(index, deriv) + end associate + end if + + !! Experimental evaporation heat flux contribution. + if (allocated(this%model%evap_flux)) then + call this%model%evap_flux%compute_deriv(t, u%tf) + associate (index => this%model%evap_flux%index, & + deriv => this%model%evap_flux%deriv) + call dm%incr_face_diag(index, this%mesh%area(index)*deriv) + end associate + endif + + !! Internal HTC interface condition contribution. + if (allocated(this%model%ic_htc)) then + call this%model%ic_htc%compute(t, u%tf) + associate (index => this%model%ic_htc%index, & + deriv => this%model%ic_htc%deriv) + if (associated(this%model%void_face)) then + do j = 1, size(index,2) !FIXME? Bad form to modify deriv? + if (any(this%model%void_face(index(:,j)))) deriv(:,j) = 0.0_r8 + end do + end if + call dm%incr_interface_flux3(index, deriv) !TODO: rename these methods + end associate + end if + + !! Internal gap radiation condition contribution. + if (allocated(this%model%ic_rad)) then + call this%model%ic_rad%compute(t, u%tf) + associate (index => this%model%ic_rad%index, & + deriv => this%model%ic_rad%deriv) + if (associated(this%model%void_face)) then + do j = 1, size(index,2) !FIXME? Bad form to modify deriv? + if (any(this%model%void_face(index(:,j)))) deriv(:,j) = 0.0_r8 + end do + end if + call dm%incr_interface_flux3(index, deriv) !TODO: rename these methods + end associate + end if + + !! Dirichlet fix-ups for void faces. + call dm%set_dir_faces(more_dir_faces) + + !! Enclosure radiation contributions to the preconditioner. + !! TODO: what about factorization coupling? Is this still correct? + if (associated(this%model%vf_rad_prob)) then + do index = 1, size(this%model%vf_rad_prob) + faces => this%model%vf_rad_prob(index)%faces + allocate(values(size(faces))) + call this%model%vf_rad_prob(index)%rhs_deriv(t, u%tf(faces), values) + where (.not.this%model%vf_rad_prob(index)%fmask) values = 0 + call dm%incr_face_diag(faces, this%mesh%area(faces) * values) + deallocate(values) + end do + end if + + !! The matrix is now complete; re-compute the preconditioner. + call this%pc%compute + + call stop_timer('ht-precon-compute') + + end subroutine compute + + + subroutine apply(this, t, u, f) + + class(ht_precon), intent(in) :: this + real(r8), intent(in) :: t + type(ht_vector), intent(inout) :: u ! data is intent(in) + type(ht_vector), intent(inout) :: f ! data is intent(inout) + + integer :: index, j, n + real(r8), pointer :: fq(:) + real(r8), allocatable :: z(:) + + call start_timer('ht-precon-apply') + + !! Precondition the radiosity components: + !! Factorization and forward Gauss-Seidel coupling methods. + if (associated(this%model%vf_rad_prob)) then + call start_timer('vf-rad-precon') + !call HTSD_model_get_face_temp_view (this%model, f, f2) + do index = 1, size(this%model%vf_rad_prob) + if (this%vfr_precon_coupling(index) == VFR_FGS .or. & + this%vfr_precon_coupling(index) == VFR_FAC) then + !TODO: call HTSD_model_get_radiosity_view (this%model, index, f, fq) + allocate(z(size(fq))) + z = fq + call this%model%vf_rad_prob(index)%precon(t, z) + if (this%vfr_precon_coupling(index) == VFR_FGS) fq = z + !! Update the heat equation face residual. + call this%model%vf_rad_prob(index)%precon_matvec1(t, z) + do j = 1, size(z) + if (this%model%vf_rad_prob(index)%fmask(j)) then + n = this%model%vf_rad_prob(index)%faces(j) + f%tf(n) = f%tf(n) + this%mesh%area(n) * z(j) + end if + end do + deallocate(z) + end if + end do + call stop_timer('vf-rad-precon') + end if + + !! Heat equation cell residual with the H/T relation residual eliminated. + !call HTSD_model_get_cell_heat_view (this%model, f, f0) + !call HTSD_model_get_cell_temp_view (this%model, f, f1) + !f1x(:this%mesh%ncell_onP) = f1 - (this%mesh%volume(:this%mesh%ncell_onP)/this%dt)*f0 + !f%tc = f%tc - (this%mesh%volume/dt)*f%hc + + if (associated(this%model%void_cell)) then + where (.not.this%model%void_cell) f%tc = f%tc - (this%mesh%volume/this%dt)*f%hc + else + f%tc = f%tc - (this%mesh%volume/this%dt)*f%hc + end if + call gather_boundary(this%mesh%cell_ip, f%tc) + + !! Heat equation face residual (with radiosity residuals optionally eliminated). + !call HTSD_model_get_face_temp_copy (this%model, f, f2x) + !call gather_boundary (this%mesh%face_ip, f2x) + call gather_boundary(this%mesh%face_ip, f%tf) + + !! Precondition the heat equation. + !call diff_precon_apply (this%hcprecon, f1x, f2x) + call this%pc%apply(f%tc, f%tf) + !call HTSD_model_set_cell_temp (this%model, f1x, f) + !call HTSD_model_set_face_temp (this%model, f2x, f) + + !! Backsubstitute to obtain the preconditioned H/T-relation residual. + !f0 = f0 + this%dHdT(:this%mesh%ncell_onP)*f1 + f%hc = f%hc + this%dHdT*f%tc + + !! Precondition the radiosity components: + !! Factorization, backward Gauss-Seidel and Jacobi coupling methods. + if (associated(this%model%vf_rad_prob)) then + call start_timer('vf-rad-precon') + !call HTSD_model_get_face_temp_view (this%model, f, f2) + !call HTSD_model_get_face_temp_view (this%model, u, Tface) + do index = 1, size(this%model%vf_rad_prob) + if (this%vfr_precon_coupling(index) == VFR_JAC .or. & + this%vfr_precon_coupling(index) == VFR_BGS .or. & + this%vfr_precon_coupling(index) == VFR_FAC) then + !! Update the radiosity residual components. + !TODO: call HTSD_model_get_radiosity_view (this%model, index, f, fq) + if (this%vfr_precon_coupling(index) /= VFR_JAC) then + allocate(z(size(fq))) + call this%model%vf_rad_prob(index)%rhs_deriv(t, u%tf(this%model%vf_rad_prob(index)%faces), z) + fq = fq + z * f%tf(this%model%vf_rad_prob(index)%faces) + deallocate(z) + end if + call this%model%vf_rad_prob(index)%precon(t, fq) + end if + end do + call stop_timer('vf-rad-precon') + end if + + call stop_timer('ht-precon-apply') + + end subroutine apply + +end module ht_precon_type diff --git a/src/truchas/physics/heat_species_transport/ht_solver_factory.F90 b/src/truchas/physics/heat_species_transport/ht_solver_factory.F90 new file mode 100644 index 000000000..fb8da81bb --- /dev/null +++ b/src/truchas/physics/heat_species_transport/ht_solver_factory.F90 @@ -0,0 +1,78 @@ +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! +!! This file is part of Truchas. 3-Clause BSD license; see the LICENSE file. +!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +#include "f90_assert.fpp" + +module ht_solver_factory + + use,intrinsic :: iso_fortran_env, only: r8 => real64 + use ht_model_type + use ht_solver_type + use matl_mesh_func_type + use parameter_list_type + implicit none + private + + public :: create_ht_solver + +contains + + function create_ht_solver(mmf, model, params, stat, errmsg) result(solver) + + use enclosure_radiation_namelist, only: er_params => params + use parallel_communication + use truchas_env, only: output_file_name + + type(matl_mesh_func), intent(in), target :: mmf + type(ht_model), intent(in), target :: model + type(parameter_list), intent(inout) :: params + type(ht_solver), pointer :: solver + integer, intent(out) :: stat + character(:), allocatable, intent(out) :: errmsg + + integer :: j, n, lun + type(parameter_list_iterator) :: piter + type(parameter_list), pointer :: plist + character(:), allocatable :: string + character(16), allocatable :: vfr_precon_coupling(:) + real(r8), allocatable :: rad_tol(:) + logical :: verbose_stepping + + call params%get('verbose-stepping', verbose_stepping) + if (verbose_stepping) then + lun = -1 + if (is_IOP) open(newunit=lun,file=output_file_name('bdf2.out'),position='rewind',action='write') + call params%set('output-unit', lun) + plist => params%sublist('norm') + call plist%set('unit', lun) + end if + +!TODO !TODO: read enclosure radiation namelists directly into sublist +!TODO if (associated(model%ht)) then +!TODO if (associated(model%ht%vf_rad_prob)) then +!TODO !TODO! Assumes model%ht%vf_rad_prob array was initialized under the same +!TODO !TODO! loop so that that array and vfr_precon_coupling are in correspondence. +!TODO !TODO! This is fragile and needs to be fixed. +!TODO piter = parameter_list_iterator(er_params) +!TODO n = piter%count() +!TODO allocate(vfr_precon_coupling(n), rad_tol(n)) +!TODO do j = 1, n +!TODO plist => piter%sublist() +!TODO call plist%get('precon-coupling-method', string, default='BACKWARD GS') +!TODO vfr_precon_coupling(j) = string +!TODO call piter%next +!TODO end do +!TODO plist => params%sublist('precon') +!TODO call plist%set('vfr-precon-coupling', vfr_precon_coupling) +!TODO end if +!TODO end if + + allocate(solver) + call solver%init(mmf, model, params, stat, errmsg) + + end function create_ht_solver + +end module ht_solver_factory diff --git a/src/truchas/physics/heat_species_transport/ht_solver_type.F90 b/src/truchas/physics/heat_species_transport/ht_solver_type.F90 new file mode 100644 index 000000000..ca9a50fe0 --- /dev/null +++ b/src/truchas/physics/heat_species_transport/ht_solver_type.F90 @@ -0,0 +1,319 @@ +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!! +!! This file is part of Truchas. 3-Clause BSD license; see the LICENSE file. +!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +#include "f90_assert.fpp" + +module ht_solver_type + + use,intrinsic :: iso_fortran_env, only: r8 => real64 + use ht_vector_type + use ht_model_type + use ht_precon_type + use ht_norm_type + use matl_mesh_func_type + use unstr_mesh_type + use index_partitioning + use ht_idaesol_model_type + use new_idaesol_type + use parameter_list_type + implicit none + private + + type, public :: ht_solver + type(matl_mesh_func), pointer :: mmf => null() + type(unstr_mesh), pointer :: mesh => null() + type(ht_idaesol_model) :: integ_model + type(idaesol) :: integ + logical :: state_is_pending = .false. + !! Pending state + real(r8) :: t, dt + type(ht_vector) :: u + type(ht_model), pointer :: model => null() + type(ht_precon), pointer :: precon => null() + type(ht_norm), pointer :: norm => null() + type(parameter_list), pointer :: ic_params => null() + contains + procedure :: init + procedure :: step + procedure :: commit_pending_state + procedure :: get_cell_heat_copy, get_cell_heat_view + procedure :: get_cell_temp_copy, get_cell_temp_view + procedure :: get_face_temp_copy, get_face_temp_view + procedure :: get_cell_temp_grad + procedure :: get_stepping_stats + procedure :: last_step_size + procedure :: last_time + procedure :: set_initial_state + procedure :: restart + final :: ht_solver_delete + end type + +contains + + !! Final subroutine for HT_SOLVER objects. + subroutine ht_solver_delete (this) + type(ht_solver), intent(inout) :: this + if (associated(this%precon)) deallocate(this%precon) + if (associated(this%norm)) deallocate(this%norm) + !! The solver owns the void cell and face mask components of the model. + if (associated(this%model)) then !TODO: SOLVER MUST BE FINALIZED BEFORE MODEL! + if (associated(this%model%void_cell)) deallocate(this%model%void_cell) + if (associated(this%model%void_face)) deallocate(this%model%void_face) + end if + if (associated(this%ic_params)) deallocate(this%ic_params) + end subroutine + +!TODO: model should hold a reference to mmf + + subroutine init(this, mmf, model, params, stat, errmsg) + + use parallel_communication, only: is_IOP + + class(ht_solver), intent(out), target :: this + type(matl_mesh_func), intent(in), target :: mmf + type(ht_model), intent(in), target :: model + type(parameter_list), intent(inout) :: params + integer, intent(out) :: stat + character(:), allocatable, intent(out) :: errmsg + + integer :: output_unit + type(parameter_list), pointer :: plist + logical :: verbose_stepping + + this%mmf => mmf + this%model => model + this%mesh => model%mesh + + allocate(this%norm) + plist => params%sublist('norm') + call this%norm%init(model, plist, stat, errmsg) + if (stat /= 0) return + + allocate(this%precon) + plist => params%sublist('precon') + call this%precon%init(model, plist, stat, errmsg) + if (stat /= 0) return + + call this%u%init(this%mesh) + + call this%integ_model%init(this%model, this%precon, this%norm) + + block + integer :: stat + character(:), allocatable :: errmsg + call this%integ%init(this%integ_model, params, stat, errmsg) + INSIST(stat == 0) + end block + + call params%get('verbose-stepping', verbose_stepping) + if (is_IOP .and. verbose_stepping) then + call params%get('output-unit', output_unit) + call this%integ%set_verbose_stepping(output_unit) + end if + + !! Grab parameters for ht_ic_solver%init + allocate(this%ic_params) + block + real(r8) :: rval + plist => params%sublist('norm') + call plist%get('abs-t-tol', rval) + call this%ic_params%set('atol-temp', 0.01_r8 * rval) + call plist%get('rel-t-tol', rval) + call this%ic_params%set('rtol-temp', 0.01_r8 * rval) + call this%ic_params%set('max-iter', 50) + call this%ic_params%set('method', 'SSOR') + plist => this%ic_params%sublist('params') + call plist%set('num-cycles', 1) + end block + + end subroutine init + + !! Get a reference to the pending/current cell temperatures + subroutine get_cell_temp_view(this, view) + class(ht_solver), intent(in), target :: this + real(r8), pointer :: view(:) + view => this%u%tc + end subroutine + + !! Get a copy of the pending/current cell temperatures + subroutine get_cell_temp_copy(this, copy) + class(ht_solver), intent(in) :: this + real(r8), intent(inout) :: copy(:) + integer :: n + n = min(size(copy), size(this%u%tc)) + copy(1:n) = this%u%tc(1:n) + end subroutine + + !! Get a reference to the pending/current face temperatures + subroutine get_face_temp_view(this, view) + class(ht_solver), intent(in), target :: this + real(r8), pointer :: view(:) + view => this%u%tf + end subroutine + + !! Get a copy of the pending/current face temperatures + subroutine get_face_temp_copy(this, copy) + class(ht_solver), intent(in) :: this + real(r8), intent(inout) :: copy(:) + integer :: n + n = min(size(copy), size(this%u%tf)) + copy(1:n) = this%u%tf(1:n) + end subroutine + + !! Get a reference to the pending/current cell enthalpies + subroutine get_cell_heat_view(this, view) + class(ht_solver), intent(in), target :: this + real(r8), pointer :: view(:) + view => this%u%hc + end subroutine + + !! Get a copy of the pending/current cell enthalpies + subroutine get_cell_heat_copy(this, copy) + class(ht_solver), intent(in) :: this + real(r8), intent(inout) :: copy(:) + integer :: n + n = min(size(copy), size(this%u%hc)) + copy(1:n) = this%u%hc(1:n) + end subroutine + + !! Compute a cell-based approximation to the gradient of the pending/current + !! cell temperatures. Note that this derived quantity is not used in the + !! discretization of the heat equation. + + subroutine get_cell_temp_grad(this, tgrad) + use mfd_disc_type + class(ht_solver), intent(inout) :: this + real(r8), intent(out) :: tgrad(:,:) + INSIST(size(tgrad,1) == 3) + INSIST(size(tgrad,2) == this%model%mesh%ncell_onP) + call gather_boundary(this%model%mesh%face_ip, this%u%tf) + call this%model%disc%compute_cell_grad(this%u%tf, tgrad) + end subroutine + + !! Return the time of the current (committed) solution + function last_time(this) result(t) + class(ht_solver), intent(in) :: this + real(r8) :: t + t = this%integ%last_time() + end function + + !! Return the size of the last successful (committed) time step + function last_step_size(this) result(h) + class(ht_solver), intent(in) :: this + real(r8) :: h + h = this%integ%last_step_size() + end function + + subroutine get_stepping_stats(this, counters) + class(ht_solver), intent(in) :: this + integer, intent(out) :: counters(:) + ASSERT(size(counters) == 6) + call this%integ%get_stepping_statistics(counters) + end subroutine + + !! Advance the current solution state by a single time step from the current + !! time to time T. STAT returns 0 if the step was successful. The advanced + !! state is only provisional, however, until COMMIT_PENDING_STATE is called + !! to commit it. If STAT returns a nonzero value the step has failed (STAT=1 + !! if the nonlinear iteration failed using a fresh preconditioner; STAT=2 if + !! if the step succeeded but the predictor error exceeded its tolerance). + !! HNEXT returns the suggested next time step size: if STAT == 0 this is the + !! estimated step size needed to keep the estimated predictor error in its + !! target range; otherwise this is a reduced step size to use when retrying + !! the time step. + + subroutine step(this, t, hnext, stat) + class(ht_solver), intent(inout) :: this + real(r8), intent(in) :: t + real(r8), intent(out) :: hnext + integer, intent(out) :: stat + call this%integ%step(t, this%u, hnext, stat) + if (stat == 0) then + this%t = t + this%state_is_pending = .true. + call this%u%gather_boundary() !TODO: Can the be made unnecessary? + else + call this%integ%get_last_state_copy(this%u) + this%state_is_pending = .false. + end if + end subroutine + + subroutine commit_pending_state(this) + class(ht_solver), intent(inout) :: this + if (this%state_is_pending) then + call this%integ%commit_state(this%t, this%u) + this%state_is_pending = .false. + end if + end subroutine + + subroutine set_initial_state(this, t, temp, dt) + + use ht_ic_solver_type + use parallel_communication, only: global_count + + class(ht_solver), intent(inout) :: this + real(r8), intent(in) :: t + real(r8), intent(in) :: temp(:) + real(r8), intent(in) :: dt + + integer :: j + real(r8), allocatable :: void_vol_frac(:) + type(ht_vector) :: udot + type(ht_ic_solver) :: ic + + INSIST(associated(this%model)) + + !! Establish the void cell and face masks in the model. + !! N.B. We assume no one else allocates and defines these arrays. + allocate(void_vol_frac(this%mesh%ncell), this%model%void_cell(this%mesh%ncell)) + call this%mmf%get_void_vol_frac(void_vol_frac) + this%model%void_cell = (void_vol_frac > 1.0_r8 - 2*epsilon(1.0_r8)) + deallocate(void_vol_frac) + if (global_count(this%model%void_cell) > 0) then + allocate(this%model%void_face(this%mesh%nface)) + this%model%void_face = .true. + do j = 1, this%mesh%ncell + if (.not.this%model%void_cell(j)) then + associate (cface => this%mesh%cface(this%mesh%xcface(j):this%mesh%xcface(j+1)-1)) + this%model%void_face(cface) = .false. + end associate + end if + end do + call gather_boundary(this%mesh%face_ip, this%model%void_face) + else + deallocate(this%model%void_cell) + this%model%void_face => null() + end if + + call udot%init(this%u) + call this%ic_params%set('dt', dt) + call ic%init(this%model, this%ic_params) + call ic%compute(t, temp, this%u, udot) + call this%integ%set_initial_state(t, this%u, udot) + + end subroutine set_initial_state + + subroutine restart(this, dt) + + use ht_ic_solver_type + + class(ht_solver), intent(inout) :: this + real(r8), intent(in) :: dt + + type(ht_vector) :: udot + type(ht_ic_solver) :: ic + + INSIST(associated(this%model)) + + call udot%init(this%u) + call this%ic_params%set('dt', dt) + call ic%init(this%model, this%ic_params) + call ic%compute_udot(this%t, this%u, udot) + call this%integ%set_initial_state(this%t, this%u, udot) + + end subroutine + +end module ht_solver_type diff --git a/src/truchas/physics/heat_species_transport/ht_vector_type.F90 b/src/truchas/physics/heat_species_transport/ht_vector_type.F90 new file mode 100644 index 000000000..fb8842817 --- /dev/null +++ b/src/truchas/physics/heat_species_transport/ht_vector_type.F90 @@ -0,0 +1,229 @@ + +module ht_vector_type + + use,intrinsic :: iso_fortran_env, only: r8 => real64 + use vector_class + use unstr_mesh_type + use parallel_communication, only: global_sum + implicit none + private + + type, extends(vector), public :: ht_vector + type(unstr_mesh), pointer :: mesh => null() + real(r8), allocatable :: hc(:), tc(:), tf(:) + contains + !! Deferred base class procedures + procedure :: clone1 + procedure :: clone2 + procedure :: copy_ + procedure :: setval + procedure :: scale + procedure :: update1_ + procedure :: update2_ + procedure :: update3_ + procedure :: update4_ + procedure :: dot_ + procedure :: norm2 => norm2_ + procedure :: checksum + !! Additional procedures specific to this type + generic :: init => init_mesh, init_mold + procedure, private :: init_mesh, init_mold + procedure :: gather_boundary => ht_gather_boundary + end type + +contains + + !! Specific subroutine for the generic INIT. Initialize a HT_VECTOR object + !! for the given unstructured MESH. The elements are initialized to 0. + + subroutine init_mesh(this, mesh) + class(ht_vector), intent(out) :: this + type(unstr_mesh), intent(in), target :: mesh + this%mesh => mesh + allocate(this%hc(mesh%ncell), this%tc(mesh%ncell), this%tf(mesh%nface)) + call setval(this, 0.0_r8) + end subroutine + + !! Specific subroutine for the generic INIT. Initialize a HT_VECTOR object + !! to be a clone of MOLD. The elements are initialized to zero. + + subroutine init_mold(this, mold) + class(ht_vector), intent(out) :: this + class(ht_vector), intent(in) :: mold + call init_mesh(this, mold%mesh) + end subroutine + + subroutine ht_gather_boundary(this) + use index_partitioning, only: gather_boundary + class(ht_vector), intent(inout) :: this + call gather_boundary(this%mesh%cell_ip, this%hc) + call gather_boundary(this%mesh%cell_ip, this%tc) + call gather_boundary(this%mesh%face_ip, this%tf) + end subroutine + + subroutine clone1(this, clone) + class(ht_vector), intent(in) :: this + class(vector), allocatable, intent(out) :: clone + type(ht_vector), allocatable :: tmp + allocate(tmp) + tmp%mesh => this%mesh + allocate(tmp%hc, mold=this%hc) + allocate(tmp%tc, mold=this%tc) + allocate(tmp%tf, mold=this%tf) + call move_alloc(tmp, clone) + end subroutine + + subroutine clone2(this, clone, n) + class(ht_vector), intent(in) :: this + class(vector), allocatable, intent(out) :: clone(:) + type(ht_vector), allocatable :: tmp(:) + integer, intent(in) :: n + integer :: j + allocate(tmp(n)) + do j = 1, n + tmp(j)%mesh => this%mesh + allocate(tmp(j)%hc, mold=this%hc) + allocate(tmp(j)%tc, mold=this%tc) + allocate(tmp(j)%tf, mold=this%tf) + end do + call move_alloc(tmp, clone) + end subroutine + + subroutine copy_(dest, src) + class(ht_vector), intent(inout) :: dest + class(vector), intent(in) :: src + select type (src) + class is (ht_vector) + dest%hc(:) = src%hc + dest%tc(:) = src%tc + dest%tf(:) = src%tf + end select + end subroutine + + subroutine setval(this, val) + class(ht_vector), intent(inout) :: this + real(r8), intent(in) :: val + this%hc = val + this%tc = val + this%tf = val + end subroutine + + subroutine scale(this, a) + class(ht_vector), intent(inout) :: this + real(r8), intent(in) :: a + this%hc = a * this%hc + this%tc = a * this%tc + this%tf = a * this%tf + end subroutine + + !! Conventional SAXPY procedure: y <-- a*x + y + subroutine update1_(this, a, x) + class(ht_vector), intent(inout) :: this + class(vector), intent(in) :: x + real(r8), intent(in) :: a + select type (x) + class is (ht_vector) + this%hc = a * x%hc + this%hc + this%tc = a * x%tc + this%tc + this%tf = a * x%tf + this%tf + end select + end subroutine + + !! SAXPY-like procedure: y <-- a*x + b*y + subroutine update2_(this, a, x, b) + class(ht_vector), intent(inout) :: this + class(vector), intent(in) :: x + real(r8), intent(in) :: a, b + select type (x) + class is (ht_vector) + this%hc = a * x%hc + b * this%hc + this%tc = a * x%tc + b * this%tc + this%tf = a * x%tf + b * this%tf + end select + end subroutine + + !! SAXPY-like procedure: z <-- a*x + b*y + z + subroutine update3_(this, a, x, b, y) + class(ht_vector), intent(inout) :: this + class(vector), intent(in) :: x, y + real(r8), intent(in) :: a, b + select type (x) + class is (ht_vector) + select type (y) + class is (ht_vector) + this%hc = a * x%hc + b * y%hc + this%hc + this%tc = a * x%tc + b * y%tc + this%tc + this%tf = a * x%tf + b * y%tf + this%tf + end select + end select + end subroutine + + !! SAXPY-like procedure: z <-- a*x + b*y + c*z + subroutine update4_(this, a, x, b, y, c) + class(ht_vector), intent(inout) :: this + class(vector), intent(in) :: x, y + real(r8), intent(in) :: a, b, c + select type (x) + class is (ht_vector) + select type (y) + class is (ht_vector) + this%hc = a * x%hc + b * y%hc + c * this%hc + this%tc = a * x%tc + b * y%tc + c * this%tc + this%tf = a * x%tf + b * y%tf + c * this%tf + end select + end select + end subroutine + + function dot_(x, y) result(dp) + class(ht_vector), intent(in) :: x + class(vector), intent(in) :: y + real(r8) :: dp + integer :: j + select type (y) + class is (ht_vector) + dp = 0 + do j = 1, x%mesh%ncell_onP + dp = dp + x%hc(j) * y%hc(j) + end do + do j = 1, x%mesh%ncell_onP + dp = dp + x%tc(j) * y%tc(j) + end do + do j = 1, x%mesh%nface_onP + dp = dp + x%tf(j) * y%tf(j) + end do + dp = global_sum(dp) + end select + end function + + function norm2_(this) + use parallel_communication, only: global_sum + class(ht_vector), intent(in) :: this + real(r8) :: norm2_ + norm2_ = norm2(this%hc(:this%mesh%ncell_onP))**2 + & + norm2(this%tc(:this%mesh%ncell_onP))**2 + & + norm2(this%tf(:this%mesh%nface_onP))**2 + norm2_ = sqrt(global_sum(norm2_)) + end function + + function checksum(this, full) result(string) + use md5_hash_type + class(ht_vector), intent(in) :: this + logical, intent(in), optional :: full ! default is FALSE + character(:), allocatable :: string + type(md5_hash) :: hash + logical :: strict + strict = .true. + if (present(full)) strict = .not.full + if (strict) then + call hash%update(this%hc(:this%mesh%ncell_onP)) + call hash%update(this%tc(:this%mesh%ncell_onP)) + call hash%update(this%tf(:this%mesh%nface_onP)) + else + call hash%update(this%hc(:this%mesh%ncell)) + call hash%update(this%tc(:this%mesh%ncell)) + call hash%update(this%tf(:this%mesh%nface)) + end if + string = hash%hexdigest() + end function + +end module ht_vector_type diff --git a/test/enthalpy/test2.py b/test/enthalpy/test2.py index c2efa42f7..ef1fe9366 100755 --- a/test/enthalpy/test2.py +++ b/test/enthalpy/test2.py @@ -18,14 +18,14 @@ def run_test(outA, outB): # test final temperature tempA = outA.field(2,"Z_TEMP") tempB = outB.field(2,"Z_TEMP") - nfail += truchas.compare_max_rel(tempA, tempB, 1e-15, "temperature", outA.time(2)) + nfail += truchas.compare_max_rel(tempA, tempB, 2e-8, "temperature", outA.time(2)) # test the enthalpy probe data filename = os.path.join(outA.directory, "probe1.dat") dataA = numpy.loadtxt(filename) filename = os.path.join(outB.directory, "probe1.dat") dataB = numpy.loadtxt(filename) - nfail += truchas.compare_max_rel(dataA[:,1], dataB[:,1], 1e-15, 'enthalpy', -1.0) + nfail += truchas.compare_max_rel(dataA[:,1], dataB[:,1], 2e-8, 'enthalpy', -1.0) truchas.report_summary(nfail) return nfail -- GitLab From d0051bca74a98629f0f0bec3200225cbeada2c1a Mon Sep 17 00:00:00 2001 From: "Neil N. Carlson" <nnc@lanl.gov> Date: Sun, 21 Feb 2021 10:10:07 -0700 Subject: [PATCH 2/6] WIP: Prototype HT solver using vector-class-based DAE integrator (#450) All test problems are now passing. --- .../enclosure_radiation/Test/test_rade.F90 | 7 +- .../enclosure_radiation_namelist.F90 | 8 +- .../enclosure_radiation/rad_problem_type.F90 | 6 + .../FHT_model_factory.F90 | 11 +- .../FHT_solver_factory.F90 | 4 +- .../HTSD_model_factory.F90 | 16 +-- .../HTSD_solver_factory.F90 | 5 +- .../diffusion_solver.F90 | 14 +-- .../ht_ic_solver_type.F90 | 76 ++++++------- .../ht_idaesol_model_type.F90 | 2 +- .../ht_model_factory.F90 | 24 ++-- .../heat_species_transport/ht_model_type.F90 | 60 +++++----- .../heat_species_transport/ht_norm_type.F90 | 40 ++++--- .../heat_species_transport/ht_precon_type.F90 | 41 ++++--- .../ht_solver_factory.F90 | 41 ++++--- .../heat_species_transport/ht_solver_type.F90 | 6 +- .../heat_species_transport/ht_vector_type.F90 | 105 ++++++++++++++---- test/vfrad-moving/test.py | 2 +- 18 files changed, 261 insertions(+), 207 deletions(-) diff --git a/src/truchas/physics/enclosure_radiation/Test/test_rade.F90 b/src/truchas/physics/enclosure_radiation/Test/test_rade.F90 index ab51d37de..2239f5761 100644 --- a/src/truchas/physics/enclosure_radiation/Test/test_rade.F90 +++ b/src/truchas/physics/enclosure_radiation/Test/test_rade.F90 @@ -14,10 +14,12 @@ program test_rade use truchas_logging_services use unstr_mesh_type use rad_problem_type + use parameter_list_type implicit none type(unstr_mesh), pointer :: mesh real(r8), allocatable :: hc_temp(:) + type(parameter_list) :: params !! Enclosures to be tested character(len=31), parameter :: ENCL_EXT_FACE = "exterior_face" @@ -81,7 +83,7 @@ contains call read_function_namelists (lun) !! Enclosure radiation namelists. - call read_enclosure_radiation_namelists(lun) + call read_enclosure_radiation_namelists(lun, params) !! The mesh must be named 'main' call enable_mesh ('main', exists) @@ -118,9 +120,6 @@ contains !! Test that radiosity solver converges within given error tolerance subroutine test_converge (mesh, encl_name, hc_temp) - use enclosure_radiation_namelist, only: params - use parameter_list_type - type(unstr_mesh), pointer, intent(in) :: mesh character(len=31), intent(in) :: encl_name real(r8), intent(in) :: hc_temp(:) diff --git a/src/truchas/physics/enclosure_radiation/enclosure_radiation_namelist.F90 b/src/truchas/physics/enclosure_radiation/enclosure_radiation_namelist.F90 index dd285a261..28807ffdb 100644 --- a/src/truchas/physics/enclosure_radiation/enclosure_radiation_namelist.F90 +++ b/src/truchas/physics/enclosure_radiation/enclosure_radiation_namelist.F90 @@ -26,15 +26,14 @@ module enclosure_radiation_namelist public :: read_enclosure_radiation_namelists - type(parameter_list), public :: params - integer, parameter :: MAX_NAME_LEN = 31, MAX_FILE_LEN = 255, MAX_FACE_BLOCK_IDS = 32 contains - subroutine read_enclosure_radiation_namelists(lun) + subroutine read_enclosure_radiation_namelists(lun, params) integer, intent(in) :: lun + type(parameter_list), intent(inout) :: params integer :: n, ios logical :: found @@ -199,9 +198,10 @@ contains end subroutine read_enclosure_radiation_namelists - subroutine read_enclosure_surface_namelists(lun) + subroutine read_enclosure_surface_namelists(lun, params) integer, intent(in) :: lun + type(parameter_list), intent(inout) :: params integer :: n, ios logical :: found diff --git a/src/truchas/physics/enclosure_radiation/rad_problem_type.F90 b/src/truchas/physics/enclosure_radiation/rad_problem_type.F90 index 349177b5c..b0695af44 100644 --- a/src/truchas/physics/enclosure_radiation/rad_problem_type.F90 +++ b/src/truchas/physics/enclosure_radiation/rad_problem_type.F90 @@ -53,6 +53,7 @@ module rad_problem_type procedure :: rhs_deriv procedure :: update_moving_vf procedure :: add_moving_vf_events + procedure :: size => prob_size end type rad_problem integer, parameter, public :: DT_POLICY_NONE = 0 @@ -1183,4 +1184,9 @@ contains end select end function + elemental integer function prob_size(this) + class(rad_problem), intent(in) :: this + prob_size = size(this%faces) + end function + end module rad_problem_type diff --git a/src/truchas/physics/heat_species_transport/FHT_model_factory.F90 b/src/truchas/physics/heat_species_transport/FHT_model_factory.F90 index ad90757a2..56203c07d 100644 --- a/src/truchas/physics/heat_species_transport/FHT_model_factory.F90 +++ b/src/truchas/physics/heat_species_transport/FHT_model_factory.F90 @@ -34,6 +34,7 @@ module FHT_model_factory use rad_problem_type use thermal_bc_factory_class use thermal_source_factory_type + use parameter_list_type implicit none private @@ -41,7 +42,7 @@ module FHT_model_factory contains - function create_FHT_model (tinit, disc, mmf, tbc_fac, tsrc_fac, stat, errmsg) result (model) + function create_FHT_model (tinit, disc, mmf, tbc_fac, tsrc_fac, er_params, stat, errmsg) result (model) use diffusion_solver_data, only: void_temperature @@ -50,6 +51,7 @@ contains type(matl_mesh_func), intent(in), target :: mmf class(thermal_bc_factory), intent(inout) :: tbc_fac type(thermal_source_factory), intent(inout) :: tsrc_fac + type(parameter_list), intent(inout) :: er_params integer, intent(out) :: stat character(*), intent(out) :: errmsg character(:), allocatable :: errmsg2 @@ -62,7 +64,7 @@ contains allocate(model) !! Initializes the VF_RAD_PROB components of MODEL. - call vf_rad_init (tinit, mesh, model, stat, errmsg) + call vf_rad_init (tinit, mesh, model, er_params, stat, errmsg) if (stat /= 0) return !! Defines the heat equation parameter components of MODEL. @@ -84,16 +86,15 @@ contains end function create_FHT_model - subroutine vf_rad_init (tinit, mesh, model, stat, errmsg) + subroutine vf_rad_init (tinit, mesh, model, params, stat, errmsg) - use enclosure_radiation_namelist, only: params use bitfield_type, only: btest use parallel_communication, only: global_any, global_all - use parameter_list_type real(r8), intent(in) :: tinit type(unstr_mesh), intent(in) :: mesh type(FHT_model), intent(inout) :: model + type(parameter_list), intent(inout) :: params integer, intent(out) :: stat character(len=*), intent(out) :: errmsg diff --git a/src/truchas/physics/heat_species_transport/FHT_solver_factory.F90 b/src/truchas/physics/heat_species_transport/FHT_solver_factory.F90 index 9c596adbb..687608504 100644 --- a/src/truchas/physics/heat_species_transport/FHT_solver_factory.F90 +++ b/src/truchas/physics/heat_species_transport/FHT_solver_factory.F90 @@ -20,7 +20,7 @@ module FHT_solver_factory contains - function create_FHT_solver(mmf, model, params) result(solver) + function create_FHT_solver(mmf, model, params, er_params) result(solver) use enclosure_radiation_namelist, only: er_params => params use parallel_communication @@ -28,7 +28,7 @@ contains type(matl_mesh_func), intent(in), target :: mmf type(FHT_model), intent(in), target :: model - type(parameter_list), intent(inout) :: params + type(parameter_list), intent(inout) :: params, er_params type(FHT_solver), pointer :: solver integer :: j, n, lun diff --git a/src/truchas/physics/heat_species_transport/HTSD_model_factory.F90 b/src/truchas/physics/heat_species_transport/HTSD_model_factory.F90 index c2e48429e..af497a768 100644 --- a/src/truchas/physics/heat_species_transport/HTSD_model_factory.F90 +++ b/src/truchas/physics/heat_species_transport/HTSD_model_factory.F90 @@ -36,6 +36,7 @@ module HTSD_model_factory use species_bc_factory_class use thermal_source_factory_type use truchas_logging_services + use parameter_list_type implicit none private @@ -43,7 +44,7 @@ module HTSD_model_factory contains - function create_HTSD_model (tinit, disc, mmf, tbc_fac, sbc_fac, tsrc_fac, stat, errmsg) result (model) + function create_HTSD_model (tinit, disc, mmf, tbc_fac, sbc_fac, tsrc_fac, er_params, stat, errmsg) result (model) use diffusion_solver_data, only: heat_eqn, num_species, void_temperature @@ -53,6 +54,7 @@ contains class(thermal_bc_factory), intent(inout) :: tbc_fac class(species_bc_factory), intent(inout) :: sbc_fac type(thermal_source_factory), intent(inout) :: tsrc_fac + type(parameter_list), intent(inout) :: er_params integer, intent(out) :: stat character(*), intent(out) :: errmsg type(HTSD_model), pointer :: model @@ -64,7 +66,7 @@ contains mesh => disc%mesh if (heat_eqn) then - htmodel => create_HT_model(tinit, mesh, mmf, tbc_fac, tsrc_fac, stat, errmsg) + htmodel => create_HT_model(tinit, mesh, mmf, tbc_fac, tsrc_fac, er_params, stat, errmsg) if (stat /= 0) return endif @@ -80,7 +82,7 @@ contains end function create_HTSD_model - function create_HT_model (tinit, mesh, mmf, bc_fac, src_fac, stat, errmsg) result (model) + function create_HT_model (tinit, mesh, mmf, bc_fac, src_fac, er_params, stat, errmsg) result (model) use rad_problem_type @@ -89,6 +91,7 @@ contains type(matl_mesh_func), intent(in), target :: mmf class(thermal_bc_factory), intent(inout) :: bc_fac type(thermal_source_factory), intent(inout) :: src_fac + type(parameter_list), intent(inout) :: er_params integer, intent(out) :: stat character(*), intent(out) :: errmsg character(:), allocatable :: errmsg2 @@ -97,7 +100,7 @@ contains allocate(model) !! Initializes the VF_RAD_PROB components of MODEL. - model%vf_rad_prob => create_vf_rad_prob(tinit, mesh, stat, errmsg) + model%vf_rad_prob => create_vf_rad_prob(tinit, mesh, er_params, stat, errmsg) if (stat /= 0) return !! Defines the equation parameter components of MODEL. @@ -345,16 +348,15 @@ contains end function create_HT_model - function create_vf_rad_prob (tinit, mesh, stat, errmsg) result (vf_rad_prob) + function create_vf_rad_prob (tinit, mesh, params, stat, errmsg) result (vf_rad_prob) use rad_problem_type - use enclosure_radiation_namelist, only: params use bitfield_type, only: btest use parallel_communication, only: global_any, global_all - use parameter_list_type real(r8), intent(in) :: tinit type(unstr_mesh), intent(in) :: mesh + type(parameter_list), intent(inout) :: params integer, intent(out) :: stat character(len=*), intent(out) :: errmsg type(rad_problem), pointer :: vf_rad_prob(:) diff --git a/src/truchas/physics/heat_species_transport/HTSD_solver_factory.F90 b/src/truchas/physics/heat_species_transport/HTSD_solver_factory.F90 index 4a0217c0f..8ed2bb503 100644 --- a/src/truchas/physics/heat_species_transport/HTSD_solver_factory.F90 +++ b/src/truchas/physics/heat_species_transport/HTSD_solver_factory.F90 @@ -20,15 +20,14 @@ module HTSD_solver_factory contains - function create_HTSD_solver(mmf, model, params) result(solver) + function create_HTSD_solver(mmf, model, params, er_params) result(solver) - use enclosure_radiation_namelist, only: er_params => params use parallel_communication use truchas_env, only: output_file_name type(matl_mesh_func), intent(in), target :: mmf type(HTSD_model), intent(in), target :: model - type(parameter_list), intent(inout) :: params + type(parameter_list), intent(inout) :: params, er_params type(HTSD_solver), pointer :: solver integer :: j, n, lun diff --git a/src/truchas/physics/heat_species_transport/diffusion_solver.F90 b/src/truchas/physics/heat_species_transport/diffusion_solver.F90 index d3b22fdb2..8e591191b 100644 --- a/src/truchas/physics/heat_species_transport/diffusion_solver.F90 +++ b/src/truchas/physics/heat_species_transport/diffusion_solver.F90 @@ -110,7 +110,7 @@ contains call read_species_bc_namelists(lun, this%species_bc_params) call read_ds_source (lun) - call read_enclosure_radiation_namelists(lun) + call read_enclosure_radiation_namelists(lun, this%er_params) call this%ds_params%get('void-temperature', void_temperature) @@ -508,12 +508,12 @@ contains call tbc_fac%init(this%mesh, stefan_boltzmann, absolute_zero, this%bc_params) call sbc_fac%init(this%mesh, this%species_bc_params) call tsrc_fac%init(this%mesh, this%thermal_source_params) - this%mod1 => create_HTSD_model(tinit, this%disc, this%mmf, tbc_fac, sbc_fac, tsrc_fac, stat, errmsg) + this%mod1 => create_HTSD_model(tinit, this%disc, this%mmf, tbc_fac, sbc_fac, tsrc_fac, this%er_params, stat, errmsg) end block if (stat /= 0) call TLS_fatal ('DS_INIT: ' // trim(errmsg)) if (this%have_heat_transfer) this%ht_source => this%mod1%ht%source if (this%have_species_diffusion) this%sd_source => this%mod1%sd%source - this%sol1 => create_HTSD_solver (this%mmf, this%mod1, this%ds_params) + this%sol1 => create_HTSD_solver (this%mmf, this%mod1, this%ds_params, this%er_params) case (SOLVER2) block @@ -521,11 +521,11 @@ contains type(thermal_source_factory) :: tsrc_fac call bc_fac%init(this%mesh, stefan_boltzmann, absolute_zero, this%bc_params) call tsrc_fac%init(this%mesh, this%thermal_source_params) - this%mod2 => create_FHT_model (tinit, this%disc, this%mmf, bc_fac, tsrc_fac, stat, errmsg) + this%mod2 => create_FHT_model (tinit, this%disc, this%mmf, bc_fac, tsrc_fac, this%er_params, stat, errmsg) end block if (stat /= 0) call TLS_fatal ('DS_INIT: ' // trim(errmsg)) this%ht_source => this%mod2%q ! we need this to set the advected heat at each step - this%sol2 => create_FHT_solver(this%mmf, this%mod2, this%ds_params) + this%sol2 => create_FHT_solver(this%mmf, this%mod2, this%ds_params, this%er_params) call move_alloc(this%hadv, this%sol2%hadv) case (SOLVER3) @@ -534,11 +534,11 @@ contains type(thermal_source_factory) :: tsrc_fac call tbc_fac%init(this%mesh, stefan_boltzmann, absolute_zero, this%bc_params) call tsrc_fac%init(this%mesh, this%thermal_source_params) - this%mod3 => create_ht_model(tinit, this%disc, this%mmf, tbc_fac, tsrc_fac, stat, errmsg2) + this%mod3 => create_ht_model(tinit, this%disc, this%mmf, tbc_fac, tsrc_fac, this%er_params, stat, errmsg2) end block if (stat /= 0) call TLS_fatal('DS_INIT: ' // trim(errmsg2)) this%ht_source => this%mod3%source - this%sol3 => create_ht_solver(this%mmf, this%mod3, this%ds_params, stat, errmsg2) + this%sol3 => create_ht_solver(this%mmf, this%mod3, this%ds_params, this%er_params, stat, errmsg2) if (stat /= 0) call TLS_fatal('DS_INIT: ' // trim(errmsg2)) case default INSIST(.false.) diff --git a/src/truchas/physics/heat_species_transport/ht_ic_solver_type.F90 b/src/truchas/physics/heat_species_transport/ht_ic_solver_type.F90 index 8d15b9b1d..51b03a386 100644 --- a/src/truchas/physics/heat_species_transport/ht_ic_solver_type.F90 +++ b/src/truchas/physics/heat_species_transport/ht_ic_solver_type.F90 @@ -97,7 +97,7 @@ contains type(ht_vector), intent(inout) :: udot ! data is intent(out) target :: u, udot - integer :: j + integer :: j, n real(r8), pointer :: state(:,:) real(r8) :: dt, Tmin, Tmax @@ -149,18 +149,6 @@ contains call udot%update(-1.0_r8, u) call udot%scale(1.0_r8/dt) -!TODO !! Finite difference approx of face temp and radiosity derivatives. -!TODO if (associated(this%model%ht)) then -!TODO call HTSD_model_get_face_temp_view (this%model, f, var) -!TODO call HTSD_model_set_face_temp (this%model, var, udot) -!TODO if (associated(this%model%vf_rad_prob)) then -!TODO do n = 1, size(this%model%vf_rad_prob) -!TODO call HTSD_model_get_radiosity_view (this%model, n, f, var) -!TODO call HTSD_model_set_radiosity (this%model, n, var, udot) -!TODO end do -!TODO end if -!TODO end if - end subroutine compute_udot !! Given a cell-based field, this auxiliary subroutine produces a face-based @@ -234,21 +222,20 @@ contains call TLS_info (' Computing consistent face temperatures and radiosities ...') -!TODO !! Solve for the radiosity components given the (approx) face temperatures. -!TODO if (associated(this%vf_rad_prob)) then -!TODO do n = 1, size(this%vf_rad_prob) -!TODO call HTSD_model_get_radiosity_view (this, n, u, var) -!TODO var = 0.0_r8 -!TODO associate (faces => this%vf_rad_prob(n)%faces) -!TODO call this%vf_rad_prob(n)%solve_radiosity (t, u%tf(faces), var, stat, num_itr, error) -!TODO end associate -!TODO if (TLS_VERBOSITY >= TLS_VERB_NOISY) then -!TODO write(string,'(4x,a,i0,a,es9.2," (",i0,")")') 'radiosity[', n, ']: |r|/|b|=', error, num_itr -!TODO call TLS_info (string) -!TODO end if -!TODO if (stat /= 0) call TLS_info (' WARNING: radiosities not converged') -!TODO end do -!TODO end if + !! Solve for the radiosity components given the (approx) face temperatures. + if (associated(this%vf_rad_prob)) then + do n = 1, size(this%vf_rad_prob) + associate (q => u%encl(n)%qrad, faces => this%vf_rad_prob(n)%faces) + q = 0.0_r8 + call this%vf_rad_prob(n)%solve_radiosity(t, u%tf(faces), q, stat, num_itr, error) + end associate + if (TLS_VERBOSITY >= TLS_VERB_NOISY) then + write(string,'(4x,a,i0,a,es9.2," (",i0,")")') 'radiosity[', n, ']: |r|/|b|=', error, num_itr + call TLS_info (string) + end if + if (stat /= 0) call TLS_info (' WARNING: radiosities not converged') + end do + end if !! Initial residual of the face temperature equations. call udot%init(u) @@ -298,20 +285,19 @@ contains dT_max = global_maxval(abs(z)) -!TODO !! Solve the radiosity components given the new face temperatures. -!TODO if (associated(this%vf_rad_prob)) then -!TODO do n = 1, size(this%vf_rad_prob) -!TODO call HTSD_model_get_radiosity_view(this, n, u, var) -!TODO associate (faces => this%vf_rad_prob(n)%faces) -!TODO call this%vf_rad_prob(n)%solve_radiosity(t, Tface(faces), var, stat, num_itr, error) -!TODO end associate -!TODO if (TLS_VERBOSITY >= TLS_VERB_NOISY) then -!TODO write(string,'(4x,a,i0,a,es9.2," (",i0,")")') 'radiosity[', n, ']: |r|/|b|=', error, num_itr -!TODO call TLS_info (string) -!TODO end if -!TODO if (stat /= 0) call TLS_info (' WARNING: radiosities not converged') -!TODO end do -!TODO end if + !! Solve the radiosity components given the new face temperatures. + if (associated(this%vf_rad_prob)) then + do n = 1, size(this%vf_rad_prob) + associate (q => u%encl(n)%qrad, faces => this%vf_rad_prob(n)%faces) + call this%vf_rad_prob(n)%solve_radiosity(t, u%tf(faces), q, stat, num_itr, error) + end associate + if (TLS_VERBOSITY >= TLS_VERB_NOISY) then + write(string,'(4x,a,i0,a,es9.2," (",i0,")")') 'radiosity[', n, ']: |r|/|b|=', error, num_itr + call TLS_info (string) + end if + if (stat /= 0) call TLS_info (' WARNING: radiosities not converged') + end do + end if !! Recompute the face temp residual. call this%compute_f(t, u, udot, f) @@ -351,10 +337,10 @@ contains type(mfd_diff_matrix), intent(inout) :: matrix - integer :: j, n, n1, n2, index + integer :: n, index integer, allocatable :: more_dir_faces(:) - real(r8) :: D(this%mesh%ncell), term - real(r8), allocatable :: values(:), values2(:,:) + real(r8) :: D(this%mesh%ncell) + real(r8), allocatable :: values(:) !! Generate list of void faces; these are treated like Dirichlet BC. if (associated(this%void_face)) then diff --git a/src/truchas/physics/heat_species_transport/ht_idaesol_model_type.F90 b/src/truchas/physics/heat_species_transport/ht_idaesol_model_type.F90 index 0708e5356..8948e10b5 100644 --- a/src/truchas/physics/heat_species_transport/ht_idaesol_model_type.F90 +++ b/src/truchas/physics/heat_species_transport/ht_idaesol_model_type.F90 @@ -50,7 +50,7 @@ contains class(vector), allocatable, intent(out) :: vec type(ht_vector), allocatable :: tmp allocate(tmp) - call tmp%init(this%model%mesh) + call this%model%init_vector(tmp) call move_alloc(tmp, vec) end subroutine diff --git a/src/truchas/physics/heat_species_transport/ht_model_factory.F90 b/src/truchas/physics/heat_species_transport/ht_model_factory.F90 index 2827abf31..08f09e69f 100644 --- a/src/truchas/physics/heat_species_transport/ht_model_factory.F90 +++ b/src/truchas/physics/heat_species_transport/ht_model_factory.F90 @@ -15,6 +15,7 @@ module ht_model_factory use matl_mesh_func_type use thermal_bc_factory_class use thermal_source_factory_type + use parameter_list_type implicit none private @@ -22,7 +23,7 @@ module ht_model_factory contains - function create_ht_model(tinit, disc, mmf, bc_fac, src_fac, stat, errmsg) result(model) + function create_ht_model(tinit, disc, mmf, bc_fac, src_fac, er_params, stat, errmsg) result(model) use diffusion_solver_data, only: void_temperature use rad_problem_type @@ -32,10 +33,9 @@ contains type(matl_mesh_func), intent(in), target :: mmf class(thermal_bc_factory), intent(inout) :: bc_fac type(thermal_source_factory), intent(inout) :: src_fac + type(parameter_list), intent(inout) :: er_params integer, intent(out) :: stat character(:), allocatable , intent(out):: errmsg - - character(128) :: errmsg2 type(ht_model), pointer :: model allocate(model) @@ -43,11 +43,8 @@ contains call model%init(disc) !! Initializes the VF_RAD_PROB components of MODEL. - model%vf_rad_prob => create_vf_rad_prob(tinit, disc%mesh, stat, errmsg2) - if (stat /= 0) then - errmsg = trim(errmsg2) - return - end if + model%vf_rad_prob => create_vf_rad_prob(tinit, disc%mesh, er_params, stat, errmsg) + if (stat /= 0) return !! Defines the equation parameter components of MODEL. call define_system_parameters(disc%mesh, mmf, model, stat, errmsg) @@ -295,34 +292,31 @@ contains end function create_HT_model - function create_vf_rad_prob (tinit, mesh, stat, errmsg) result (vf_rad_prob) + function create_vf_rad_prob(tinit, mesh, params, stat, errmsg) result (vf_rad_prob) use rad_problem_type - use enclosure_radiation_namelist, only: params use bitfield_type, only: btest use parallel_communication, only: global_any, global_all - use parameter_list_type real(r8), intent(in) :: tinit type(unstr_mesh), intent(in) :: mesh + type(parameter_list), intent(inout) :: params integer, intent(out) :: stat - character(len=*), intent(out) :: errmsg + character(:), allocatable, intent(out) :: errmsg type(rad_problem), pointer :: vf_rad_prob(:) integer :: n, j logical, allocatable :: mask(:) - character(len=31), allocatable :: encl_name(:) type(parameter_list_iterator) :: piter type(parameter_list), pointer :: plist stat = 0 - errmsg = '' !! Initialize the enclosure radiation (view factor) problems, if any. piter = parameter_list_iterator(params, sublists_only=.true.) n = piter%count() if (n > 0) then - allocate(vf_rad_prob(n), encl_name(n)) + allocate(vf_rad_prob(n)) do j = 1, n plist => piter%sublist() call vf_rad_prob(j)%init (mesh, piter%name(), plist, tinit) diff --git a/src/truchas/physics/heat_species_transport/ht_model_type.F90 b/src/truchas/physics/heat_species_transport/ht_model_type.F90 index 489d0d51f..0119493cb 100644 --- a/src/truchas/physics/heat_species_transport/ht_model_type.F90 +++ b/src/truchas/physics/heat_species_transport/ht_model_type.F90 @@ -49,6 +49,7 @@ module ht_model_type type(rad_problem), pointer :: vf_rad_prob(:) => null() contains procedure :: init + procedure :: init_vector procedure :: compute_f procedure :: update_moving_vf procedure :: add_moving_vf_events @@ -63,6 +64,16 @@ contains if (associated(this%vf_rad_prob)) deallocate(this%vf_rad_prob) end subroutine + subroutine init_vector(this, vec) + class(ht_model), intent(in) :: this + type(ht_vector), intent(out) :: vec + if (associated(this%vf_rad_prob)) then + call vec%init(this%mesh, this%vf_rad_prob%size()) + else + call vec%init(this%mesh) + end if + end subroutine + subroutine init(this, disc) class(ht_model), intent(out), target :: this type(mfd_disc), intent(in), target :: disc @@ -77,7 +88,7 @@ contains call params%get('tofh-tol', eps, default=0.0_r8) call params%get('tofh-max-try', max_try, default=50) call params%get('tofh-delta', delta, default=1.0_r8) - call this%T_of_H%init (this%H_of_T, eps=eps, max_try=max_try, delta=delta) + call this%T_of_H%init(this%H_of_T, eps=eps, max_try=max_try, delta=delta) end block end subroutine @@ -93,11 +104,9 @@ contains target :: u integer :: j, n, n1, n2 - real(r8) :: term real(r8), pointer :: qrad(:) real(r8), dimension(this%mesh%ncell) :: value - real(r8), allocatable :: Tdir(:), flux(:) - integer, pointer :: faces(:) + real(r8), allocatable :: Tdir(:) logical, allocatable :: void_link(:) real(r8), pointer :: state(:,:) @@ -242,30 +251,25 @@ contains !!!! RESIDUALS OF THE ENCLOSURE RADIATION SYSTEMS !!!!!!!!!!!!!!!!!!!!!!!!!!!! -!TODO if (associated(this%vf_rad_prob)) then -!TODO do n = 1, size(this%vf_rad_prob) -!TODO call HTSD_model_get_radiosity_view (this, n, u, qrad) -!TODO faces => this%vf_rad_prob(n)%faces -!TODO !! Radiative heat flux contribution to the heat conduction face residual. -!TODO allocate(flux(size(faces))) -!TODO call this%vf_rad_prob(n)%heat_flux (t, qrad, Tface(faces), flux) -!TODO do j = 1, size(faces) -!TODO if (this%vf_rad_prob(n)%fmask(j)) & -!TODO Fface(faces(j)) = Fface(faces(j)) + this%mesh%area(faces(j)) * flux(j) -!TODO end do -!TODO deallocate(flux) -!TODO !! Residual of the algebraic radiosity system. -!TODO call HTSD_model_get_radiosity_view (this, n, f, fptr) -!TODO call this%vf_rad_prob(n)%residual (t, qrad, Tface(faces), fptr) -!TODO fptr = -fptr -!TODO end do -!TODO end if - - !TODO: is this necessary? Off-process values are not needed, but may be - ! used in dummy vector operations, and we don't want fp exceptions. - call gather_boundary(this%mesh%cell_ip, f%hc) - call gather_boundary(this%mesh%cell_ip, f%tc) - call gather_boundary(this%mesh%face_ip, f%tf) + if (associated(this%vf_rad_prob)) then + do n = 1, size(this%vf_rad_prob) + associate (q => u%encl(n)%qrad, r => f%encl(n)%qrad, faces => this%vf_rad_prob(n)%faces) + !! Radiative heat flux contribution to the heat conduction face residual. + call this%vf_rad_prob(n)%heat_flux(t, q, u%tf(faces), r) + do j = 1, size(faces) + if (this%vf_rad_prob(n)%fmask(j)) & + f%tf(faces(j)) = f%tf(faces(j)) + this%mesh%area(faces(j)) * r(j) + end do + !! Residual of the algebraic radiosity system. + call this%vf_rad_prob(n)%residual(t, q, u%tf(faces), r) + r = -r + end associate + end do + end if + + !TODO: is this necessary? Off-process values are not needed, but may be + ! used in dummy vector operations, and we don't want fp exceptions. + call f%gather_boundary() call stop_timer('ht-function') diff --git a/src/truchas/physics/heat_species_transport/ht_norm_type.F90 b/src/truchas/physics/heat_species_transport/ht_norm_type.F90 index a9e902c52..b0e3f85b8 100644 --- a/src/truchas/physics/heat_species_transport/ht_norm_type.F90 +++ b/src/truchas/physics/heat_species_transport/ht_norm_type.F90 @@ -107,11 +107,9 @@ contains type(ht_vector), intent(in) :: u, du real(r8), intent(out) :: du_norm - !integer :: n - !real(r8) :: qerror - !real(r8), pointer :: qrad(:) - !integer, pointer :: faces(:) - !real(r8), allocatable :: res(:), rhs(:) + integer :: n + real(r8) :: qerror + real(r8), allocatable :: res(:) integer :: ncell, nface ncell = this%model%mesh%ncell_onP @@ -128,22 +126,22 @@ contains !TODO! integrator thinks it's computing the norm of the correction du but !TODO! in this case it is the actual residual norm. We need more general !TODO! norms in the integrator. -!TODO if (associated(this%model%ht%vf_rad_prob)) then -!TODO !if (is_IOP) write(*,'(a)',advance='no')'ER error: ||res||/||rhs||=' -!TODO do n = 1, size(this%model%ht%vf_rad_prob) -!TODO faces => this%model%ht%vf_rad_prob(n)%faces -!TODO allocate(res(size(faces)), rhs(size(faces))) -!TODO call HTSD_model_get_radiosity_view (this%model, n, u, qrad) -!TODO call HTSD_model_get_face_temp_view (this%model, u, temp) -!TODO call this%model%ht%vf_rad_prob(n)%residual (t, qrad, temp(faces), res) -!TODO call this%model%ht%vf_rad_prob(n)%rhs (t, temp(faces), rhs) -!TODO qerror = sqrt(global_sum(res**2)) / sqrt(global_sum(rhs**2)) -!TODO !if (is_IOP) write(*,'(e10.3)',advance='no') qerror -!TODO ht_du_norm = max(ht_du_norm, qerror/1.0d-3) -!TODO deallocate(res, rhs) -!TODO end do -!TODO !if (is_IOP) write(*,*) -!TODO end if + if (associated(this%model%vf_rad_prob)) then + !if (is_IOP) write(*,'(a)',advance='no')'ER error: ||res||/||rhs||=' + do n = 1, size(this%model%vf_rad_prob) + associate (q => u%encl(n)%qrad, faces => this%model%vf_rad_prob(n)%faces) + allocate(res(size(faces))) + call this%model%vf_rad_prob(n)%residual(t, q, u%tf(faces), res) + qerror = sqrt(global_sum(norm2(res)**2)) + call this%model%vf_rad_prob(n)%rhs(t, u%tf(faces), res) + qerror = qerror / sqrt(global_sum(norm2(res)**2)) + !if (is_IOP) write(*,'(e10.3)',advance='no') qerror + du_norm = max(du_norm, qerror/1.0d-3) + deallocate(res) + end associate + end do + !if (is_IOP) write(*,*) + end if contains diff --git a/src/truchas/physics/heat_species_transport/ht_precon_type.F90 b/src/truchas/physics/heat_species_transport/ht_precon_type.F90 index 0dc4d563c..5da930029 100644 --- a/src/truchas/physics/heat_species_transport/ht_precon_type.F90 +++ b/src/truchas/physics/heat_species_transport/ht_precon_type.F90 @@ -103,7 +103,6 @@ contains real(r8) :: D(this%mesh%ncell), A(this%mesh%ncell), term real(r8), allocatable :: values(:), values2(:,:) type(mfd_diff_matrix), pointer :: dm - integer, pointer :: faces(:) => null() ASSERT(dt > 0.0_r8) @@ -221,12 +220,13 @@ contains !! TODO: what about factorization coupling? Is this still correct? if (associated(this%model%vf_rad_prob)) then do index = 1, size(this%model%vf_rad_prob) - faces => this%model%vf_rad_prob(index)%faces - allocate(values(size(faces))) - call this%model%vf_rad_prob(index)%rhs_deriv(t, u%tf(faces), values) - where (.not.this%model%vf_rad_prob(index)%fmask) values = 0 - call dm%incr_face_diag(faces, this%mesh%area(faces) * values) - deallocate(values) + associate (faces => this%model%vf_rad_prob(index)%faces) + allocate(values(size(faces))) + call this%model%vf_rad_prob(index)%rhs_deriv(t, u%tf(faces), values) + where (.not.this%model%vf_rad_prob(index)%fmask) values = 0 + call dm%incr_face_diag(faces, this%mesh%area(faces) * values) + deallocate(values) + end associate end do end if @@ -246,7 +246,6 @@ contains type(ht_vector), intent(inout) :: f ! data is intent(inout) integer :: index, j, n - real(r8), pointer :: fq(:) real(r8), allocatable :: z(:) call start_timer('ht-precon-apply') @@ -259,11 +258,9 @@ contains do index = 1, size(this%model%vf_rad_prob) if (this%vfr_precon_coupling(index) == VFR_FGS .or. & this%vfr_precon_coupling(index) == VFR_FAC) then - !TODO: call HTSD_model_get_radiosity_view (this%model, index, f, fq) - allocate(z(size(fq))) - z = fq + z = f%encl(index)%qrad call this%model%vf_rad_prob(index)%precon(t, z) - if (this%vfr_precon_coupling(index) == VFR_FGS) fq = z + if (this%vfr_precon_coupling(index) == VFR_FGS) f%encl(index)%qrad = z !! Update the heat equation face residual. call this%model%vf_rad_prob(index)%precon_matvec1(t, z) do j = 1, size(z) @@ -316,15 +313,17 @@ contains if (this%vfr_precon_coupling(index) == VFR_JAC .or. & this%vfr_precon_coupling(index) == VFR_BGS .or. & this%vfr_precon_coupling(index) == VFR_FAC) then - !! Update the radiosity residual components. - !TODO: call HTSD_model_get_radiosity_view (this%model, index, f, fq) - if (this%vfr_precon_coupling(index) /= VFR_JAC) then - allocate(z(size(fq))) - call this%model%vf_rad_prob(index)%rhs_deriv(t, u%tf(this%model%vf_rad_prob(index)%faces), z) - fq = fq + z * f%tf(this%model%vf_rad_prob(index)%faces) - deallocate(z) - end if - call this%model%vf_rad_prob(index)%precon(t, fq) + associate (fq => f%encl(index)%qrad, faces => this%model%vf_rad_prob(index)%faces) + !! Update the radiosity residual components. + !TODO: call HTSD_model_get_radiosity_view (this%model, index, f, fq) + if (this%vfr_precon_coupling(index) /= VFR_JAC) then + allocate(z(size(fq))) + call this%model%vf_rad_prob(index)%rhs_deriv(t, u%tf(faces), z) + fq = fq + z * f%tf(faces) + deallocate(z) + end if + call this%model%vf_rad_prob(index)%precon(t, fq) + end associate end if end do call stop_timer('vf-rad-precon') diff --git a/src/truchas/physics/heat_species_transport/ht_solver_factory.F90 b/src/truchas/physics/heat_species_transport/ht_solver_factory.F90 index fb8da81bb..f0d9340fc 100644 --- a/src/truchas/physics/heat_species_transport/ht_solver_factory.F90 +++ b/src/truchas/physics/heat_species_transport/ht_solver_factory.F90 @@ -20,15 +20,14 @@ module ht_solver_factory contains - function create_ht_solver(mmf, model, params, stat, errmsg) result(solver) + function create_ht_solver(mmf, model, params, er_params, stat, errmsg) result(solver) - use enclosure_radiation_namelist, only: er_params => params use parallel_communication use truchas_env, only: output_file_name type(matl_mesh_func), intent(in), target :: mmf type(ht_model), intent(in), target :: model - type(parameter_list), intent(inout) :: params + type(parameter_list), intent(inout) :: params, er_params type(ht_solver), pointer :: solver integer, intent(out) :: stat character(:), allocatable, intent(out) :: errmsg @@ -50,25 +49,23 @@ contains call plist%set('unit', lun) end if -!TODO !TODO: read enclosure radiation namelists directly into sublist -!TODO if (associated(model%ht)) then -!TODO if (associated(model%ht%vf_rad_prob)) then -!TODO !TODO! Assumes model%ht%vf_rad_prob array was initialized under the same -!TODO !TODO! loop so that that array and vfr_precon_coupling are in correspondence. -!TODO !TODO! This is fragile and needs to be fixed. -!TODO piter = parameter_list_iterator(er_params) -!TODO n = piter%count() -!TODO allocate(vfr_precon_coupling(n), rad_tol(n)) -!TODO do j = 1, n -!TODO plist => piter%sublist() -!TODO call plist%get('precon-coupling-method', string, default='BACKWARD GS') -!TODO vfr_precon_coupling(j) = string -!TODO call piter%next -!TODO end do -!TODO plist => params%sublist('precon') -!TODO call plist%set('vfr-precon-coupling', vfr_precon_coupling) -!TODO end if -!TODO end if + !TODO: read enclosure radiation namelists directly into sublist + if (associated(model%vf_rad_prob)) then + !TODO! Assumes model%vf_rad_prob array was initialized under the same + !TODO! loop so that that array and vfr_precon_coupling are in correspondence. + !TODO! This is fragile and needs to be fixed. + piter = parameter_list_iterator(er_params) + n = piter%count() + allocate(vfr_precon_coupling(n), rad_tol(n)) + do j = 1, n + plist => piter%sublist() + call plist%get('precon-coupling-method', string, default='BACKWARD GS') + vfr_precon_coupling(j) = string + call piter%next + end do + plist => params%sublist('precon') + call plist%set('vfr-precon-coupling', vfr_precon_coupling) + end if allocate(solver) call solver%init(mmf, model, params, stat, errmsg) diff --git a/src/truchas/physics/heat_species_transport/ht_solver_type.F90 b/src/truchas/physics/heat_species_transport/ht_solver_type.F90 index ca9a50fe0..c24e7f5c8 100644 --- a/src/truchas/physics/heat_species_transport/ht_solver_type.F90 +++ b/src/truchas/physics/heat_species_transport/ht_solver_type.F90 @@ -97,7 +97,11 @@ contains call this%precon%init(model, plist, stat, errmsg) if (stat /= 0) return - call this%u%init(this%mesh) + if (associated(this%model%vf_rad_prob)) then + call this%u%init(this%mesh, this%model%vf_rad_prob%size()) + else + call this%u%init(this%mesh) + end if call this%integ_model%init(this%model, this%precon, this%norm) diff --git a/src/truchas/physics/heat_species_transport/ht_vector_type.F90 b/src/truchas/physics/heat_species_transport/ht_vector_type.F90 index fb8842817..233d18b5a 100644 --- a/src/truchas/physics/heat_species_transport/ht_vector_type.F90 +++ b/src/truchas/physics/heat_species_transport/ht_vector_type.F90 @@ -8,9 +8,14 @@ module ht_vector_type implicit none private + type :: encl_vector + real(r8), allocatable :: qrad(:) + end type + type, extends(vector), public :: ht_vector type(unstr_mesh), pointer :: mesh => null() real(r8), allocatable :: hc(:), tc(:), tf(:) + type(encl_vector), allocatable :: encl(:) contains !! Deferred base class procedures procedure :: clone1 @@ -36,12 +41,20 @@ contains !! Specific subroutine for the generic INIT. Initialize a HT_VECTOR object !! for the given unstructured MESH. The elements are initialized to 0. - subroutine init_mesh(this, mesh) + subroutine init_mesh(this, mesh, encl_size) class(ht_vector), intent(out) :: this type(unstr_mesh), intent(in), target :: mesh + integer, intent(in), optional :: encl_size(:) + integer :: n this%mesh => mesh allocate(this%hc(mesh%ncell), this%tc(mesh%ncell), this%tf(mesh%nface)) - call setval(this, 0.0_r8) + if (present(encl_size)) then + allocate(this%encl(size(encl_size))) + do n = 1, size(this%encl) + allocate(this%encl(n)%qrad(encl_size(n))) + end do + end if + call this%setval(0.0_r8) end subroutine !! Specific subroutine for the generic INIT. Initialize a HT_VECTOR object @@ -50,7 +63,12 @@ contains subroutine init_mold(this, mold) class(ht_vector), intent(out) :: this class(ht_vector), intent(in) :: mold - call init_mesh(this, mold%mesh) + integer :: n + if (allocated(mold%encl)) then + call init_mesh(this, mold%mesh, [(size(mold%encl(n)%qrad), n=1, size(mold%encl))]) + else + call init_mesh(this, mold%mesh) + end if end subroutine subroutine ht_gather_boundary(this) @@ -64,13 +82,7 @@ contains subroutine clone1(this, clone) class(ht_vector), intent(in) :: this class(vector), allocatable, intent(out) :: clone - type(ht_vector), allocatable :: tmp - allocate(tmp) - tmp%mesh => this%mesh - allocate(tmp%hc, mold=this%hc) - allocate(tmp%tc, mold=this%tc) - allocate(tmp%tf, mold=this%tf) - call move_alloc(tmp, clone) + allocate(clone, source=this) ! easy, but with unwanted copy of data too end subroutine subroutine clone2(this, clone, n) @@ -78,42 +90,52 @@ contains class(vector), allocatable, intent(out) :: clone(:) type(ht_vector), allocatable :: tmp(:) integer, intent(in) :: n - integer :: j - allocate(tmp(n)) - do j = 1, n - tmp(j)%mesh => this%mesh - allocate(tmp(j)%hc, mold=this%hc) - allocate(tmp(j)%tc, mold=this%tc) - allocate(tmp(j)%tf, mold=this%tf) - end do - call move_alloc(tmp, clone) + allocate(clone(n), source=this) ! easy, but with unwanted copy of data too end subroutine subroutine copy_(dest, src) class(ht_vector), intent(inout) :: dest class(vector), intent(in) :: src + integer :: n select type (src) class is (ht_vector) dest%hc(:) = src%hc dest%tc(:) = src%tc dest%tf(:) = src%tf + if (allocated(dest%encl)) then + do n = 1, size(dest%encl) + dest%encl(n)%qrad(:) = src%encl(n)%qrad + end do + end if end select end subroutine subroutine setval(this, val) class(ht_vector), intent(inout) :: this real(r8), intent(in) :: val + integer :: n this%hc = val this%tc = val this%tf = val + if (allocated(this%encl)) then + do n = 1, size(this%encl) + this%encl(n)%qrad = val + end do + end if end subroutine subroutine scale(this, a) class(ht_vector), intent(inout) :: this real(r8), intent(in) :: a + integer :: n this%hc = a * this%hc this%tc = a * this%tc this%tf = a * this%tf + if (allocated(this%encl)) then + do n = 1, size(this%encl) + this%encl(n)%qrad = a * this%encl(n)%qrad + end do + end if end subroutine !! Conventional SAXPY procedure: y <-- a*x + y @@ -121,11 +143,17 @@ contains class(ht_vector), intent(inout) :: this class(vector), intent(in) :: x real(r8), intent(in) :: a + integer :: n select type (x) class is (ht_vector) this%hc = a * x%hc + this%hc this%tc = a * x%tc + this%tc this%tf = a * x%tf + this%tf + if (allocated(this%encl)) then + do n = 1, size(this%encl) + this%encl(n)%qrad = a * x%encl(n)%qrad + this%encl(n)%qrad + end do + end if end select end subroutine @@ -134,11 +162,17 @@ contains class(ht_vector), intent(inout) :: this class(vector), intent(in) :: x real(r8), intent(in) :: a, b + integer :: n select type (x) class is (ht_vector) this%hc = a * x%hc + b * this%hc this%tc = a * x%tc + b * this%tc this%tf = a * x%tf + b * this%tf + if (allocated(this%encl)) then + do n = 1, size(this%encl) + this%encl(n)%qrad = a * x%encl(n)%qrad + b * this%encl(n)%qrad + end do + end if end select end subroutine @@ -147,6 +181,7 @@ contains class(ht_vector), intent(inout) :: this class(vector), intent(in) :: x, y real(r8), intent(in) :: a, b + integer :: n select type (x) class is (ht_vector) select type (y) @@ -154,6 +189,11 @@ contains this%hc = a * x%hc + b * y%hc + this%hc this%tc = a * x%tc + b * y%tc + this%tc this%tf = a * x%tf + b * y%tf + this%tf + if (allocated(this%encl)) then + do n = 1, size(this%encl) + this%encl(n)%qrad = a * x%encl(n)%qrad + b * y%encl(n)%qrad + this%encl(n)%qrad + end do + end if end select end select end subroutine @@ -163,6 +203,7 @@ contains class(ht_vector), intent(inout) :: this class(vector), intent(in) :: x, y real(r8), intent(in) :: a, b, c + integer :: n select type (x) class is (ht_vector) select type (y) @@ -170,6 +211,11 @@ contains this%hc = a * x%hc + b * y%hc + c * this%hc this%tc = a * x%tc + b * y%tc + c * this%tc this%tf = a * x%tf + b * y%tf + c * this%tf + if (allocated(this%encl)) then + do n = 1, size(this%encl) + this%encl(n)%qrad = a * x%encl(n)%qrad + b * y%encl(n)%qrad + c * this%encl(n)%qrad + end do + end if end select end select end subroutine @@ -178,7 +224,7 @@ contains class(ht_vector), intent(in) :: x class(vector), intent(in) :: y real(r8) :: dp - integer :: j + integer :: j, n select type (y) class is (ht_vector) dp = 0 @@ -191,6 +237,13 @@ contains do j = 1, x%mesh%nface_onP dp = dp + x%tf(j) * y%tf(j) end do + if (allocated(x%encl)) then + do n = 1, size(x%encl) + do j = 1, size(x%encl(n)%qrad) + dp = dp + x%encl(n)%qrad(j) * y%encl(n)%qrad(j) + end do + end do + end if dp = global_sum(dp) end select end function @@ -199,9 +252,15 @@ contains use parallel_communication, only: global_sum class(ht_vector), intent(in) :: this real(r8) :: norm2_ + integer :: n norm2_ = norm2(this%hc(:this%mesh%ncell_onP))**2 + & norm2(this%tc(:this%mesh%ncell_onP))**2 + & norm2(this%tf(:this%mesh%nface_onP))**2 + if (allocated(this%encl)) then + do n = 1, size(this%encl) + norm2_ = norm2_ + norm2(this%encl(n)%qrad)**2 + end do + end if norm2_ = sqrt(global_sum(norm2_)) end function @@ -210,6 +269,7 @@ contains class(ht_vector), intent(in) :: this logical, intent(in), optional :: full ! default is FALSE character(:), allocatable :: string + integer :: n type(md5_hash) :: hash logical :: strict strict = .true. @@ -223,6 +283,11 @@ contains call hash%update(this%tc(:this%mesh%ncell)) call hash%update(this%tf(:this%mesh%nface)) end if + if (allocated(this%encl)) then + do n = 1, size(this%encl) + call hash%update(this%encl(n)%qrad) + end do + end if string = hash%hexdigest() end function diff --git a/test/vfrad-moving/test.py b/test/vfrad-moving/test.py index 517baf4c4..bea49bd88 100755 --- a/test/vfrad-moving/test.py +++ b/test/vfrad-moving/test.py @@ -11,7 +11,7 @@ def run_test(tenv): dataA = numpy.loadtxt(filename) filename = os.path.join(tenv._input_dir, "golden.dat") dataB = numpy.loadtxt(filename) - nfail += truchas.compare_max_rel(dataA[:,1], dataB[:,1], 5e-6, 'temperature', -1.0) + nfail += truchas.compare_max_rel(dataA[:,1], dataB[:,1], 7e-5, 'temperature', -1.0) return nfail if __name__=="__main__": -- GitLab From e14574fe56f945ef90ac25e3c5d38487e9c5b0db Mon Sep 17 00:00:00 2001 From: "Neil N. Carlson" <neil.n.carlson@gmail.com> Date: Wed, 11 May 2022 15:29:24 -0600 Subject: [PATCH 3/6] Updates required by rebasing --- .../enclosure_radiation_namelist.F90 | 2 +- .../ht_ic_solver_type.F90 | 7 +++---- .../ht_model_factory.F90 | 21 ++++++++++++++----- .../heat_species_transport/ht_model_type.F90 | 21 ++++++++++++++----- .../heat_species_transport/ht_precon_type.F90 | 11 +++++----- .../heat_species_transport/ht_solver_type.F90 | 7 +++---- .../heat_species_transport/ht_vector_type.F90 | 11 +++++----- 7 files changed, 49 insertions(+), 31 deletions(-) diff --git a/src/truchas/physics/enclosure_radiation/enclosure_radiation_namelist.F90 b/src/truchas/physics/enclosure_radiation/enclosure_radiation_namelist.F90 index 28807ffdb..9f2e5985e 100644 --- a/src/truchas/physics/enclosure_radiation/enclosure_radiation_namelist.F90 +++ b/src/truchas/physics/enclosure_radiation/enclosure_radiation_namelist.F90 @@ -190,7 +190,7 @@ contains end do if (n > 0) then - call read_enclosure_surface_namelists(lun) + call read_enclosure_surface_namelists(lun, params) else call TLS_info(' none found') end if diff --git a/src/truchas/physics/heat_species_transport/ht_ic_solver_type.F90 b/src/truchas/physics/heat_species_transport/ht_ic_solver_type.F90 index 51b03a386..7592032c8 100644 --- a/src/truchas/physics/heat_species_transport/ht_ic_solver_type.F90 +++ b/src/truchas/physics/heat_species_transport/ht_ic_solver_type.F90 @@ -10,7 +10,6 @@ module ht_ic_solver_type use,intrinsic :: iso_fortran_env, only: r8 => real64 use unstr_mesh_type - use index_partitioning use ht_vector_type use ht_model_type use truchas_logging_services @@ -65,7 +64,7 @@ contains u%tc(:this%mesh%ncell_onP) = temp(:this%mesh%ncell_onP) if (associated(this%model%void_cell)) & where (this%model%void_cell) u%tc = this%model%void_temp - call gather_boundary(this%mesh%cell_ip, u%tc) + call this%mesh%cell_imap%gather_offp(u%tc) !! Solve for consistent face temperatures and radiosities. !! Averaging adjacent cell temps provides a cheap initial guess. @@ -180,8 +179,8 @@ contains scale(cface) = scale(cface) + 1 end associate end do - call gather_boundary(mesh%face_ip, uface) - call gather_boundary(mesh%face_ip, scale) + call mesh%face_imap%gather_offp(uface) + call mesh%face_imap%gather_offp(scale) where (scale == 0) uface = 0.0_r8 diff --git a/src/truchas/physics/heat_species_transport/ht_model_factory.F90 b/src/truchas/physics/heat_species_transport/ht_model_factory.F90 index 08f09e69f..ed6ca725d 100644 --- a/src/truchas/physics/heat_species_transport/ht_model_factory.F90 +++ b/src/truchas/physics/heat_species_transport/ht_model_factory.F90 @@ -16,6 +16,7 @@ module ht_model_factory use thermal_bc_factory_class use thermal_source_factory_type use parameter_list_type + use truchas_logging_services implicit none private @@ -115,7 +116,6 @@ contains use evaporation_namelist, only: evap_params => params use evap_heat_flux_type - use index_partitioning, only: gather_boundary use bitfield_type use parallel_communication, only: global_all, global_any, global_count use string_utilities, only: i_to_c @@ -143,7 +143,7 @@ contains if (allocated(model%ic_htc)) then mask(model%ic_htc%index(1,:)) = .true. mask(model%ic_htc%index(2,:)) = .true. - call gather_boundary(mesh%face_ip, mask) + call mesh%face_imap%gather_offp(mask) end if !! Define the gap radiation interface conditions; @@ -153,7 +153,7 @@ contains if (allocated(model%ic_rad)) then mask(model%ic_rad%index(1,:)) = .true. mask(model%ic_rad%index(2,:)) = .true. - call gather_boundary(mesh%face_ip, mask) + call mesh%face_imap%gather_offp(mask) end if !! Flux-type boundary conditions. These may be superimposed. @@ -172,6 +172,17 @@ contains fmask(model%bc_flux%index) = .true. ! mark the simple flux faces end if + !! Define the oriented flux boundary conditions. + call bc_fac%alloc_vflux_bc(model%bc_vflux, stat, errmsg) + if (stat /= 0) return + if (allocated(model%bc_vflux)) then + if (global_any(mask(model%bc_vflux%index))) & + call TLS_info(' NOTE: oriented-flux condition is superimposed with interface conditions') + do j = 1, size(model%bc_vflux%index) ! index not necessarily 1-1 + fmask(model%bc_vflux%index(j)) = .true. ! mark the oriented flux faces + end do + end if + !! Define the external HTC boundary conditions. call bc_fac%alloc_htc_bc(model%bc_htc, stat, errmsg) if (stat /= 0) return @@ -193,8 +204,8 @@ contains rmask(model%vf_rad_prob(j)%faces) = .true. !fmask(model%vf_rad_prob(j)%faces) = .true. end do - call gather_boundary (mesh%face_ip, rmask) - !call gather_boundary (mesh%face_ip, fmask) + call mesh%face_imap%gather_offp(rmask) + !call mesh%face_imap%gather_offp(fmask) end if !! Define the (simple) radiation boundary conditions. diff --git a/src/truchas/physics/heat_species_transport/ht_model_type.F90 b/src/truchas/physics/heat_species_transport/ht_model_type.F90 index 0119493cb..4a5b789cf 100644 --- a/src/truchas/physics/heat_species_transport/ht_model_type.F90 +++ b/src/truchas/physics/heat_species_transport/ht_model_type.F90 @@ -15,11 +15,12 @@ module ht_model_type use prop_mesh_func_type use source_mesh_function use scalar_mesh_func_class + use bndry_vfunc_class use bndry_func1_class use bndry_func2_class use intfc_func2_class use rad_problem_type - use index_partitioning +! use index_partitioning use TofH_type use truchas_timers use parameter_list_type @@ -40,6 +41,7 @@ module ht_model_type !! Boundary condition data class(bndry_func1), allocatable :: bc_dir ! Dirichlet class(bndry_func1), allocatable :: bc_flux ! simple flux + class(bndry_vfunc), allocatable :: bc_vflux ! oriented flux class(bndry_func2), allocatable :: bc_htc ! external HTC class(bndry_func2), allocatable :: bc_rad ! simple radiation class(intfc_func2), allocatable :: ic_htc ! internal HTC @@ -112,8 +114,8 @@ contains call start_timer('ht-function') - call gather_boundary(this%mesh%cell_ip, u%tc) - call gather_boundary(this%mesh%face_ip, u%tf) + call this%mesh%cell_imap%gather_offp(u%tc) + call this%mesh%face_imap%gather_offp(u%tf) !TODO: The existing prop_mesh_func%compute_value function expects a rank-2 ! state array. This is a workaround until prop_mesh_func is redesigned. @@ -147,7 +149,7 @@ contains if (associated(this%void_cell)) where (this%void_cell) value = 0.0_r8 call this%disc%apply_diff(value, u%tc, u%tf, f%tc, f%tf) call smf_eval(this%source, t, value) - call gather_boundary(this%mesh%cell_ip, udot%hc) + call this%mesh%cell_imap%gather_offp(udot%hc) f%tc = f%tc + this%mesh%volume*(udot%hc - value) !! Additional heat source @@ -176,6 +178,15 @@ contains end do end if + !! Oriented flux BC contribution. + if (allocated(this%bc_vflux)) then + call this%bc_vflux%compute(t) + do j = 1, size(this%bc_vflux%index) + n = this%bc_vflux%index(j) + f%tf(n) = f%tf(n) + dot_product(this%mesh%normal(:,n), this%bc_vflux%value(:,j)) + end do + end if + !! External HTC flux contribution. if (allocated(this%bc_htc)) then call this%bc_htc%compute(t, u%tf) @@ -269,7 +280,7 @@ contains !TODO: is this necessary? Off-process values are not needed, but may be ! used in dummy vector operations, and we don't want fp exceptions. - call f%gather_boundary() + call f%gather_offp call stop_timer('ht-function') diff --git a/src/truchas/physics/heat_species_transport/ht_precon_type.F90 b/src/truchas/physics/heat_species_transport/ht_precon_type.F90 index 5da930029..92e2d1f92 100644 --- a/src/truchas/physics/heat_species_transport/ht_precon_type.F90 +++ b/src/truchas/physics/heat_species_transport/ht_precon_type.F90 @@ -15,7 +15,6 @@ module ht_precon_type use unstr_mesh_type use mfd_diff_precon_type use mfd_diff_matrix_type - use index_partitioning use truchas_timers implicit none private @@ -127,8 +126,8 @@ contains ! state array. This is a workaround until prop_mesh_func is redesigned. state(1:this%mesh%ncell,1:1) => u%tc - call gather_boundary(this%mesh%cell_ip, u%tc) - call gather_boundary(this%mesh%face_ip, u%tf) + call this%mesh%cell_imap%gather_offp(u%tc) + call this%mesh%face_imap%gather_offp(u%tf) this%dt = dt @@ -286,12 +285,12 @@ contains else f%tc = f%tc - (this%mesh%volume/this%dt)*f%hc end if - call gather_boundary(this%mesh%cell_ip, f%tc) + call this%mesh%cell_imap%gather_offp(f%tc) !! Heat equation face residual (with radiosity residuals optionally eliminated). !call HTSD_model_get_face_temp_copy (this%model, f, f2x) - !call gather_boundary (this%mesh%face_ip, f2x) - call gather_boundary(this%mesh%face_ip, f%tf) + !this%mesh%face_imap%gather_offp(f2x) + call this%mesh%face_imap%gather_offp(f%tf) !! Precondition the heat equation. !call diff_precon_apply (this%hcprecon, f1x, f2x) diff --git a/src/truchas/physics/heat_species_transport/ht_solver_type.F90 b/src/truchas/physics/heat_species_transport/ht_solver_type.F90 index c24e7f5c8..0bc6b2560 100644 --- a/src/truchas/physics/heat_species_transport/ht_solver_type.F90 +++ b/src/truchas/physics/heat_species_transport/ht_solver_type.F90 @@ -15,7 +15,6 @@ module ht_solver_type use ht_norm_type use matl_mesh_func_type use unstr_mesh_type - use index_partitioning use ht_idaesol_model_type use new_idaesol_type use parameter_list_type @@ -193,7 +192,7 @@ contains real(r8), intent(out) :: tgrad(:,:) INSIST(size(tgrad,1) == 3) INSIST(size(tgrad,2) == this%model%mesh%ncell_onP) - call gather_boundary(this%model%mesh%face_ip, this%u%tf) + call this%model%mesh%face_imap%gather_offp(this%u%tf) call this%model%disc%compute_cell_grad(this%u%tf, tgrad) end subroutine @@ -238,7 +237,7 @@ contains if (stat == 0) then this%t = t this%state_is_pending = .true. - call this%u%gather_boundary() !TODO: Can the be made unnecessary? + call this%u%gather_offp !TODO: Can the be made unnecessary? else call this%integ%get_last_state_copy(this%u) this%state_is_pending = .false. @@ -286,7 +285,7 @@ contains end associate end if end do - call gather_boundary(this%mesh%face_ip, this%model%void_face) + call this%mesh%face_imap%gather_offp(this%model%void_face) else deallocate(this%model%void_cell) this%model%void_face => null() diff --git a/src/truchas/physics/heat_species_transport/ht_vector_type.F90 b/src/truchas/physics/heat_species_transport/ht_vector_type.F90 index 233d18b5a..1987dde4e 100644 --- a/src/truchas/physics/heat_species_transport/ht_vector_type.F90 +++ b/src/truchas/physics/heat_species_transport/ht_vector_type.F90 @@ -33,7 +33,7 @@ module ht_vector_type !! Additional procedures specific to this type generic :: init => init_mesh, init_mold procedure, private :: init_mesh, init_mold - procedure :: gather_boundary => ht_gather_boundary + procedure :: gather_offp end type contains @@ -71,12 +71,11 @@ contains end if end subroutine - subroutine ht_gather_boundary(this) - use index_partitioning, only: gather_boundary + subroutine gather_offp(this) class(ht_vector), intent(inout) :: this - call gather_boundary(this%mesh%cell_ip, this%hc) - call gather_boundary(this%mesh%cell_ip, this%tc) - call gather_boundary(this%mesh%face_ip, this%tf) + call this%mesh%cell_imap%gather_offp(this%hc) + call this%mesh%cell_imap%gather_offp(this%tc) + call this%mesh%face_imap%gather_offp(this%tf) end subroutine subroutine clone1(this, clone) -- GitLab From 0cac0a246a6b09d960f4b784a506d78eedfa5cb6 Mon Sep 17 00:00:00 2001 From: "Neil N. Carlson" <neil.n.carlson@gmail.com> Date: Fri, 10 Nov 2023 06:59:51 -0700 Subject: [PATCH 4/6] Updates required to rebase branch on 22.11 --- .../enclosure_radiation/Test/test_rade.F90 | 7 ++++--- .../enclosure_radiation_namelist.F90 | 10 +++++----- .../heat_species_transport/FHT_model_factory.F90 | 11 +++++------ .../FHT_solver_factory.F90 | 4 ++-- .../HTSD_model_factory.F90 | 16 +++++++--------- .../HTSD_solver_factory.F90 | 5 +++-- .../heat_species_transport/diffusion_solver.F90 | 14 +++++++------- .../heat_species_transport/ht_model_factory.F90 | 5 ++--- .../heat_species_transport/ht_model_type.F90 | 7 +++---- .../heat_species_transport/ht_solver_factory.F90 | 5 +++-- 10 files changed, 41 insertions(+), 43 deletions(-) diff --git a/src/truchas/physics/enclosure_radiation/Test/test_rade.F90 b/src/truchas/physics/enclosure_radiation/Test/test_rade.F90 index 2239f5761..ab51d37de 100644 --- a/src/truchas/physics/enclosure_radiation/Test/test_rade.F90 +++ b/src/truchas/physics/enclosure_radiation/Test/test_rade.F90 @@ -14,12 +14,10 @@ program test_rade use truchas_logging_services use unstr_mesh_type use rad_problem_type - use parameter_list_type implicit none type(unstr_mesh), pointer :: mesh real(r8), allocatable :: hc_temp(:) - type(parameter_list) :: params !! Enclosures to be tested character(len=31), parameter :: ENCL_EXT_FACE = "exterior_face" @@ -83,7 +81,7 @@ contains call read_function_namelists (lun) !! Enclosure radiation namelists. - call read_enclosure_radiation_namelists(lun, params) + call read_enclosure_radiation_namelists(lun) !! The mesh must be named 'main' call enable_mesh ('main', exists) @@ -120,6 +118,9 @@ contains !! Test that radiosity solver converges within given error tolerance subroutine test_converge (mesh, encl_name, hc_temp) + use enclosure_radiation_namelist, only: params + use parameter_list_type + type(unstr_mesh), pointer, intent(in) :: mesh character(len=31), intent(in) :: encl_name real(r8), intent(in) :: hc_temp(:) diff --git a/src/truchas/physics/enclosure_radiation/enclosure_radiation_namelist.F90 b/src/truchas/physics/enclosure_radiation/enclosure_radiation_namelist.F90 index 9f2e5985e..dd285a261 100644 --- a/src/truchas/physics/enclosure_radiation/enclosure_radiation_namelist.F90 +++ b/src/truchas/physics/enclosure_radiation/enclosure_radiation_namelist.F90 @@ -26,14 +26,15 @@ module enclosure_radiation_namelist public :: read_enclosure_radiation_namelists + type(parameter_list), public :: params + integer, parameter :: MAX_NAME_LEN = 31, MAX_FILE_LEN = 255, MAX_FACE_BLOCK_IDS = 32 contains - subroutine read_enclosure_radiation_namelists(lun, params) + subroutine read_enclosure_radiation_namelists(lun) integer, intent(in) :: lun - type(parameter_list), intent(inout) :: params integer :: n, ios logical :: found @@ -190,7 +191,7 @@ contains end do if (n > 0) then - call read_enclosure_surface_namelists(lun, params) + call read_enclosure_surface_namelists(lun) else call TLS_info(' none found') end if @@ -198,10 +199,9 @@ contains end subroutine read_enclosure_radiation_namelists - subroutine read_enclosure_surface_namelists(lun, params) + subroutine read_enclosure_surface_namelists(lun) integer, intent(in) :: lun - type(parameter_list), intent(inout) :: params integer :: n, ios logical :: found diff --git a/src/truchas/physics/heat_species_transport/FHT_model_factory.F90 b/src/truchas/physics/heat_species_transport/FHT_model_factory.F90 index 56203c07d..ad90757a2 100644 --- a/src/truchas/physics/heat_species_transport/FHT_model_factory.F90 +++ b/src/truchas/physics/heat_species_transport/FHT_model_factory.F90 @@ -34,7 +34,6 @@ module FHT_model_factory use rad_problem_type use thermal_bc_factory_class use thermal_source_factory_type - use parameter_list_type implicit none private @@ -42,7 +41,7 @@ module FHT_model_factory contains - function create_FHT_model (tinit, disc, mmf, tbc_fac, tsrc_fac, er_params, stat, errmsg) result (model) + function create_FHT_model (tinit, disc, mmf, tbc_fac, tsrc_fac, stat, errmsg) result (model) use diffusion_solver_data, only: void_temperature @@ -51,7 +50,6 @@ contains type(matl_mesh_func), intent(in), target :: mmf class(thermal_bc_factory), intent(inout) :: tbc_fac type(thermal_source_factory), intent(inout) :: tsrc_fac - type(parameter_list), intent(inout) :: er_params integer, intent(out) :: stat character(*), intent(out) :: errmsg character(:), allocatable :: errmsg2 @@ -64,7 +62,7 @@ contains allocate(model) !! Initializes the VF_RAD_PROB components of MODEL. - call vf_rad_init (tinit, mesh, model, er_params, stat, errmsg) + call vf_rad_init (tinit, mesh, model, stat, errmsg) if (stat /= 0) return !! Defines the heat equation parameter components of MODEL. @@ -86,15 +84,16 @@ contains end function create_FHT_model - subroutine vf_rad_init (tinit, mesh, model, params, stat, errmsg) + subroutine vf_rad_init (tinit, mesh, model, stat, errmsg) + use enclosure_radiation_namelist, only: params use bitfield_type, only: btest use parallel_communication, only: global_any, global_all + use parameter_list_type real(r8), intent(in) :: tinit type(unstr_mesh), intent(in) :: mesh type(FHT_model), intent(inout) :: model - type(parameter_list), intent(inout) :: params integer, intent(out) :: stat character(len=*), intent(out) :: errmsg diff --git a/src/truchas/physics/heat_species_transport/FHT_solver_factory.F90 b/src/truchas/physics/heat_species_transport/FHT_solver_factory.F90 index 687608504..9c596adbb 100644 --- a/src/truchas/physics/heat_species_transport/FHT_solver_factory.F90 +++ b/src/truchas/physics/heat_species_transport/FHT_solver_factory.F90 @@ -20,7 +20,7 @@ module FHT_solver_factory contains - function create_FHT_solver(mmf, model, params, er_params) result(solver) + function create_FHT_solver(mmf, model, params) result(solver) use enclosure_radiation_namelist, only: er_params => params use parallel_communication @@ -28,7 +28,7 @@ contains type(matl_mesh_func), intent(in), target :: mmf type(FHT_model), intent(in), target :: model - type(parameter_list), intent(inout) :: params, er_params + type(parameter_list), intent(inout) :: params type(FHT_solver), pointer :: solver integer :: j, n, lun diff --git a/src/truchas/physics/heat_species_transport/HTSD_model_factory.F90 b/src/truchas/physics/heat_species_transport/HTSD_model_factory.F90 index af497a768..c2e48429e 100644 --- a/src/truchas/physics/heat_species_transport/HTSD_model_factory.F90 +++ b/src/truchas/physics/heat_species_transport/HTSD_model_factory.F90 @@ -36,7 +36,6 @@ module HTSD_model_factory use species_bc_factory_class use thermal_source_factory_type use truchas_logging_services - use parameter_list_type implicit none private @@ -44,7 +43,7 @@ module HTSD_model_factory contains - function create_HTSD_model (tinit, disc, mmf, tbc_fac, sbc_fac, tsrc_fac, er_params, stat, errmsg) result (model) + function create_HTSD_model (tinit, disc, mmf, tbc_fac, sbc_fac, tsrc_fac, stat, errmsg) result (model) use diffusion_solver_data, only: heat_eqn, num_species, void_temperature @@ -54,7 +53,6 @@ contains class(thermal_bc_factory), intent(inout) :: tbc_fac class(species_bc_factory), intent(inout) :: sbc_fac type(thermal_source_factory), intent(inout) :: tsrc_fac - type(parameter_list), intent(inout) :: er_params integer, intent(out) :: stat character(*), intent(out) :: errmsg type(HTSD_model), pointer :: model @@ -66,7 +64,7 @@ contains mesh => disc%mesh if (heat_eqn) then - htmodel => create_HT_model(tinit, mesh, mmf, tbc_fac, tsrc_fac, er_params, stat, errmsg) + htmodel => create_HT_model(tinit, mesh, mmf, tbc_fac, tsrc_fac, stat, errmsg) if (stat /= 0) return endif @@ -82,7 +80,7 @@ contains end function create_HTSD_model - function create_HT_model (tinit, mesh, mmf, bc_fac, src_fac, er_params, stat, errmsg) result (model) + function create_HT_model (tinit, mesh, mmf, bc_fac, src_fac, stat, errmsg) result (model) use rad_problem_type @@ -91,7 +89,6 @@ contains type(matl_mesh_func), intent(in), target :: mmf class(thermal_bc_factory), intent(inout) :: bc_fac type(thermal_source_factory), intent(inout) :: src_fac - type(parameter_list), intent(inout) :: er_params integer, intent(out) :: stat character(*), intent(out) :: errmsg character(:), allocatable :: errmsg2 @@ -100,7 +97,7 @@ contains allocate(model) !! Initializes the VF_RAD_PROB components of MODEL. - model%vf_rad_prob => create_vf_rad_prob(tinit, mesh, er_params, stat, errmsg) + model%vf_rad_prob => create_vf_rad_prob(tinit, mesh, stat, errmsg) if (stat /= 0) return !! Defines the equation parameter components of MODEL. @@ -348,15 +345,16 @@ contains end function create_HT_model - function create_vf_rad_prob (tinit, mesh, params, stat, errmsg) result (vf_rad_prob) + function create_vf_rad_prob (tinit, mesh, stat, errmsg) result (vf_rad_prob) use rad_problem_type + use enclosure_radiation_namelist, only: params use bitfield_type, only: btest use parallel_communication, only: global_any, global_all + use parameter_list_type real(r8), intent(in) :: tinit type(unstr_mesh), intent(in) :: mesh - type(parameter_list), intent(inout) :: params integer, intent(out) :: stat character(len=*), intent(out) :: errmsg type(rad_problem), pointer :: vf_rad_prob(:) diff --git a/src/truchas/physics/heat_species_transport/HTSD_solver_factory.F90 b/src/truchas/physics/heat_species_transport/HTSD_solver_factory.F90 index 8ed2bb503..4a0217c0f 100644 --- a/src/truchas/physics/heat_species_transport/HTSD_solver_factory.F90 +++ b/src/truchas/physics/heat_species_transport/HTSD_solver_factory.F90 @@ -20,14 +20,15 @@ module HTSD_solver_factory contains - function create_HTSD_solver(mmf, model, params, er_params) result(solver) + function create_HTSD_solver(mmf, model, params) result(solver) + use enclosure_radiation_namelist, only: er_params => params use parallel_communication use truchas_env, only: output_file_name type(matl_mesh_func), intent(in), target :: mmf type(HTSD_model), intent(in), target :: model - type(parameter_list), intent(inout) :: params, er_params + type(parameter_list), intent(inout) :: params type(HTSD_solver), pointer :: solver integer :: j, n, lun diff --git a/src/truchas/physics/heat_species_transport/diffusion_solver.F90 b/src/truchas/physics/heat_species_transport/diffusion_solver.F90 index 8e591191b..d3b22fdb2 100644 --- a/src/truchas/physics/heat_species_transport/diffusion_solver.F90 +++ b/src/truchas/physics/heat_species_transport/diffusion_solver.F90 @@ -110,7 +110,7 @@ contains call read_species_bc_namelists(lun, this%species_bc_params) call read_ds_source (lun) - call read_enclosure_radiation_namelists(lun, this%er_params) + call read_enclosure_radiation_namelists(lun) call this%ds_params%get('void-temperature', void_temperature) @@ -508,12 +508,12 @@ contains call tbc_fac%init(this%mesh, stefan_boltzmann, absolute_zero, this%bc_params) call sbc_fac%init(this%mesh, this%species_bc_params) call tsrc_fac%init(this%mesh, this%thermal_source_params) - this%mod1 => create_HTSD_model(tinit, this%disc, this%mmf, tbc_fac, sbc_fac, tsrc_fac, this%er_params, stat, errmsg) + this%mod1 => create_HTSD_model(tinit, this%disc, this%mmf, tbc_fac, sbc_fac, tsrc_fac, stat, errmsg) end block if (stat /= 0) call TLS_fatal ('DS_INIT: ' // trim(errmsg)) if (this%have_heat_transfer) this%ht_source => this%mod1%ht%source if (this%have_species_diffusion) this%sd_source => this%mod1%sd%source - this%sol1 => create_HTSD_solver (this%mmf, this%mod1, this%ds_params, this%er_params) + this%sol1 => create_HTSD_solver (this%mmf, this%mod1, this%ds_params) case (SOLVER2) block @@ -521,11 +521,11 @@ contains type(thermal_source_factory) :: tsrc_fac call bc_fac%init(this%mesh, stefan_boltzmann, absolute_zero, this%bc_params) call tsrc_fac%init(this%mesh, this%thermal_source_params) - this%mod2 => create_FHT_model (tinit, this%disc, this%mmf, bc_fac, tsrc_fac, this%er_params, stat, errmsg) + this%mod2 => create_FHT_model (tinit, this%disc, this%mmf, bc_fac, tsrc_fac, stat, errmsg) end block if (stat /= 0) call TLS_fatal ('DS_INIT: ' // trim(errmsg)) this%ht_source => this%mod2%q ! we need this to set the advected heat at each step - this%sol2 => create_FHT_solver(this%mmf, this%mod2, this%ds_params, this%er_params) + this%sol2 => create_FHT_solver(this%mmf, this%mod2, this%ds_params) call move_alloc(this%hadv, this%sol2%hadv) case (SOLVER3) @@ -534,11 +534,11 @@ contains type(thermal_source_factory) :: tsrc_fac call tbc_fac%init(this%mesh, stefan_boltzmann, absolute_zero, this%bc_params) call tsrc_fac%init(this%mesh, this%thermal_source_params) - this%mod3 => create_ht_model(tinit, this%disc, this%mmf, tbc_fac, tsrc_fac, this%er_params, stat, errmsg2) + this%mod3 => create_ht_model(tinit, this%disc, this%mmf, tbc_fac, tsrc_fac, stat, errmsg2) end block if (stat /= 0) call TLS_fatal('DS_INIT: ' // trim(errmsg2)) this%ht_source => this%mod3%source - this%sol3 => create_ht_solver(this%mmf, this%mod3, this%ds_params, this%er_params, stat, errmsg2) + this%sol3 => create_ht_solver(this%mmf, this%mod3, this%ds_params, stat, errmsg2) if (stat /= 0) call TLS_fatal('DS_INIT: ' // trim(errmsg2)) case default INSIST(.false.) diff --git a/src/truchas/physics/heat_species_transport/ht_model_factory.F90 b/src/truchas/physics/heat_species_transport/ht_model_factory.F90 index ed6ca725d..b943388cd 100644 --- a/src/truchas/physics/heat_species_transport/ht_model_factory.F90 +++ b/src/truchas/physics/heat_species_transport/ht_model_factory.F90 @@ -24,8 +24,9 @@ module ht_model_factory contains - function create_ht_model(tinit, disc, mmf, bc_fac, src_fac, er_params, stat, errmsg) result(model) + function create_ht_model(tinit, disc, mmf, bc_fac, src_fac, stat, errmsg) result(model) + use enclosure_radiation_namelist, only: er_params => params use diffusion_solver_data, only: void_temperature use rad_problem_type @@ -34,7 +35,6 @@ contains type(matl_mesh_func), intent(in), target :: mmf class(thermal_bc_factory), intent(inout) :: bc_fac type(thermal_source_factory), intent(inout) :: src_fac - type(parameter_list), intent(inout) :: er_params integer, intent(out) :: stat character(:), allocatable , intent(out):: errmsg type(ht_model), pointer :: model @@ -119,7 +119,6 @@ contains use bitfield_type use parallel_communication, only: global_all, global_any, global_count use string_utilities, only: i_to_c - use f08_intrinsics, only: findloc type(unstr_mesh), intent(in), target :: mesh type(HT_model), intent(inout) :: model diff --git a/src/truchas/physics/heat_species_transport/ht_model_type.F90 b/src/truchas/physics/heat_species_transport/ht_model_type.F90 index 4a5b789cf..a610f8cef 100644 --- a/src/truchas/physics/heat_species_transport/ht_model_type.F90 +++ b/src/truchas/physics/heat_species_transport/ht_model_type.F90 @@ -15,7 +15,6 @@ module ht_model_type use prop_mesh_func_type use source_mesh_function use scalar_mesh_func_class - use bndry_vfunc_class use bndry_func1_class use bndry_func2_class use intfc_func2_class @@ -41,7 +40,7 @@ module ht_model_type !! Boundary condition data class(bndry_func1), allocatable :: bc_dir ! Dirichlet class(bndry_func1), allocatable :: bc_flux ! simple flux - class(bndry_vfunc), allocatable :: bc_vflux ! oriented flux + class(bndry_func2), allocatable :: bc_vflux ! oriented flux class(bndry_func2), allocatable :: bc_htc ! external HTC class(bndry_func2), allocatable :: bc_rad ! simple radiation class(intfc_func2), allocatable :: ic_htc ! internal HTC @@ -180,10 +179,10 @@ contains !! Oriented flux BC contribution. if (allocated(this%bc_vflux)) then - call this%bc_vflux%compute(t) + call this%bc_vflux%compute(t, u%tf) do j = 1, size(this%bc_vflux%index) n = this%bc_vflux%index(j) - f%tf(n) = f%tf(n) + dot_product(this%mesh%normal(:,n), this%bc_vflux%value(:,j)) + f%tf(n) = f%tf(n) + this%bc_vflux%value(j) end do end if diff --git a/src/truchas/physics/heat_species_transport/ht_solver_factory.F90 b/src/truchas/physics/heat_species_transport/ht_solver_factory.F90 index f0d9340fc..2b4b0307d 100644 --- a/src/truchas/physics/heat_species_transport/ht_solver_factory.F90 +++ b/src/truchas/physics/heat_species_transport/ht_solver_factory.F90 @@ -20,14 +20,15 @@ module ht_solver_factory contains - function create_ht_solver(mmf, model, params, er_params, stat, errmsg) result(solver) + function create_ht_solver(mmf, model, params, stat, errmsg) result(solver) + use enclosure_radiation_namelist, only: er_params => params use parallel_communication use truchas_env, only: output_file_name type(matl_mesh_func), intent(in), target :: mmf type(ht_model), intent(in), target :: model - type(parameter_list), intent(inout) :: params, er_params + type(parameter_list), intent(inout) :: params type(ht_solver), pointer :: solver integer, intent(out) :: stat character(:), allocatable, intent(out) :: errmsg -- GitLab From d1be91a902b9874c0b8be837a3ceecafc7fe877b Mon Sep 17 00:00:00 2001 From: "Neil N. Carlson" <neil.n.carlson@gmail.com> Date: Wed, 1 May 2024 07:23:48 -0600 Subject: [PATCH 5/6] Updates required to rebase on 318bfa7 --- .../physics/heat_species_transport/ht_model_factory.F90 | 2 +- .../physics/heat_species_transport/ht_model_type.F90 | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/truchas/physics/heat_species_transport/ht_model_factory.F90 b/src/truchas/physics/heat_species_transport/ht_model_factory.F90 index b943388cd..ea3c82616 100644 --- a/src/truchas/physics/heat_species_transport/ht_model_factory.F90 +++ b/src/truchas/physics/heat_species_transport/ht_model_factory.F90 @@ -104,7 +104,7 @@ contains call define_external_source (mesh, 'temperature', model%source) !! Additional heat sources - call src_fac%alloc_source_func1(model%src, stat, errmsg) + call src_fac%alloc_source_funcs(model%src, stat, errmsg) if (stat /= 0) return stat = 0 diff --git a/src/truchas/physics/heat_species_transport/ht_model_type.F90 b/src/truchas/physics/heat_species_transport/ht_model_type.F90 index a610f8cef..57227ea23 100644 --- a/src/truchas/physics/heat_species_transport/ht_model_type.F90 +++ b/src/truchas/physics/heat_species_transport/ht_model_type.F90 @@ -14,7 +14,7 @@ module ht_model_type use mfd_disc_type use prop_mesh_func_type use source_mesh_function - use scalar_mesh_func_class + use scalar_mesh_multifunc_type use bndry_func1_class use bndry_func2_class use intfc_func2_class @@ -36,7 +36,7 @@ module ht_model_type type(prop_mesh_func) :: H_of_T ! enthalpy as a function of temperature type(TofH) :: T_of_H type(source_mf) :: source ! external heat source - class(scalar_mesh_func), allocatable :: src ! another external heat source + type(scalar_mesh_multifunc), allocatable :: src ! another external heat source !! Boundary condition data class(bndry_func1), allocatable :: bc_dir ! Dirichlet class(bndry_func1), allocatable :: bc_flux ! simple flux @@ -153,7 +153,7 @@ contains !! Additional heat source if (allocated(this%src)) then - call this%src%compute(t) + call this%src%compute(t, u%tc) f%tc = f%tc - this%mesh%volume*this%src%value end if -- GitLab From 7ca71dee46107c2cd77724440d9fbf476cd8a594 Mon Sep 17 00:00:00 2001 From: "Neil N. Carlson" <neil.n.carlson@gmail.com> Date: Fri, 10 May 2024 17:28:34 -0600 Subject: [PATCH 6/6] Updates required to rebase branch on 9ead673 --- src/truchas/ode/new_idaesol_type.F90 | 12 ++++++------ .../physics/heat_species_transport/ht_norm_type.F90 | 10 +++++----- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/truchas/ode/new_idaesol_type.F90 b/src/truchas/ode/new_idaesol_type.F90 index 224786dd9..d6c20d344 100644 --- a/src/truchas/ode/new_idaesol_type.F90 +++ b/src/truchas/ode/new_idaesol_type.F90 @@ -168,8 +168,8 @@ contains this%model => model - context = 'processing ' // params%name() // ': ' - call params%get('nlk-max-iter', this%mitr, default=5, stat=stat, errmsg=errmsg) + context = 'processing ' // params%path() // ': ' + call params%get('nlk-max-iter', this%mitr, stat, errmsg, default=5) if (stat /= 0) then errmsg = context//errmsg return @@ -180,7 +180,7 @@ contains return end if - call params%get('nlk-tol', this%ntol, default=0.1_r8, stat=stat, errmsg=errmsg) + call params%get('nlk-tol', this%ntol, stat, errmsg, default=0.1_r8) if (stat /= 0) then errmsg = context//errmsg return @@ -191,7 +191,7 @@ contains return end if - call params%get('nlk-max-vec', maxv, default=this%mitr-1, stat=stat, errmsg=errmsg) + call params%get('nlk-max-vec', maxv, stat, errmsg, default=this%mitr-1) if (stat /= 0) then errmsg = context//errmsg return @@ -203,7 +203,7 @@ contains end if maxv = min(maxv, this%mitr-1) - call params%get('nlk-vec-tol', vtol, default=0.01_r8, stat=stat, errmsg=errmsg) + call params%get('nlk-vec-tol', vtol, stat, errmsg, default=0.01_r8) if (stat /= 0) then errmsg = context//errmsg return @@ -214,7 +214,7 @@ contains return end if - call params%get('pc-freq', this%pc_freq, default=huge(this%pc_freq), stat=stat, errmsg=errmsg) + call params%get('pc-freq', this%pc_freq, stat, errmsg, default=huge(this%pc_freq)) if (stat /= 0) then errmsg = context//errmsg return diff --git a/src/truchas/physics/heat_species_transport/ht_norm_type.F90 b/src/truchas/physics/heat_species_transport/ht_norm_type.F90 index b0e3f85b8..59c075168 100644 --- a/src/truchas/physics/heat_species_transport/ht_norm_type.F90 +++ b/src/truchas/physics/heat_species_transport/ht_norm_type.F90 @@ -44,8 +44,8 @@ contains this%model => model - context = 'processing ' // params%name() // ': ' - call params%get('abs-t-tol', this%temp_atol, default=0.0_r8, stat=stat, errmsg=errmsg) + context = 'processing ' // params%path() // ': ' + call params%get('abs-t-tol', this%temp_atol, stat, errmsg, default=0.0_r8) if (stat /= 0) then errmsg = context // errmsg return @@ -55,7 +55,7 @@ contains return end if - call params%get('rel-t-tol', this%temp_rtol, default=0.0_r8, stat=stat, errmsg=errmsg) + call params%get('rel-t-tol', this%temp_rtol, stat, errmsg, default=0.0_r8) if (stat /= 0) then errmsg = context // errmsg return @@ -71,7 +71,7 @@ contains return end if - call params%get('abs-h-tol', this%enth_atol, default=0.0_r8, stat=stat, errmsg=errmsg) + call params%get('abs-h-tol', this%enth_atol, stat, errmsg, default=0.0_r8) if (stat /= 0) then errmsg = context // errmsg return @@ -81,7 +81,7 @@ contains return end if - call params%get('rel-h-tol', this%enth_rtol, default=this%temp_rtol, stat=stat, errmsg=errmsg) + call params%get('rel-h-tol', this%enth_rtol, stat, errmsg, default=this%temp_rtol) if (stat /= 0) then errmsg = context // errmsg return -- GitLab