From 61ed99382170f44d5046eadb65880dcb7201c03c Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 20 Sep 2021 08:44:38 -0300 Subject: [PATCH] =?UTF-8?q?EVOLUTION:=20Fix=20decay=20rate=20of=20=CF=88?= =?UTF-8?q?=20field.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This scalar potential tracks the development of div(B). Its decay depends on the cell size, so make it calculated per refinement level. Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 92 ++++++++++++++++++++++++++++--------------- 1 file changed, 60 insertions(+), 32 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 329ba12..3619ab5 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -69,7 +69,6 @@ module evolution ! coefficients controlling the decay of scalar potential ѱ ! real(kind=8) , save :: glm_alpha = 5.0d-01 - real(kind=8) , save :: decay = 7.788007830714049d-01 ! time variables ! @@ -98,6 +97,10 @@ module evolution ! real(kind=8), dimension(3), save :: errs +! GLM variables +! + real(kind=8), dimension(:), allocatable, save :: adecay, aglm + ! by default everything is private ! private @@ -140,7 +143,8 @@ module evolution ! include external procedures ! - use parameters, only : get_parameter + use coordinates, only : toplev, adx, ady, adz + use parameters , only : get_parameter ! local variables are not implicit by default ! @@ -281,20 +285,19 @@ module evolution end select -! proceed if everything is fine -! - if (status == 0) then - -! calculate the decay factor for magnetic field divergence scalar source term -! - decay = exp(- glm_alpha * cfl) - - end if ! status - ! reset recent error history ! errs = 1.0d+00 +! GLM coefficients for all refinement levels +! + allocate(adecay(toplev), aglm(toplev)) +#if NDIMS == 3 + aglm(:) = - glm_alpha / min(adx(:), ady(:), adz(:)) +#else /* NDIMS == 3 */ + aglm(:) = - glm_alpha / min(adx(:), ady(:)) +#endif /* NDIMS == 3 */ + #ifdef PROFILE ! stop accounting time for module initialization/finalization ! @@ -345,6 +348,9 @@ module evolution call start_timer(imi) #endif /* PROFILE */ + if (allocated(adecay)) deallocate(adecay) + if (allocated(aglm)) deallocate(aglm) + ! reset the status flag ! status = 0 @@ -682,7 +688,7 @@ module evolution use blocks , only : block_data, list_data use boundaries, only : boundary_fluxes - use equations , only : ibp + use equations , only : ibp, cmax use sources , only : update_sources implicit none @@ -722,10 +728,12 @@ module evolution end do if (ibp > 0) then + adecay(:) = exp(aglm(:) * cmax * dt) + pdata => list_data do while (associated(pdata)) - pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) + pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:) pdata => pdata%next end do @@ -761,8 +769,8 @@ module evolution use blocks , only : block_data, list_data use boundaries, only : boundary_fluxes - use equations , only : ibp - use sources , only : update_sources + use equations , only : ibp, cmax + use sources , only : update_sources implicit none @@ -833,10 +841,12 @@ module evolution end do if (ibp > 0) then + adecay(:) = exp(aglm(:) * cmax * dt) + pdata => list_data do while (associated(pdata)) - pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) + pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:) pdata => pdata%next end do @@ -876,7 +886,7 @@ module evolution use blocks , only : block_data, list_data, get_dblocks use boundaries , only : boundary_fluxes use coordinates, only : nc => ncells, nb, ne - use equations , only : errors, nf, ibp + use equations , only : errors, nf, ibp, cmax #ifdef MPI use mpitools , only : reduce_maximum #endif /* MPI */ @@ -1100,13 +1110,21 @@ module evolution pdata%u => pdata%uu(:,:,:,:,1) -! update ψ with its source term -! - if (ibp > 0) pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) - pdata => pdata%next end do + if (ibp > 0) then + adecay(:) = exp(aglm(:) * cmax * dt) + + pdata => list_data + do while (associated(pdata)) + + pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:) + + pdata => pdata%next + end do + end if + call update_variables(tm, dtm) #ifdef PROFILE @@ -1138,7 +1156,7 @@ module evolution use blocks , only : block_data, list_data use boundaries, only : boundary_fluxes - use equations , only : ibp + use equations , only : ibp, cmax use sources , only : update_sources implicit none @@ -1243,10 +1261,12 @@ module evolution end do if (ibp > 0) then + adecay(:) = exp(aglm(:) * cmax * dt) + pdata => list_data do while (associated(pdata)) - pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) + pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:) pdata => pdata%next end do @@ -1284,7 +1304,7 @@ module evolution use blocks , only : block_data, list_data use boundaries, only : boundary_fluxes - use equations , only : ibp + use equations , only : ibp, cmax use sources , only : update_sources implicit none @@ -1407,10 +1427,12 @@ module evolution end do if (ibp > 0) then + adecay(:) = exp(aglm(:) * cmax * dt) + pdata => list_data do while (associated(pdata)) - pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) + pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:) pdata => pdata%next end do @@ -1448,7 +1470,7 @@ module evolution use blocks , only : block_data, list_data use boundaries, only : boundary_fluxes - use equations , only : ibp + use equations , only : ibp, cmax use sources , only : update_sources implicit none @@ -1615,10 +1637,12 @@ module evolution end do if (ibp > 0) then + adecay(:) = exp(aglm(:) * cmax * dt) + pdata => list_data do while (associated(pdata)) - pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) + pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:) pdata => pdata%next end do @@ -1656,7 +1680,7 @@ module evolution use blocks , only : block_data, list_data, get_dblocks use boundaries , only : boundary_fluxes use coordinates, only : nc => ncells, nb, ne - use equations , only : errors, nf, ibp + use equations , only : errors, nf, ibp, cmax #ifdef MPI use mpitools , only : reduce_maximum #endif /* MPI */ @@ -1983,10 +2007,12 @@ module evolution call update_variables(tm, dtm) if (ibp > 0) then + adecay(:) = exp(aglm(:) * cmax * dt) + pdata => list_data do while (associated(pdata)) - pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) + pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:) pdata => pdata%next end do @@ -2021,7 +2047,7 @@ module evolution use blocks , only : block_data, list_data use boundaries, only : boundary_fluxes - use equations , only : ibp + use equations , only : ibp, cmax use sources , only : update_sources implicit none @@ -2157,10 +2183,12 @@ module evolution end do if (ibp > 0) then + adecay(:) = exp(aglm(:) * cmax * dt) + pdata => list_data do while (associated(pdata)) - pdata%u(ibp,:,:,:) = decay * pdata%u(ibp,:,:,:) + pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:) pdata => pdata%next end do