From 454c3f53121888d40f58146e3d885954a3857bcf Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 19 Dec 2013 21:18:21 -0200 Subject: [PATCH 01/10] PROBLEMS: Implement reconnection problem. Signed-off-by: Grzegorz Kowal --- src/problems.F90 | 156 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 156 insertions(+) diff --git a/src/problems.F90 b/src/problems.F90 index 766d768..c811ede 100644 --- a/src/problems.F90 +++ b/src/problems.F90 @@ -126,6 +126,11 @@ module problems case("blast") call setup_problem_blast(pdata) +! specific problems +! + case("reconnection") + call setup_problem_reconnection(pdata) + case default call print_error("problems::init_problem()" & , "Setup subroutime is not implemented for this problem!") @@ -430,6 +435,157 @@ module problems !------------------------------------------------------------------------------- ! end subroutine setup_problem_blast +! +!=============================================================================== +! +! subroutine SETUP_PROBLEM_RECONNECTION: +! ------------------------------------- +! +! Subroutine sets the initial conditions for the reconnection problem. +! +! Arguments: +! +! pdata - pointer to the datablock structure of the currently initialized +! block; +! +!=============================================================================== +! + subroutine setup_problem_reconnection(pdata) + +! include external procedures and variables +! + use blocks , only : block_data + use coordinates, only : im, jm, km + use coordinates, only : ay + use equations , only : prim2cons + use equations , only : nv + use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp + use parameters , only : get_parameter_real + use random , only : randomn + +! local variables are not implicit by default +! + implicit none + +! input arguments +! + type(block_data), pointer, intent(inout) :: pdata + +! default parameter values +! + real(kind=8), save :: dens = 1.00d+00 + real(kind=8), save :: pres = 1.00d+00 + real(kind=8), save :: bamp = 1.00d+00 + real(kind=8), save :: bgui = 0.00d+00 + real(kind=8), save :: vper = 1.00d-02 + real(kind=8), save :: ycut = 1.00d-01 + +! local saved parameters +! + logical , save :: first = .true. + +! local variables +! + integer :: i, j, k + +! local arrays +! + real(kind=8), dimension(nv,jm) :: q, u + real(kind=8), dimension(jm) :: y +! +!------------------------------------------------------------------------------- +! +! prepare problem constants during the first subroutine call +! + if (first) then + +! get problem parameters +! + call get_parameter_real("dens" , dens) + call get_parameter_real("pres" , pres) + call get_parameter_real("brecon", bamp) + call get_parameter_real("bguide", bgui) + call get_parameter_real("vper" , vper) + call get_parameter_real("ycut" , ycut) + +! reset the first execution flag +! + first = .false. + + end if ! first call + +! prepare block coordinates +! + y(1:jm) = pdata%meta%ymin + ay(pdata%meta%level,1:jm) + +! set the density and pressure +! + q(idn,:) = dens + if (ipr > 0) q(ipr,:) = pres + +! reset velocity components +! + q(ivx,:) = 0.0d+00 + q(ivy,:) = 0.0d+00 + q(ivz,:) = 0.0d+00 + +! if magnetic field is present, set it to be uniform with the desired strength +! and orientation +! + if (ibx > 0) then + + q(ibx,:) = 0.0d+00 + q(iby,:) = 0.0d+00 + q(ibz,:) = bgui + q(ibp,:) = 0.0d+00 + + end if + +! iterate over all positions in the XZ plane +! + do k = 1, km + do i = 1, im + +! sweep along the Y coordinate +! + do j = 1, jm + +! set the random velocity field in a layer near current sheet +! + if (abs(y(j)) <= ycut) then + + q(ivx,j) = vper * randomn() + q(ivy,j) = vper * randomn() +#if NDIMS == 3 + q(ivz,j) = vper * randomn() +#endif /* NDIMS == 3 */ + + end if + +! set antiparallel magnetic field component +! + if (ibx > 0) q(ibx,j) = sign(bamp, y(j)) + + end do ! j = 1, jm + +! convert the primitive variables to conservative ones +! + call prim2cons(jm, q(1:nv,1:jm), u(1:nv,1:jm)) + +! copy the conserved variables to the current block +! + pdata%u(1:nv,i,1:jm,k) = u(1:nv,1:jm) + +! copy the primitive variables to the current block +! + pdata%q(1:nv,i,1:jm,k) = q(1:nv,1:jm) + + end do ! i = 1, im + end do ! k = 1, km + +!------------------------------------------------------------------------------- +! + end subroutine setup_problem_reconnection !=============================================================================== ! From 1721758285296e8b901e3784757b6b653c61f3dd Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 8 Aug 2014 11:39:00 -0300 Subject: [PATCH 02/10] PROBLEMS: Control x-size of the perturbation region for reconnection. Signed-off-by: Grzegorz Kowal --- src/problems.F90 | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/problems.F90 b/src/problems.F90 index b74d205..95d4d4e 100644 --- a/src/problems.F90 +++ b/src/problems.F90 @@ -1269,7 +1269,7 @@ module problems ! use blocks , only : block_data use coordinates, only : im, jm, km - use coordinates, only : ay + use coordinates, only : ax, ay use equations , only : prim2cons use equations , only : nv use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp @@ -1291,6 +1291,7 @@ module problems real(kind=8), save :: bamp = 1.00d+00 real(kind=8), save :: bgui = 0.00d+00 real(kind=8), save :: vper = 1.00d-02 + real(kind=8), save :: xcut = 4.00d-01 real(kind=8), save :: ycut = 1.00d-01 ! local saved parameters @@ -1304,6 +1305,7 @@ module problems ! local arrays ! real(kind=8), dimension(nv,jm) :: q, u + real(kind=8), dimension(im) :: x real(kind=8), dimension(jm) :: y ! !------------------------------------------------------------------------------- @@ -1325,6 +1327,7 @@ module problems call get_parameter_real("brecon", bamp) call get_parameter_real("bguide", bgui) call get_parameter_real("vper" , vper) + call get_parameter_real("xcut" , xcut) call get_parameter_real("ycut" , ycut) ! reset the first execution flag @@ -1335,6 +1338,7 @@ module problems ! prepare block coordinates ! + x(1:im) = pdata%meta%xmin + ax(pdata%meta%level,1:im) y(1:jm) = pdata%meta%ymin + ay(pdata%meta%level,1:jm) ! set the density and pressure @@ -1373,7 +1377,7 @@ module problems ! set the random velocity field in a layer near current sheet ! - if (abs(y(j)) <= ycut) then + if (abs(x(i)) <= xcut .and. abs(y(j)) <= ycut) then q(ivx,j) = vper * randomn() q(ivy,j) = vper * randomn() From 68ff67c2439f6fc13f6483774ec0c73037d0d6bc Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 8 Aug 2014 11:39:43 -0300 Subject: [PATCH 03/10] PROBLEMS: Change the y-size of perturbation region in reconnection. Signed-off-by: Grzegorz Kowal --- src/problems.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/problems.F90 b/src/problems.F90 index 95d4d4e..165f276 100644 --- a/src/problems.F90 +++ b/src/problems.F90 @@ -1292,7 +1292,7 @@ module problems real(kind=8), save :: bgui = 0.00d+00 real(kind=8), save :: vper = 1.00d-02 real(kind=8), save :: xcut = 4.00d-01 - real(kind=8), save :: ycut = 1.00d-01 + real(kind=8), save :: ycut = 1.00d-02 ! local saved parameters ! From d2b84c8accd94c548e158e095e1a5145c2b81f71 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 8 Aug 2014 11:58:29 -0300 Subject: [PATCH 04/10] PROBLEMS: Rewrite the reconnection problem. Reset all variables in the loops, avoiding to have them set in the region which shouldn't be touched. Signed-off-by: Grzegorz Kowal --- src/problems.F90 | 75 ++++++++++++++++++++++++------------------------ 1 file changed, 37 insertions(+), 38 deletions(-) diff --git a/src/problems.F90 b/src/problems.F90 index 165f276..78a7c9f 100644 --- a/src/problems.F90 +++ b/src/problems.F90 @@ -1341,57 +1341,56 @@ module problems x(1:im) = pdata%meta%xmin + ax(pdata%meta%level,1:im) y(1:jm) = pdata%meta%ymin + ay(pdata%meta%level,1:jm) -! set the density and pressure -! - q(idn,:) = dens - if (ipr > 0) q(ipr,:) = pres - -! reset velocity components -! - q(ivx,:) = 0.0d+00 - q(ivy,:) = 0.0d+00 - q(ivz,:) = 0.0d+00 - -! if magnetic field is present, set it to be uniform with the desired strength -! and orientation -! - if (ibx > 0) then - -! set magnetic field components -! - q(ibx,:) = 0.0d+00 - q(iby,:) = 0.0d+00 - q(ibz,:) = bgui - q(ibp,:) = 0.0d+00 - - end if - ! iterate over all positions in the XZ plane ! do k = 1, km do i = 1, im -! sweep along the Y coordinate +! set the uniform density and pressure ! - do j = 1, jm + q(idn,:) = dens + if (ipr > 0) q(ipr,:) = pres -! set the random velocity field in a layer near current sheet +! if magnetic field is present, set its initial configuration ! - if (abs(x(i)) <= xcut .and. abs(y(j)) <= ycut) then + if (ibx > 0) then - q(ivx,j) = vper * randomn() - q(ivy,j) = vper * randomn() -#if NDIMS == 3 - q(ivz,j) = vper * randomn() -#endif /* NDIMS == 3 */ - - end if +! reset magnetic field components +! + q(ibx,:) = 0.0d+00 + q(iby,:) = 0.0d+00 + q(ibz,:) = bgui + q(ibp,:) = 0.0d+00 ! set antiparallel magnetic field component ! - if (ibx > 0) q(ibx,j) = sign(bamp, y(j)) + do j = 1, jm + q(ibx,j) = sign(bamp, y(j)) + end do ! j = 1, jm - end do ! j = 1, jm + end if ! ibx > 0 + +! reset velocity components +! + q(ivx,:) = 0.0d+00 + q(ivy,:) = 0.0d+00 + q(ivz,:) = 0.0d+00 + +! set the random velocity field in a layer near current sheet +! + if (abs(x(i)) <= xcut) then + do j = 1, jm + if (abs(y(j)) <= ycut) then + + q(ivx,j) = vper * randomn() + q(ivy,j) = vper * randomn() +#if NDIMS == 3 + q(ivz,j) = vper * randomn() +#endif /* NDIMS == 3 */ + + end if ! |y| < ycut + end do ! j = 1, jm + end if ! |x| < xcut ! convert the primitive variables to conservative ones ! From e5536c37c4da18455bc435e964562c4b04733511 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 20 Aug 2014 18:17:39 -0300 Subject: [PATCH 05/10] COORDINATES: Add decay factors, used in boundaries. Decay factors determines the factor by which the variable is decayed at the boundaries. Signed-off-by: Grzegorz Kowal --- src/coordinates.F90 | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/src/coordinates.F90 b/src/coordinates.F90 index 9281ec2..5fd6978 100644 --- a/src/coordinates.F90 +++ b/src/coordinates.F90 @@ -88,6 +88,10 @@ module coordinates real(kind=8), save :: vol = 1.0d+00 real(kind=8), save :: voli = 1.0d+00 +! the characteristic decay length +! + real(kind=8), save :: ldec = 1.0d-03 + ! the block coordinates for all levels of refinement ! real(kind=8), dimension(:,:), allocatable, save :: ax , ay , az @@ -95,6 +99,10 @@ module coordinates real(kind=8), dimension(: ), allocatable, save :: adxi, adyi, adzi real(kind=8), dimension(: ), allocatable, save :: advol +! the decay factors for all levels +! + real(kind=8), dimension(: ), allocatable, save :: adec + ! by default everything is private ! public @@ -261,6 +269,7 @@ module coordinates allocate(adyi (toplev)) allocate(adzi (toplev)) allocate(advol(toplev)) + allocate(adec (toplev)) ! reset all coordinate variables to initial values ! @@ -323,6 +332,16 @@ module coordinates end do ! l = 1, toplev +! get the characteristic decay scale +! + call get_parameter_real("ldecay", ldec) + +! calculate the decay factors +! + do l = 1, toplev + adec(l) = exp(- adx(l) / ldec) + end do ! l = 1, toplev + ! print general information about the level resolutions ! if (verbose) then @@ -403,6 +422,7 @@ module coordinates if (allocated(adyi) ) deallocate(adyi) if (allocated(adzi) ) deallocate(adzi) if (allocated(advol)) deallocate(advol) + if (allocated(adec) ) deallocate(adec) !------------------------------------------------------------------------------- ! From b855dbb87fe369b39a018b3e85b937984842e414 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 26 Aug 2014 07:33:19 -0300 Subject: [PATCH 06/10] PROBLEMS: Rewrite reconnection setup. Add magnetic field perturbation calculated from the vector potential. This is used in the studies of laminar and fast reconnection problems. Signed-off-by: Grzegorz Kowal --- src/makefile | 2 +- src/problems.F90 | 121 +++++++++++++++++++++++++++++++++++++---------- 2 files changed, 97 insertions(+), 26 deletions(-) diff --git a/src/makefile b/src/makefile index f680933..d81ae3b 100644 --- a/src/makefile +++ b/src/makefile @@ -151,7 +151,7 @@ mpitools.o : mpitools.F90 timers.o operators.o : operators.F90 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 + error.o operators.o parameters.o random.o timers.o random.o : random.F90 mpitools.o parameters.o refinement.o : refinement.F90 blocks.o coordinates.o equations.o \ parameters.o timers.o diff --git a/src/problems.F90 b/src/problems.F90 index 78a7c9f..961beb9 100644 --- a/src/problems.F90 +++ b/src/problems.F90 @@ -1268,11 +1268,14 @@ module problems ! include external procedures and variables ! use blocks , only : block_data + use constants , only : pi2 use coordinates, only : im, jm, km - use coordinates, only : ax, ay + use coordinates, only : ax, ay, adx, ady, adz use equations , only : prim2cons use equations , only : nv use equations , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp + use equations , only : csnd2 + use operators , only : curl use parameters , only : get_parameter_real use random , only : randomn @@ -1289,10 +1292,14 @@ module problems real(kind=8), save :: dens = 1.00d+00 real(kind=8), save :: pres = 1.00d+00 real(kind=8), save :: bamp = 1.00d+00 + real(kind=8), save :: bper = 0.00d+00 real(kind=8), save :: bgui = 0.00d+00 real(kind=8), save :: vper = 1.00d-02 real(kind=8), save :: xcut = 4.00d-01 real(kind=8), save :: ycut = 1.00d-02 + real(kind=8), save :: yth = 1.00d-16 + real(kind=8), save :: pth = 1.00d-02 + real(kind=8), save :: pmag = 5.00d-01 ! local saved parameters ! @@ -1301,12 +1308,14 @@ module problems ! local variables ! integer :: i, j, k + real(kind=8) :: yt, yp ! local arrays ! real(kind=8), dimension(nv,jm) :: q, u real(kind=8), dimension(im) :: x - real(kind=8), dimension(jm) :: y + real(kind=8), dimension(jm) :: y, pm + real(kind=8), dimension(3) :: dh ! !------------------------------------------------------------------------------- ! @@ -1322,13 +1331,20 @@ module problems ! get problem parameters ! - call get_parameter_real("dens" , dens) - call get_parameter_real("pres" , pres) - call get_parameter_real("brecon", bamp) - call get_parameter_real("bguide", bgui) - call get_parameter_real("vper" , vper) - call get_parameter_real("xcut" , xcut) - call get_parameter_real("ycut" , ycut) + call get_parameter_real("dens", dens) + call get_parameter_real("pres", pres) + call get_parameter_real("bamp", bamp) + call get_parameter_real("bper", bper) + call get_parameter_real("bgui", bgui) + call get_parameter_real("vper", vper) + call get_parameter_real("xcut", xcut) + call get_parameter_real("ycut", ycut) + call get_parameter_real("yth" , yth ) + call get_parameter_real("pth" , pth ) + +! calculate the maximum magnetic pressure +! + pmag = 0.5d+00 * (bamp**2 + bgui**2) ! reset the first execution flag ! @@ -1341,40 +1357,95 @@ module problems x(1:im) = pdata%meta%xmin + ax(pdata%meta%level,1:im) y(1:jm) = pdata%meta%ymin + ay(pdata%meta%level,1:jm) +! prepare cell sizes +! + dh(1) = adx(pdata%meta%level) + dh(2) = ady(pdata%meta%level) + dh(3) = adz(pdata%meta%level) + +! calculate the perturbation of magnetic field +! + if (bper /= 0.0d+00) then + +! initiate the vector potential (we use velocity components to store vector +! potential temporarily, and we store magnetic field perturbation in U) +! + do k = 1, km + do j = 1, jm + do i = 1, im + yt = abs(y(j) / yth) + yt = log(exp(yt) + exp(- yt)) + yp = y(j) / pth + + pdata%q(ivx,i,j,k) = 0.0d+00 + pdata%q(ivy,i,j,k) = 0.0d+00 + pdata%q(ivz,i,j,k) = bper * cos(pi2 * x(i)) * exp(- yp * yp) / pi2 + end do ! i = 1, im + end do ! j = 1, jm + end do ! k = 1, km + +! calculate magnetic field components from vector potential +! + call curl(dh(1:3), pdata%q(ivx:ivz,1:im,1:jm,1:km) & + , pdata%q(ibx:ibz,1:im,1:jm,1:km)) + + else + +! reset magnetic field components +! + pdata%q(ibx:ibz,:,:,:) = 0.0d+00 + + end if ! bper /= 0.0 + ! iterate over all positions in the XZ plane ! do k = 1, km do i = 1, im -! set the uniform density and pressure -! - q(idn,:) = dens - if (ipr > 0) q(ipr,:) = pres - ! if magnetic field is present, set its initial configuration ! if (ibx > 0) then -! reset magnetic field components -! - q(ibx,:) = 0.0d+00 - q(iby,:) = 0.0d+00 - q(ibz,:) = bgui - q(ibp,:) = 0.0d+00 - ! set antiparallel magnetic field component ! do j = 1, jm - q(ibx,j) = sign(bamp, y(j)) + q(ibx,j) = bamp * tanh(y(j) / yth) end do ! j = 1, jm +! set tangential magnetic field components +! + q(iby,1:jm) = 0.0d+00 + q(ibz,1:jm) = bgui + q(ibp,1:jm) = 0.0d+00 + +! calculate local magnetic pressure +! + pm(1:jm) = 0.5d+00 * sum(q(ibx:ibz,:) * q(ibx:ibz,:), 1) + +! add magnetic field perturbation +! + if (bper /= 0.0d+00) then + q(ibx,1:jm) = q(ibx,1:jm) + pdata%q(ibx,i,1:jm,k) + q(iby,1:jm) = q(iby,1:jm) + pdata%q(iby,i,1:jm,k) + q(ibz,1:jm) = q(ibz,1:jm) + pdata%q(ibz,i,1:jm,k) + end if ! bper /= 0.0 + end if ! ibx > 0 +! set the uniform density and pressure +! + if (ipr > 0) then + q(idn,1:jm) = dens + q(ipr,1:jm) = pres + (pmag - pm(1:jm)) + else + q(idn,1:jm) = dens + (pmag - pm(1:jm)) / csnd2 + end if + ! reset velocity components ! - q(ivx,:) = 0.0d+00 - q(ivy,:) = 0.0d+00 - q(ivz,:) = 0.0d+00 + q(ivx,1:jm) = 0.0d+00 + q(ivy,1:jm) = 0.0d+00 + q(ivz,1:jm) = 0.0d+00 ! set the random velocity field in a layer near current sheet ! From 02352a271c9659e5826c5157da257a11a30aec92 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 29 Aug 2014 09:11:18 -0300 Subject: [PATCH 07/10] INTEGRALS: Add reconnection rate calculation. The reconnection rate is calculated by averagin the inflow velocity at the bottom and top boundaries. It is stored in reconnection_??.dat files. Signed-off-by: Grzegorz Kowal --- src/integrals.F90 | 42 ++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 40 insertions(+), 2 deletions(-) diff --git a/src/integrals.F90 b/src/integrals.F90 index 17d2f99..1ce6804 100644 --- a/src/integrals.F90 +++ b/src/integrals.F90 @@ -55,6 +55,7 @@ module integrals ! integer(kind=4), save :: funit = 7 integer(kind=4), save :: sunit = 8 + integer(kind=4), save :: runit = 9 integer(kind=4), save :: iintd = 1 ! by default everything is private @@ -193,6 +194,30 @@ module integrals , 'mean(Malf)', 'min(Malf)', 'max(Malf)' write(sunit,"('#')") +! generate the reconnection file name +! + write(fname, "('reconnection_',i2.2,'.dat')") irun + +! create a new statistics file +! +#ifdef GNU + open (newunit = runit, file = fname, form = 'formatted' & + , status = 'replace') +#endif /* GNU */ +#ifdef INTEL + open (newunit = runit, file = fname, form = 'formatted' & + , status = 'replace', buffered = 'yes') +#endif /* INTEL */ +#ifdef IBM + open (unit = runit, file = fname, form = 'formatted' & + , status = 'replace') +#endif /* IBM */ + +! write the integral file header +! + write(runit,"('#',a8,3(1x,a18))") 'step', 'time', 'Vrec_l', 'Vrec_u' + write(runit,"('#')") + end if ! master #ifdef PROFILE @@ -238,6 +263,7 @@ module integrals if (master) then close(funit) close(sunit) + close(runit) end if #ifdef PROFILE @@ -267,7 +293,8 @@ module integrals ! use blocks , only : block_meta, block_data, list_data use coordinates , only : in, jn, kn, ib, jb, kb, ie, je, ke - use coordinates , only : advol, voli + use coordinates , only : adx, adz, advol, voli + use coordinates , only : ymin, ymax, xlen, zlen, yarea use equations , only : idn, ipr, ivx, ivy, ivz, ibx, iby, ibz, ibp use equations , only : ien, imx, imy, imz use equations , only : gamma, csnd @@ -286,7 +313,7 @@ module integrals ! local variables ! integer :: iret - real(kind=8) :: dvol, dvolh + real(kind=8) :: dvol, dvolh, dxz ! local pointers ! @@ -352,6 +379,7 @@ module integrals ! dvol = advol(pdata%meta%level) dvolh = 0.5d+00 * dvol + dxz = adx(pdata%meta%level) * adz(pdata%meta%level) / yarea ! sum up density and momenta components ! @@ -444,6 +472,15 @@ module integrals mxarr(7) = max(mxarr(7), maxval(tmp(:,:,:))) end if +! get the inflow speed +! + if (pdata%meta%ymin == ymin) then + inarr(10) = inarr(10) + sum(pdata%q(ivy,ib:ie,jb,kb:ke)) * dxz + end if + if (pdata%meta%ymax == ymax) then + inarr(11) = inarr(11) + sum(pdata%q(ivy,ib:ie,je,kb:ke)) * dxz + end if + ! associate the pointer with the next block on the list ! pdata => pdata%next @@ -482,6 +519,7 @@ module integrals , avarr(5), mnarr(5), mxarr(5) & , avarr(6), mnarr(6), mxarr(6) & , avarr(7), mnarr(7), mxarr(7) + write(runit,"(i9, 3(1x,1e18.8))") step, time, inarr(10:11) end if #ifdef PROFILE From b2d256aef8d845cd0125fb3625ec32b98996aa45 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 27 Sep 2014 14:13:38 -0300 Subject: [PATCH 08/10] INTEGRALS: Extend the exponent to 3 digits in reconnection rate. Signed-off-by: Grzegorz Kowal --- src/integrals.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/integrals.F90 b/src/integrals.F90 index 61be9a5..75aecf2 100644 --- a/src/integrals.F90 +++ b/src/integrals.F90 @@ -519,7 +519,7 @@ module integrals , avarr(5), mnarr(5), mxarr(5) & , avarr(6), mnarr(6), mxarr(6) & , avarr(7), mnarr(7), mxarr(7) - write(runit,"(i9, 3(1x,1e18.8))") step, time, inarr(10:11) + write(runit,"(i9, 3(1x,1e18.8e3))") step, time, inarr(10:11) end if #ifdef PROFILE From eee6117f6594419b8f4ef1269286ae3cb1eed670 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 13 Oct 2014 19:24:24 -0300 Subject: [PATCH 09/10] BOUNDARIES: Implement vertical reconnection boundary conditions. These conditions keep decays the antiparallel and guide magnetic components to their initial values. The same is done for density and pressure. Also the vertical (Y) magnetic field component is calculated from the div(B)=0 condition. Signed-off-by: Grzegorz Kowal --- src/boundaries.F90 | 201 +++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 194 insertions(+), 7 deletions(-) diff --git a/src/boundaries.F90 b/src/boundaries.F90 index be2f7d2..3b9d2d7 100644 --- a/src/boundaries.F90 +++ b/src/boundaries.F90 @@ -47,10 +47,11 @@ module boundaries ! parameters corresponding to the boundary type ! - integer, parameter :: bnd_periodic = 0 - integer, parameter :: bnd_open = 1 - integer, parameter :: bnd_outflow = 2 - integer, parameter :: bnd_reflective = 3 + integer, parameter :: bnd_periodic = 0 + integer, parameter :: bnd_open = 1 + integer, parameter :: bnd_outflow = 2 + integer, parameter :: bnd_reflective = 3 + integer, parameter :: bnd_reconnection = 4 ! variable to store boundary type flags ! @@ -175,6 +176,8 @@ module boundaries bnd_type(2,1) = bnd_outflow case("reflective", "reflecting", "reflect") bnd_type(2,1) = bnd_reflective + case("reconnection", "recon", "rec") + bnd_type(2,1) = bnd_reconnection case default bnd_type(2,1) = bnd_periodic end select @@ -186,6 +189,8 @@ module boundaries bnd_type(2,2) = bnd_outflow case("reflective", "reflecting", "reflect") bnd_type(2,2) = bnd_reflective + case("reconnection", "recon", "rec") + bnd_type(2,2) = bnd_reconnection case default bnd_type(2,2) = bnd_periodic end select @@ -1141,6 +1146,7 @@ module boundaries ! if (.not. associated(pmeta%edges(i,j,m)%ptr)) & call block_boundary_specific(i, j, k, n & + , pmeta%level & , pmeta%data%q(1:nv,1:im,1:jm,1:km)) end do ! i = 1, sides @@ -1169,6 +1175,7 @@ module boundaries ! if (.not. associated(pmeta%faces(i,j,k,n)%ptr)) & call block_boundary_specific(i, j, k, n & + , pmeta%level & , pmeta%data%q(1:nv,1:im,1:jm,1:km)) end do ! i = 1, sides @@ -6022,20 +6029,23 @@ module boundaries ! ! nc - the edge direction; ! ic, jc, kc - the corner position; +! lv - the block level; ! qn - the variable array; ! !=============================================================================== ! - subroutine block_boundary_specific(ic, jc, kc, nc, qn) + subroutine block_boundary_specific(ic, jc, kc, nc, lv, qn) ! import external procedures and variables ! use coordinates , only : im , jm , km , ng use coordinates , only : ib , jb , kb , ie , je , ke use coordinates , only : ibl, jbl, kbl, ieu, jeu, keu + use coordinates , only : ady, adxi, adzi use equations , only : nv - use equations , only : idn, ivx, ivy, ivz, ibx, iby, ibz, ibp + use equations , only : idn, ipr, ivx, ivy, ivz, ibx, iby, ibz, ibp use error , only : print_error, print_warning + use parameters , only : get_parameter_real ! local variables are not implicit by default ! @@ -6044,9 +6054,21 @@ module boundaries ! subroutine arguments ! integer , intent(in) :: ic, jc, kc - integer , intent(in) :: nc + integer , intent(in) :: nc, lv real(kind=8), dimension(1:nv,1:im,1:jm,1:km), intent(inout) :: qn +! default parameter values +! + real(kind=8), save :: dens = 1.00d+00 + real(kind=8), save :: pres = 1.00d+00 + real(kind=8), save :: bamp = 1.00d+00 + real(kind=8), save :: bgui = 0.00d+00 + real(kind=8), save :: blim = 1.00d+00 + +! local saved parameters +! + logical , save :: first = .true. + ! local variables ! integer :: i , j , k @@ -6054,9 +6076,38 @@ module boundaries integer :: iu, ju, ku integer :: is, js, ks integer :: it, jt, kt + integer :: im1, ip1 + integer :: jm2, jm1, jp1, jp2 +#if NDIMS == 3 + integer :: km1, kp1 +#endif /* NDIMS == 3 */ + real(kind=8) :: dyx, dyz + real(kind=8) :: fl, fr ! !------------------------------------------------------------------------------- ! +! prepare problem constants during the first subroutine call +! + if (first) then + +! get problem parameters +! + call get_parameter_real("dens" , dens) + call get_parameter_real("pres" , pres) + call get_parameter_real("bamp" , bamp) + call get_parameter_real("bgui" , bgui) + call get_parameter_real("blimit", blim) + +! upper limit for blim +! + blim = max(blim, ng * ady(1)) + +! reset the first execution flag +! + first = .false. + + end if ! first call + ! apply specific boundaries depending on the direction ! select case(nc) @@ -6212,6 +6263,142 @@ module boundaries end do ! j = jeu, jm end if +! "reconnection" boundary conditions +! + case(bnd_reconnection) + +! process case with magnetic field, otherwise revert to standard outflow +! + if (ibx > 0) then + +! get the cell size ratios +! + dyx = ady(lv) * adxi(lv) + dyz = ady(lv) * adzi(lv) + +! process left and right side boundary separatelly +! + if (jc == 1) then + +! iterate over left-side ghost layers +! + do j = jbl, 1, -1 + +! calculate neighbor cell indices +! + jp1 = min(jm, j + 1) + jp2 = min(jm, j + 2) + +! calculate variable decay coefficients +! + fr = (ady(lv) * (jb - j - 5.0d-01)) / blim + fl = 1.0d+00 - fr + +! iterate over boundary layer +! + do k = kl, ku +#if NDIMS == 3 + km1 = max( 1, k - 1) + kp1 = min(km, k + 1) +#endif /* NDIMS == 3 */ + do i = il, iu + im1 = max( 1, i - 1) + ip1 = min(im, i + 1) + +! make normal derivatives zero +! + qn(1:nv,i,j,k) = qn(1:nv,i,jb,k) + +! decay density and pressure to their limits +! + qn(idn,i,j,k) = fl * qn(idn,i,jb,k) + fr * dens + if (ipr > 0) qn(ipr,i,j,k) = fl * qn(ipr,i,jb,k) + fr * pres + +! decay magnetic field to its limit +! + qn(ibx,i,j,k) = fl * qn(ibx,i,jb,k) - fr * bamp + qn(ibz,i,j,k) = fl * qn(ibz,i,jb,k) + fr * bgui + +! update By from div(B)=0 +! + qn(iby,i,j,k) = qn(iby,i,jp2,k) & + + (qn(ibx,ip1,jp1,k) - qn(ibx,im1,jp1,k)) * dyx +#if NDIMS == 3 + qn(iby,i,j,k) = qn(iby,i,j ,k) & + + (qn(ibz,i,jp1,kp1) - qn(ibz,i,jp1,km1)) * dyz +#endif /* NDIMS == 3 */ + qn(ibp,i,j,k) = 0.0d+00 + end do ! i = il, iu + end do ! k = kl, ku + end do ! j = jbl, 1, -1 + else ! jc = 1 + +! iterate over right-side ghost layers +! + do j = jeu, jm + +! calculate neighbor cell indices +! + jm1 = max( 1, j - 1) + jm2 = max( 1, j - 2) + +! calculate variable decay coefficients +! + fr = (ady(lv) * (j - je - 5.0d-01)) / blim + fl = 1.0d+00 - fr + +! iterate over boundary layer +! + do k = kl, ku +#if NDIMS == 3 + km1 = max( 1, k - 1) + kp1 = min(km, k + 1) +#endif /* NDIMS == 3 */ + do i = il, iu + im1 = max( 1, i - 1) + ip1 = min(im, i + 1) + +! make normal derivatives zero +! + qn(1:nv,i,j,k) = qn(1:nv,i,je,k) + +! decay density and pressure to their limits +! + qn(idn,i,j,k) = fl * qn(idn,i,je,k) + fr * dens + if (ipr > 0) qn(ipr,i,j,k) = fl * qn(ipr,i,je,k) + fr * pres + +! decay magnetic field to its limit +! + qn(ibx,i,j,k) = fl * qn(ibx,i,je,k) + fr * bamp + qn(ibz,i,j,k) = fl * qn(ibz,i,je,k) + fr * bgui + +! update By from div(B)=0 +! + qn(iby,i,j,k) = qn(iby,i,jm2,k) & + + (qn(ibx,im1,jm1,k) - qn(ibx,ip1,jm1,k)) * dyx +#if NDIMS == 3 + qn(iby,i,j,k) = qn(iby,i,j ,k) & + + (qn(ibz,i,jm1,km1) - qn(ibz,i,jm1,kp1)) * dyz +#endif /* NDIMS == 3 */ + qn(ibp,i,j,k) = 0.0d+00 + end do ! i = il, iu + end do ! k = kl, ku + end do ! j = jeu, jm + end if ! jc = 1 + else ! ibx > 0 + if (jc == 1) then + do j = jbl, 1, -1 + qn(1:nv,il:iu,j,kl:ku) = qn(1:nv,il:iu,jb,kl:ku) + qn(ivy ,il:iu,j,kl:ku) = min(0.0d+00, qn(ivy,il:iu,jb,kl:ku)) + end do ! j = jbl, 1, -1 + else + do j = jeu, jm + qn(1:nv,il:iu,j,kl:ku) = qn(1:nv,il:iu,je,kl:ku) + qn(ivy ,il:iu,j,kl:ku) = max(0.0d+00, qn(ivy,il:iu,je,kl:ku)) + end do ! j = jeu, jm + end if + end if ! ibx > 0 + ! "reflective" boundary conditions ! case(bnd_reflective) From 77a668215efb180f326d5fe6f06a16d8dedef49d Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 13 Oct 2014 19:34:24 -0300 Subject: [PATCH 10/10] BOUNDARIES: Implement horizontal reconnection boundary conditions. These boundaries force J=0 condition near the current sheet and calculate Bx from div(B)=0 condition. They also prevent the inflow. Signed-off-by: Grzegorz Kowal --- src/boundaries.F90 | 188 +++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 183 insertions(+), 5 deletions(-) diff --git a/src/boundaries.F90 b/src/boundaries.F90 index 3b9d2d7..2b0e1fd 100644 --- a/src/boundaries.F90 +++ b/src/boundaries.F90 @@ -154,6 +154,8 @@ module boundaries bnd_type(1,1) = bnd_outflow case("reflective", "reflecting", "reflect") bnd_type(1,1) = bnd_reflective + case("reconnection", "recon", "rec") + bnd_type(1,1) = bnd_reconnection case default bnd_type(1,1) = bnd_periodic end select @@ -165,6 +167,8 @@ module boundaries bnd_type(1,2) = bnd_outflow case("reflective", "reflecting", "reflect") bnd_type(1,2) = bnd_reflective + case("reconnection", "recon", "rec") + bnd_type(1,2) = bnd_reconnection case default bnd_type(1,2) = bnd_periodic end select @@ -6041,7 +6045,7 @@ module boundaries use coordinates , only : im , jm , km , ng use coordinates , only : ib , jb , kb , ie , je , ke use coordinates , only : ibl, jbl, kbl, ieu, jeu, keu - use coordinates , only : ady, adxi, adzi + use coordinates , only : adx, ady, adxi, adyi, adzi use equations , only : nv use equations , only : idn, ipr, ivx, ivy, ivz, ibx, iby, ibz, ibp use error , only : print_error, print_warning @@ -6076,13 +6080,13 @@ module boundaries integer :: iu, ju, ku integer :: is, js, ks integer :: it, jt, kt - integer :: im1, ip1 + integer :: im2, im1, ip1, ip2 integer :: jm2, jm1, jp1, jp2 #if NDIMS == 3 - integer :: km1, kp1 + integer :: km2, km1, kp1, kp2 #endif /* NDIMS == 3 */ - real(kind=8) :: dyx, dyz - real(kind=8) :: fl, fr + real(kind=8) :: dxy, dxz, dyx, dyz + real(kind=8) :: fl, fr, ds ! !------------------------------------------------------------------------------- ! @@ -6169,6 +6173,180 @@ module boundaries end do ! i = ieu, im end if +! "reconnection" boundary conditions +! + case(bnd_reconnection) + +! process case with magnetic field, otherwise revert to standard outflow +! + if (ibx > 0) then + +! get the cell size ratios +! + dxy = adx(lv) * adyi(lv) + dxz = adx(lv) * adzi(lv) + +! process left and right side boundary separatelly +! + if (ic == 1) then + +! iterate over left-side ghost layers +! + do i = ibl, 1, -1 + +! calculate neighbor cell indices +! + ip1 = min(im, i + 1) + ip2 = min(im, i + 2) + +! iterate over boundary layer +! + do k = kl, ku +#if NDIMS == 3 + km2 = max( 1, k - 2) + km1 = max( 1, k - 1) + kp1 = min(km, k + 1) + kp2 = min(km, k + 2) +#endif /* NDIMS == 3 */ + do j = jl, ju + jm2 = max( 1, j - 2) + jm1 = max( 1, j - 1) + jp1 = min(jm, j + 1) + jp2 = min(jm, j + 2) + +! make normal derivatives zero +! + qn(1:nv,i,j,k) = qn(1:nv,ib,j,k) + +! prevent the inflow +! + qn(ivx,i,j,k) = min(0.0d+00, qn(ivx,ib,j,k)) + +! prevent from creating strong reconnection at the boundary (only near +! the plane of Bx polarity change) +! +#if NDIMS == 3 +! detect the current sheet +! + ds = min(qn(ibx,ib,jm2,k) * qn(ibx,ib,jp2,k) & + , qn(ibx,ib,j,km2) * qn(ibx,ib,j,kp2)) + +! apply curl-free condition in the vicinity of current sheet +! + if (ds < 0.0d+00) then + qn(iby,i,j,k) = qn(iby,ip2,j,k) & + + (qn(ibx,ip1,jm1,k) - qn(ibx,ip1,jp1,k)) * dxy + qn(ibz,i,j,k) = qn(ibz,ip2,j,k) & + + (qn(ibx,ip1,j,km1) - qn(ibx,ip1,j,kp1)) * dxz + end if +#else /* NDIMS == 3 */ +! apply curl-free condition in the vicinity of current sheet +! + if (qn(ibx,ib,jm2,k) * qn(ibx,ib,jp2,k) < 0.0d+00) then + qn(iby,i,j,k) = qn(iby,ip2,j,k) & + + (qn(ibx,ip1,jm1,k) - qn(ibx,ip1,jp1,k)) * dxy + end if +#endif /* NDIMS == 3 */ + +! update Bx from div(B)=0 +! + qn(ibx,i,j,k) = qn(ibx,ip2,j,k) & + + (qn(iby,ip1,jp1,k) - qn(iby,ip1,jm1,k)) * dxy +#if NDIMS == 3 + qn(ibx,i,j,k) = qn(ibx,i ,j,k) & + + (qn(ibz,ip1,j,kp1) - qn(ibz,ip1,j,km1)) * dxz +#endif /* NDIMS == 3 */ + qn(ibp,i,j,k) = 0.0d+00 + end do ! j = jl, ju + end do ! k = kl, ku + end do ! i = ibl, 1, -1 + else ! ic == 1 + +! iterate over right-side ghost layers +! + do i = ieu, im + +! calculate neighbor cell indices +! + im1 = max( 1, i - 1) + im2 = max( 1, i - 2) + +! iterate over boundary layer +! + do k = kl, ku +#if NDIMS == 3 + km1 = max( 1, k - 1) + kp1 = min(km, k + 1) + km2 = max( 1, k - 2) + kp2 = min(km, k + 2) +#endif /* NDIMS == 3 */ + do j = jl, ju + jm1 = max( 1, j - 1) + jp1 = min(jm, j + 1) + jm2 = max( 1, j - 2) + jp2 = min(jm, j + 2) + +! make normal derivatives zero +! + qn(1:nv,i,j,k) = qn(1:nv,ie,j,k) + +! prevent the inflow +! + qn(ivx,i,j,k) = max(0.0d+00, qn(ivx,ie,j,k)) + +! prevent from creating strong reconnection at the boundary (only near +! the plane of Bx polarity change) +! +#if NDIMS == 3 +! detect the current sheet +! + ds = min(qn(ibx,ie,jm2,k) * qn(ibx,ie,jp2,k) & + , qn(ibx,ie,j,km2) * qn(ibx,ie,j,kp2)) + +! apply curl-free condition in the vicinity of current sheet +! + if (ds < 0.0d+00) then + qn(iby,i,j,k) = qn(iby,im2,j,k) & + + (qn(ibx,im1,jp1,k) - qn(ibx,im1,jm1,k)) * dxy + qn(ibz,i,j,k) = qn(ibz,im2,j,k) & + + (qn(ibx,im1,j,kp1) - qn(ibx,im1,j,km1)) * dxz + end if +#else /* NDIMS == 3 */ +! apply curl-free condition in the vicinity of current sheet +! + if (qn(ibx,ie,jm2,k) * qn(ibx,ie,jp2,k) < 0.0d+00) then + qn(iby,i,j,k) = qn(iby,im2,j,k) & + + (qn(ibx,im1,jp1,k) - qn(ibx,im1,jm1,k)) * dxy + end if +#endif /* NDIMS == 3 */ + +! update Bx from div(B)=0 +! + qn(ibx,i,j,k) = qn(ibx,im2,j,k) & + + (qn(iby,im1,jm1,k) - qn(iby,im1,jp1,k)) * dxy +#if NDIMS == 3 + qn(ibx,i,j,k) = qn(ibx,i ,j,k) & + + (qn(ibz,im1,j,km1) - qn(ibz,im1,j,kp1)) * dxz +#endif /* NDIMS == 3 */ + qn(ibp,i,j,k) = 0.0d+00 + end do ! j = jl, ju + end do ! k = kl, ku + end do ! i = ieu, im + end if ! ic == 1 + else ! ibx > 0 + if (ic == 1) then + do i = ibl, 1, -1 + qn(1:nv,i,jl:ju,kl:ku) = qn(1:nv,ib,jl:ju,kl:ku) + qn(ivx ,i,jl:ju,kl:ku) = min(0.0d+00, qn(ivx,ib,jl:ju,kl:ku)) + end do ! i = ibl, 1, -1 + else + do i = ieu, im + qn(1:nv,i,jl:ju,kl:ku) = qn(1:nv,ie,jl:ju,kl:ku) + qn(ivx ,i,jl:ju,kl:ku) = max(0.0d+00, qn(ivx,ie,jl:ju,kl:ku)) + end do ! i = ieu, im + end if + end if ! ibx > 0 + ! "reflective" boundary conditions ! case(bnd_reflective)