! Initial condition (density, magnetic field, velocity)
! for a particular configuration of a magnetic tube.
!
! 04-aug-10/simon: created from trefoil.f90
!
!** AUTOMATIC CPARAM.INC GENERATION ****************************
! Declare (for generation of cparam.inc) the number of f array
! variables and auxiliary variables added by this module
!
! CPARAM logical, parameter :: linitial_condition = .true.
!
!***************************************************************
module InitialCondition
!
use Cparam
use Cdata
use General, only: keep_compiler_quiet
use Messages
!
implicit none
!
include '../initial_condition.h'
!
! ampl = amplitude of the magnetic field
! width_ring = width of the flux tube
! C = offset from the center of the knot, should be > 1
! D = extention coefficient for the z direction
! [xyz]scale = scale in each dimension
! [xyz]shift = shift in each dimension in a 2*pi box
! twist = B_phi/B_0
real :: ampl=1.0,width_ring=0.3, C=2.0, D=2.0
real :: xscale = 1.0, yscale = 1.0, zscale = 1.0
real :: xshift = 0.2, yshift = 0.5, zshift = 0.0
real :: twist = 0
integer :: n_foil = 3
character (len=labellen) :: prof='constant'
!
namelist /initial_condition_pars/ &
ampl,width_ring,prof,C,D,xscale,yscale,zscale,xshift,yshift,zshift,n_foil,twist
!
contains
!***********************************************************************
subroutine register_initial_condition()
!
! Configure pre-initialised (i.e. before parameter read) variables
! which should be know to be able to evaluate
!
! 07-oct-09/wlad: coded
!
! Identify CVS/SVN version information.
!
if (lroot) call svn_id( &
"$Id$")
!
endsubroutine register_initial_condition
!***********************************************************************
subroutine initial_condition_uu(f)
!
! Initialize the velocity field.
!
! 07-may-09/wlad: coded
!
real, dimension (mx,my,mz,mfarray), intent(inout) :: f
!
call keep_compiler_quiet(f)
!
endsubroutine initial_condition_uu
!***********************************************************************
subroutine initial_condition_lnrho(f)
!
! Initialize logarithmic density. init_lnrho
! will take care of converting it to linear
! density if you use ldensity_nolog
!
! 07-may-09/wlad: coded
!
real, dimension (mx,my,mz,mfarray), intent(inout) :: f
!
call keep_compiler_quiet(f)
!
endsubroutine initial_condition_lnrho
!***********************************************************************
subroutine initial_condition_aa(f)
!
! Initialize the magnetic vector potential.
!
! 04-aug-10/simon: coded
!
! Magnetic flux ring which has the form of a n-foil knot.
! NB: Curerently this works only for one CPU. For multi CPU usage
! initialize on one CPU and do the remeshing.
!
! Created 2010-08-04 by Simon Candelaresi (Iomsn)
!
use Mpicomm, only: stop_it
use Poisson
use Sub
real, dimension (mx,my,mz,mfarray) :: f
real :: knot_param, circle_param, circle_radius
real :: delta_knot_param, delta_circle_param, delta_circle_radius
real :: curvature_radius, curvature_factor
real, dimension(3) :: knot_pos, circle_pos, tangent, normal, ee_phi
integer :: domain_width, domain_depth, domain_height
integer :: l, j, ju
! The next 2 variables are used for the uncurling.
real, dimension (nx,ny,nz,3) :: jj, tmpJ ! This is phi for poisson.f90
!
! initialize the magnetic flux tube
!
domain_width = l2-l1; domain_depth = m2-m1; domain_height = n2-n1
!
! Calculate the minimum step size of the curve parameters
! to avoid discretation issues, like mesh points without magnetic field
!
delta_knot_param = min(dx, dy, dz)/10 * 3./n_foil
! delta_knot_param = min(dx, dy, dz)/2
delta_circle_param = delta_knot_param/(width_ring/2.)
delta_circle_radius = delta_circle_param
!
knot_param = 0.
!
! loop which moves along the n-foil knot
!
do
if (knot_param .gt. 2.*pi) exit
! Position along the central spine.
knot_pos(1) = (C+sin(knot_param*n_foil))*sin(knot_param*(n_foil-1))
knot_pos(2) = (C+sin(knot_param*n_foil))*cos(knot_param*(n_foil-1))
knot_pos(3) = D*cos(knot_param*n_foil)
! Tangent to the trajectory and direction of the magnetic field.
tangent(1) = n_foil*cos(knot_param*n_foil)*sin(knot_param*(n_foil-1))+&
(n_foil-1)*(C+sin(knot_param*n_foil))*cos(knot_param*(n_foil-1))
tangent(2) = n_foil*cos(knot_param*n_foil)*cos(knot_param*(n_foil-1))-&
(n_foil-1)*(C+sin(knot_param*n_foil))*sin(knot_param*(n_foil-1))
tangent(3) = -D*n_foil*sin(knot_param*n_foil)
tangent = tangent / sqrt(tangent(1)**2+tangent(2)**2+tangent(3)**2)
! Normal component using the length of the curve as parametrization.
normal(1) = (n_foil*sin(knot_param*(n_foil - 1.0d0))*cos(knot_param*n_foil) + (C + &
sin(knot_param*n_foil))*(n_foil - 1)*cos(knot_param*(n_foil - &
1.0d0)))*(-D**2*n_foil**3*sin(knot_param*n_foil)*cos(knot_param* &
n_foil) - 1.0d0/2.0d0*(n_foil*sin(knot_param*(n_foil - 1.0d0))* &
cos(knot_param*n_foil) + (C + sin(knot_param*n_foil))*(n_foil - 1 &
)*cos(knot_param*(n_foil - 1.0d0)))*(-2*n_foil**2*sin(knot_param* &
n_foil)*sin(knot_param*(n_foil - 1.0d0)) + 4*n_foil*(n_foil - 1)* &
cos(knot_param*n_foil)*cos(knot_param*(n_foil - 1.0d0)) - 2*(C + &
sin(knot_param*n_foil))*(n_foil - 1)**2*sin(knot_param*(n_foil - &
1.0d0))) - 1.0d0/2.0d0*(n_foil*cos(knot_param*n_foil)*cos( &
knot_param*(n_foil - 1.0d0)) - (C + sin(knot_param*n_foil))*( &
n_foil - 1)*sin(knot_param*(n_foil - 1.0d0)))*(-2*n_foil**2*sin( &
knot_param*n_foil)*cos(knot_param*(n_foil - 1.0d0)) - 4*n_foil*( &
n_foil - 1)*sin(knot_param*(n_foil - 1.0d0))*cos(knot_param* &
n_foil) - 2*(C + sin(knot_param*n_foil))*(n_foil - 1)**2*cos( &
knot_param*(n_foil - 1.0d0))))/(D**2*n_foil**2*sin(knot_param* &
n_foil)**2 + (n_foil*sin(knot_param*(n_foil - 1.0d0))*cos( &
knot_param*n_foil) + (C + sin(knot_param*n_foil))*(n_foil - 1)* &
cos(knot_param*(n_foil - 1.0d0)))**2 + (n_foil*cos(knot_param* &
n_foil)*cos(knot_param*(n_foil - 1.0d0)) - (C + sin(knot_param* &
n_foil))*(n_foil - 1)*sin(knot_param*(n_foil - 1.0d0)))**2)**( &
3.0d0/2.0d0) + (-n_foil**2*sin(knot_param*n_foil)*sin(knot_param* &
(n_foil - 1.0d0)) + 2*n_foil*(n_foil - 1)*cos(knot_param*n_foil)* &
cos(knot_param*(n_foil - 1.0d0)) - (C + sin(knot_param*n_foil))*( &
n_foil - 1)**2*sin(knot_param*(n_foil - 1.0d0)))/sqrt(D**2*n_foil &
**2*sin(knot_param*n_foil)**2 + (n_foil*sin(knot_param*(n_foil - &
1.0d0))*cos(knot_param*n_foil) + (C + sin(knot_param*n_foil))*( &
n_foil - 1)*cos(knot_param*(n_foil - 1.0d0)))**2 + (n_foil*cos( &
knot_param*n_foil)*cos(knot_param*(n_foil - 1.0d0)) - (C + sin( &
knot_param*n_foil))*(n_foil - 1)*sin(knot_param*(n_foil - 1.0d0 &
)))**2)
normal(2) = (n_foil*cos(knot_param*n_foil)*cos(knot_param*(n_foil - 1.0d0)) - (C + &
sin(knot_param*n_foil))*(n_foil - 1)*sin(knot_param*(n_foil - &
1.0d0)))*(-D**2*n_foil**3*sin(knot_param*n_foil)*cos(knot_param* &
n_foil) - 1.0d0/2.0d0*(n_foil*sin(knot_param*(n_foil - 1.0d0))* &
cos(knot_param*n_foil) + (C + sin(knot_param*n_foil))*(n_foil - 1 &
)*cos(knot_param*(n_foil - 1.0d0)))*(-2*n_foil**2*sin(knot_param* &
n_foil)*sin(knot_param*(n_foil - 1.0d0)) + 4*n_foil*(n_foil - 1)* &
cos(knot_param*n_foil)*cos(knot_param*(n_foil - 1.0d0)) - 2*(C + &
sin(knot_param*n_foil))*(n_foil - 1)**2*sin(knot_param*(n_foil - &
1.0d0))) - 1.0d0/2.0d0*(n_foil*cos(knot_param*n_foil)*cos( &
knot_param*(n_foil - 1.0d0)) - (C + sin(knot_param*n_foil))*( &
n_foil - 1)*sin(knot_param*(n_foil - 1.0d0)))*(-2*n_foil**2*sin( &
knot_param*n_foil)*cos(knot_param*(n_foil - 1.0d0)) - 4*n_foil*( &
n_foil - 1)*sin(knot_param*(n_foil - 1.0d0))*cos(knot_param* &
n_foil) - 2*(C + sin(knot_param*n_foil))*(n_foil - 1)**2*cos( &
knot_param*(n_foil - 1.0d0))))/(D**2*n_foil**2*sin(knot_param* &
n_foil)**2 + (n_foil*sin(knot_param*(n_foil - 1.0d0))*cos( &
knot_param*n_foil) + (C + sin(knot_param*n_foil))*(n_foil - 1)* &
cos(knot_param*(n_foil - 1.0d0)))**2 + (n_foil*cos(knot_param* &
n_foil)*cos(knot_param*(n_foil - 1.0d0)) - (C + sin(knot_param* &
n_foil))*(n_foil - 1)*sin(knot_param*(n_foil - 1.0d0)))**2)**( &
3.0d0/2.0d0) + (-n_foil**2*sin(knot_param*n_foil)*cos(knot_param* &
(n_foil - 1.0d0)) - 2*n_foil*(n_foil - 1)*sin(knot_param*(n_foil &
- 1.0d0))*cos(knot_param*n_foil) - (C + sin(knot_param*n_foil))*( &
n_foil - 1)**2*cos(knot_param*(n_foil - 1.0d0)))/sqrt(D**2*n_foil &
**2*sin(knot_param*n_foil)**2 + (n_foil*sin(knot_param*(n_foil - &
1.0d0))*cos(knot_param*n_foil) + (C + sin(knot_param*n_foil))*( &
n_foil - 1)*cos(knot_param*(n_foil - 1.0d0)))**2 + (n_foil*cos( &
knot_param*n_foil)*cos(knot_param*(n_foil - 1.0d0)) - (C + sin( &
knot_param*n_foil))*(n_foil - 1)*sin(knot_param*(n_foil - 1.0d0 &
)))**2)
normal(3) = -D*n_foil**2*cos(knot_param*n_foil)/sqrt(D**2*n_foil**2*sin(knot_param* &
n_foil)**2 + (n_foil*sin(knot_param*(n_foil - 1.0d0))*cos( &
knot_param*n_foil) + (C + sin(knot_param*n_foil))*(n_foil - 1)* &
cos(knot_param*(n_foil - 1.0d0)))**2 + (n_foil*cos(knot_param* &
n_foil)*cos(knot_param*(n_foil - 1.0d0)) - (C + sin(knot_param* &
n_foil))*(n_foil - 1)*sin(knot_param*(n_foil - 1.0d0)))**2) - D* &
n_foil*(-D**2*n_foil**3*sin(knot_param*n_foil)*cos(knot_param* &
n_foil) - 1.0d0/2.0d0*(n_foil*sin(knot_param*(n_foil - 1.0d0))* &
cos(knot_param*n_foil) + (C + sin(knot_param*n_foil))*(n_foil - 1 &
)*cos(knot_param*(n_foil - 1.0d0)))*(-2*n_foil**2*sin(knot_param* &
n_foil)*sin(knot_param*(n_foil - 1.0d0)) + 4*n_foil*(n_foil - 1)* &
cos(knot_param*n_foil)*cos(knot_param*(n_foil - 1.0d0)) - 2*(C + &
sin(knot_param*n_foil))*(n_foil - 1)**2*sin(knot_param*(n_foil - &
1.0d0))) - 1.0d0/2.0d0*(n_foil*cos(knot_param*n_foil)*cos( &
knot_param*(n_foil - 1.0d0)) - (C + sin(knot_param*n_foil))*( &
n_foil - 1)*sin(knot_param*(n_foil - 1.0d0)))*(-2*n_foil**2*sin( &
knot_param*n_foil)*cos(knot_param*(n_foil - 1.0d0)) - 4*n_foil*( &
n_foil - 1)*sin(knot_param*(n_foil - 1.0d0))*cos(knot_param* &
n_foil) - 2*(C + sin(knot_param*n_foil))*(n_foil - 1)**2*cos( &
knot_param*(n_foil - 1.0d0))))*sin(knot_param*n_foil)/(D**2* &
n_foil**2*sin(knot_param*n_foil)**2 + (n_foil*sin(knot_param*( &
n_foil - 1.0d0))*cos(knot_param*n_foil) + (C + sin(knot_param* &
n_foil))*(n_foil - 1)*cos(knot_param*(n_foil - 1.0d0)))**2 + ( &
n_foil*cos(knot_param*n_foil)*cos(knot_param*(n_foil - 1.0d0)) - &
(C + sin(knot_param*n_foil))*(n_foil - 1)*sin(knot_param*(n_foil &
- 1.0d0)))**2)**(3.0d0/2.0d0)
curvature_radius = 1/sqrt(normal(1)**2+normal(2)**2+normal(3)**2)
normal = normal*curvature_radius
! !
! ! Find vector which is orthonormal to tangent vector.
! !
! if (abs(tangent(1)) .le. 0.5) then
! normal(1) = tangent(1)**2 - 1.0
! normal(2) = tangent(2)*tangent(1)
! normal(3) = tangent(3)*tangent(1)
! elseif (abs(tangent(2)) .le. 0.5) then
! normal(1) = tangent(1)*tangent(2)
! normal(2) = tangent(2)**2 - 1.0
! normal(3) = tangent(3)*tangent(2)
! else
! normal(1) = tangent(1)*tangent(3)
! normal(2) = tangent(2)*tangent(3)
! normal(3) = tangent(3)**2 - 1.0
! endif
! !
! ! ! normalize the normal vector
! !
! normal = normal / sqrt(normal(1)**2+normal(2)**2+normal(3)**2)
circle_radius = 0.
!
! loop which changes the circle's radius
!
do
if (circle_radius .gt. width_ring/2.) exit
circle_param = 0.
!
! loop which goes around the circle
!
do
if (circle_param .gt. 2.*pi) exit
circle_pos(1) = knot_pos(1) + circle_radius * &
((tangent(1)*tangent(1)*(1-cos(circle_param))+cos(circle_param))*normal(1) + &
(tangent(1)*tangent(2)*(1-cos(circle_param))-tangent(3)*sin(circle_param))*normal(2) + &
(tangent(1)*tangent(3)*(1-cos(circle_param))+tangent(2)*sin(circle_param))*normal(3))
circle_pos(2) = knot_pos(2) + circle_radius * &
((tangent(1)*tangent(2)*(1-cos(circle_param))+tangent(3)*sin(circle_param))*normal(1) + &
(tangent(2)*tangent(2)*(1-cos(circle_param))+cos(circle_param))*normal(2) + &
(tangent(2)*tangent(3)*(1-cos(circle_param))-tangent(1)*sin(circle_param))*normal(3))
circle_pos(3) = knot_pos(3) + circle_radius * &
((tangent(1)*tangent(3)*(1-cos(circle_param))-tangent(2)*sin(circle_param))*normal(1) + &
(tangent(2)*tangent(3)*(1-cos(circle_param))+tangent(1)*sin(circle_param))*normal(2) + &
(tangent(3)*tangent(3)*(1-cos(circle_param))+cos(circle_param))*normal(3))
!
! Add the azimuthal field in case of twist > 0.
!
if (sum(abs(circle_pos-knot_pos)) > 0) then
ee_phi(1) = tangent(2)*(circle_pos(3)-knot_pos(3)) - tangent(3)*(circle_pos(2)-knot_pos(2))
ee_phi(2) = tangent(3)*(circle_pos(1)-knot_pos(1)) - tangent(1)*(circle_pos(3)-knot_pos(3))
ee_phi(3) = tangent(1)*(circle_pos(2)-knot_pos(2)) - tangent(2)*(circle_pos(1)-knot_pos(1))
ee_phi = ee_phi/sqrt(ee_phi(1)**2 + ee_phi(2)**2 + ee_phi(3)**2)
else
ee_phi = [0, 0, 0]
endif
!
! Correct for the loop curvature with fields inside being weakend and outside being strengthend.
!
! curvature_factor = 1-sum(normal*(circle_pos-knot_pos))/curvature_radius
curvature_factor = 1
!
! Find the corresponding mesh point to this position.
!
l = nint((circle_pos(1)*xscale+xshift)/(2*pi)*nxgrid + nxgrid/2.) + 1 - nx*ipx
m = nint((circle_pos(2)*yscale+yshift)/(2*pi)*nygrid + nygrid/2.) + 1 - ny*ipy
n = nint((circle_pos(3)*zscale+zshift)/(2*pi)*nzgrid + nzgrid/2.) + 1 - nz*ipz
!
! Write the magnetic field B.
! Note that B is written in the f-array where A is stored.
! This is corrected further in the code.
!
if ((l > mx .or. m > my .or. n > mz .or. l < 1 .or. m < 1 .or. n < 1) .eqv. .false.) then
f(l,m,n,iax:iaz) = (tangent+ee_phi*twist*sqrt(sum((circle_pos-knot_pos)**2))) * &
ampl*curvature_factor
endif
circle_param = circle_param + delta_circle_param
enddo
circle_radius = circle_radius + delta_circle_radius
enddo
knot_param = knot_param + delta_knot_param
enddo
!
! Transform the magnetic field into a vector potential
!
! Compute curl(B) = J for the Poisson solver
do m=m1,m2
do n=n1,n2
call curl(f,iaa,jj(:,m-nghost,n-nghost,:))
enddo
enddo
tmpJ = -jj
! Use the Poisson solver to solve \nabla^2 A = -J for A
do j=1,3
call inverse_laplacian(tmpJ(:,:,:,j))
enddo
! Overwrite the f-array with the correct vector potential A
do j=1,3
ju=iaa-1+j
f(l1:l2,m1:m2,n1:n2,ju) = tmpJ(:,:,:,j)
enddo
!
endsubroutine initial_condition_aa
!***********************************************************************
subroutine read_initial_condition_pars(iostat)
!
use File_io, only: parallel_unit
!
integer, intent(out) :: iostat
!
read(parallel_unit, NML=initial_condition_pars, IOSTAT=iostat)
!
endsubroutine read_initial_condition_pars
!***********************************************************************
subroutine write_initial_condition_pars(unit)
!
integer, intent(in) :: unit
!
write(unit, NML=initial_condition_pars)
!
endsubroutine write_initial_condition_pars
!***********************************************************************
!********************************************************************
!************ DO NOT DELETE THE FOLLOWING **************
!********************************************************************
!** This is an automatically generated include file that creates **
!** copies dummy routines from noinitial_condition.f90 for any **
!** InitialCondition routines not implemented in this file **
!** **
include '../initial_condition_dummies.inc'
!********************************************************************
endmodule InitialCondition