<p><b>cnewman@lanl.gov</b> 2010-08-05 15:04:06 -0600 (Thu, 05 Aug 2010)</p><p>Added trilinos interface files in implicit branch<br>
</p><hr noshade><pre><font color="gray">Modified: branches/implicit/Makefile
===================================================================
--- branches/implicit/Makefile 2010-08-05 20:02:04 UTC (rev 462)
+++ branches/implicit/Makefile 2010-08-05 21:04:06 UTC (rev 463)
@@ -12,7 +12,17 @@
endif
#########################
+TRILINOS_HOME = /usr/projects/climate/cnewman/trilinos-10.2.2-Source/GCC_MPI_DBG
+#########################
+# Section for Trilinos TPL
+#########################
+ifdef TRILINOS_HOME
+ TRILINOS_DEFINE = -DHAVE_TRILINOS
+endif
+#########################
+
+
dummy:
@( echo "try one of:"; \
echo " make xlf"; \
@@ -144,7 +154,31 @@
endif
#########################
+#########################
+# Section for Trilinos TPL
+#########################
+ifdef TRILINOS_HOME
+ include $(TRILINOS_HOME)/include/Makefile.export.Epetra
+ MPI_HOME=/opt/OpenMPI/openmpi-1.2.8-gcc/ib
+ MPIINC = -I$(MPI_HOME)/include
+ MPILIB = -L$(MPI_HOME)/lib64 -lmpi_cxx
+# LIBS += -lstdc++ #cn added -lstdc++ for c++ capability
+ ifdef TRILINOS_INC_PATH
+ FCINCLUDES += -I$(TRILINOS_INC_PATH)
+ else
+ CPPINCLUDES += -I$(TRILINOS_HOME)/include $(MPIINC)
+ #CPPFLAGS += $(NOX_CXXFLAGS) $(NOX_DEFS) $(NOX_CPPFLAGS) -I$(TRILINOS_HOME)/include $(CPPINCLUDES)
+ endif
+ ifdef TRILINOS_LIB_PATH
+ LIBS += -L$(TRILINOS_LIB_PATH) -lzoltan
+ else
+ LIBS += -L$(TRILINOS_HOME)/lib $(EPETRA_LIBRARIES) $(MPILIB) $(EPETRA_TPL_LIBRARIES) -lstdc++
+ endif
+endif
+#########################
+
+
ifdef CORE
all: mpas_main
Modified: branches/implicit/src/core_sw/Makefile
===================================================================
--- branches/implicit/src/core_sw/Makefile 2010-08-05 20:02:04 UTC (rev 462)
+++ branches/implicit/src/core_sw/Makefile 2010-08-05 21:04:06 UTC (rev 463)
@@ -1,12 +1,12 @@
.SUFFIXES: .F .o
-#IMPOBJS=module_implicit_integration.o #cn
+#IMPOBJS=module_implicit_integration.o nox_driver.o #cn
OBJS = module_test_cases.o \
module_time_integration_support.o \
module_implicit_integration.o \
module_time_integration.o \
- mpas_interface.o
+ mpas_interface.o nox_driver.o
all: core_sw
@@ -17,9 +17,9 @@
module_time_integration.o:
-module_implicit_integration.o:
+module_implicit_integration.o:
-mpas_interface.o: module_time_integration_support.o module_implicit_integration.o module_test_cases.o module_time_integration.o
+mpas_interface.o: module_time_integration_support.o module_implicit_integration.o module_test_cases.o module_time_integration.o nox_driver.o
clean:
$(RM) *.o *.mod *.f90 libdycore.a
@@ -28,3 +28,6 @@
$(RM) $@ $*.mod
$(CPP) $(CPPFLAGS) $(CPPINCLUDES) $< > $*.f90
$(FC) $(FFLAGS) -c $*.f90 $(FCINCLUDES) -I../framework -I../operators
+
+.C.o:
+ g++ -c $(CXXFLAGS) $(CPPFLAGS) $*.C $(CPPINCLUDES)
Modified: branches/implicit/src/core_sw/module_implicit_integration.F
===================================================================
--- branches/implicit/src/core_sw/module_implicit_integration.F 2010-08-05 20:02:04 UTC (rev 462)
+++ branches/implicit/src/core_sw/module_implicit_integration.F 2010-08-05 21:04:06 UTC (rev 463)
@@ -1,335 +1,407 @@
module implicit_integration
- use vector_reconstruction
- use grid_types
- use configure
- use constants
- use dmpar
- use time_integration_support
+ use vector_reconstruction
+ use grid_types
+ use configure
+ use constants
+ use dmpar
+ use time_integration_support
+ use, intrinsic :: iso_c_binding, only : c_double, c_int, c_ptr, &
+ c_loc, c_funptr, c_funloc, c_f_pointer
- type (domain_type), save ,pointer:: dom
- real (kind=RKIND),save:: d_t
+ type (domain_type), save ,pointer:: dom
+ real (kind=RKIND),save:: d_t
+ real(kind=RKIND), allocatable,dimension(:)::narray
- contains
+ type, public :: objectType
+ sequence
- subroutine implicit_timestep(dt)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Advance model state forward in time by the specified time step
- !
- ! Input: domain - current model state in time level 1 (e.g., time_levs(1)state%h(:,:))
- ! plus grid meta-data
- ! Output: domain - upon exit, time level 2 (e.g., time_levs(2)%state%h(:,:)) contains
- ! model state advanced forward in time by dt seconds
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ integer :: ndum
+ real (kind=RKIND) :: rdum
- implicit none
+ end type objectType
- type (block_type), pointer :: block
+ type, public :: preconType
+ sequence
-!cn this code may be better off in an init subroutine
- real (kind=RKIND), intent(in) :: dt
- d_t=dt
-!cn end this code may be better off in an init subroutine
+ integer :: ndum
+ real (kind=RKIND) :: rdum
- if (trim(config_time_integration) == 'NRK4') then
- call do_nrk4()
- else
- write(0,*) 'Unknown time integration option '//trim(config_time_integration)
- stop
- end if
+ end type preconType
- block => dom % blocklist
- do while (associated(block))
- block % time_levs(2) % state % xtime % scalar = block % time_levs(1) % state % xtime % scalar + d_t
- block => block % next
- end do
+ type(objectType) ,target :: state_object
+ type(objectType) ,pointer :: fptr=>NULL()
+ type(c_ptr) :: c_ptr_to_object
+ type(preconType) ,target :: pre_object
+ type(preconType) ,pointer :: pptr=>NULL()
+ type(c_ptr) :: c_ptr_to_pre
- end subroutine implicit_timestep
+contains
- subroutine do_nrk4()
+ subroutine implicit_timestep(dt)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Advance model state forward in time by the specified time step
+ !
+ ! Input: domain - current model state in time level 1 (e.g., time_levs(1)state%h(:,:))
+ ! plus grid meta-data
+ ! Output: domain - upon exit, time level 2 (e.g., time_levs(2)%state%h(:,:)) contains
+ ! model state advanced forward in time by dt seconds
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- real(kind=RKIND), allocatable,dimension(:)::narray
- allocate(narray(dom%blocklist%mesh%nvertlevelssolve*dom%blocklist%mesh%nedgessolve &
- +dom%blocklist%mesh%nvertlevelssolve*dom%blocklist%mesh%ncellssolve &
- +dom%blocklist%mesh%ntracerssolve*dom%blocklist%mesh%nvertlevelssolve*dom%blocklist%mesh%ncellssolve))
+ implicit none
- call pack(narray, dom%blocklist % intermediate_step(2) % u % array(:,:), &
- dom%blocklist % intermediate_step(2) % h % array(:,:), &
- dom%blocklist % intermediate_step(2) % tracers % array(:,:,:))
+ type (block_type), pointer :: block
- call nrk4()
- end subroutine do_nrk4
+ !cn this code may be better off in an init subroutine
+ real (kind=RKIND), intent(in) :: dt
+ d_t=dt
+ !cn end this code may be better off in an init subroutine
+ if (trim(config_time_integration) == 'NRK4') then
+ call do_nrk4()
+ else
+ write(0,*) 'Unknown time integration option '//trim(config_time_integration)
+ stop
+ end if
+
+ block => dom % blocklist
+ do while (associated(block))
+ block % time_levs(2) % state % xtime % scalar = block % time_levs(1) % state % xtime % scalar + d_t
+ block => block % next
+ end do
+
+ end subroutine implicit_timestep
+
+ subroutine do_nrk4()
+ implicit none
+ integer :: iCell, k
+ type (block_type), pointer :: block
+ integer, parameter :: PROVIS = 1
+ integer, parameter :: TEND = 2
+
+ block => dom % blocklist
+ do while (associated(block))
+ block % time_levs(2) % state % u % array(:,:) = block % time_levs(1) % state % u % array(:,:)
+ block % time_levs(2) % state % h % array(:,:) = block % time_levs(1) % state % h % array(:,:)
+ do iCell=1,block % mesh % nCells ! couple tracers to h
+ do k=1,block % mesh % nVertLevels
+ block % time_levs(2) % state % tracers % array(:,k,iCell) = block % time_levs(1) % state % tracers % array(:,k,iCell) &
+ * block % time_levs(1) % state % h % array(k,iCell)
+ end do
+ end do
+
+ call copy_state(block % time_levs(1) % state, block % intermediate_step(PROVIS))
+
+ block => block % next
+ end do
+
+ !need to get this right--this will be values at start of timestep
+ call pack(narray, dom%blocklist % time_levs(1) % state %u % array(:,:), &
+ dom%blocklist % time_levs(1) % state % h % array(:,:), &
+ dom%blocklist % time_levs(1) % state %tracers % array(:,:,:))
+
+ call nrk4(narray)
+
+ !need to get this right--this will be values at end of timestep
+ call unpack(narray, dom%blocklist % time_levs(2) % state %u % array(:,:), &
+ dom%blocklist % time_levs(2) % state %h % array(:,:), &
+ dom%blocklist % time_levs(2) % state %tracers % array(:,:,:))
+ !
+ ! A little clean up at the end: decouple new scalar fields and compute diagnostics for new state
+ !
+ block => dom % blocklist
+ do while (associated(block))
+ do iCell=1,block % mesh % nCells
+ do k=1,block % mesh % nVertLevels
+ block % time_levs(2) % state % tracers % array(:,k,iCell) = &
+ block % time_levs(2) % state % tracers % array(:,k,iCell) &
+ / block % time_levs(2) % state % h % array(k,iCell)
+ end do
+ end do
+
+ if (config_test_case == 1) then ! For case 1, wind field should be fixed
+ block % time_levs(2) % state % u % array(:,:) = block % time_levs(1) % state % u % array(:,:)
+ end if
+
+ call compute_solve_diagnostics(d_t, block % time_levs(2) % state, block % mesh)
+
+ call reconstruct(block % time_levs(2) % state, block % mesh)
+
+ block => block % next
+ end do
+ end subroutine do_nrk4
+
subroutine implicit_init(domain)
- implicit none
- type (domain_type), intent(inout) ,target:: domain
- dom=>domain
- end subroutine implicit_init
+ implicit none
+ type (domain_type), intent(inout) ,target:: domain
+ !real (kind=RKIND),allocatable,dimension(:) :: xstate!cn
+ interface
- subroutine nrk4()
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Advance model state forward in time by the specified time step using
- ! 4th order Runge-Kutta
- !
- ! Input: domain - current model state in time level 1 (e.g., time_levs(1)state%h(:,:))
- ! plus grid meta-data
- ! Output: domain - upon exit, time level 2 (e.g., time_levs(2)%state%h(:,:)) contains
- ! model state advanced forward in time by dt seconds
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ subroutine nox_init(vectorSize,vector,comm,v_container,p_container) &
+ bind(C,name='noxinit')
+ use, intrinsic :: iso_c_binding
+ integer(c_int) :: vectorSize,comm
+ real(c_double), dimension(*) :: vector
+ type(c_ptr) :: v_container
+ type(c_ptr) :: p_container !precon ptr
+ end subroutine nox_init
- implicit none
+ end interface
- integer :: iCell, k
- type (block_type), pointer :: block
+ dom=>domain
+ allocate(narray(array_size()))
- integer, parameter :: PROVIS = 1
- integer, parameter :: TEND = 2
- integer :: rk_step
+ call nox_init(size(narray), narray, 1, c_ptr_to_object, c_ptr_to_pre)
- real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights
+ end subroutine implicit_init
- !real (kind=RKIND), pointer :: uu(:,:)!cn
+ subroutine nrk4(noxarray)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! Advance model state forward in time by the specified time step using
+ ! 4th order Runge-Kutta
+ !
+ ! Input: domain - current model state in time level 1 (e.g., time_levs(1)state%h(:,:))
+ ! plus grid meta-data
+ ! Output: domain - upon exit, time level 2 (e.g., time_levs(2)%state%h(:,:)) contains
+ ! model state advanced forward in time by dt seconds
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !call unpack(narray, dom%blocklist % intermediate_step(TEND) % u % array(:,:), &
- ! dom%blocklist % intermediate_step(TEND) % h % array(:,:), &
- ! dom%blocklist % intermediate_step(TEND) % tracers % array(:,:,:))
-
- !
- ! Initialize time_levs(2) with state at current time
- ! Initialize first RK state
- ! Couple tracers time_levs(2) with h in time-levels
- ! Initialize RK weights
- !
- !block => domain % blocklist
- block => dom % blocklist
- do while (associated(block))
- block % time_levs(2) % state % u % array(:,:) = block % time_levs(1) % state % u % array(:,:)
- block % time_levs(2) % state % h % array(:,:) = block % time_levs(1) % state % h % array(:,:)
- do iCell=1,block % mesh % nCells ! couple tracers to h
- do k=1,block % mesh % nVertLevels
- block % time_levs(2) % state % tracers % array(:,k,iCell) = block % time_levs(1) % state % tracers % array(:,k,iCell) &
- * block % time_levs(1) % state % h % array(k,iCell)
- end do
- end do
+ implicit none
- call copy_state(block % time_levs(1) % state, block % intermediate_step(PROVIS))
+ real(kind=RKIND), dimension(:),intent(inout)::noxarray
+ integer :: iCell, k
+ type (block_type), pointer :: block
- block => block % next
- end do
+ integer, parameter :: PROVIS = 1
+ integer, parameter :: TEND = 2
+ integer :: rk_step
- rk_weights(1) = d_t/6.
- rk_weights(2) = d_t/3.
- rk_weights(3) = d_t/3.
- rk_weights(4) = d_t/6.
+ real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights
- rk_substep_weights(1) = d_t/2.
- rk_substep_weights(2) = d_t/2.
- rk_substep_weights(3) = d_t
- rk_substep_weights(4) = 0.
+ call unpack(noxarray, dom%blocklist % time_levs(1) % state %u % array(:,:), &
+ dom%blocklist % time_levs(1) % state %h % array(:,:), &
+ dom%blocklist % time_levs(1) % state %tracers % array(:,:,:))
+ !real (kind=RKIND), pointer :: uu(:,:)!cn
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! BEGIN RK loop
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- do rk_step = 1, 4
+ !
+ ! Initialize time_levs(2) with state at current time
+ ! Initialize first RK state
+ ! Couple tracers time_levs(2) with h in time-levels
+ ! Initialize RK weights
+ !
+ !block => domain % blocklist
-! --- update halos for diagnostic variables
+ rk_weights(1) = d_t/6.
+ rk_weights(2) = d_t/3.
+ rk_weights(3) = d_t/3.
+ rk_weights(4) = d_t/6.
- block => dom % blocklist
- do while (associated(block))
- call dmpar_exch_halo_field2dReal(dom % dminfo, block % intermediate_step(PROVIS) % pv_edge % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nEdges, &
- block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
- block => block % next
- end do
+ rk_substep_weights(1) = d_t/2.
+ rk_substep_weights(2) = d_t/2.
+ rk_substep_weights(3) = d_t
+ rk_substep_weights(4) = 0.
-! --- compute tendencies
- block => dom % blocklist
- do while (associated(block))
- call compute_tend(block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
- call compute_scalar_tend(block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
- call enforce_boundaryEdge(block % intermediate_step(TEND), block % mesh)
- block => block % next
- end do
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! BEGIN RK loop
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ do rk_step = 1, 4
-! --- update halos for prognostic variables
+ ! --- update halos for diagnostic variables
- block => dom % blocklist
- do while (associated(block))
- call dmpar_exch_halo_field2dReal(dom % dminfo, block % intermediate_step(TEND) % u % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nEdges, &
- block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
- call dmpar_exch_halo_field2dReal(dom % dminfo, block % intermediate_step(TEND) % h % array(:,:), &
- block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- call dmpar_exch_halo_field3dReal(dom % dminfo, block % intermediate_step(TEND) % tracers % array(:,:,:), &
- block % mesh % nTracers, block % mesh % nVertLevels, block % mesh % nCells, &
- block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
- block => block % next
- end do
+ block => dom % blocklist
+ do while (associated(block))
+ call dmpar_exch_halo_field2dReal(dom % dminfo, block % intermediate_step(PROVIS) % pv_edge % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ block => block % next
+ end do
-! --- compute next substep state
+ ! --- compute tendencies
- if (rk_step < 4) then
-!cn this is k1, k2, k3, k4
- block => dom % blocklist
- do while (associated(block))
- !uu=>block % intermediate_step(PROVIS) % u % array(:,:) !cn experimenting here
- !uu = block % time_levs(1) % state % u % array(:,:) &
- ! + rk_substep_weights(rk_step) * block % intermediate_step(TEND) % u % array(:,:)
- block % intermediate_step(PROVIS) % u % array(:,:) = block % time_levs(1) % state % u % array(:,:) &
- + rk_substep_weights(rk_step) * block % intermediate_step(TEND) % u % array(:,:)
- block % intermediate_step(PROVIS) % h % array(:,:) = block % time_levs(1) % state % h % array(:,:) &
- + rk_substep_weights(rk_step) * block % intermediate_step(TEND) % h % array(:,:)
- do iCell=1,block % mesh % nCells
- do k=1,block % mesh % nVertLevels
- block % intermediate_step(PROVIS) % tracers % array(:,k,iCell) = ( &
- block % time_levs(1) % state % h % array(k,iCell) * &
- block % time_levs(1) % state % tracers % array(:,k,iCell) &
- + rk_substep_weights(rk_step) * block % intermediate_step(TEND) % tracers % array(:,k,iCell) &
- ) / block % intermediate_step(PROVIS) % h % array(k,iCell)
- end do
- end do
- if (config_test_case == 1) then ! For case 1, wind field should be fixed
- block % intermediate_step(PROVIS) % u % array(:,:) = block % time_levs(1) % state % u % array(:,:)
- end if
- call compute_solve_diagnostics(d_t, block % intermediate_step(PROVIS), block % mesh)
- block => block % next
- end do
- end if
+ block => dom % blocklist
+ do while (associated(block))
+ call compute_tend(block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
+ call compute_scalar_tend(block % intermediate_step(TEND), block % intermediate_step(PROVIS), block % mesh)
+ call enforce_boundaryEdge(block % intermediate_step(TEND), block % mesh)
+ block => block % next
+ end do
-!--- accumulate update (for RK4)
+ ! --- update halos for prognostic variables
-!cn this is dt/6 (k1+2 k2+2 k3+k4)
- block => dom % blocklist
- do while (associated(block))
- block % time_levs(2) % state % u % array(:,:) = block % time_levs(2) % state % u % array(:,:) &
- + rk_weights(rk_step) * block % intermediate_step(TEND) % u % array(:,:)
- block % time_levs(2) % state % h % array(:,:) = block % time_levs(2) % state % h % array(:,:) &
- + rk_weights(rk_step) * block % intermediate_step(TEND) % h % array(:,:)
- do iCell=1,block % mesh % nCells
- do k=1,block % mesh % nVertLevels
- block % time_levs(2) % state % tracers % array(:,k,iCell) = &
- block % time_levs(2) % state % tracers % array(:,k,iCell) &
- + rk_weights(rk_step) * block % intermediate_step(TEND) % tracers % array(:,k,iCell)
- end do
- end do
- block => block % next
- end do
+ block => dom % blocklist
+ do while (associated(block))
+ call dmpar_exch_halo_field2dReal(dom % dminfo, block % intermediate_step(TEND) % u % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nEdges, &
+ block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
+ call dmpar_exch_halo_field2dReal(dom % dminfo, block % intermediate_step(TEND) % h % array(:,:), &
+ block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ call dmpar_exch_halo_field3dReal(dom % dminfo, block % intermediate_step(TEND) % tracers % array(:,:,:), &
+ block % mesh % nTracers, block % mesh % nVertLevels, block % mesh % nCells, &
+ block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
+ block => block % next
+ end do
- end do
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! END RK loop
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! --- compute next substep state
- !
- ! A little clean up at the end: decouple new scalar fields and compute diagnostics for new state
- !
- block => dom % blocklist
- do while (associated(block))
- do iCell=1,block % mesh % nCells
- do k=1,block % mesh % nVertLevels
- block % time_levs(2) % state % tracers % array(:,k,iCell) = &
- block % time_levs(2) % state % tracers % array(:,k,iCell) &
- / block % time_levs(2) % state % h % array(k,iCell)
- end do
- end do
+ if (rk_step < 4) then
+ !cn this is k1, k2, k3, k4
+ block => dom % blocklist
+ do while (associated(block))
+ !uu=>block % intermediate_step(PROVIS) % u % array(:,:) !cn experimenting here
+ !uu = block % time_levs(1) % state % u % array(:,:) &
+ ! + rk_substep_weights(rk_step) * block % intermediate_step(TEND) % u % array(:,:)
+ block % intermediate_step(PROVIS) % u % array(:,:) = block % time_levs(1) % state % u % array(:,:) &
+ + rk_substep_weights(rk_step) * block % intermediate_step(TEND) % u % array(:,:)
+ block % intermediate_step(PROVIS) % h % array(:,:) = block % time_levs(1) % state % h % array(:,:) &
+ + rk_substep_weights(rk_step) * block % intermediate_step(TEND) % h % array(:,:)
+ do iCell=1,block % mesh % nCells
+ do k=1,block % mesh % nVertLevels
+ block % intermediate_step(PROVIS) % tracers % array(:,k,iCell) = ( &
+ block % time_levs(1) % state % h % array(k,iCell) * &
+ block % time_levs(1) % state % tracers % array(:,k,iCell) &
+ + rk_substep_weights(rk_step) * block % intermediate_step(TEND) % tracers % array(:,k,iCell) &
+ ) / block % intermediate_step(PROVIS) % h % array(k,iCell)
+ end do
+ end do
+ if (config_test_case == 1) then ! For case 1, wind field should be fixed
+ block % intermediate_step(PROVIS) % u % array(:,:) = block % time_levs(1) % state % u % array(:,:)
+ end if
+ call compute_solve_diagnostics(d_t, block % intermediate_step(PROVIS), block % mesh)
+ block => block % next
+ end do
+ end if
- if (config_test_case == 1) then ! For case 1, wind field should be fixed
- block % time_levs(2) % state % u % array(:,:) = block % time_levs(1) % state % u % array(:,:)
- end if
+ !--- accumulate update (for RK4)
- call compute_solve_diagnostics(d_t, block % time_levs(2) % state, block % mesh)
+ !cn this is dt/6 (k1+2 k2+2 k3+k4)
+ block => dom % blocklist
+ do while (associated(block))
+ block % time_levs(2) % state % u % array(:,:) = block % time_levs(2) % state % u % array(:,:) &
+ + rk_weights(rk_step) * block % intermediate_step(TEND) % u % array(:,:)
+ block % time_levs(2) % state % h % array(:,:) = block % time_levs(2) % state % h % array(:,:) &
+ + rk_weights(rk_step) * block % intermediate_step(TEND) % h % array(:,:)
+ do iCell=1,block % mesh % nCells
+ do k=1,block % mesh % nVertLevels
+ block % time_levs(2) % state % tracers % array(:,k,iCell) = &
+ block % time_levs(2) % state % tracers % array(:,k,iCell) &
+ + rk_weights(rk_step) * block % intermediate_step(TEND) % tracers % array(:,k,iCell)
+ end do
+ end do
+ block => block % next
+ end do
- call reconstruct(block % time_levs(2) % state, block % mesh)
+ end do
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ ! END RK loop
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ call pack(noxarray, dom%blocklist % time_levs(2) % state %u % array(:,:), &
+ dom%blocklist % time_levs(2) % state %h % array(:,:), &
+ dom%blocklist % time_levs(2) % state %tracers % array(:,:,:))
+ end subroutine nrk4
- block => block % next
- end do
+ subroutine pack(noxarray, u, h, tracers)
+ !this is acrude copy at this point, what about multiple blocks?
+ !can probably clean it up a bit also
+ !could we just pass state instead of u,h,tracers?
- end subroutine nrk4
+ !this will put all variables into one long array; might be better if noxarray was a pointer that lives on the mod....
+ !u is nvertlevelssolve, nedgessolve
+ !h is nvertlevelssolve, ncellssolve
+ !t is ntracerssolve,nvertlevelssolve, ncellssolve is it ntracers?
+ ! => noxarray is nvertlevelssolve,*nedgessolve+nvertlevelssolve,*ncellssolve+ntracerssolve*nvertlevelssolve*ncellssolve
+ implicit none
+ real (kind=RKIND), dimension(:),intent(out):: noxarray
+ real (kind=RKIND), dimension(:,:),intent(in):: u
+ real (kind=RKIND), dimension(:,:),intent(in):: h
+ real (kind=RKIND), dimension(:,:,:),intent(in):: tracers
+ integer :: nedge, nvert, ncell, ntracer
+ integer:: iedge, ivert, icell, itracer
+ integer :: n
+ nedge = dom%blocklist%mesh%nedgessolve
+ nvert = dom%blocklist%mesh%nvertlevelssolve
+ ncell = dom%blocklist%mesh%ncellssolve
+ ntracer = dom%blocklist%mesh%ntracerssolve
+ do iedge = 1, nedge
+ do ivert = 1, nvert
+ n=ivert+nvert*(iedge-1)
+ noxarray(n)=u(ivert, iedge)
+ end do
+ end do
+ do icell = 1, ncell
+ do ivert = 1, nvert
+ n=ivert+nvert*(icell-1)+nvert*nedge
+ noxarray(n)=h(ivert, icell)
+ end do
+ end do
+ do icell = 1, ncell !cn this loop is unverified right now
+ do ivert = 1, nvert
+ do itracer = 1, ntracer
+ n=itracer+ntracer*(ivert-1)+ntracer*nvert*(icell-1)+nvert*nedge+nvert*ncell
+ noxarray(n)=tracers(itracer, ivert, icell)
+ end do
+ end do
+ end do
+ end subroutine pack
- subroutine pack(noxarray, u, h, tracers)
-!this is acrude copy at this point, what about multiple blocks?
+ subroutine unpack(noxarray, u, h, tracers)
+ ! needs to be adjusted for multiple blocks
+ implicit none
+ real (kind=RKIND), dimension(:),intent(in):: noxarray
+ real (kind=RKIND), dimension(:,:),intent(out):: u
+ real (kind=RKIND), dimension(:,:),intent(out):: h
+ real (kind=RKIND), dimension(:,:,:),intent(out):: tracers
+ integer :: nedge, nvert, ncell, ntracer
+ integer:: iedge, ivert, icell, itracer
+ integer :: n
+ nedge = dom%blocklist%mesh%nedgessolve
+ nvert = dom%blocklist%mesh%nvertlevelssolve
+ ncell = dom%blocklist%mesh%ncellssolve
+ ntracer = dom%blocklist%mesh%ntracerssolve
+ do iedge = 1, nedge
+ do ivert = 1, nvert
+ n=ivert+nvert*(iedge-1)
+ u(ivert, iedge)=noxarray(n)
+ end do
+ end do
+ do icell = 1, ncell
+ do ivert = 1, nvert
+ n=ivert+nvert*(icell-1)+nvert*nedge
+ h(ivert, icell)=noxarray(n)
+ end do
+ end do
+ do icell = 1, ncell !cn this loop is unverified right now
+ do ivert = 1, nvert
+ do itracer = 1, ntracer
+ n=itracer+ntracer*(ivert-1)+ntracer*nvert*(icell-1)+nvert*nedge+nvert*ncell
+ tracers(itracer, ivert, icell)=noxarray(n)
+ end do
+ end do
+ end do
-!this will put all variables into one long array; might be better if noxarray was a pointer that lives on the mod....
-!u is nvertlevelssolve, nedgessolve
-!h is nvertlevelssolve, ncellssolve
-!t is ntracerssolve,nvertlevelssolve, ncellssolve is it ntracers?
-! => noxarray is nvertlevelssolve,*nedgessolve+nvertlevelssolve,*ncellssolve+ntracerssolve*nvertlevelssolve*ncellssolve
- implicit none
- real (kind=RKIND), dimension(:),intent(out):: noxarray
- real (kind=RKIND), dimension(:,:),intent(in):: u
- real (kind=RKIND), dimension(:,:),intent(in):: h
- real (kind=RKIND), dimension(:,:,:),intent(in):: tracers
- integer :: nedge, nvert, ncell, ntracer
- integer:: iedge, ivert, icell, itracer
- integer :: n
-! dom%blocklist%mesh%nvertlevelssolve*dom%blocklist%mesh%nedgessolve &
-! +dom%blocklist%mesh%nvertlevelssolve*dom%blocklist%mesh%ncellssolve &
-! +dom%blocklist%mesh%ntracerssolve*dom%blocklist%mesh%nvertlevelssolve*dom%blocklist%mesh%ncellssolve
- nedge = dom%blocklist%mesh%nedgessolve
- nvert = dom%blocklist%mesh%nvertlevelssolve
- ncell = dom%blocklist%mesh%ncellssolve
- ntracer = dom%blocklist%mesh%ntracerssolve
- do iedge = 1, nedge
- do ivert = 1, nvert
- n=ivert+nvert*(iedge-1)
- noxarray(n)=u(ivert, iedge)
- end do
- end do
- do icell = 1, ncell
- do ivert = 1, nvert
- n=ivert+nvert*(icell-1)+nvert*nedge
- noxarray(n)=h(ivert, icell)
- end do
- end do
- do icell = 1, ncell !cn this loop is unverified right now
- do ivert = 1, nvert
- do itracer = 1, ntracer
- n=itracer+ntracer*(ivert-1)+ntracer*nvert*(icell-1)+nvert*nedge+nvert*ncell
- noxarray(n)=tracers(itracer, ivert, icell)
- end do
- end do
- end do
- end subroutine pack
+ end subroutine unpack
- subroutine unpack(noxarray, u, h, tracers)
- implicit none
- real (kind=RKIND), dimension(:),intent(in):: noxarray
- real (kind=RKIND), dimension(:,:),intent(out):: u
- real (kind=RKIND), dimension(:,:),intent(out):: h
- real (kind=RKIND), dimension(:,:,:),intent(out):: tracers
- integer :: nedge, nvert, ncell, ntracer
- integer:: iedge, ivert, icell, itracer
- integer :: n
- nedge = dom%blocklist%mesh%nedgessolve
- nvert = dom%blocklist%mesh%nvertlevelssolve
- ncell = dom%blocklist%mesh%ncellssolve
- ntracer = dom%blocklist%mesh%ntracerssolve
- do iedge = 1, nedge
- do ivert = 1, nvert
- n=ivert+nvert*(iedge-1)
- u(ivert, iedge)=noxarray(n)
- end do
- end do
- do icell = 1, ncell
- do ivert = 1, nvert
- n=ivert+nvert*(icell-1)+nvert*nedge
- h(ivert, icell)=noxarray(n)
- end do
- end do
- do icell = 1, ncell !cn this loop is unverified right now
- do ivert = 1, nvert
- do itracer = 1, ntracer
- n=itracer+ntracer*(ivert-1)+ntracer*nvert*(icell-1)+nvert*nedge+nvert*ncell
- tracers(itracer, ivert, icell)=noxarray(n)
- end do
- end do
- end do
-
- end subroutine unpack
+ integer function array_size()
+ ! needs to be adjusted for multiple blocks
+ implicit none
+ array_size = dom%blocklist%mesh%nvertlevelssolve*dom%blocklist%mesh%nedgessolve &
+ +dom%blocklist%mesh%nvertlevelssolve*dom%blocklist%mesh%ncellssolve &
+ +dom%blocklist%mesh%ntracerssolve*dom%blocklist%mesh%nvertlevelssolve*dom%blocklist%mesh%ncellssolve
+ end function array_size
+ subroutine Residual(xstate, fx, nelemd, c_ptr_to_object, &
+ timeScalingFactor) bind(C,name='calc_f')
+
+ real (c_double), intent(in) :: xstate(nelemd)
+ integer(c_int), intent(in) ,value :: nelemd
+ real (c_double), intent(out) :: fx(nelemd)
+ type(objectType), pointer :: fptr=>NULL()
+ type(c_ptr) :: c_ptr_to_object
+ real (c_double), intent(in),value :: timeScalingFactor
+ end subroutine Residual
+
end module implicit_integration
Added: branches/implicit/src/core_sw/nox_driver.C
===================================================================
--- branches/implicit/src/core_sw/nox_driver.C (rev 0)
+++ branches/implicit/src/core_sw/nox_driver.C 2010-08-05 21:04:06 UTC (rev 463)
@@ -0,0 +1,127 @@
+
+#include "NOX.H"
+#include "NOX_Epetra.H"
+
+#ifdef HAVE_MPI
+#include "Epetra_MpiComm.h"
+#else
+#include "Epetra_SerialComm.h"
+#endif
+
+#include "nox_interface.h"
+
+extern "C" {
+ // input: state vector x, output vector f, number of local elements nelem,
+ // C-wrapped fortran object
+ void calc_f(double *, double *, int, void *, double);
+ // input: some output vector vv (not sure), input vector f, number of local elements nelem,
+ // state vector, C-wrapped Fortran precond object
+ //void precon(double *, double *, int, double*, void *);
+ // matrix size and number of nonzeros
+ // input: quality flag
+ // output: nrows, nnz
+ //void matrix_size(int, int*, int*);
+ // make matrix graph
+ // input: quality flag, nrows, nnz, allocated arrays
+ // output: row pointer[nrows+1], column indices[nnz] (CRS graph)
+ //void make_graph(int, int, int, int*, int*);
+ // void fil_matrix(...)
+ //void get_quality(int*);
+ //void set_quality(int*);
+}
+
+static Teuchos::RCP<Problem_Interface> interface;
+
+extern "C" {
+
+void noxinit(int* nelems, double* statevector, int* mpi_comm_ignored,
+ void* blackbox_res, void* blackbox_prec)
+{
+
+ void (*residualFunction)(double *, double *, int, void *, double) = calc_f;
+ //void (*precFunction)(double *, double *, int, double*, void *) = precon;
+ //void (*matrixSizeFunction)(int, int*, int*) = matrix_size;
+ // void (*makeGraphFunction)(int, int, int, int*, int*) = make_graph;
+ // Create a communicator for Epetra objects
+#ifdef HAVE_MPI
+ Epetra_MpiComm Comm(MPI_COMM_WORLD);
+#else
+ Epetra_SerialComm Comm;
+#endif
+ int MyPID = Comm.MyPID();
+ // Create parameter (options) list
+ Teuchos::RCP<Teuchos::ParameterList> paramList =
+ Teuchos::rcp(new Teuchos::ParameterList);
+ Teuchos::ParameterList& nlParams = paramList->sublist("NOX");
+ Teuchos::ParameterList& nlDir = nlParams.sublist("Direction");
+ Teuchos::ParameterList& nlNewton = nlDir.sublist("Newton");
+ Teuchos::ParameterList& lsParams = nlNewton.sublist("Linear Solver");
+ Teuchos::ParameterList& nlPrintParams = nlParams.sublist("Printing");
+ nlPrintParams.set("MyPID", MyPID);
+ if (!nlPrintParams.isParameter("Output Information"))
+ nlPrintParams.set("Output Information",
+ NOX::Utils::OuterIteration +
+// NOX::Utils::OuterIterationStatusTest +
+ NOX::Utils::InnerIteration +
+ // NOX::Utils::Details +
+ NOX::Utils::LinearSolverDetails +
+ NOX::Utils::Warning
+// NOX::Utils::StepperIteration +
+// NOX::Utils::StepperDetails +
+// NOX::Utils::StepperParameters
+ );
+
+ string prec_type=lsParams.get("Preconditioner","None");
+
+ // Sublist for line search
+ Teuchos::ParameterList& searchParams = nlParams.sublist("Line Search");
+ searchParams.set("Method", "Full Step");
+ //searchParams.set("Method", "Interval Halving");
+ //searchParams.set("Method", "Polynomial");
+ //searchParams.set("Method", "NonlinearCG");
+ //searchParams.set("Method", "Quadratic");
+ //searchParams.set("Method", "More'-Thuente");
+
+ // Sublist for direction
+ Teuchos::ParameterList& dirParams = nlParams.sublist("Direction");
+// dirParams.set("Method", "Modified-Newton");
+// Teuchos::ParameterList& newtonParams = dirParams.sublist("Modified-Newton");
+// newtonParams.set("Max Age of Jacobian", 2);
+ dirParams.set("Method", "Newton");
+ Teuchos::ParameterList& newtonParams = dirParams.sublist("Newton");
+ //newtonParams.set("Forcing Term Method", "Constant");
+ //newtonParams.set("Forcing Term Method", "Type 1");
+ newtonParams.set("Forcing Term Method", "Type 2");
+ //newtonParams.set("Forcing Term Minimum Tolerance", 1.0e-4);
+ //newtonParams.set("Forcing Term Maximum Tolerance", 0.1);
+ //dirParams.set("Method", "Steepest Descent");
+ //Teuchos::ParameterList& sdParams = dirParams.sublist("Steepest Descent");
+ //sdParams.set("Scaling Type", "None");
+ //sdParams.set("Scaling Type", "2-Norm");
+ //sdParams.set("Scaling Type", "Quadratic Model Min");
+ //dirParams.set("Method", "NonlinearCG");
+ //Teuchos::ParameterList& nlcgParams = dirParams.sublist("Nonlinear CG");
+ //nlcgParams.set("Restart Frequency", 2000);
+ //nlcgParams.set("Precondition", "On");
+ //nlcgParams.set("Orthogonalize", "Polak-Ribiere");
+ //nlcgParams.set("Orthogonalize", "Fletcher-Reeves");
+
+
+ lsParams.set("Aztec Solver", "GMRES");
+// lsParams.set("Max Iterations", 800);
+ lsParams.set("Max Iterations", 1600);
+ lsParams.set("Tolerance", 1e-4);
+ // lsParams.set("Tolerance", 1e-2);
+ lsParams.set("Preconditioner", "None");
+ //lsParams.set("Preconditioner", "Ifpack");
+ lsParams.set("Size of Krylov Subspace", 200);
+ lsParams.set("Max Age Of Prec", 5);
+
+ // Create the interface between the test problem and the nonlinear solver
+ // This is created by the user using inheritance of the abstract base class:
+ // NOX_Epetra_Interface
+ //interface = Teuchos::rcp(new Problem_Interface(*nelems, statevector, Comm,
+ // blackbox_res, residualFunction));
+}
+
+}
Added: branches/implicit/src/core_sw/nox_interface.C
===================================================================
--- branches/implicit/src/core_sw/nox_interface.C (rev 0)
+++ branches/implicit/src/core_sw/nox_interface.C 2010-08-05 21:04:06 UTC (rev 463)
@@ -0,0 +1,27 @@
+#include nox_interface.h
+
+
+Problem_Interface::Problem_Interface(*nelems, statevector, Comm,
+ blackbox_res, residualFunction)
+{ }
+ ~Problem_Interface()
+{ }
+
+ //! Compute and return F
+ bool Problem_Interface::computeF(const Epetra_Vector& x, Epetra_Vector& FVec,
+ FillType flag = Residual)
+{ }
+
+ //! Compute an explicit Jacobian
+ bool Problem_Interface::computeJacobian(const Epetra_Vector& x, Epetra_Operator& Jac)
+{ }
+
+ //! Compute the Epetra_RowMatrix M, that will be used by the Aztec preconditioner instead of the Jacobian. This is used when there is no explicit Jacobian present (i.e. Matrix-Free Newton-Krylov). This MUST BE an Epetra_RowMatrix since the Aztec preconditioners need to know the sparsity pattern of the matrix. Returns true if computation was successful.
+ bool Problem_Interface::computePrecMatrix(const Epetra_Vector& x, Epetra_RowMatrix& M)
+{ }
+
+ //! Computes a user supplied preconditioner based on input vector x. Returns true if computation was successful.
+ bool Problem_Interface::computePreconditioner(const Epetra_Vector& x, Epetra_Operator& M,
+ Teuchos::ParameterList* precParams = 0)
+{ }
+
Added: branches/implicit/src/core_sw/nox_interface.h
===================================================================
--- branches/implicit/src/core_sw/nox_interface.h (rev 0)
+++ branches/implicit/src/core_sw/nox_interface.h 2010-08-05 21:04:06 UTC (rev 463)
@@ -0,0 +1,43 @@
+
+#ifndef Problem_Interface_H
+#define Problem_Interface_H
+
+// Interface to the NLS_PetraGroup to provide for
+// residual and matrix fill routines.
+
+// ---------- Standard Includes ----------
+#include <iostream>
+#include "Epetra_Vector.h"
+#include "Epetra_Operator.h"
+#include "Epetra_RowMatrix.h"
+#include "NOX_Epetra_Interface_Required.H" // base class
+#include "NOX_Epetra_Interface_Jacobian.H" // base class
+#include "NOX_Epetra_Interface_Preconditioner.H" // base class
+
+class Problem_Interface : public NOX::Epetra::Interface::Required,
+ public NOX::Epetra::Interface::Jacobian,
+ public NOX::Epetra::Interface::Preconditioner
+{
+public:
+ Problem_Interface(int nelems, double* statevector,
+ const Epetra_Comm& comm_,
+ void* blackbox_res, void (*residualFunction)(double *, double *, int, void *, double));
+ ~Problem_Interface();
+
+ //! Compute and return F
+ bool computeF(const Epetra_Vector& x, Epetra_Vector& FVec,
+ FillType flag = Residual);
+
+ //! Compute an explicit Jacobian
+ bool computeJacobian(const Epetra_Vector& x, Epetra_Operator& Jac);
+
+ //! Compute the Epetra_RowMatrix M, that will be used by the Aztec preconditioner instead of the Jacobian. This is used when there is no explicit Jacobian present (i.e. Matrix-Free Newton-Krylov). This MUST BE an Epetra_RowMatrix since the Aztec preconditioners need to know the sparsity pattern of the matrix. Returns true if computation was successful.
+ bool computePrecMatrix(const Epetra_Vector& x, Epetra_RowMatrix& M);
+
+ //! Computes a user supplied preconditioner based on input vector x. Returns true if computation was successful.
+ bool computePreconditioner(const Epetra_Vector& x, Epetra_Operator& M,
+ Teuchos::ParameterList* precParams = 0);
+
+};
+
+#endif
</font>
</pre>