From b5d1896d30a37a282d47ac7ffd023ee7d0adb131 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 13 Sep 2014 09:47:07 -0300 Subject: [PATCH 01/12] EVOLUTION: Slightly rewrite new_time_step(). Signed-off-by: Grzegorz Kowal --- src/evolution.F90 | 54 +++++++++++++++++++++++------------------------ 1 file changed, 26 insertions(+), 28 deletions(-) diff --git a/src/evolution.F90 b/src/evolution.F90 index f07d0ed..d2562f0 100644 --- a/src/evolution.F90 +++ b/src/evolution.F90 @@ -334,18 +334,17 @@ module evolution ! implicit none -! input variables +! subroutine arguments ! - real(kind=8), intent(in) :: dtnext + real(kind=8), intent(in) :: dtnext ! local pointers ! - type(block_data), pointer :: pblock + type(block_data), pointer :: pdata ! local variables ! - integer :: iret - integer(kind=4) :: lev + integer :: iret, mlev real(kind=8) :: cm, dx_min ! local parameters @@ -357,50 +356,53 @@ module evolution ! reset the maximum speed, and the highest level ! cmax = eps - lev = 1 + mlev = 1 + +! assign pdata with the first block on the data block list +! + pdata => list_data ! iterate over all data blocks in order to find the maximum speed among them -! and the highest level which is required to obtain the spatial step +! and the highest level which is required to obtain the minimum spacial step ! - pblock => list_data - do while (associated(pblock)) + do while (associated(pdata)) -! find the maximum level occupied by blocks (can be smaller than toplev) +! update the maximum level ! - lev = max(lev, pblock%meta%level) + mlev = max(mlev, pdata%meta%level) -! obtain the maximum speed for the current block +! get the maximum speed for the current block ! - cm = maxspeed(pblock%q(:,:,:,:)) + cm = maxspeed(pdata%q(:,:,:,:)) -! compare global and local maximum speeds +! update the maximum speed ! cmax = max(cmax, cm) -! assiociate the pointer with the next block +! assign pdata to the next block ! - pblock => pblock%next + pdata => pdata%next - end do + end do ! over data blocks #ifdef MPI -! find maximum speed in the system from all processors +! reduce maximum speed and level over all processors ! call reduce_maximum_real (cmax, iret) - call reduce_maximum_integer(lev , iret) + call reduce_maximum_integer(mlev, iret) #endif /* MPI */ -! calculate squared cmax +! calculate the squared cmax ! cmax2 = cmax * cmax -! find the smallest spatial step +! find the smallest spacial step ! #if NDIMS == 2 - dx_min = min(adx(lev), ady(lev)) + dx_min = min(adx(mlev), ady(mlev)) #endif /* NDIMS == 2 */ #if NDIMS == 3 - dx_min = min(adx(lev), ady(lev), adz(lev)) + dx_min = min(adx(mlev), ady(mlev), adz(mlev)) #endif /* NDIMS == 3 */ ! calculate the new time step @@ -408,13 +410,9 @@ module evolution dtn = cfl * dx_min / max(cmax & + 2.0d+00 * max(viscosity, resistivity) / dx_min, eps) -! substitute the new time step -! - dt = dtn - ! round the time ! - if (dtnext > 0.0d+00) dt = min(dt, dtnext) + if (dtnext > 0.0d+00) dt = min(dtn, dtnext) !------------------------------------------------------------------------------- ! From 2c359fd5f28133cab1f4eafd4a4eb54adea60d0a Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 13 Sep 2014 09:53:40 -0300 Subject: [PATCH 02/12] EVOLUTION: Slightly rewrite evolve_euler(). Signed-off-by: Grzegorz Kowal --- src/evolution.F90 | 44 +++++++++++++++++++++++++++----------------- 1 file changed, 27 insertions(+), 17 deletions(-) diff --git a/src/evolution.F90 b/src/evolution.F90 index d2562f0..6ce8887 100644 --- a/src/evolution.F90 +++ b/src/evolution.F90 @@ -432,6 +432,12 @@ module evolution ! Subroutine advances the solution by one time step using the 1st order ! Euler integration method. ! +! References: +! +! [1] Press, W. H, Teukolsky, S. A., Vetterling, W. T., Flannery, B. P., +! "Numerical Recipes in Fortran", +! Cambridge University Press, Cambridge, 1992 +! !=============================================================================== ! subroutine evolve_euler() @@ -453,7 +459,7 @@ module evolution ! local pointers ! - type(block_data), pointer :: pblock + type(block_data), pointer :: pdata ! local arrays ! @@ -461,40 +467,44 @@ module evolution ! !------------------------------------------------------------------------------- ! -! update fluxes for the first step of the RK2 integration +! update fluxes ! call update_fluxes() -! update the solution using numerical fluxes stored in the data blocks +! assign pdata with the first block on the data block list ! - pblock => list_data - do while (associated(pblock)) + pdata => list_data -! calculate variable increment for the current block +! iterate over all data blocks ! - call update_increment(pblock, du(:,:,:,:)) + do while (associated(pdata)) -! add source terms +! calculate the variable increment ! - call update_sources(pblock, du(:,:,:,:)) + call update_increment(pdata, du(1:nv,1:im,1:jm,1:km)) -! update the solution for the fluid variables +! add the source terms ! - pblock%u0(1:nv,:,:,:) = pblock%u0(1:nv,:,:,:) + dt * du(1:nv,:,:,:) + call update_sources(pdata, du(1:nv,1:im,1:jm,1:km)) + +! update the solution +! + pdata%u0(1:nv,1:im,1:jm,1:km) = pdata%u0(1:nv,1:im,1:jm,1:km) & + + dt * du(1:nv,1:im,1:jm,1:km) ! update the conservative variable pointer ! - pblock%u => pblock%u0 + pdata%u => pdata%u0 -! update ψ by its source term +! update ψ with its source term ! - if (ibp > 0) pblock%u(ibp,:,:,:) = decay * pblock%u(ibp,:,:,:) + if (ibp > 0) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) -! assign pointer to the next block +! assign pdata to the next block ! - pblock => pblock%next + pdata => pdata%next - end do + end do ! over data blocks ! update primitive variables ! From 984fec49ed10b956363ea1f6f5878b2ad9904292 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 13 Sep 2014 10:00:16 -0300 Subject: [PATCH 03/12] EVOLUTION: Slightly rewrite evolve_rk2(). Signed-off-by: Grzegorz Kowal --- src/evolution.F90 | 74 +++++++++++++++++++++++++++-------------------- 1 file changed, 42 insertions(+), 32 deletions(-) diff --git a/src/evolution.F90 b/src/evolution.F90 index 6ce8887..2dbaa96 100644 --- a/src/evolution.F90 +++ b/src/evolution.F90 @@ -498,7 +498,8 @@ module evolution ! update ψ with its source term ! - if (ibp > 0) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) + if (ibp > 0) pdata%u(ibp,1:im,1:jm,1:km) = & + decay * pdata%u(ibp,1:im,1:jm,1:km) ! assign pdata to the next block ! @@ -553,7 +554,7 @@ module evolution ! local pointers ! - type(block_data), pointer :: pblock + type(block_data), pointer :: pdata ! local arrays ! @@ -561,36 +562,40 @@ module evolution ! !------------------------------------------------------------------------------- ! -! update fluxes for the first step of the RK2 integration +! update fluxes ! call update_fluxes() -! update the solution using numerical fluxes stored in the data blocks +! assign pdata with the first block on the data block list ! - pblock => list_data - do while (associated(pblock)) + pdata => list_data -! calculate variable increment for the current block +! iterate over all data blocks ! - call update_increment(pblock, du(:,:,:,:)) + do while (associated(pdata)) -! add source terms +! calculate the variable increment ! - call update_sources(pblock, du(:,:,:,:)) + call update_increment(pdata, du(1:nv,1:im,1:jm,1:km)) -! update the solution for the fluid variables +! add the source terms ! - pblock%u1(1:nv,:,:,:) = pblock%u0(1:nv,:,:,:) + dt * du(1:nv,:,:,:) + call update_sources(pdata, du(1:nv,1:im,1:jm,1:km)) + +! update the intermediate solution +! + pdata%u1(1:nv,1:im,1:jm,1:km) = pdata%u0(1:nv,1:im,1:jm,1:km) & + + dt * du(1:nv,1:im,1:jm,1:km) ! update the conservative variable pointer ! - pblock%u => pblock%u1 + pdata%u => pdata%u1 -! assign pointer to the next block +! assign pdata to the next block ! - pblock => pblock%next + pdata => pdata%next - end do + end do ! over data blocks ! update primitive variables ! @@ -600,41 +605,46 @@ module evolution ! call boundary_variables() -! update fluxes for the second step of the RK2 integration +! update fluxes from the intermediate stage ! call update_fluxes() -! update the solution using numerical fluxes stored in the data blocks +! assign pdata with the first block on the data block list ! - pblock => list_data - do while (associated(pblock)) + pdata => list_data -! calculate variable increment for the current block +! iterate over all data blocks ! - call update_increment(pblock, du(:,:,:,:)) + do while (associated(pdata)) -! add source terms +! calculate the variable increment ! - call update_sources(pblock, du(:,:,:,:)) + call update_increment(pdata, du(1:nv,1:im,1:jm,1:km)) -! update the solution for the fluid variables +! add the source terms ! - pblock%u0(1:nv,:,:,:) = 0.5d+00 * (pblock%u0(1:nv,:,:,:) & - + pblock%u1(1:nv,:,:,:) + dt * du(1:nv,:,:,:)) + call update_sources(pdata, du(1:nv,1:im,1:jm,1:km)) + +! update the final solution +! + pdata%u0(1:nv,1:im,1:jm,1:km) = 0.5d+00 * (pdata%u0(1:nv,1:im,1:jm,1:km) & + + pdata%u1(1:nv,1:im,1:jm,1:km) & + + dt * du(1:nv,1:im,1:jm,1:km)) ! update the conservative variable pointer ! - pblock%u => pblock%u0 + pdata%u => pdata%u0 -! update ψ by its source term +! update ψ with its source term ! - if (ibp > 0) pblock%u(ibp,:,:,:) = decay * pblock%u(ibp,:,:,:) + if (ibp > 0) pdata%u(ibp,1:im,1:jm,1:km) = & + decay * pdata%u(ibp,1:im,1:jm,1:km) ! assign pointer to the next block ! - pblock => pblock%next + pdata => pdata%next - end do + end do ! over data blocks ! update primitive variables ! From f11cdf70975164f36c2bbb175ad3fdcdb4275a98 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 13 Sep 2014 10:10:54 -0300 Subject: [PATCH 04/12] EVOLUTION: Slightly rewrite evolve_rk3(). Signed-off-by: Grzegorz Kowal --- src/evolution.F90 | 113 ++++++++++++++++++++++++++-------------------- 1 file changed, 63 insertions(+), 50 deletions(-) diff --git a/src/evolution.F90 b/src/evolution.F90 index 2dbaa96..f47a25b 100644 --- a/src/evolution.F90 +++ b/src/evolution.F90 @@ -640,7 +640,7 @@ module evolution if (ibp > 0) pdata%u(ibp,1:im,1:jm,1:km) = & decay * pdata%u(ibp,1:im,1:jm,1:km) -! assign pointer to the next block +! assign pdata to the next block ! pdata => pdata%next @@ -896,7 +896,7 @@ module evolution ! local pointers ! - type(block_data), pointer :: pblock + type(block_data), pointer :: pdata ! local variables ! @@ -915,40 +915,44 @@ module evolution ! !! 1st substep of integration !! -! prepare fractional time step +! prepare the fractional time step ! ds = dt -! update fluxes for the first step of the RK2 integration +! update fluxes ! call update_fluxes() -! update the solution using numerical fluxes stored in the data blocks +! assign pdata with the first block on the data block list ! - pblock => list_data - do while (associated(pblock)) + pdata => list_data -! calculate variable increment for the current block +! iterate over all data blocks ! - call update_increment(pblock, du(:,:,:,:)) + do while (associated(pdata)) -! add source terms +! calculate the variable increment ! - call update_sources(pblock, du(:,:,:,:)) + call update_increment(pdata, du(1:nv,1:im,1:jm,1:km)) -! update the solution for the fluid variables +! add the source terms ! - pblock%u1(1:nv,:,:,:) = pblock%u0(1:nv,:,:,:) + ds * du(1:nv,:,:,:) + call update_sources(pdata, du(1:nv,1:im,1:jm,1:km)) + +! update the first intermediate solution +! + pdata%u1(1:nv,1:im,1:jm,1:km) = pdata%u0(1:nv,1:im,1:jm,1:km) & + + ds * du(1:nv,1:im,1:jm,1:km) ! update the conservative variable pointer ! - pblock%u => pblock%u1 + pdata%u => pdata%u1 -! assign pointer to the next block +! assign pdata to the next block ! - pblock => pblock%next + pdata => pdata%next - end do + end do ! over data blocks ! update primitive variables ! @@ -960,37 +964,41 @@ module evolution !! 2nd substep of integration !! -! prepare fractional time step +! prepare the fractional time step ! ds = f22 * dt -! update fluxes for the first step of the RK2 integration +! update fluxes ! call update_fluxes() -! update the solution using numerical fluxes stored in the data blocks +! assign pdata with the first block on the data block list ! - pblock => list_data - do while (associated(pblock)) + pdata => list_data -! calculate variable increment for the current block +! iterate over all data blocks ! - call update_increment(pblock, du(:,:,:,:)) + do while (associated(pdata)) -! add source terms +! calculate the variable increment ! - call update_sources(pblock, du(:,:,:,:)) + call update_increment(pdata, du(1:nv,1:im,1:jm,1:km)) -! update the solution for the fluid variables +! add the source terms ! - pblock%u1(1:nv,:,:,:) = f21 * pblock%u0(1:nv,:,:,:) & - + f22 * pblock%u1(1:nv,:,:,:) + ds * du(1:nv,:,:,:) + call update_sources(pdata, du(1:nv,1:im,1:jm,1:km)) -! assign pointer to the next block +! update the second intermediate solution ! - pblock => pblock%next + pdata%u1(1:nv,1:im,1:jm,1:km) = f21 * pdata%u0(1:nv,1:im,1:jm,1:km) & + + f22 * pdata%u1(1:nv,1:im,1:jm,1:km) & + + ds * du(1:nv,1:im,1:jm,1:km) - end do +! assign pdata to the next block +! + pdata => pdata%next + + end do ! over data blocks ! update primitive variables ! @@ -1002,45 +1010,50 @@ module evolution !! 3rd substep of integration !! -! prepare fractional time step +! prepare the fractional time step ! ds = f32 * dt -! update fluxes for the second step of the RK2 integration +! update fluxes ! call update_fluxes() -! update the solution using numerical fluxes stored in the data blocks +! assign pdata with the first block on the data block list ! - pblock => list_data - do while (associated(pblock)) + pdata => list_data -! calculate variable increment for the current block +! iterate over all data blocks ! - call update_increment(pblock, du(:,:,:,:)) + do while (associated(pdata)) -! add source terms +! calculate the variable increment ! - call update_sources(pblock, du(:,:,:,:)) + call update_increment(pdata, du(1:nv,1:im,1:jm,1:km)) -! update the solution for the fluid variables +! add the source terms ! - pblock%u0(1:nv,:,:,:) = f31 * pblock%u0(1:nv,:,:,:) & - + f32 * pblock%u1(1:nv,:,:,:) + ds * du(1:nv,:,:,:) + call update_sources(pdata, du(1:nv,1:im,1:jm,1:km)) + +! update the final solution +! + pdata%u0(1:nv,1:im,1:jm,1:km) = f31 * pdata%u0(1:nv,1:im,1:jm,1:km) & + + f32 * pdata%u1(1:nv,1:im,1:jm,1:km) & + + ds * du(1:nv,1:im,1:jm,1:km) ! update the conservative variable pointer ! - pblock%u => pblock%u0 + pdata%u => pdata%u0 -! update ψ by its source term +! update ψ with its source term ! - if (ibp > 0) pblock%u(ibp,:,:,:) = decay * pblock%u(ibp,:,:,:) + if (ibp > 0) pdata%u(ibp,1:im,1:jm,1:km) = & + decay * pdata%u(ibp,1:im,1:jm,1:km) -! assign pointer to the next block +! assign pdata to the next block ! - pblock => pblock%next + pdata => pdata%next - end do + end do ! over data blocks ! update primitive variables ! From 8f86be865fa53094a9653be109fcba3d9b461ef5 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 13 Sep 2014 10:15:34 -0300 Subject: [PATCH 05/12] EVOLUTION: Slightly rewrite update_fluxes(). Signed-off-by: Grzegorz Kowal --- src/evolution.F90 | 30 ++++++++++++++++++------------ 1 file changed, 18 insertions(+), 12 deletions(-) diff --git a/src/evolution.F90 b/src/evolution.F90 index f47a25b..5a6cbf6 100644 --- a/src/evolution.F90 +++ b/src/evolution.F90 @@ -1641,6 +1641,8 @@ module evolution ! use blocks , only : block_data, list_data use coordinates , only : adx, ady, adz + use coordinates , only : im, jm, km + use equations , only : nv ! local variables are not implicit by default ! @@ -1648,7 +1650,7 @@ module evolution ! local pointers ! - type(block_data), pointer :: pblock + type(block_data), pointer :: pdata ! local vectors ! @@ -1660,31 +1662,35 @@ module evolution ! !------------------------------------------------------------------------------- ! +! assign pdata with the first block on the data block list +! + pdata => list_data + ! iterate over all data blocks ! - pblock => list_data - do while (associated(pblock)) + do while (associated(pdata)) ! obtain dx, dy, and dz for the current block ! - dx(1) = adx(pblock%meta%level) - dx(2) = ady(pblock%meta%level) - dx(3) = adz(pblock%meta%level) + dx(1) = adx(pdata%meta%level) + dx(2) = ady(pdata%meta%level) + dx(3) = adz(pdata%meta%level) -! update the flux for the current block +! update fluxes for the current block ! do n = 1, NDIMS - call update_flux(n, dx(n), pblock%q(:,:,:,:), pblock%f(n,:,:,:,:)) + call update_flux(n, dx(n), pdata%q(1:nv,1:im,1:jm,1:km) & + , pdata%f(n,1:nv,1:im,1:jm,1:km)) end do -! assign pointer to the next block +! assign pdata to the next block ! - pblock => pblock%next + pdata => pdata%next - end do + end do ! over data blocks ! correct the numerical fluxes of the blocks which have neighbours at higher -! level +! levels ! call boundary_fluxes() From 2960aea9c4a38b75a1b36a1a811a60762e73870f Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 13 Sep 2014 10:17:02 -0300 Subject: [PATCH 06/12] EVOLUTION: Update some comments. Signed-off-by: Grzegorz Kowal --- src/evolution.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/evolution.F90 b/src/evolution.F90 index 5a6cbf6..d732177 100644 --- a/src/evolution.F90 +++ b/src/evolution.F90 @@ -837,7 +837,7 @@ module evolution ! pdata%u => pdata%u0 -! update ψ by its source term +! update ψ with its source term ! if (ibp > 0) pdata%u(ibp,1:im,1:jm,1:km) = & decay * pdata%u(ibp,1:im,1:jm,1:km) @@ -1292,7 +1292,7 @@ module evolution ! pdata%u => pdata%u0 -! update ψ by its source term +! update ψ with its source term ! if (ibp > 0) pdata%u(ibp,1:im,1:jm,1:km) = & decay * pdata%u(ibp,1:im,1:jm,1:km) @@ -1595,7 +1595,7 @@ module evolution + a55 * pdata%u0(1:nv,1:im,1:jm,1:km) & + ds * du(1:nv,1:im,1:jm,1:km) -! update ψ by its source term +! update ψ with its source term ! if (ibp > 0) pdata%u(ibp,1:im,1:jm,1:km) = & decay * pdata%u(ibp,1:im,1:jm,1:km) From d1e78eb0b7977c3a61bf806ae511633b4739791c Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 13 Sep 2014 10:19:25 -0300 Subject: [PATCH 07/12] MESH: Replace pblock with pmeta in some subroutines. Signed-off-by: Grzegorz Kowal --- src/mesh.F90 | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/mesh.F90 b/src/mesh.F90 index 2c3524f..6537f28 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -1031,11 +1031,11 @@ module mesh ! ! Arguments: ! -! pblock - the input meta block; +! pmeta - the input meta block; ! !=============================================================================== ! - subroutine prolong_block(pblock) + subroutine prolong_block(pmeta) ! import external procedures and variables ! @@ -1051,7 +1051,7 @@ module mesh ! input arguments ! - type(block_meta), pointer, intent(inout) :: pblock + type(block_meta), pointer, intent(inout) :: pmeta ! local variables ! @@ -1083,7 +1083,7 @@ module mesh ! assign the pdata pointer ! - pdata => pblock%data + pdata => pmeta%data ! prepare dimensions ! @@ -1175,7 +1175,7 @@ module mesh ! assign pointer to the current child ! - pchild => pblock%child(p)%ptr + pchild => pmeta%child(p)%ptr ! obtain the position of child in the parent block ! @@ -1236,11 +1236,11 @@ module mesh ! ! Arguments: ! -! pblock - the input meta block; +! pmeta - the input meta block; ! !=============================================================================== ! - subroutine restrict_block(pblock) + subroutine restrict_block(pmeta) ! import external procedures and variables ! @@ -1256,7 +1256,7 @@ module mesh ! subroutine arguments ! - type(block_meta), pointer, intent(inout) :: pblock + type(block_meta), pointer, intent(inout) :: pmeta ! local variables ! @@ -1282,7 +1282,7 @@ module mesh ! assign the parent data pointer ! - pparent => pblock%data + pparent => pmeta%data ! iterate over all children ! @@ -1290,7 +1290,7 @@ module mesh ! assign a pointer to the current child ! - pchild => pblock%child(p)%ptr%data + pchild => pmeta%child(p)%ptr%data ! obtain the child position in the parent block ! From f89d2133864252b6bc7ef1e7695ff82e940fdedd Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sun, 14 Sep 2014 20:38:53 -0300 Subject: [PATCH 08/12] EVOLUTION: Add support for profiling. Signed-off-by: Grzegorz Kowal --- src/evolution.F90 | 180 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 180 insertions(+) diff --git a/src/evolution.F90 b/src/evolution.F90 index d732177..ab6be92 100644 --- a/src/evolution.F90 +++ b/src/evolution.F90 @@ -31,10 +31,22 @@ ! module evolution +#ifdef PROFILE +! import external subroutines +! + use timers, only : set_timer, start_timer, stop_timer +#endif /* PROFILE */ + ! module variables are not implicit by default ! implicit none +#ifdef PROFILE +! timer indices +! + integer , save :: imi, ima, imt, imu, imf, imn, imv +#endif /* PROFILE */ + ! pointer to the temporal integration subroutine ! procedure(evolve_euler), pointer, save :: evolve => null() @@ -116,6 +128,22 @@ module evolution ! !------------------------------------------------------------------------------- ! +#ifdef PROFILE +! set timer descriptions +! + call set_timer('evolution:: initialization' , imi) + call set_timer('evolution:: solution advance', ima) + call set_timer('evolution:: new time step' , imt) + call set_timer('evolution:: solution update' , imu) + call set_timer('evolution:: flux update' , imf) + call set_timer('evolution:: increment update', imn) + call set_timer('evolution:: variable update' , imv) + +! start accounting time for module initialization/finalization +! + call start_timer(imi) +#endif /* PROFILE */ + ! get the integration method and the value of the CFL coefficient ! call get_parameter_string ("time_advance", integration) @@ -185,6 +213,12 @@ module evolution end if +#ifdef PROFILE +! stop accounting time for module initialization/finalization +! + call stop_timer(imi) +#endif /* PROFILE */ + !------------------------------------------------------------------------------- ! end subroutine initialize_evolution @@ -213,9 +247,23 @@ module evolution integer, intent(inout) :: iret ! !------------------------------------------------------------------------------- +! +#ifdef PROFILE +! start accounting time for module initialization/finalization +! + call start_timer(imi) +#endif /* PROFILE */ + +! nullify pointer to integration subroutine ! nullify(evolve) +#ifdef PROFILE +! stop accounting time for module initialization/finalization +! + call stop_timer(imi) +#endif /* PROFILE */ + !------------------------------------------------------------------------------- ! end subroutine finalize_evolution @@ -256,6 +304,12 @@ module evolution ! !------------------------------------------------------------------------------- ! +#ifdef PROFILE +! start accounting time for solution advance +! + call start_timer(ima) +#endif /* PROFILE */ + ! find new time step ! call new_time_step(dtnext) @@ -296,6 +350,12 @@ module evolution end if ! toplev > 1 +#ifdef PROFILE +! stop accounting time for solution advance +! + call stop_timer(ima) +#endif /* PROFILE */ + !------------------------------------------------------------------------------- ! end subroutine advance @@ -353,6 +413,12 @@ module evolution ! !------------------------------------------------------------------------------- ! +#ifdef PROFILE +! start accounting time for new time step estimation +! + call start_timer(imt) +#endif /* PROFILE */ + ! reset the maximum speed, and the highest level ! cmax = eps @@ -414,6 +480,12 @@ module evolution ! if (dtnext > 0.0d+00) dt = min(dtn, dtnext) +#ifdef PROFILE +! stop accounting time for new time step estimation +! + call stop_timer(imt) +#endif /* PROFILE */ + !------------------------------------------------------------------------------- ! end subroutine new_time_step @@ -467,6 +539,12 @@ module evolution ! !------------------------------------------------------------------------------- ! +#ifdef PROFILE +! start accounting time for one step update +! + call start_timer(imu) +#endif /* PROFILE */ + ! update fluxes ! call update_fluxes() @@ -515,6 +593,12 @@ module evolution ! call boundary_variables() +#ifdef PROFILE +! stop accounting time for one step update +! + call stop_timer(imu) +#endif /* PROFILE */ + !------------------------------------------------------------------------------- ! end subroutine evolve_euler @@ -562,6 +646,12 @@ module evolution ! !------------------------------------------------------------------------------- ! +#ifdef PROFILE +! start accounting time for one step update +! + call start_timer(imu) +#endif /* PROFILE */ + ! update fluxes ! call update_fluxes() @@ -654,6 +744,12 @@ module evolution ! call boundary_variables() +#ifdef PROFILE +! stop accounting time for one step update +! + call stop_timer(imu) +#endif /* PROFILE */ + !------------------------------------------------------------------------------- ! end subroutine evolve_rk2 @@ -714,6 +810,12 @@ module evolution ! !------------------------------------------------------------------------------- ! +#ifdef PROFILE +! start accounting time for one step update +! + call start_timer(imu) +#endif /* PROFILE */ + ! prepare things which don't change later ! if (first) then @@ -856,6 +958,12 @@ module evolution ! call boundary_variables() +#ifdef PROFILE +! stop accounting time for one step update +! + call stop_timer(imu) +#endif /* PROFILE */ + !------------------------------------------------------------------------------- ! end subroutine evolve_ssprk2 @@ -913,6 +1021,12 @@ module evolution ! !------------------------------------------------------------------------------- ! +#ifdef PROFILE +! start accounting time for one step update +! + call start_timer(imu) +#endif /* PROFILE */ + !! 1st substep of integration !! ! prepare the fractional time step @@ -1063,6 +1177,12 @@ module evolution ! call boundary_variables() +#ifdef PROFILE +! stop accounting time for one step update +! + call stop_timer(imu) +#endif /* PROFILE */ + !------------------------------------------------------------------------------- ! end subroutine evolve_rk3 @@ -1121,6 +1241,12 @@ module evolution ! !------------------------------------------------------------------------------- ! +#ifdef PROFILE +! start accounting time for one step update +! + call start_timer(imu) +#endif /* PROFILE */ + != 1st step: U(1) = U(n) + 1/2 dt F[U(n)] ! ! calculate the fractional time step @@ -1311,6 +1437,12 @@ module evolution ! call boundary_variables() +#ifdef PROFILE +! stop accounting time for one step update +! + call stop_timer(imu) +#endif /* PROFILE */ + !------------------------------------------------------------------------------- ! end subroutine evolve_ssprk34 @@ -1377,6 +1509,12 @@ module evolution ! !------------------------------------------------------------------------------- ! +#ifdef PROFILE +! start accounting time for one step update +! + call start_timer(imu) +#endif /* PROFILE */ + != 1st step: U(1) = U(n) + b1 dt F[U(n)] ! ! calculate the fractional time step @@ -1614,6 +1752,12 @@ module evolution ! call boundary_variables() +#ifdef PROFILE +! stop accounting time for one step update +! + call stop_timer(imu) +#endif /* PROFILE */ + !------------------------------------------------------------------------------- ! end subroutine evolve_ssprk35 @@ -1662,6 +1806,12 @@ module evolution ! !------------------------------------------------------------------------------- ! +#ifdef PROFILE +! start accounting time for fluxe update +! + call start_timer(imf) +#endif /* PROFILE */ + ! assign pdata with the first block on the data block list ! pdata => list_data @@ -1694,6 +1844,12 @@ module evolution ! call boundary_fluxes() +#ifdef PROFILE +! stop accounting time for flux update +! + call stop_timer(imf) +#endif /* PROFILE */ + !------------------------------------------------------------------------------- ! end subroutine update_fluxes @@ -1738,6 +1894,12 @@ module evolution ! !------------------------------------------------------------------------------- ! +#ifdef PROFILE +! start accounting time for increment update +! + call start_timer(imn) +#endif /* PROFILE */ + ! reset the increment array du ! du(:,:,:,:) = 0.0d+00 @@ -1773,6 +1935,12 @@ module evolution end do #endif /* NDIMS == 3 */ +#ifdef PROFILE +! stop accounting time for increment update +! + call stop_timer(imn) +#endif /* PROFILE */ + !------------------------------------------------------------------------------- ! end subroutine update_increment @@ -1811,6 +1979,12 @@ module evolution ! !------------------------------------------------------------------------------- ! +#ifdef PROFILE +! start accounting time for variable update +! + call start_timer(imv) +#endif /* PROFILE */ + ! associate the pointer with the first block on the data block list ! pdata => list_data @@ -1837,6 +2011,12 @@ module evolution end do +#ifdef PROFILE +! stop accounting time for variable update +! + call stop_timer(imv) +#endif /* PROFILE */ + !------------------------------------------------------------------------------- ! end subroutine update_variables From 7c3a621096821a1723ea16f9456e51745bed79f0 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 17 Sep 2014 08:34:47 -0300 Subject: [PATCH 09/12] MESH: Redistribute all variable arrays in redistribute_blocks(). Signed-off-by: Grzegorz Kowal --- src/mesh.F90 | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/mesh.F90 b/src/mesh.F90 index 6537f28..6eb6758 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -878,7 +878,7 @@ module mesh ! local buffer for data block exchange ! - real(kind=8) , dimension(2,nv,im,jm,km) :: rbuf + real(kind=8) , dimension(3,nv,im,jm,km) :: rbuf #endif /* MPI */ !------------------------------------------------------------------------------- @@ -931,8 +931,9 @@ module mesh ! copy data to buffer ! - rbuf(1,:,:,:,:) = pmeta%data%u(:,:,:,:) - rbuf(2,:,:,:,:) = pmeta%data%q(:,:,:,:) + rbuf(1,:,:,:,:) = pmeta%data%q (:,:,:,:) + rbuf(2,:,:,:,:) = pmeta%data%u0(:,:,:,:) + rbuf(3,:,:,:,:) = pmeta%data%u1(:,:,:,:) ! send data ! @@ -961,8 +962,9 @@ module mesh ! coppy the buffer to data block ! - pmeta%data%u(:,:,:,:) = rbuf(1,:,:,:,:) - pmeta%data%q(:,:,:,:) = rbuf(2,:,:,:,:) + pmeta%data%q (:,:,:,:) = rbuf(1,:,:,:,:) + pmeta%data%u0(:,:,:,:) = rbuf(2,:,:,:,:) + pmeta%data%u1(:,:,:,:) = rbuf(3,:,:,:,:) end if ! nproc == n From e00dd584e40cf0b31f83dcc24928926fc6d54621 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 19 Sep 2014 07:17:28 -0300 Subject: [PATCH 10/12] OPERATORS: Add gradient operator. Signed-off-by: Grzegorz Kowal --- src/operators.F90 | 53 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) diff --git a/src/operators.F90 b/src/operators.F90 index b94b4ef..98920a2 100644 --- a/src/operators.F90 +++ b/src/operators.F90 @@ -187,6 +187,59 @@ module operators ! !=============================================================================== ! +! subroutine GRADIENT: +! ------------------- +! +! Subroutine calculates the cell centered gradient of the input scalar field. +! +! grad(U) = ∇ U = [ ∂x U, ∂y U, ∂z U ] +! +! Arguments: +! +! dh - the spacial intervals in all direction; +! u - the input scalar field; +! v - the output gradient vector field; +! +!=============================================================================== +! + subroutine gradient(dh, u, v) + +! local variables are not implicit by default +! + implicit none + +! input and output variables +! + real(kind=8), dimension(3) , intent(in) :: dh + real(kind=8), dimension(:,:,:) , intent(in) :: u + real(kind=8), dimension(:,:,:,:), intent(out) :: v + +! local variables +! + integer :: dir +! +!------------------------------------------------------------------------------- +! +! reset the output array +! + v(:,:,:,:) = 0.0d+00 + +! iterate over directions and update divergence with directional derivatives +! + do dir = 1, NDIMS + +! calculate contribution from the Y derivative of the Y component +! + call derivative_1st(dir, dh(dir), u(:,:,:), v(dir,:,:,:)) + + end do ! directions + +!------------------------------------------------------------------------------- +! + end subroutine gradient +! +!=============================================================================== +! ! subroutine CURL: ! --------------- ! From f9d3ed87dfb16ec147e48023110994d568656b0f Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 19 Sep 2014 07:34:04 -0300 Subject: [PATCH 11/12] OPERATORS: Add support for module profiling. Signed-off-by: Grzegorz Kowal --- src/makefile | 2 +- src/operators.F90 | 124 ++++++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 121 insertions(+), 5 deletions(-) diff --git a/src/makefile b/src/makefile index f680933..ae6e6e8 100644 --- a/src/makefile +++ b/src/makefile @@ -148,7 +148,7 @@ mesh.o : mesh.F90 blocks.o coordinates.o domains.o equations.o \ error.o interpolations.o mpitools.o problems.o refinement.o \ timers.o mpitools.o : mpitools.F90 timers.o -operators.o : operators.F90 +operators.o : operators.F90 timers.o parameters.o : parameters.F90 mpitools.o problems.o : problems.F90 blocks.o constants.o coordinates.o equations.o \ error.o parameters.o random.o timers.o diff --git a/src/operators.F90 b/src/operators.F90 index 98920a2..917e0a0 100644 --- a/src/operators.F90 +++ b/src/operators.F90 @@ -30,10 +30,22 @@ ! module operators +#ifdef PROFILE +! import external subroutines +! + use timers, only : set_timer, start_timer, stop_timer +#endif /* PROFILE */ + ! module variables are not implicit by default ! implicit none +#ifdef PROFILE +! timer indices +! + integer , save :: imi, imd, img, imc, iml, im1, im2 +#endif /* PROFILE */ + ! by default everything is public ! private @@ -41,7 +53,7 @@ module operators ! declare public subroutines ! public :: initialize_operators, finalize_operators - public :: divergence, curl, laplace + public :: divergence, gradient, curl, laplace !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! @@ -80,6 +92,27 @@ module operators ! !------------------------------------------------------------------------------- ! +#ifdef PROFILE +! set timer descriptions +! + call set_timer('operators:: initialization', imi) + call set_timer('operators:: divergence' , imd) + call set_timer('operators:: gradient' , img) + call set_timer('operators:: curl' , imc) + call set_timer('operators:: laplace' , iml) + call set_timer('operators:: 1st derivative', im1) + call set_timer('operators:: 2nd derivative', im2) + +! start accounting time for the module initialization/finalization +! + call start_timer(imi) +#endif /* PROFILE */ + +#ifdef PROFILE +! stop accounting time for the module initialization/finalization +! + call stop_timer(imi) +#endif /* PROFILE */ !------------------------------------------------------------------------------- ! @@ -110,6 +143,17 @@ module operators ! !------------------------------------------------------------------------------- ! +#ifdef PROFILE +! start accounting time for the module initialization/finalization +! + call start_timer(imi) +#endif /* PROFILE */ + +#ifdef PROFILE +! stop accounting time for the module initialization/finalization +! + call stop_timer(imi) +#endif /* PROFILE */ !------------------------------------------------------------------------------- ! @@ -155,6 +199,12 @@ module operators ! !------------------------------------------------------------------------------- ! +#ifdef PROFILE +! start accounting time for the divergence operator calculation +! + call start_timer(imd) +#endif /* PROFILE */ + ! allocate temporary array ! allocate(w(size(u,2), size(u,3), size(u,4))) @@ -167,7 +217,7 @@ module operators ! do dir = 1, NDIMS -! calculate contribution from the Y derivative of the Y component +! calculate derivative along the current direction ! call derivative_1st(dir, dh(dir), u(dir,:,:,:), w(:,:,:)) @@ -181,6 +231,12 @@ module operators ! deallocate(w) +#ifdef PROFILE +! stop accounting time for the divergence operator calculation +! + call stop_timer(imd) +#endif /* PROFILE */ + !------------------------------------------------------------------------------- ! end subroutine divergence @@ -220,20 +276,32 @@ module operators ! !------------------------------------------------------------------------------- ! +#ifdef PROFILE +! start accounting time for the gradient operator calculation +! + call start_timer(img) +#endif /* PROFILE */ + ! reset the output array ! v(:,:,:,:) = 0.0d+00 -! iterate over directions and update divergence with directional derivatives +! iterate over directions and calculate gradient components ! do dir = 1, NDIMS -! calculate contribution from the Y derivative of the Y component +! calculate derivative along the current direction ! call derivative_1st(dir, dh(dir), u(:,:,:), v(dir,:,:,:)) end do ! directions +#ifdef PROFILE +! stop accounting time for the gradient operator calculation +! + call stop_timer(img) +#endif /* PROFILE */ + !------------------------------------------------------------------------------- ! end subroutine gradient @@ -274,6 +342,12 @@ module operators ! !------------------------------------------------------------------------------- ! +#ifdef PROFILE +! start accounting time for the rotation operator calculation +! + call start_timer(imc) +#endif /* PROFILE */ + ! allocate temporary array ! allocate(w(size(u,2), size(u,3), size(u,4))) @@ -348,6 +422,12 @@ module operators ! deallocate(w) +#ifdef PROFILE +! stop accounting time for the rotation operator calculation +! + call stop_timer(imc) +#endif /* PROFILE */ + !------------------------------------------------------------------------------- ! end subroutine curl @@ -387,6 +467,12 @@ module operators ! !------------------------------------------------------------------------------- ! +#ifdef PROFILE +! start accounting time for the laplace operator calculation +! + call start_timer(iml) +#endif /* PROFILE */ + ! allocate temporary array ! allocate(w(size(u,1), size(u,2), size(u,3))) @@ -421,6 +507,12 @@ module operators ! deallocate(w) +#ifdef PROFILE +! stop accounting time for the laplace operator calculation +! + call stop_timer(iml) +#endif /* PROFILE */ + !------------------------------------------------------------------------------- ! end subroutine laplace @@ -473,6 +565,12 @@ module operators if (dir < 1 .or. dir > NDIMS .or. dh == 0.0d+00) return #endif /* DEBUG */ +#ifdef PROFILE +! start accounting time for the 1st derivative calculation +! + call start_timer(im1) +#endif /* PROFILE */ + ! prepare index limits ! m0 = size(u, dir) @@ -511,6 +609,12 @@ module operators end select +#ifdef PROFILE +! stop accounting time for the 1st derivative calculation +! + call stop_timer(im1) +#endif /* PROFILE */ + !------------------------------------------------------------------------------- ! end subroutine derivative_1st @@ -557,6 +661,12 @@ module operators if (dir < 1 .or. dir > NDIMS .or. dh == 0.0d+00) return #endif /* DEBUG */ +#ifdef PROFILE +! start accounting time for the 2nd derivative calculation +! + call start_timer(im2) +#endif /* PROFILE */ + ! prepare index limits ! m0 = size(u, dir) @@ -595,6 +705,12 @@ module operators end select +#ifdef PROFILE +! stop accounting time for the 2nd derivative calculation +! + call stop_timer(im2) +#endif /* PROFILE */ + !------------------------------------------------------------------------------- ! end subroutine derivative_2nd From 310961ca7eb93075b2629a9f63a53e53ceffca0d Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 19 Sep 2014 08:04:48 -0300 Subject: [PATCH 12/12] SOURCES: Implement GLM-MHD source terms (EGLM and HEGLM). Signed-off-by: Grzegorz Kowal --- src/sources.F90 | 129 +++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 112 insertions(+), 17 deletions(-) diff --git a/src/sources.F90 b/src/sources.F90 index 753888f..bcd9b7e 100644 --- a/src/sources.F90 +++ b/src/sources.F90 @@ -45,6 +45,10 @@ module sources integer, save :: imi, imu #endif /* PROFILE */ +! GLM-MHD source terms type (1 - EGLM, 2 - HEGLM) +! + integer , save :: glm_type = 0 + ! gravitational acceleration coefficient ! real(kind=8), save :: gpoint = 0.0d+00 @@ -95,7 +99,7 @@ module sources ! include external procedures and variables ! - use parameters , only : get_parameter_real + use parameters , only : get_parameter_string, get_parameter_real ! local variables are not implicit by default ! @@ -105,6 +109,10 @@ module sources ! logical, intent(in) :: verbose integer, intent(inout) :: iret + +! local variables +! + character(len=8) :: tglm = "none" ! !------------------------------------------------------------------------------- ! @@ -119,6 +127,23 @@ module sources call start_timer(imi) #endif /* PROFILE */ +! get the type of the GLM source terms +! + call get_parameter_string("glm_source_terms", tglm) + +! set the glm_type variable to correct value +! + select case(trim(tglm)) + case("eglm", "EGLM") + glm_type = 1 + + case("heglm", "HEGLM") + glm_type = 2 + + case default + glm_type = 0 + end select + ! get acceleration coefficient ! call get_parameter_real("gpoint" , gpoint ) @@ -135,6 +160,7 @@ module sources ! if (verbose) then + write (*,"(4x,a,1x,a) ") "glm source terms =", trim(tglm) write (*,"(4x,a,1x,1e9.2)") "point mass constant =", gpoint write (*,"(4x,a,1x,1e9.2)") "viscosity =", viscosity write (*,"(4x,a,1x,1e9.2)") "resistivity =", resistivity @@ -215,8 +241,8 @@ module sources use coordinates , only : ax, ay, az, adx, ady, adz, adxi, adyi, adzi use equations , only : nv, inx, iny, inz use equations , only : idn, ivx, ivy, ivz, imx, imy, imz, ien - use equations , only : ibx, iby, ibz - use operators , only : laplace, curl + use equations , only : ibx, iby, ibz, ibp + use operators , only : divergence, gradient, laplace, curl ! local variables are not implicit by default ! @@ -241,6 +267,7 @@ module sources real(kind=8), dimension(im) :: x real(kind=8), dimension(jm) :: y real(kind=8), dimension(km) :: z + real(kind=8), dimension(im,jm,km) :: db real(kind=8), dimension(3,im,jm,km) :: jc ! !------------------------------------------------------------------------------- @@ -440,9 +467,9 @@ module sources end if ! viscosity is not zero -! proceed only if the resistivity coefficient is not zero +!=== add magnetic field related source terms === ! - if (resistivity > 0.0d+00 .and. ibx > 0) then + if (ibx > 0) then ! prepare coordinate increments ! @@ -450,43 +477,111 @@ module sources dh(2) = ady(pdata%meta%level) dh(3) = adz(pdata%meta%level) +! add the EGLM-MHD source terms +! + if (glm_type > 0) then + +! calculate the magnetic field divergence +! + call divergence(dh(:), pdata%q(ibx:ibz,:,:,:), db(:,:,:)) + +! update the momentum component increments, i.e. +! d/dt (ρv) + ∇.F = - (∇.B)B +! + du(imx,:,:,:) = du(imx,:,:,:) - db(:,:,:) * pdata%q(ibx,:,:,:) + du(imy,:,:,:) = du(imy,:,:,:) - db(:,:,:) * pdata%q(iby,:,:,:) + du(imz,:,:,:) = du(imz,:,:,:) - db(:,:,:) * pdata%q(ibz,:,:,:) + +! update the energy equation +! + if (ien > 0 .and. ibp > 0) then + +! calculate the gradient of divergence potential +! + call gradient(dh(:), pdata%q(ibp,:,:,:), jc(inx:inz,:,:,:)) + +! add the divergence potential source term to the energy equation, i.e. +! d/dt E + ∇.F = - B.(∇ψ) +! + du(ien,:,:,:) = du(ien,:,:,:) & + - sum(pdata%q(ibx:ibz,:,:,:) * jc(inx:inz,:,:,:), 1) + + end if ! ien > 0 + +! add the HEGLM-MHD source terms +! + if (glm_type > 1) then + +! update magnetic field component increments, i.e. +! d/dt B + ∇.F = - (∇.B)v +! + du(ibx,:,:,:) = du(ibx,:,:,:) - db(:,:,:) * pdata%q(ivx,:,:,:) + du(iby,:,:,:) = du(iby,:,:,:) - db(:,:,:) * pdata%q(ivy,:,:,:) + du(ibz,:,:,:) = du(ibz,:,:,:) - db(:,:,:) * pdata%q(ivz,:,:,:) + +! update the energy equation +! + if (ien > 0) then + +! calculate scalar product of velocity and magnetic field +! + jc(inx,:,:,:) = sum(pdata%q(ivx:ivz,:,:,:) & + * pdata%q(ibx:ibz,:,:,:), 1) + +! add the divergence potential source term to the energy equation, i.e. +! d/dt E + ∇.F = - (∇.B) (v.B) +! + du(ien,:,:,:) = du(ien,:,:,:) - db(:,:,:) * jc(inx,:,:,:) + + end if ! ien > 0 + + end if ! glm_type > 1 + + end if ! glmtype > 0 + +! proceed only if the resistivity coefficient is not zero +! + if (resistivity > 0.0d+00) then + ! calculate the Laplace operator of B, i.e. Δ(B) ! - call laplace(dh(:), pdata%q(ibx,:,:,:), jc(inx,:,:,:)) - call laplace(dh(:), pdata%q(iby,:,:,:), jc(iny,:,:,:)) - call laplace(dh(:), pdata%q(ibz,:,:,:), jc(inz,:,:,:)) + call laplace(dh(:), pdata%q(ibx,:,:,:), jc(inx,:,:,:)) + call laplace(dh(:), pdata%q(iby,:,:,:), jc(iny,:,:,:)) + call laplace(dh(:), pdata%q(ibz,:,:,:), jc(inz,:,:,:)) ! update magnetic field component increments ! - du(ibx,:,:,:) = du(ibx,:,:,:) + resistivity * jc(inx,:,:,:) - du(iby,:,:,:) = du(iby,:,:,:) + resistivity * jc(iny,:,:,:) - du(ibz,:,:,:) = du(ibz,:,:,:) + resistivity * jc(inz,:,:,:) + du(ibx,:,:,:) = du(ibx,:,:,:) + resistivity * jc(inx,:,:,:) + du(iby,:,:,:) = du(iby,:,:,:) + resistivity * jc(iny,:,:,:) + du(ibz,:,:,:) = du(ibz,:,:,:) + resistivity * jc(inz,:,:,:) ! update energy equation ! - if (ien > 0) then + if (ien > 0) then ! add the first resistive source term to the energy equation, i.e. ! d/dt E + ∇.F = η B.[Δ(B)] ! - du(ien,:,:,:) = du(ien,:,:,:) & + du(ien,:,:,:) = du(ien,:,:,:) & + resistivity * (pdata%q(ibx,:,:,:) * jc(inx,:,:,:) & + pdata%q(iby,:,:,:) * jc(iny,:,:,:) & + pdata%q(ibz,:,:,:) * jc(inz,:,:,:)) ! calculate current density J = ∇xB ! - call curl(dh(:), pdata%q(ibx:ibz,:,:,:), jc(inx:inz,:,:,:)) + call curl(dh(:), pdata%q(ibx:ibz,:,:,:), jc(inx:inz,:,:,:)) ! add the second resistive source term to the energy equation, i.e. ! d/dt E + ∇.F = η J² ! - du(ien,:,:,:) = du(ien,:,:,:) & + du(ien,:,:,:) = du(ien,:,:,:) & + resistivity * sum(jc(:,:,:,:) * jc(:,:,:,:), 1) - end if ! energy equation present + end if ! energy equation present - end if ! resistivity is not zero + end if ! resistivity is not zero + + end if ! ibx > 0 #ifdef PROFILE ! stop accounting time for source terms