diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml new file mode 100644 index 0000000..75d7c92 --- /dev/null +++ b/.github/workflows/test.yaml @@ -0,0 +1,55 @@ +name: Test +on: [push] + +jobs: + standard_tests: + name: Tests with ${{ matrix.config.name }} + timeout-minutes: 60 + runs-on: ubuntu-latest + container: + image: oi4ai/bout3d:db-outer-${{ matrix.config.mode }} + strategy: + fail-fast: true + matrix: + config: + - name: debugging + mode: debug + - name: optimisations + mode: opt + steps: + - name: Job information + run: | + echo Build: ${{ matrix.config.name }}, ${{ matrix.config.os }} + echo CMake options: ${{ matrix.config.cmake_options }} + cat /etc/os-release + + - name: Setup dependencies + run: | + sudo dnf install -y git-lfs + + - name: Checkout hermes-2 + run: | + cd /home/boutuser/ + git clone ${{ github.server_url }}/${{ github.repository }} -b ${{ github.ref_name }} hermes-2 + cd hermes-2 + git checkout ${{ github.sha }} + git submodule update --init --recursive + + - name: Build + run: | + cd /home/boutuser/hermes-2 + cmake --version + cmake -S . -B build -Dbout++_DIR=/home/boutuser/BOUT-dev/build + make -C build -j 2 + + - name: Run + run: | + cd /home/boutuser/hermes-2 + # Run only shortly + sed -e 's/nout = .*/nout = 1/' -i */BOUT.inp + sed -e 's/timestep = 100/timestep = 1/' -i */BOUT.inp + # Change solver type - beuler is quite heavy on setup, but we want + # to run only short + sed -e 's/type = beuler/type = pvode/' -i */BOUT.inp + build/hermes-2 -d 1-lowres + build/hermes-2 -d 3-medium diff --git a/.gitignore b/.gitignore index 38425c0..2bc2e53 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,13 @@ *.o *.nc BOUT.log.* - +.BOUT.pid.* +BOUT.settings +/atomicpp.a +/doc/hermes-2-manual.aux +/doc/hermes-2-manual.log +/doc/hermes-2-manual.out +/doc/hermes-2-manual.pdf +/hermes-2 +/core* +/Makefile diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..c81be56 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "grids"] + path = grids + url = https://github.com/dschwoerer/hermes-grids.git diff --git a/1-lowres/BOUT.inp b/1-lowres/BOUT.inp index 2a1e604..f2659b1 100644 --- a/1-lowres/BOUT.inp +++ b/1-lowres/BOUT.inp @@ -1,23 +1,23 @@ -NOUT = 100 -TIMESTEP = 10 +nout = 10 +timestep = 100 -TwistShift = false # use twist-shift condition? +twistshift = false # use twist-shift condition? ballooning = false # ballooning transformation? shiftinitial = false -ShiftWithoutTwist=false +ShiftWithoutTwist = false -MYG=1 -MXG=2 +MYG = 1 +MXG = 2 -grid = "rotating-ellipse-36x8x32-ix20.fci.nc" +grid = "grids/rotating-ellipse-36x8x32-ix20.fci.nc" -NXPE=1 +NXPE = 1 dump_format = nc -[restart] -init_missing=true # initialize missing variables on restart? +[restart_files] +init_missing = true # initialize missing variables on restart? [mesh] symmetricGlobalX = true @@ -30,19 +30,19 @@ type = fci ################################################## # derivative methods -[ddx] +[mesh:ddx] first = C2 second = C2 upwind = W3 -[ddy] +[mesh:ddy] first = C2 second = C2 upwind = W3 -[ddz] +[mesh:ddz] first = C2 second = C2 @@ -55,10 +55,11 @@ upwind = W3 # Note: If evolving neutrals, need preconditioning # type = cvode +type = beuler # use_precon = true -ATOL = 1.0e-10 # absolute tolerance -RTOL = 1.0e-5 # relative tolerance +atol = 1.0e-10 # absolute tolerance +rtol = 1.0e-5 # relative tolerance mxstep = 1000000 # Maximum internal steps per output cvode_max_order = 2 @@ -72,7 +73,7 @@ nonuniform = true # Electrostatic potential solver [phiSolver] -type = petsc # Needed if Boussinesq = false +type = petsc # Needed if Boussinesq = false ksptype = gmres # Linear iterative method pctype = lu # Preconditioner. Direct "lu" or "ilu"; iterative "jacobi", "sor" @@ -115,7 +116,8 @@ nonuniform = true # general settings for the model [input] -transform_from_field_aligned=False +transform_from_field_aligned = false +error_on_unused_options = false [Hermes] loadmetric = false @@ -140,7 +142,7 @@ resistivity_boundary_width = 0 # Numerical dissipation vepsi_dissipation = false # Parallel dissipation on Ve-Vi -vort_dissipation = false +vort_dissipation = true # false # Changed 2021-10-25 numdiff = 0.01 hyperpar = -1 @@ -171,16 +173,18 @@ split_n0_psi = false j_diamag = false # Diamagnetic current: Vort <-> Pe j_par = false # Parallel current: Vort <-> Psi -evolve_te = false +evolve_te = true evolve_ti = false +evolve_vort = false + pe_par = true # Parallel pressure gradient: Pe <-> Psi resistivity = true # Resistivity: Psi -> Pe thermal_flux = false thermal_force = false electron_ion_transfer = false electron_viscosity = false -ion_viscosity = false # Ion parallel viscosity +ion_viscosity = false # Ion parallel viscosity thermal_conduction = false frecycle = 0.0 # Neutral gas recycling fraction @@ -199,8 +203,8 @@ drift_wave = false # Transport coefficients classical_diffusion = false # Collisional diffusion -anomalous_D = -1 # Anomalous density diffusion [m^2/s] -anomalous_chi = -1 # Anomalous thermal diffusion [m^2/s] +anomalous_D = 0.2 # Anomalous density diffusion [m^2/s] +anomalous_chi = 0.6 # Anomalous thermal diffusion [m^2/s] anomalous_nu = -1 # Anomalous viscosity poloidal_flows = false @@ -208,7 +212,7 @@ poloidal_flows = false magnetic_drift = true ion_velocity = true -ion_neutral = 0.0 +ion_neutral = false neutral_friction = false # Friction between plasma and neutrals boussinesq = true # Use Boussinesq approximation @@ -229,6 +233,7 @@ sheath_ydown = true sheath_allow_supersonic = false sheath_gamma_e = 4 # Electron sheath heat transmission sheath_gamma_i = 2.5 # Ion sheath heat transmission +parallel_sheaths = true neutral_gamma = 0.0 startprofiles = false @@ -254,18 +259,21 @@ AA = 1 # Atomic mass. 1 = Hydrogen, 2 = Deuterium type = none # Neutral model: none, diffusion2d, recycling, fullvelocity, mixed viscosity = 1 # Dynamic viscosity bulk = 0 # Bulk (volume) viscosity -conduction = 1 +conduction = 1 neutral_gamma = 0.0 nn_floor = 1e-2 # Floor applied when calculating Vn = NVn / Nn low_n_equilibriate = -1e-4 # If n < nn_floor, evolve Tn and Vn towards plasma values -[All] +[all] scale = 0.0 bndry_all = neumann_o2 bndry_xin = neumann_o2 bndry_xout = neumann_o2 +# Maybe not documented? +bndry_par_all = parallel_neumann + [Ne] # Electron density scale = 1 @@ -275,9 +283,11 @@ function = 0.1 + 0.01*gauss(x,0.21) #function = 1 + 0.1*gauss(x-0.5, 0.14)*gauss(z-pi/2.,0.14/2.) + 0.1*gauss(x-0.5, 0.14)*gauss(z-pi,0.14/2.) + 0.1*gauss(x-0.5, 0.14)*gauss(z,0.14/2.) + 0.1*gauss(x-0.5, 0.14)*gauss(z-3*pi/2.,0.14/2.) source = 2.8e3*gauss(x, 0.21) +bndry_xin = neumann_o2 +#dirichlet(0.1) [Vort] -scale=0 +scale = 0 function = sin(7*z)*sin(3*x-z)*gauss(-x,0.21)*(1+sin(y))#mixmode(y) bndry_all = dirichlet_o2 @@ -291,20 +301,22 @@ scale = 1 function = 0.1 + 0.01*cos(z)*(sin(8*pi*x - 30*z)+sin(7*z)+0.5*sin(13*z - 2*pi*x))*gauss(-x,0.21) #function = 1 + 1.2*(Ne:function-1) source = Ne:source +# bndry_par_all = parallel_free [Pi] scale = 1 #function = 1 + 1.2*(Ne:function-1) function = 0.1 + 0.01*cos(z)*(sin(8*pi*x - 30*z)+sin(7*z)+0.5*sin(13*z - 2*pi*x))*gauss(-x,0.21) source = Ne:source +# bndry_xin = dirichlet(1) [Ve] -parallel_bndry_yup=1 -parallel_bndry_ydown=-1 +parallel_bndry_yup = 1 +parallel_bndry_ydown = -1 [Jpar] -parallel_bndry_yup=1 -parallel_bndry_ydown=-1 +parallel_bndry_yup = 1 +parallel_bndry_ydown = -1 [phi] # Radial boundaries determined by Laplacian inversion @@ -348,4 +360,3 @@ bndry_all = neumann_o2 [NVn] bndry_all = dirichlet_o2 - diff --git a/3-medium/BOUT.inp b/3-medium/BOUT.inp new file mode 100644 index 0000000..61e23bf --- /dev/null +++ b/3-medium/BOUT.inp @@ -0,0 +1,362 @@ +nout = 1 +timestep = 1e-3 + +twistshift = false # use twist-shift condition? +ballooning = false # ballooning transformation? +shiftinitial = false +ShiftWithoutTwist=false + +MYG=1 +MXG=2 + +grid = "grids/W7X-conf4-132x16x512.fci.f2.nc" + +NXPE=1 + +dump_format = nc + +[restart_files] +init_missing=true # initialize missing variables on restart? + +[mesh] +symmetricGlobalX = true + +extrapolate_y = false # Extrapolate metrics into Y boundary cells? + +[mesh:paralleltransform] +type = fci + +################################################## +# derivative methods + +[ddx] + +first = C2 +second = C2 +upwind = W3 + +[ddy] + +first = C2 +second = C2 +upwind = W3 + +[ddz] + +first = C2 +second = C2 +upwind = W3 + +################################################### +# Time-integration solver + +[solver] + +# Note: If evolving neutrals, need preconditioning +# type = cvode +type = beuler +# use_precon = true + +atol = 1.0e-10 # absolute tolerance +rtol = 1.0e-5 # relative tolerance +mxstep = 1000000 # Maximum internal steps per output + +cvode_max_order = 2 +cvode_stability_limit_detection = true + +[laplace] # This is used for Delp2 operator +all_terms = true +nonuniform = true + +################################################## +# Electrostatic potential solver + +[phiSolver] +type = petsc # Needed if Boussinesq = false +ksptype = gmres # Linear iterative method +pctype = lu # Preconditioner. Direct "lu" or "ilu"; iterative "jacobi", "sor" + +# Set package to perform factorisation for direct solves +# "petsc" built-in solver only serial +# "superlu", "superlu_dist", "mumps", "cusparse" +factor_package = superlu_dist + +inner_boundary_flags = 2 + +all_terms = true #false +nonuniform = true # NOTE: Necessary to avoid numerical instability + +[laplacexy] # 2D solver in X-Y +pctype = sor # Preconditioner + +atol = 1e-12 +rtol = 1e-8 + +core_bndry_dirichlet = false +pf_bndry_dirichlet = true +y_bndry_dirichlet = false + +include_y_derivs = true + +[laplacexz] +type = petsc +inner_boundary_flags = 32 +outer_boundary_flags = 0 + + +[aparSolver] +type = petsc +inner_boundary_flags = 0 +outer_boundary_flags = 0 + +all_terms = true +nonuniform = true + +# general settings for the model +[input] + +transform_from_field_aligned = false +error_on_unused_options = false + +[Hermes] +loadmetric = false +newXZsolver = true + +####################### +# Output variables +output_ddt = false # Output time derivatives +verbose = true # Output additional fields + +####################### +# radial buffers +radial_buffers = false +radial_inner_width = 4 +radial_outer_width = 4 +radial_buffer_D = 0.1 + +resistivity_boundary = 0 +resistivity_boundary_width = 0 + +####################### +# Numerical dissipation + +vepsi_dissipation = true # Parallel dissipation on Ve-Vi +vort_dissipation = true # false # Changed 2021-10-25 +numdiff = 0.01 +hyperpar = -1 + +ne_num_diff = 1e-4 +vi_num_diff = 1e-4 +ve_num_diff = -1 +ve_num_hyper = 1e-1 + +x_hyper_viscos = 1e-3 +z_hyper_viscos = 1e-3 + +# Flux limiters +kappa_limit_alpha = -1#0.2 # SOLPS style heat flux limiter +eta_limit_alpha = -1#0.5 # SOLPS style viscosity limiter + +####################### +# Electric field and Ohm's law +electromagnetic = true # Electromagnetic? Otherwise electrostatic +FiniteElMass = false # Finite electron mass? + +parallel_flow = true + +# Electrostatic potential +split_n0 = false # Solve n=0 separately? +split_n0_psi = false + +# NOTE: all currents switched off for fluid run +j_diamag = false # Diamagnetic current: Vort <-> Pe +j_par = false # Parallel current: Vort <-> Psi + +evolve_te = true +evolve_ti = true + +evolve_vort = true + +pe_par = true # Parallel pressure gradient: Pe <-> Psi +resistivity = true # Resistivity: Psi -> Pe +thermal_flux = true +thermal_force = true +electron_ion_transfer = true +electron_viscosity = false +ion_viscosity = false # Ion parallel viscosity +thermal_conduction = true + +frecycle = 0.9 # Neutral gas recycling fraction + +carbon_fraction = 0.0 + +excitation = false # Hydrogen neutral excitation radiation + +## Settings for 2D parallel closures +sinks = false +sink_invlpar = 0.2 # 5m parallel connection length +sheath_closure = false +drift_wave = false + +####################### +# Transport coefficients +classical_diffusion = true # Collisional diffusion + +anomalous_D = 1 # Anomalous density diffusion [m^2/s] +anomalous_chi = 3 # Anomalous thermal diffusion [m^2/s] +anomalous_nu = 1 # Anomalous viscosity + +poloidal_flows = false + +magnetic_drift = true +ion_velocity = true + +ion_neutral = true +neutral_friction = true # Friction between plasma and neutrals + +boussinesq = true # Use Boussinesq approximation + +# Radial boundary fluxes +ne_bndry_flux = true +pe_bndry_flux = true +vort_bndry_flux = true + +ramp_mesh = false +ramp_timescale = 1e4 + +####################### +# Plasma sheath +sheath_model = 0 # 0 = Bohm, 1 = Loizu, 2 = Bohm + free density +sheath_yup = true +sheath_ydown = true +sheath_allow_supersonic = false +sheath_gamma_e = 4 # Electron sheath heat transmission +sheath_gamma_i = 2.5 # Ion sheath heat transmission +parallel_sheaths = true +neutral_gamma = 0.0 + +startprofiles = false + +core_sources = false # Only sources in the core +adapt_source = false # Feedback on profiles (PI controller) +energy_source = false # Constant power per particle (at fixed x). False -> Constant power per volume +# source_p = 1e-2 # Rate based on current error (1/timescale) +# source_i = 1e-6 # Integral of error (1/timescale^2) +source_vary_g11 = false # Vary source in poloidal angle to better match radial transport + +staggered = false + +# Normalisation parameters + +# Normalisation factors +Nnorm = 1e20 +Tnorm = 100 +Bnorm = 1 +AA = 1 # Atomic mass. 1 = Hydrogen, 2 = Deuterium + +[neutral] +type = diffusion2d # Neutral model: none, diffusion2d, recycling, fullvelocity, mixed +viscosity = 1 # Dynamic viscosity +bulk = 0 # Bulk (volume) viscosity +conduction = 1 +neutral_gamma = 0.0 + +nn_floor = 1e-2 # Floor applied when calculating Vn = NVn / Nn +low_n_equilibriate = -1e-4 # If n < nn_floor, evolve Tn and Vn towards plasma values + +[all] +scale = 0.0 + +bndry_all = neumann_o2 +bndry_xin = neumann_o2 +bndry_xout = neumann_o2 +# Maybe not documented? +bndry_par_all = parallel_neumann + + +[Ne] # Electron density +scale = 1 + +#function = 0.1 + 0.01*sin(8*pi*x - 9*z)*gauss(x, 0.21)#1 + 0.27*gauss(x, 0.42)#*gauss(z-3*pi/2.,0.21/4.)#*gauss(y-pi,0.26) +function = 0.1 + 0.01*gauss(x,0.21) +#function = 1 + 0.1*gauss(x-0.5, 0.14)*gauss(z-pi/2.,0.14/2.) + 0.1*gauss(x-0.5, 0.14)*gauss(z-pi,0.14/2.) + 0.1*gauss(x-0.5, 0.14)*gauss(z,0.14/2.) + 0.1*gauss(x-0.5, 0.14)*gauss(z-3*pi/2.,0.14/2.) + +source = 2.8e3*gauss(x, 0.21) +bndry_xin = neumann_o2 +#dirichlet(0.1) + +[Vort] +scale=0 +function = sin(7*z)*sin(3*x-z)*gauss(-x,0.21)*(1+sin(y))#mixmode(y) +bndry_all = dirichlet_o2 + +[VePsi] # Ve + 0.5*beta_e*mi_me*psi +bndry_core = dirichlet_o2#zerolaplace +bndry_pf = dirichlet_o2 +bndry_xout = dirichlet_o2 + +[Pe] # Electron pressure +scale = 1 +function = 0.1 + 0.01*cos(z)*(sin(8*pi*x - 30*z)+sin(7*z)+0.5*sin(13*z - 2*pi*x))*gauss(-x,0.21) +#function = 1 + 1.2*(Ne:function-1) +source = Ne:source +# bndry_par_all = parallel_free + +[Pi] +scale = 1 +#function = 1 + 1.2*(Ne:function-1) +function = 0.1 + 0.01*cos(z)*(sin(8*pi*x - 30*z)+sin(7*z)+0.5*sin(13*z - 2*pi*x))*gauss(-x,0.21) +source = Ne:source +# bndry_xin = dirichlet(1) + +[Ve] +parallel_bndry_yup=1 +parallel_bndry_ydown=-1 + +[Jpar] +parallel_bndry_yup=1 +parallel_bndry_ydown=-1 + +[phi] +# Radial boundaries determined by Laplacian inversion +bndry_xin = none +bndry_xout = none + +bndry_all = dirichlet_o2 + +[Nn] + +scale = 5e-2 +function = 1 + +[Vn] + +scale = 0 +function = 0 + +bndry_all = dirichlet_o2 + +[Vn_x] +scale = 0 +function = 0 +bndry_all = dirichlet_o2 + +[Vn_y] +scale = 0 +function = 0 +bndry_all = dirichlet_o2 + +[Vn_z] +scale = 0 +function = 0 +bndry_all = dirichlet_o2 + +[Pn] + +scale = 1e-5 +function = 1 +bndry_all = neumann_o2 + +[NVn] +bndry_all = dirichlet_o2 + diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..0a3e527 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,8 @@ +cmake_minimum_required(VERSION 3.9...3.12) + +project(hermes-2 LANGUAGES CXX) +set(CMAKE_CXX_STANDARD 14) +find_package(bout++ REQUIRED) +add_executable(hermes-2 hermes-2.cxx div_ops.cxx loadmetric.cxx radiation.cxx neutral-model.cxx diffusion2d.cxx recycling.cxx full-velocity.cxx mixed.cxx atomicpp/ImpuritySpecies.cxx atomicpp/Prad.cxx atomicpp/RateCoefficient.cxx atomicpp/sharedFunctions.cxx) +target_link_libraries(hermes-2 PRIVATE bout++::bout++) + diff --git a/atomicpp/makefile b/atomicpp/makefile index 2371b47..198ccff 100644 --- a/atomicpp/makefile +++ b/atomicpp/makefile @@ -1,4 +1,4 @@ -BOUT_TOP=../../.. +BOUT_TOP ?= ../../.. SOURCEC = ImpuritySpecies.cxx Prad.cxx RateCoefficient.cxx sharedFunctions.cxx SOURCEH = $(SOURCEC:%.cxx=%.hxx) json.hxx diff --git a/diffusion2d.cxx b/diffusion2d.cxx index fb112ca..da32206 100644 --- a/diffusion2d.cxx +++ b/diffusion2d.cxx @@ -18,84 +18,83 @@ Diffusion2D::Diffusion2D(Solver *solver, Mesh*, Options &options) : NeutralModel SAVE_REPEAT(Dnn); Lmax = options["Lmax"].doc("Maximum mean free path [m]").withDefault(1.0); - - // Set Laplacian inversion to null - inv = 0; -} - -Diffusion2D::~Diffusion2D() { - if(inv) { - delete &inv; - } } void Diffusion2D::update(const Field3D &Ne, const Field3D &Te, const Field3D &UNUSED(Ti), const Field3D &UNUSED(Vi)) { mesh->communicate(Nn, Pn); - - Nn = floor(Nn, 1e-8); - Field3D Tn = Pn / Nn; - Tn = floor(Tn, 0.01/Tnorm); - - Field3D Nelim = floor(Ne, 1e-19); // Smaller limit for rate coefficients - + // Calculate atomic processes - for(int i=0;iLocalNx;i++) - for(int j=0;jLocalNy;j++) - for(int k=0;kLocalNz;k++) { - // Charge exchange frequency, normalised to ion cyclotron frequency - BoutReal sigma_cx = Nelim(i,j,k)*Nnorm*hydrogen.chargeExchange(Te(i,j,k)*Tnorm)/Fnorm; - - // Ionisation frequency, normalised to ion cyclotron frequency - BoutReal sigma_iz = Nelim(i,j,k)*Nnorm*Nn(i,j,k) * hydrogen.ionisation(Te(i,j,k)*Tnorm)/Fnorm; - - // Neutral thermal velocity - BoutReal vth_n = sqrt(Tn(i,j,k)); // Normalised to Cs0 - - // Neutral-neutral mean free path - BoutReal a0 = PI*SQ(5.29e-11); - BoutReal lambda_nn = 1. / (Nnorm*Nn(i,j,k)*a0); // meters - if(lambda_nn > Lmax) { - // Limit maximum mean free path - lambda_nn = Lmax; - } - - lambda_nn /= Lnorm; // Normalised length to Lnorm - // Neutral-Neutral collision rate, normalised to ion cyclotron frequency - BoutReal sigma_nn = vth_n / lambda_nn; - - // Total neutral collision frequency, normalised to ion cyclotron frequency - BoutReal sigma = sigma_cx + sigma_nn + sigma_iz; - - // Neutral gas diffusion - Dnn(i,j,k) = SQ(vth_n) / sigma; - - // Rates - BoutReal R_rc = SQ(Nelim(i,j,k)) * hydrogen.recombination(Nelim(i,j,k)*Nnorm, Te(i,j,k)*Tnorm) * Nnorm / Fnorm; // Time normalisation - BoutReal R_iz = Nelim(i,j,k)*Nn(i,j,k) * hydrogen.ionisation(Te(i,j,k)*Tnorm) * Nnorm / Fnorm; // Time normalisation - BoutReal R_cx = sigma_cx * Nn(i,j,k); - - // Plasma sink / neutral source - S(i,j,k) = R_rc - R_iz; - - // Power transfer from plasma to neutrals - - Qi(i,j,k) = R_cx * (3./2)*(Te(i,j,k) - Tn(i,j,k)); - - // Power transfer due to ionisation and recombination - Qi(i,j,k) += (3./2)*(Te(i,j,k)*R_rc - Tn(i,j,k)*R_iz); - - // Ion-neutral friction - Fperp(i,j,k) = R_cx // Charge-Exchange - + R_rc; // Recombination - - // Radiated power from plasma - // Factor of 1.09 so that recombination becomes an energy source at 5.25eV - Rp(i,j,k) = (1.09*Te(i,j,k) - 13.6/Tnorm)*R_rc - + (Eionize/Tnorm) * R_iz; // Ionisation energy - - } - + BOUT_FOR(i, Ne.getRegion("RGN_ALL")) { + Nn[i] = std::max(Nn[i], 1e-8); + BoutReal Tn = Pn[i] / Nn[i]; + Tn = std::max(Tn, 0.01 / Tnorm); + + BoutReal Nelim = + std::max(Ne[i], 1e-19); // Smaller limit for rate coefficients + + // Charge exchange frequency, normalised to ion cyclotron frequency + BoutReal sigma_cx = + Nelim * Nnorm * hydrogen.chargeExchange(Te[i] * Tnorm) / Fnorm; + + // Ionisation frequency, normalised to ion cyclotron frequency + BoutReal sigma_iz = + Nelim * Nnorm * Nn[i] * hydrogen.ionisation(Te[i] * Tnorm) / Fnorm; + + // Neutral thermal velocity + BoutReal vth_n = sqrt(Tn); // Normalised to Cs0 + + // Neutral-neutral mean free path + BoutReal a0 = PI * SQ(5.29e-11); + BoutReal lambda_nn = 1. / (Nnorm * Nn[i] * a0); // meters + if (lambda_nn > Lmax) { + // Limit maximum mean free path + lambda_nn = Lmax; + } + + lambda_nn /= Lnorm; // Normalised length to Lnorm + // Neutral-Neutral collision rate, normalised to ion cyclotron frequency + BoutReal sigma_nn = vth_n / lambda_nn; + + // Total neutral collision frequency, normalised to ion cyclotron frequency + BoutReal sigma = sigma_cx + sigma_nn + sigma_iz; + + // Neutral gas diffusion + Dnn[i] = SQ(vth_n) / sigma; + + // Rates + BoutReal R_rc = SQ(Nelim) * + hydrogen.recombination(Nelim * Nnorm, Te[i] * Tnorm) * + Nnorm / Fnorm; // Time normalisation + BoutReal R_iz = Nelim * Nn[i] * hydrogen.ionisation(Te[i] * Tnorm) * Nnorm / + Fnorm; // Time normalisation + BoutReal R_cx = sigma_cx * Nn[i]; + + // Plasma sink / neutral source + S[i] = R_rc - R_iz; + + // Power transfer from plasma to neutrals + + Qi[i] = R_cx * (3. / 2) * (Te[i] - Tn); + + // Power transfer due to ionisation and recombination + Qi[i] += (3. / 2) * (Te[i] * R_rc - Tn * R_iz); + + // Ion-neutral friction + Fperp[i] = R_cx // Charge-Exchange + + R_rc; // Recombination + + // Radiated power from plasma + // Factor of 1.09 so that recombination becomes an energy source at 5.25eV + Rp[i] = (1.09 * Te[i] - 13.6 / Tnorm) * R_rc + + (Eionize / Tnorm) * R_iz; // Ionisation energy + } + mesh->communicate(Dnn, Fperp); + Dnn.applyParallelBoundary("parallel_neumann"); + Fperp.applyParallelBoundary("parallel_neumann"); + Nn.applyParallelBoundary(); + Pn.applyParallelBoundary(); + // Neutral density ddt(Nn) = + S diff --git a/diffusion2d.hxx b/diffusion2d.hxx index e5e3ffa..2607596 100644 --- a/diffusion2d.hxx +++ b/diffusion2d.hxx @@ -8,11 +8,14 @@ class Diffusion2D : public NeutralModel { public: Diffusion2D(Solver *solver, Mesh *mesh, Options &options); - ~Diffusion2D(); + ~Diffusion2D() = default; void update(const Field3D &Ne, const Field3D &Te, const Field3D &Ti, const Field3D &Vi); void precon(BoutReal t, BoutReal gamma, BoutReal delta); + + virtual Field3D getDensity() override { return Nn; } + private: Field3D Nn; // Neutral gas density (evolving) Field3D Pn; // Neutral gas pressure (evolving) diff --git a/grids b/grids new file mode 160000 index 0000000..a90eacc --- /dev/null +++ b/grids @@ -0,0 +1 @@ +Subproject commit a90eacc32305599818e896590d37dae3f34a0a47 diff --git a/hermes-2.cxx b/hermes-2.cxx index ef0c773..674ef32 100644 --- a/hermes-2.cxx +++ b/hermes-2.cxx @@ -64,7 +64,7 @@ namespace FV { Coordinates *coord = f_in.getCoordinates(); Field3D result{zeroFrom(f)}; - + // Only need one guard cell, so no need to communicate fluxes // Instead calculate in guard cells to preserve fluxes int ys = mesh->ystart-1; @@ -206,7 +206,7 @@ namespace FV { } } -BoutReal floor(const BoutReal &var, const BoutReal &f) { +BoutReal floor(BoutReal var, BoutReal f) { if (var < f) return f; return var; @@ -230,67 +230,173 @@ const Field3D ceil(const Field3D &var, BoutReal f, REGION rgn = RGN_ALL) { // Square function for vectors Field3D SQ(const Vector3D &v) { return v * v; } -Field3D floor_all(const Field3D &f, BoutReal minval) { - Field3D result = floor(f, minval); - result.splitParallelSlices(); - result.yup() = floor(f.yup(), minval); - result.ydown() = floor(f.ydown(), minval); - return result; +void setRegions(Field3D &f) { + f.yup().setRegion("RGN_YPAR_+1"); + f.ydown().setRegion("RGN_YPAR_-1"); } -Field3D copy_all(const Field3D &f) { - Field3D result = copy(f); - result.splitParallelSlices(); - result.yup() = copy(f.yup()); - result.ydown() = copy(f.ydown()); - return result; -} - -Field3D div_all(const Field3D &num, const Field3D &den) { - Field3D result = num / den; - result.splitParallelSlices(); - result.yup() = num.yup() / den.yup(); - result.ydown() = num.ydown() / den.ydown(); - return result; -} - -Field3D div_all(const Field3D &num, BoutReal den) { - Field3D result = num / den; - result.splitParallelSlices(); - result.yup() = num.yup() / den; - result.ydown() = num.ydown() / den; - return result; -} - -Field3D mul_all(const Field3D &a, const Field3D &b) { - Field3D result = a * b; - result.splitParallelSlices(); - result.yup() = a.yup() * b.yup(); - result.ydown() = a.ydown() * b.ydown(); - return result; +const Field3D &yup(const Field3D &f) { return f.yup(); } +BoutReal yup(BoutReal f) { return f; }; +const Field3D &ydown(const Field3D &f) { return f.ydown(); } +BoutReal ydown(BoutReal f) { return f; }; +const BoutReal yup(BoutReal f, Ind3D i) { return f; }; +const BoutReal ydown(BoutReal f, Ind3D i) { return f; }; +// const BoutReal& yup(const Field3D &f, Ind3D i) { return f.yup()[i.yp()]; } +// const BoutReal& ydown(const Field3D &f, Ind3D i) { return f.ydown()[i.ym()]; +// } BoutReal& yup(Field3D &f, Ind3D i) { return f.yup()[i.yp()]; } BoutReal& +// ydown(Field3D &f, Ind3D i) { return f.ydown()[i.ym()]; } +const BoutReal &yup(const Field3D &f, Ind3D i) { return f.yup()[i]; } +const BoutReal &ydown(const Field3D &f, Ind3D i) { return f.ydown()[i]; } +BoutReal &yup(Field3D &f, Ind3D i) { return f.yup()[i]; } +BoutReal &ydown(Field3D &f, Ind3D i) { return f.ydown()[i]; } +const BoutReal &_get(const Field3D &f, Ind3D i) { return f[i]; } +BoutReal &_get(Field3D &f, Ind3D i) { return f[i]; } +BoutReal _get(BoutReal f, Ind3D i) { return f; }; +BoutReal copy(BoutReal f) { return f; }; + +void alloc_all(Field3D &f) { + f.allocate(); + f.splitParallelSlices(); + f.yup().allocate(); + f.ydown().allocate(); + setRegions(f); } -Field3D sub_all(const Field3D &a, const Field3D &b) { - Field3D result = a - b; - result.splitParallelSlices(); - result.yup() = a.yup() - b.yup(); - result.ydown() = a.ydown() - b.ydown(); - return result; +#define GET_ALL(name) \ + auto *name##a = &name[Ind3D(0)]; \ + auto *name##b = &name.yup()[Ind3D(0)]; \ + auto *name##c = &name.ydown()[Ind3D(0)]; + +#define DO_ALL(op, name) \ + template Field3D name##_all(const A &a, const B &b) { \ + Field3D result; \ + alloc_all(result); \ + BOUT_FOR(i, result.getRegion("RGN_ALL")) { name##_all(result, a, b, i); } \ + setRegions(result); \ + return result; \ + } \ + template \ + void name##_all(Field3D &result, const A &a, const B &b, Ind3D i) { \ + result[i] = op(_get(a, i), _get(b, i)); \ + yup(result, i) = op(yup(a, i), yup(b, i)); \ + ydown(result, i) = op(ydown(a, i), ydown(b, i)); \ + } \ + template void name##_all(Field3D &result, const B &b, Ind3D i) { \ + result[i] = op(result[i], _get(b, i)); \ + yup(result, i) = op(yup(result, i), yup(b, i)); \ + ydown(result, i) = op(ydown(result, i), ydown(b, i)); \ + } + +DO_ALL(floor, floor) +DO_ALL(pow, pow) + +#undef DO_ALL +#define DO_ALL(op, name) \ + template Field3D name##_all(const A &a, const B &b) { \ + Field3D result; \ + alloc_all(result); \ + BOUT_FOR(i, result.getRegion("RGN_ALL")) { name##_all(result, a, b, i); } \ + checkData(result, "RGN_ALL"); \ + setRegions(result); \ + return result; \ + } \ + Field3D name##_all(const Field3D &a, const Field3D &b) { \ + Field3D result; \ + alloc_all(result); \ + const int n = result.getNx() * result.getNy() * result.getNz(); \ + GET_ALL(result); \ + GET_ALL(a); \ + GET_ALL(b); \ + BOUT_OMP(omp parallel for simd) \ + for (int i = 0; i < n; ++i) { \ + resulta[i] = aa[i] op ba[i]; \ + resultb[i] = ab[i] op bb[i]; \ + resultc[i] = ac[i] op bc[i]; \ + } \ + setRegions(result); \ + return result; \ + } \ + Field3D name##_all(const Field3D &a, BoutReal b) { \ + Field3D result; \ + alloc_all(result); \ + const int n = result.getNx() * result.getNy() * result.getNz(); \ + GET_ALL(result); \ + GET_ALL(a); \ + BOUT_OMP(omp parallel for simd) \ + for (int i = 0; i < n; ++i) { \ + resulta[i] = aa[i] op b; \ + resultb[i] = ab[i] op b; \ + resultc[i] = ac[i] op b; \ + } \ + setRegions(result); \ + return result; \ + } \ + template \ + void name##_all(Field3D &result, const A &a, const B &b, Ind3D i) { \ + result[i] = _get(a, i) op _get(b, i); \ + yup(result, i) = yup(a, i) op yup(b, i); \ + ydown(result, i) = ydown(a, i) op ydown(b, i); \ + } + +// void div_all(Field3D & result, const Field3D & a, const Field3D & b, Ind3D i) +// { +// result[i] = a[i] / b[i]; +// result.yup()[i.yp()] = a.yup()[i.yp()] / b.yup()[i.yp()]; +// result.ydown()[i.ym()] = a.ydown()[i.ym()] / b.ydown()[i.ym()]; +// } + +//#include "mul_all.cxx" +DO_ALL(*, mul) +DO_ALL(/, div) +DO_ALL(+, add) +DO_ALL(-, sub) + +#undef DO_ALL +#define DO_ALL(op) \ + inline void op##_all(Field3D &result, const Field3D &a, Ind3D i) { \ + result[i] = op(a[i]); \ + yup(result, i) = op(yup(a, i)); \ + ydown(result, i) = op(ydown(a, i)); \ + } \ + inline Field3D op##_all(const Field3D &a) { \ + Field3D result; \ + alloc_all(result); \ + BOUT_FOR(i, result.getRegion("RGN_ALL")) { op##_all(result, a, i); } \ + checkData(result, "RGN_ALL"); \ + setRegions(result); \ + return result; \ + } + +DO_ALL(sqrt) +DO_ALL(SQ) +DO_ALL(copy) +DO_ALL(exp) +DO_ALL(log) + +#undef DO_ALL + +void set_all(Field3D &f, BoutReal val) { + alloc_all(f); + BOUT_FOR(i, f.getRegion("RGN_ALL")) { + f[i] = val; + f.yup()[i] = val; + f.ydown()[i] = val; + } } +void zero_all(Field3D &f) { set_all(f, 0); } -Field3D add_all(const Field3D &a, const Field3D &b) { - Field3D result = a + b; - result.splitParallelSlices(); - result.yup() = a.yup() + b.yup(); - result.ydown() = a.ydown() + b.ydown(); - return result; +void check_all(Field3D &f) { + checkData(f); + checkData(f.yup()); + checkData(f.ydown()); } -void zero_all(Field3D &f) { - f = 0.0; - f.splitParallelSlices(); - f.yup() = 0.0; - f.ydown() = 0.0; +void ASSERT_CLOSE_ALL(const Field3D &a, const Field3D &b) { + BOUT_FOR(i, a.getRegion("RGN_NOY")) { + ASSERT0(std::abs(a[i] - b[i]) < 1e-10); + ASSERT0(std::abs(yup(a, i) - yup(b, i)) < 1e-10); + ASSERT0(std::abs(ydown(a, i) - ydown(b, i)) < 1e-10); + } } /// Modifies and returns the first argument, taking the boundary from second argument @@ -300,6 +406,45 @@ Field3D withBoundary(Field3D &&f, const Field3D &bndry) { return f; } +void apply_flux_bcs(Field3D &f1, Field3D &f2) { + // Boundary conditions on fluxes: + // * Zero gradient on the targets + // * Zero value on upstream conditions + const auto mesh = f1.getMesh(); + for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + f1(r.ind, mesh->ystart - 1, jz) = f1(r.ind, mesh->ystart, jz); + f2(r.ind, mesh->ystart - 1, jz) = f2(r.ind, mesh->ystart, jz); + } + } + for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { + for (int jz = 0; jz < mesh->LocalNz; jz++) { + f1(r.ind, mesh->yend + 1, jz) = f1(r.ind, mesh->yend, jz); + f2(r.ind, mesh->yend + 1, jz) = f2(r.ind, mesh->yend, jz); + } + } + const int ny = mesh->LocalNy; + const int nz = mesh->LocalNz; + for (const auto &bndry_par : mesh->getBoundariesPar(BoundaryParType::xout)) { + for (bndry_par->first(); !bndry_par->isDone(); bndry_par->next()) { + const Ind3D i{(bndry_par->x * ny + bndry_par->y) * nz + bndry_par->z, ny, + nz}; + const auto inext = i.yp(bndry_par->dir); + f1.ynext(bndry_par->dir)[inext] = f1[i]; + f2.ynext(bndry_par->dir)[inext] = f2[i]; + } + } + for (const auto &bndry_par : mesh->getBoundariesPar(BoundaryParType::xin)) { + for (bndry_par->first(); !bndry_par->isDone(); bndry_par->next()) { + const Ind3D i{(bndry_par->x * ny + bndry_par->y) * nz + bndry_par->z, ny, + nz}; + const auto inext = i.yp(bndry_par->dir); + f1.ynext(bndry_par->dir)[inext] = -f1[i]; + f2.ynext(bndry_par->dir)[inext] = -f2[i]; + } + } +} + int Hermes::init(bool restarting) { auto& opt = Options::root(); @@ -498,7 +643,7 @@ int Hermes::init(bool restarting) { OPTION(optsc, AA, 2.0); // Ion mass (2 = Deuterium) - output.write("Normalisation Te=%e, Ne=%e, B=%e\n", Tnorm, Nnorm, Bnorm); + output.write("Normalisation Te={:e}, Ne={:e}, B={:e}\n", Tnorm, Nnorm, Bnorm); SAVE_ONCE(Tnorm, Nnorm, Bnorm, AA); // Save Cs0 = sqrt(qe * Tnorm / (AA * Mp)); // Reference sound speed [m/s] @@ -512,7 +657,7 @@ int Hermes::init(bool restarting) { output.write("\tmi_me={}, beta_e={}\n", mi_me, beta_e); SAVE_ONCE(mi_me, beta_e, me_mi); - output.write("\t Cs=%e, rho_s=%e, Omega_ci=%e\n", Cs0, rho_s0, Omega_ci); + output.write("\t Cs={:e}, rho_s={:e}, Omega_ci={:e}\n", Cs0, rho_s0, Omega_ci); SAVE_ONCE(Cs0, rho_s0, Omega_ci); // Collision times @@ -523,30 +668,34 @@ int Hermes::init(bool restarting) { tau_i0 = sqrt(AA) / (4.78e-8 * (Nnorm / 1e6) * lambda_ii * pow(Tnorm, -3. / 2)); - output.write("\ttau_e0=%e, tau_i0=%e\n", tau_e0, tau_i0); + output.write("\ttau_e0={:e}, tau_i0={:e}\n", tau_e0, tau_i0); if (anomalous_D > 0.0) { // Normalise anomalous_D /= rho_s0 * rho_s0 * Omega_ci; // m^2/s - output.write("\tnormalised anomalous D_perp = %e\n", anomalous_D); + output.write("\tnormalised anomalous D_perp = {:e}\n", anomalous_D); a_d3d = anomalous_D; mesh->communicate(a_d3d); - + a_d3d.yup() = anomalous_D; + a_d3d.ydown() = anomalous_D; } if (anomalous_chi > 0.0) { // Normalise anomalous_chi /= rho_s0 * rho_s0 * Omega_ci; // m^2/s - output.write("\tnormalised anomalous chi_perp = %e\n", anomalous_chi); + output.write("\tnormalised anomalous chi_perp = {:e}\n", anomalous_chi); a_chi3d = anomalous_chi; mesh->communicate(a_chi3d); + a_chi3d.yup() = anomalous_D; + a_chi3d.ydown() = anomalous_D; } if (anomalous_nu > 0.0) { // Normalise anomalous_nu /= rho_s0 * rho_s0 * Omega_ci; // m^2/s - output.write("\tnormalised anomalous nu_perp = %e\n", anomalous_nu); + output.write("\tnormalised anomalous nu_perp = {:e}\n", anomalous_nu); a_nu3d = anomalous_nu; mesh->communicate(a_nu3d); - + a_nu3d.yup() = anomalous_D; + a_nu3d.ydown() = anomalous_D; } if (ramp_mesh) { @@ -570,13 +719,15 @@ int Hermes::init(bool restarting) { alpha_dw = fact.create2D("Hermes:alpha_dw"); SAVE_ONCE(alpha_dw); } + } else { + optsc["sink_invlpar"].setConditionallyUsed(); } // Get switches from each variable section auto& optne = opt["Ne"]; NeSource = optne["source"].doc("Source term in ddt(Ne)").withDefault(Field3D{0.0}); NeSource /= Omega_ci; - Sn = DC(NeSource); + Sn = NeSource; // Inflowing density carries momentum OPTION(optne, density_inflow, false); @@ -584,12 +735,12 @@ int Hermes::init(bool restarting) { auto& optpe = opt["Pe"]; PeSource = optpe["source"].withDefault(Field3D{0.0}); PeSource /= Omega_ci; - Spe = DC(PeSource); + Spe = PeSource; auto& optpi = opt["Pi"]; PiSource = optpi["source"].withDefault(Field3D{0.0}); PiSource /= Omega_ci; - Spi = DC(PiSource); + Spi = PiSource; OPTION(optsc, core_sources, false); if (core_sources) { @@ -752,21 +903,17 @@ int Hermes::init(bool restarting) { coord->geometry(); // Calculate other metrics } - sqrtB = sqrt(coord->Bxy); // B^(1/2) - B32 = sqrtB * coord->Bxy; // B^(3/2) - - SAVE_ONCE(B32); - if (Options::root()["mesh:paralleltransform"]["type"].as() == "fci") { fci_transform = true; }else{ fci_transform = false; } - if(fci_transform){ poloidal_flows = false; mesh->get(Bxyz, "B",1.0); + mesh->get(coord->Bxy, "Bxy", 1.0); Bxyz /= Bnorm; + coord->Bxy /= Bnorm; SAVE_ONCE(Bxyz); Bxyz.applyBoundary("neumann"); ASSERT1( min(Bxyz) > 0.0 ); @@ -776,16 +923,32 @@ int Hermes::init(bool restarting) { // Note: A Neumann condition simplifies boundary conditions on fluxes // where the condition e.g. on J should be on flux (J/B) Bxyz.applyParallelBoundary("parallel_neumann"); - + coord->dz.applyParallelBoundary("parallel_neumann"); + coord->dy.applyParallelBoundary("parallel_neumann"); + coord->J.applyParallelBoundary("parallel_neumann"); + coord->g_22.applyParallelBoundary("parallel_neumann"); + coord->g_23.applyParallelBoundary("parallel_neumann"); + coord->g23.applyParallelBoundary("parallel_neumann"); + // coord->Bxy.applyParallelBoundary("parallel_neumann"); + auto logBxy = log(coord->Bxy); + logBxy.applyBoundary("neumann"); + mesh->communicate(logBxy); + logBxy.applyParallelBoundary("parallel_neumann"); + printf("Setting from log"); + coord->Bxy = exp_all(logBxy); + + bout::checkPositive(coord->Bxy, "f", "RGN_NOCORNERS"); + bout::checkPositive(coord->Bxy.yup(), "fyup", "RGN_YPAR_+1"); + bout::checkPositive(coord->Bxy.ydown(), "fdown", "RGN_YPAR_-1"); logB = log(Bxyz); bracket_factor = sqrt(coord->g_22) / (coord->J * Bxyz); SAVE_ONCE(bracket_factor); - }else{ - bracket_factor = sqrt(coord->g_22) / (coord->J * coord->Bxy); - SAVE_ONCE(bracket_factor); - } - + + B12 = sqrt_all(coord->Bxy); // B^(1/2) + B32 = mul_all(B12, coord->Bxy); // B^(3/2) + B42 = SQ_all(coord->Bxy); + ///////////////////////////////////////////////////////// // Neutral models @@ -795,6 +958,8 @@ int Hermes::init(bool restarting) { // Set normalisations if (neutrals) { neutrals->setNormalisation(Tnorm, Nnorm, Bnorm, rho_s0, Omega_ci); + } else { + opt["neutral"].setConditionallyUsed(); } ///////////////////////////////////////////////////////// @@ -961,7 +1126,15 @@ int Hermes::init(bool restarting) { ////////////////////////////////////////////////////////////// // Electromagnetic fields - + + opt["phiSolver"].setConditionallyUsed(); + opt["aparSolver"].setConditionallyUsed(); + opt["laplace"].setConditionallyUsed(); + opt["laplacexy"].setConditionallyUsed(); + opt["laplacexz"].setConditionallyUsed(); + optsc["newXZsolver"].setConditionallyUsed(); + optsc["split_n0"].setConditionallyUsed(); + optsc["split_n0_psi"].setConditionallyUsed(); if (j_par | j_diamag) { // Only needed if there are any currents SAVE_REPEAT(phi); @@ -995,7 +1168,7 @@ int Hermes::init(bool restarting) { // Test new LaplaceXZ solver newSolver = LaplaceXZ::create(bout::globals::mesh); // Set coefficients for Boussinesq solve - newSolver->setCoefs(Field3D(1.0), Field3D(0.0)); + newSolver->setCoefs(1. / SQ(coord->Bxy), Field3D(0.0)); } else { // Use older Laplacian solver phiSolver = Laplacian::create(&opt["phiSolver"]); @@ -1019,12 +1192,12 @@ int Hermes::init(bool restarting) { if (!restarting) { // Start by setting to the sheath current = 0 boundary value - Field3D Nelim = floor(Ne, 1e-5); - Telim = floor(Pe / Nelim, 0.1 / Tnorm); - Tilim = floor(Pi / Nelim, 0.1 / Tnorm); + Ne = floor(Ne, 1e-5); + Te = floor(Pe / Ne, 0.1 / Tnorm); + Ti = floor(Pi / Ne, 0.1 / Tnorm); phi.setBoundaryTo(DC( - (log(0.5 * sqrt(mi_me / PI)) + log(sqrt(Telim / (Telim + Tilim)))) * Telim)); + (log(0.5 * sqrt(mi_me / PI)) + log(sqrt(Te / (Te + Ti)))) * Te)); } // Set the last update time to -1, so it will reset @@ -1057,7 +1230,7 @@ int Hermes::init(bool restarting) { kappa_ipar = 0.0; Dn = 0.0; - SAVE_REPEAT(Telim, Tilim); + SAVE_REPEAT(Te, Ti); if (verbose) { // Save additional fields @@ -1089,13 +1262,11 @@ int Hermes::init(bool restarting) { // Preconditioner setPrecon((preconfunc)&Hermes::precon); - SAVE_REPEAT(a,b,c,d); if (evolve_te && evolve_ti && parallel_sheaths){ SAVE_REPEAT(sheath_dpe, sheath_dpi); } // Magnetic field in boundary auto& Bxy = mesh->getCoordinates()->Bxy; - mesh->communicate(Bxy); for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { @@ -1110,6 +1281,26 @@ int Hermes::init(bool restarting) { } } + opt["Jpar"].setConditionallyUsed(); + opt["Ve"].setConditionallyUsed(); + opt["Pn"].setConditionallyUsed(); + opt["Nn"].setConditionallyUsed(); + opt["NVn"].setConditionallyUsed(); + opt["Pe"].setConditionallyUsed(); + opt["Pi"].setConditionallyUsed(); + opt["Vn"].setConditionallyUsed(); + opt["Vn_x"].setConditionallyUsed(); + opt["Vn_y"].setConditionallyUsed(); + opt["Vn_z"].setConditionallyUsed(); + opt["phi"].setConditionallyUsed(); + opt["phiSolver"].setConditionallyUsed(); + opt["Vort"].setConditionallyUsed(); + opt["VePsi"].setConditionallyUsed(); + optsc["neutral_gamma"].setConditionallyUsed(); + + alloc_all(Te); + alloc_all(Ti); + alloc_all(Vi); return 0; } @@ -1119,11 +1310,6 @@ int Hermes::rhs(BoutReal t) { } Coordinates *coord = mesh->getCoordinates(); - a = 0; - b = 0; - c = 0; - d = 0; - f = 0; if (!evolve_plasma) { Ne = 0.0; Pe = 0.0; @@ -1138,27 +1324,59 @@ int Hermes::rhs(BoutReal t) { // Note: Parallel slices are not calculated because parallel derivatives // are calculated using field aligned quantities mesh->communicate(EvolvingVars); + Ne.applyParallelBoundary(); + Pe.applyParallelBoundary(); + Pi.applyParallelBoundary(); + NVi.applyParallelBoundary(); - Field3D Nelim = floor_all(Ne, 1e-5); - - if (!evolve_te) { - Pe = copy_all(Nelim); // Fixed electron temperature - } - - Te = div_all(Pe, Nelim); - Vi = div_all(NVi, Nelim); + // Are there any currents? If not, then there is no source + // for vorticity, phi = 0 and jpar = 0 + bool currents = j_par | j_diamag; + + // Local sound speed. Used for parallel advection operator + // Assumes isothermal electrons, adiabatic ions + // The factor scale_num_cs can be used to test sensitity + Field3D sound_speed; + sound_speed.allocate(); + + alloc_all(Te); + alloc_all(Ti); + alloc_all(Vi); + alloc_all(Pi); + BOUT_FOR(i, Ne.getRegion("RGN_ALL")) { + // Field3D Ne = floor_all(Ne, 1e-5); + floor_all(Ne, 1e-5, i); + + if (!evolve_te) { + copy_all(Pe, Ne, i); // Fixed electron temperature + } + + div_all(Te, Pe, Ne, i); + // ASSERT0(Te[i] > 1e-10); + /// printf("%f\n", Te[i]); + div_all(Vi, NVi, Ne, i); + + floor_all(Te, 0.1 / Tnorm, i); + // ASSERT0(Te[i] > 1e-10); + + mul_all(Pe, Te, Ne, i); - Telim = floor_all(Te, 0.1 / Tnorm); + if (!evolve_ti) { + copy_all(Pi, Ne, i); // Fixed ion temperature + } - Field3D Pelim = mul_all(Telim, Nelim); + div_all(Ti, Pi, Ne, i); + floor_all(Ti, 0.1 / Tnorm, i); + mul_all(Pi, Ti, Ne, i); + div_all(Te, Pe, Ne, i); + // ASSERT0(Te[i] > 1e-10); - if (!evolve_ti) { - Pi = copy_all(Nelim); // Fixed ion temperature + sound_speed[i] = scale_num_cs * sqrt(Te[i] + Ti[i] * (5. / 3)); + if (floor_num_cs > 0.0) { + // Apply a floor function to the sound speed + sound_speed[i] = floor(sound_speed[i], floor_num_cs); + } } - - Ti = div_all(Pi, Nelim); - Tilim = floor_all(Ti, 0.1 / Tnorm); - Field3D Pilim = mul_all(Tilim, Nelim); // Set radial boundary conditions on Te, Ti, Vi // @@ -1183,8 +1401,8 @@ int Hermes::rhs(BoutReal t) { if (ti_bndry < 0.1 / Tnorm) ti_bndry = 0.1 / Tnorm; - Telim(1, j, k) = 2. * te_bndry - Telim(2, j, k); - Tilim(1, j, k) = 2. * ti_bndry - Tilim(2, j, k); + Te(1, j, k) = 2. * te_bndry - Te(2, j, k); + Ti(1, j, k) = 2. * ti_bndry - Ti(2, j, k); } } } @@ -1210,33 +1428,19 @@ int Hermes::rhs(BoutReal t) { if (ti_bndry < 0.1 / Tnorm) ti_bndry = 0.1 / Tnorm; - Telim(n - 1, j, k) = 2. * te_bndry - Telim(n - 2, j, k); - Tilim(n - 1, j, k) = 2. * ti_bndry - Tilim(n - 2, j, k); + Te(n - 1, j, k) = 2. * te_bndry - Te(n - 2, j, k); + Ti(n - 1, j, k) = 2. * ti_bndry - Ti(n - 2, j, k); } } } - - // Are there any currents? If not, then there is no source - // for vorticity, phi = 0 and jpar = 0 - bool currents = j_par | j_diamag; - - // Local sound speed. Used for parallel advection operator - // Assumes isothermal electrons, adiabatic ions - // The factor scale_num_cs can be used to test sensitity - Field3D sound_speed = scale_num_cs * sqrt(Telim + Tilim * (5. / 3)); - if (floor_num_cs > 0.0) { - // Apply a floor function to the sound speed - sound_speed = floor(sound_speed, floor_num_cs); - } sound_speed.applyBoundary("neumann"); - + ////////////////////////////////////////////////////////////// // Calculate electrostatic potential phi // // TRACE("Electrostatic potential"); - if (!currents) { // Disabling electric fields // phi = 0.0; // Already set in initialisation @@ -1249,7 +1453,7 @@ int Hermes::rhs(BoutReal t) { // and only phi is considered, not phi + Pi which is handled in Boussinesq solves Field2D phi_boundary2d; Field3D phi_boundary3d; - + if (phi_boundary_relax) { // Update the boundary regions by relaxing towards zero gradient // on a given timescale. @@ -1330,7 +1534,7 @@ int Hermes::rhs(BoutReal t) { // Sheath multiplier Te -> phi (2.84522 for Deuterium if Ti = 0) phi_boundary2d = - DC((log(0.5 * sqrt(mi_me / PI)) + log(sqrt(Telim / (Telim + Tilim)))) * Telim); + DC((log(0.5 * sqrt(mi_me / PI)) + log(sqrt(Te / (Te + Ti)))) * Te); phi_boundary3d = phi_boundary2d; } @@ -1446,6 +1650,7 @@ int Hermes::rhs(BoutReal t) { } } phi.applyBoundary(t); +#warning no parallel phi boundaries needed? mesh->communicate(phi); } @@ -1480,6 +1685,8 @@ int Hermes::rhs(BoutReal t) { Ve.applyBoundary(t); mesh->communicate(Ve, psi); + // Ve.applyParallelBoundary(); + // psi.applyParallelBoundary(); Jpar = mul_all(Ne, sub_all(Vi, Ve)); Jpar.applyBoundary(t); @@ -1489,19 +1696,15 @@ int Hermes::rhs(BoutReal t) { // No Ve term in VePsi, only electromagnetic term psi = div_all(VePsi, 0.5 * mi_me * beta_e); - // Ve = (NVi - Delp2(psi)) / Nelim; - if(fci_transform){ - Field3D atmp = 1.0; - mesh->communicate(atmp,psi); - Jpar = FV::Div_a_Laplace_perp(atmp, psi); - } else { - Jpar = FV::Div_a_Laplace_perp(1.0, psi); - } + // Ve = (NVi - Delp2(psi)) / Ne; + Field3D one; + set_all(one, 1.0); + Jpar = FV::Div_a_Laplace_perp(one, psi); mesh->communicate(Jpar); Jpar.applyBoundary(t); - Ve = div_all(sub_all(NVi, Jpar), Nelim); + Ve = div_all(sub_all(NVi, Jpar), Ne); } // psi -= psi.DC(); // Remove toroidal average, only keep fluctuations @@ -1514,28 +1717,20 @@ int Hermes::rhs(BoutReal t) { } else { // Zero electron mass and electrostatic. // Special case where Ohm's law has no time-derivatives - mesh->communicate(phi,Pe); + // mesh->communicate(phi,Pe); - if(fci_transform){ - tau_e = (Cs0 / rho_s0) * tau_e0 * pow(Telim, 1.5) / Nelim; - nu = resistivity_multiply / (1.96 * tau_e * mi_me); + tau_e = (Cs0 / rho_s0) * tau_e0 * pow(Te, 1.5) / Ne; + nu = resistivity_multiply / (1.96 * tau_e * mi_me); - Ve = Vi + (Grad_parP(phi) - Grad_parP(Pe) / Ne) / nu; - }else{ - Ve = Vi + (Grad_parP(phi) - Grad_parP(Pe) / Ne) / nu; - } + Ve = Vi + (Grad_parP(phi) - Grad_parP(Pe) / Ne) / nu; if (thermal_force) { - if(fci_transform) {mesh->communicate(Te);} Ve -= 0.71 * Grad_parP(Te) / nu; } } Ve.applyBoundary(t); - // Communicate auxilliary variables - mesh->communicate(Ve); Field3D neve = mul_all(Ne,Ve); - mesh->communicate(neve); Jpar = sub_all(NVi, neve); } // Ve -= Jpar0 / Ne; // Equilibrium current density @@ -1549,207 +1744,26 @@ int Hermes::rhs(BoutReal t) { TRACE("Sheath boundaries"); - if (sheath_ydown) { - switch (sheath_model) { - case 0: { // Normal Bohm sheath - for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - // Zero-gradient density - BoutReal nesheath = floor(Ne(r.ind, mesh->ystart, jz), 0.0); - - // Temperature at the sheath entrance - BoutReal tesheath = floor(Te(r.ind, mesh->ystart, jz), 0.0); - BoutReal tisheath = floor(Ti(r.ind, mesh->ystart, jz), 0.0); - - // Zero-gradient potential - BoutReal phisheath = phi(r.ind, mesh->ystart, jz); - - // Ion velocity goes to the sound speed. Note negative since out of the domain - BoutReal visheath = -sqrt(tesheath + tisheath); - - if (sheath_allow_supersonic && (Vi(r.ind, mesh->ystart, jz) < visheath)) { - // If plasma is faster, go to plasma velocity - visheath = Vi(r.ind, mesh->ystart, jz); - } - - // Sheath current - // Note that phi/Te >= 0.0 since for phi < 0 - // vesheath is the electron saturation current - BoutReal phi_te = - floor(phisheath / Telim(r.ind, mesh->ystart, jz), 0.0); - BoutReal vesheath = - -sqrt(tesheath) * (sqrt(mi_me) / (2. * sqrt(PI))) * exp(-phi_te); - - // J = n*(Vi - Ve) - BoutReal jsheath = nesheath * (visheath - vesheath); - if (nesheath < 1e-10) { - vesheath = visheath; - jsheath = 0.0; - } - - // Apply boundary condition half-way between cells - // Neumann conditions - Ne.ydown()(r.ind, mesh->ystart - 1, jz) = nesheath; - phi.ydown()(r.ind, mesh->ystart - 1, jz) = phisheath; - Vort.ydown()(r.ind, mesh->ystart - 1, jz) = Vort(r.ind, mesh->ystart, jz); - - // Here zero-gradient Te, heat flux applied later - Te.ydown()(r.ind, mesh->ystart - 1, jz) = Te(r.ind, mesh->ystart, jz); - Ti.ydown()(r.ind, mesh->ystart - 1, jz) = Ti(r.ind, mesh->ystart, jz); - - Pe.ydown()(r.ind, mesh->ystart - 1, jz) = Pe(r.ind, mesh->ystart, jz); - Pelim.ydown()(r.ind, mesh->ystart - 1, jz) = Pelim(r.ind, mesh->ystart, jz); - Pi.ydown()(r.ind, mesh->ystart - 1, jz) = Pi(r.ind, mesh->ystart, jz); - Pilim.ydown()(r.ind, mesh->ystart - 1, jz) = Pilim(r.ind, mesh->ystart, jz); - - // Dirichlet conditions - Vi.ydown()(r.ind, mesh->ystart - 1, jz) = 2. * visheath - Vi(r.ind, mesh->ystart, jz); - Ve.ydown()(r.ind, mesh->ystart - 1, jz) = 2. * vesheath - Ve(r.ind, mesh->ystart, jz); - Jpar.ydown()(r.ind, mesh->ystart - 1, jz) = - 2. * jsheath - Jpar(r.ind, mesh->ystart, jz); - NVi.ydown()(r.ind, mesh->ystart - 1, jz) = - 2. * nesheath * visheath - NVi(r.ind, mesh->ystart, jz); - } - } - break; - } - case 2: { // Bohm sheath with free density - for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - // Zero-gradient density - BoutReal nesheath = 0.5 * (3. * Ne(r.ind, mesh->ystart, jz) - - Ne.yup()(r.ind, mesh->ystart + 1, jz)); - if (nesheath < 0.0) - nesheath = 0.0; - - // Temperature at the sheath entrance - BoutReal tesheath = floor(Te(r.ind, mesh->ystart, jz), 0.0); - BoutReal tisheath = floor(Ti(r.ind, mesh->ystart, jz), 0.0); - - // Zero-gradient potential - BoutReal phisheath = phi(r.ind, mesh->ystart, jz); - - // Ion velocity goes to the sound speed - BoutReal visheath = - -sqrt(tesheath + tisheath); // Sound speed outwards - - if (sheath_allow_supersonic && (Vi(r.ind, mesh->ystart, jz) < visheath)) { - // If plasma is faster, go to plasma velocity - visheath = Vi(r.ind, mesh->ystart, jz); - } - - // Sheath current - BoutReal phi_te = - floor(phisheath / Telim(r.ind, mesh->ystart, jz), 0.0); - BoutReal vesheath = - -sqrt(tesheath) * (sqrt(mi_me) / (2. * sqrt(PI))) * exp(-phi_te); - // J = n*(Vi - Ve) - BoutReal jsheath = nesheath * (visheath - vesheath); - if (nesheath < 1e-10) { - vesheath = visheath; - jsheath = 0.0; - } - - // Apply boundary condition half-way between cells - - // Neumann conditions - phi.ydown()(r.ind, mesh->ystart - 1, jz) = phisheath; - Vort.ydown()(r.ind, mesh->ystart - 1, jz) = Vort(r.ind, mesh->ystart, jz); - - // Here zero-gradient Te, heat flux applied later - Te.ydown()(r.ind, mesh->ystart - 1, jz) = Te(r.ind, mesh->ystart, jz); - Ti.ydown()(r.ind, mesh->ystart - 1, jz) = Ti(r.ind, mesh->ystart, jz); - - // Dirichlet conditions - Ne.ydown()(r.ind, mesh->ystart - 1, jz) = 2. * nesheath - Ne(r.ind, mesh->ystart, jz); - Pe.ydown()(r.ind, mesh->ystart - 1, jz) = - 2. * nesheath * tesheath - Pe(r.ind, mesh->ystart, jz); - Pi.ydown()(r.ind, mesh->ystart - 1, jz) = - 2. * nesheath * tisheath - Pi(r.ind, mesh->ystart, jz); - - Vi.ydown()(r.ind, mesh->ystart - 1, jz) = 2. * visheath - Vi(r.ind, mesh->ystart, jz); - Ve.ydown()(r.ind, mesh->ystart - 1, jz) = 2. * vesheath - Ve(r.ind, mesh->ystart, jz); - Jpar.ydown()(r.ind, mesh->ystart - 1, jz) = 2. * jsheath - Jpar(r.ind, mesh->ystart, jz); - NVi.ydown()(r.ind, mesh->ystart - 1, jz) = - 2. * nesheath * visheath - NVi(r.ind, mesh->ystart, jz); - } - } - break; - } - case 3: { // Insulating Bohm sheath with free density - for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - // Free density - BoutReal nesheath = 0.5 * (3. * Ne(r.ind, mesh->ystart, jz) - - Ne.yup()(r.ind, mesh->ystart + 1, jz)); - if (nesheath < 0.0) - nesheath = 0.0; - - // Temperature at the sheath entrance - BoutReal tesheath = floor(Te(r.ind, mesh->ystart, jz), 0.0); - BoutReal tisheath = floor(Ti(r.ind, mesh->ystart, jz), 0.0); - - // Zero-gradient potential - // NOTE: This should probably not be zero-gradient, - // since insulator could have electric field across it - BoutReal phisheath = phi(r.ind, mesh->ystart, jz); - - // Ion velocity goes to the sound speed - BoutReal visheath = - -sqrt(tesheath + tisheath); // Sound speed outwards - - if (sheath_allow_supersonic && (Vi(r.ind, mesh->ystart, jz) < visheath)) { - // If plasma is faster, go to plasma velocity - visheath = Vi(r.ind, mesh->ystart, jz); - } - - // Sheath current set to zero, as insulating boundary - BoutReal vesheath = visheath; - - // Apply boundary condition half-way between cells - Ne(r.ind, mesh->ystart - 1, jz) = 2. * nesheath - Ne(r.ind, mesh->ystart, jz); - phi(r.ind, mesh->ystart - 1, jz) = - 2. * phisheath - phi(r.ind, mesh->ystart, jz); - Vort(r.ind, mesh->ystart - 1, jz) = Vort(r.ind, mesh->ystart, jz); - - // Here zero-gradient Te, heat flux applied later - Te.ydown()(r.ind, mesh->ystart - 1, jz) = Te(r.ind, mesh->ystart, jz); - Ti.ydown()(r.ind, mesh->ystart - 1, jz) = Ti(r.ind, mesh->ystart, jz); - - Pe.ydown()(r.ind, mesh->ystart - 1, jz) = Pe(r.ind, mesh->ystart, jz); - Pi.ydown()(r.ind, mesh->ystart - 1, jz) = Pi(r.ind, mesh->ystart, jz); - - // Dirichlet conditions - Vi.ydown()(r.ind, mesh->ystart - 1, jz) = 2. * visheath - Vi(r.ind, mesh->ystart, jz); - Ve.ydown()(r.ind, mesh->ystart - 1, jz) = 2. * vesheath - Ve(r.ind, mesh->ystart, jz); - Jpar.ydown()(r.ind, mesh->ystart - 1, jz) = -Jpar(r.ind, mesh->ystart, jz); - NVi.ydown()(r.ind, mesh->ystart - 1, jz) = - 2. * nesheath * visheath - NVi(r.ind, mesh->ystart, jz); - } - } - break; - } - } - } else { + { // sheath_ydown == false for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { // Zero-gradient Te Te.ydown()(r.ind, mesh->ystart - 1, jz) = Te(r.ind, mesh->ystart, jz); - Telim.ydown()(r.ind, mesh->ystart - 1, jz) = Telim(r.ind, mesh->ystart, jz); + Te.ydown()(r.ind, mesh->ystart - 1, jz) = Te(r.ind, mesh->ystart, jz); Ti.ydown()(r.ind, mesh->ystart - 1, jz) = Ti(r.ind, mesh->ystart, jz); - Tilim.ydown()(r.ind, mesh->ystart - 1, jz) = Tilim(r.ind, mesh->ystart, jz); + Ti.ydown()(r.ind, mesh->ystart - 1, jz) = Ti(r.ind, mesh->ystart, jz); Pe.ydown()(r.ind, mesh->ystart - 1, jz) = Ne.ydown()(r.ind, mesh->ystart - 1, jz) * Te.ydown()(r.ind, mesh->ystart - 1, jz); - Pelim.ydown()(r.ind, mesh->ystart - 1, jz) = - Nelim.ydown()(r.ind, mesh->ystart - 1, jz) * Telim.ydown()(r.ind, mesh->ystart - 1, jz); + Pe.ydown()(r.ind, mesh->ystart - 1, jz) = + Ne.ydown()(r.ind, mesh->ystart - 1, jz) * Te.ydown()(r.ind, mesh->ystart - 1, jz); Pi.ydown()(r.ind, mesh->ystart - 1, jz) = Ne.ydown()(r.ind, mesh->ystart - 1, jz) * Ti.ydown()(r.ind, mesh->ystart - 1, jz); - Pilim.ydown()(r.ind, mesh->ystart - 1, jz) = - Nelim.ydown()(r.ind, mesh->ystart - 1, jz) * Tilim.ydown()(r.ind, mesh->ystart - 1, jz); + Pi.ydown()(r.ind, mesh->ystart - 1, jz) = + Ne.ydown()(r.ind, mesh->ystart - 1, jz) * Ti.ydown()(r.ind, mesh->ystart - 1, jz); // No flows Vi.ydown()(r.ind, mesh->ystart - 1, jz) = -Vi(r.ind, mesh->ystart, jz); @@ -1760,200 +1774,26 @@ int Hermes::rhs(BoutReal t) { } } - if (sheath_yup) { - switch (sheath_model) { - case 0: { // Normal Bohm sheath - for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - // Zero-gradient density - BoutReal nesheath = floor(Ne(r.ind, mesh->yend, jz), 0.0); - - // Temperature at the sheath entrance - BoutReal tesheath = floor(Te(r.ind, mesh->yend, jz), 0.0); - BoutReal tisheath = floor(Ti(r.ind, mesh->yend, jz), 0.0); - - // Zero-gradient potential - BoutReal phisheath = phi(r.ind, mesh->yend, jz); - - // Ion velocity goes to the sound speed - BoutReal visheath = sqrt(tesheath + tisheath); // Sound speed outwards - - if (sheath_allow_supersonic && (Vi(r.ind, mesh->yend, jz) > visheath)) { - // If plasma is faster, go to plasma velocity - visheath = Vi(r.ind, mesh->yend, jz); - } - - // Sheath current - // Note that phi/Te >= 0.0 since for phi < 0 - // vesheath is the electron saturation current - BoutReal phi_te = - floor(phisheath / Telim(r.ind, mesh->yend, jz), 0.0); - BoutReal vesheath = - sqrt(tesheath) * (sqrt(mi_me) / (2. * sqrt(PI))) * exp(-phi_te); - // J = n*(Vi - Ve) - BoutReal jsheath = nesheath * (visheath - vesheath); - if (nesheath < 1e-10) { - vesheath = visheath; - jsheath = 0.0; - } - - // Apply boundary condition half-way between cells - // Neumann conditions - Ne.yup()(r.ind, mesh->yend + 1, jz) = nesheath; - phi.yup()(r.ind, mesh->yend + 1, jz) = phisheath; - Vort.yup()(r.ind, mesh->yend + 1, jz) = Vort(r.ind, mesh->yend, jz); - - // Here zero-gradient Te, heat flux applied later - Te.yup()(r.ind, mesh->yend + 1, jz) = Te(r.ind, mesh->yend, jz); - Ti.yup()(r.ind, mesh->yend + 1, jz) = Ti(r.ind, mesh->yend, jz); - - Pe.yup()(r.ind, mesh->yend + 1, jz) = Pe(r.ind, mesh->yend, jz); - Pelim.yup()(r.ind, mesh->yend + 1, jz) = Pelim(r.ind, mesh->yend, jz); - Pi.yup()(r.ind, mesh->yend + 1, jz) = Pi(r.ind, mesh->yend, jz); - Pilim.yup()(r.ind, mesh->yend + 1, jz) = Pilim(r.ind, mesh->yend, jz); - - // Dirichlet conditions - Vi.yup()(r.ind, mesh->yend + 1, jz) = 2. * visheath - Vi(r.ind, mesh->yend, jz); - Ve.yup()(r.ind, mesh->yend + 1, jz) = 2. * vesheath - Ve(r.ind, mesh->yend, jz); - Jpar.yup()(r.ind, mesh->yend + 1, jz) = 2. * jsheath - Jpar(r.ind, mesh->yend, jz); - NVi.yup()(r.ind, mesh->yend + 1, jz) = - 2. * nesheath * visheath - NVi(r.ind, mesh->yend, jz); - } - } - break; - } - case 2: { // Bohm sheath with free density - for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - // Zero-gradient density - BoutReal nesheath = 0.5 * (3. * Ne(r.ind, mesh->yend, jz) - - Ne.ydown()(r.ind, mesh->yend - 1, jz)); - if (nesheath < 0.0) - nesheath = 0.0; - - // Temperature at the sheath entrance - BoutReal tesheath = floor(Te(r.ind, mesh->yend, jz), 0.0); - BoutReal tisheath = floor(Ti(r.ind, mesh->yend, jz), 0.0); - - // Zero-gradient potential - BoutReal phisheath = phi(r.ind, mesh->yend, jz); - - // Ion velocity goes to the sound speed - BoutReal visheath = sqrt(tesheath + tisheath); // Sound speed outwards - - if (sheath_allow_supersonic && (Vi(r.ind, mesh->yend, jz) > visheath)) { - // If plasma is faster, go to plasma velocity - visheath = Vi(r.ind, mesh->yend, jz); - } - - // Sheath current - BoutReal phi_te = - floor(phisheath / Telim(r.ind, mesh->yend, jz), 0.0); - BoutReal vesheath = - sqrt(tesheath) * (sqrt(mi_me) / (2. * sqrt(PI))) * exp(-phi_te); - // J = n*(Vi - Ve) - BoutReal jsheath = nesheath * (visheath - vesheath); - if (nesheath < 1e-10) { - vesheath = visheath; - jsheath = 0.0; - } - - // Apply boundary condition half-way between cells - // Neumann conditions - phi.yup()(r.ind, mesh->yend + 1, jz) = phisheath; - Vort.yup()(r.ind, mesh->yend + 1, jz) = Vort(r.ind, mesh->yend, jz); - - // Here zero-gradient Te, heat flux applied later - Te.yup()(r.ind, mesh->yend + 1, jz) = tesheath; - Ti.yup()(r.ind, mesh->yend + 1, jz) = Ti(r.ind, mesh->yend, jz); - - // Dirichlet conditions - Ne.yup()(r.ind, mesh->yend + 1, jz) = 2. * nesheath - Ne(r.ind, mesh->yend, jz); - Pe.yup()(r.ind, mesh->yend + 1, jz) = - 2. * nesheath * tesheath - Pe(r.ind, mesh->yend, jz); - Pi.yup()(r.ind, mesh->yend + 1, jz) = - 2. * nesheath * tisheath - Pi(r.ind, mesh->yend, jz); - - Vi.yup()(r.ind, mesh->yend + 1, jz) = 2. * visheath - Vi(r.ind, mesh->yend, jz); - Ve.yup()(r.ind, mesh->yend + 1, jz) = 2. * vesheath - Ve(r.ind, mesh->yend, jz); - Jpar.yup()(r.ind, mesh->yend + 1, jz) = 2. * jsheath - Jpar(r.ind, mesh->yend, jz); - NVi.yup()(r.ind, mesh->yend + 1, jz) = - 2. * nesheath * visheath - NVi(r.ind, mesh->yend, jz); - } - } - break; - } - case 3: { // Insulating Bohm sheath with free density - for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - // Zero-gradient density - BoutReal nesheath = 0.5 * (3. * Ne(r.ind, mesh->yend, jz) - - Ne.ydown()(r.ind, mesh->yend - 1, jz)); - if (nesheath < 0.0) { - nesheath = 0.0; - } - - // Temperature at the sheath entrance - BoutReal tesheath = floor(Te(r.ind, mesh->yend, jz), 0.0); - BoutReal tisheath = floor(Ti(r.ind, mesh->yend, jz), 0.0); - - // Zero-gradient potential - BoutReal phisheath = phi(r.ind, mesh->yend, jz); - - // Ion velocity goes to the sound speed - BoutReal visheath = sqrt(tesheath + tisheath); // Sound speed outwards - - if (sheath_allow_supersonic && (Vi(r.ind, mesh->yend, jz) > visheath)) { - // If plasma is faster, go to plasma velocity - visheath = Vi(r.ind, mesh->yend, jz); - } - - // Zero sheath current - BoutReal vesheath = visheath; - - // Apply boundary condition half-way between cells - Ne.yup()(r.ind, mesh->yend + 1, jz) = 2. * nesheath - Ne(r.ind, mesh->yend, jz); - phi.yup()(r.ind, mesh->yend + 1, jz) = 2. * phisheath - phi(r.ind, mesh->yend, jz); - Vort.yup()(r.ind, mesh->yend + 1, jz) = Vort(r.ind, mesh->yend, jz); - - // Here zero-gradient Te, heat flux applied later - Te.yup()(r.ind, mesh->yend + 1, jz) = Te(r.ind, mesh->yend, jz); - Ti.yup()(r.ind, mesh->yend + 1, jz) = Ti(r.ind, mesh->yend, jz); - - Pe.yup()(r.ind, mesh->yend + 1, jz) = Pe(r.ind, mesh->yend, jz); - Pi.yup()(r.ind, mesh->yend + 1, jz) = Pi(r.ind, mesh->yend, jz); - - // Dirichlet conditions - Vi.yup()(r.ind, mesh->yend + 1, jz) = 2. * visheath - Vi(r.ind, mesh->yend, jz); - Ve.yup()(r.ind, mesh->yend + 1, jz) = 2. * vesheath - Ve(r.ind, mesh->yend, jz); - Jpar.yup()(r.ind, mesh->yend + 1, jz) = -Jpar(r.ind, mesh->yend, jz); - NVi.yup()(r.ind, mesh->yend + 1, jz) = - 2. * nesheath * visheath - NVi(r.ind, mesh->yend, jz); - } - } - break; - } - } - } else { + { // sheath_yup == false for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { for (int jz = 0; jz < mesh->LocalNz; jz++) { // Zero-gradient Te Te.yup()(r.ind, mesh->yend + 1, jz) = Te(r.ind, mesh->yend, jz); - Telim.yup()(r.ind, mesh->yend + 1, jz) = Telim(r.ind, mesh->yend, jz); + Te.yup()(r.ind, mesh->yend + 1, jz) = Te(r.ind, mesh->yend, jz); Ti.yup()(r.ind, mesh->yend + 1, jz) = Ti(r.ind, mesh->yend, jz); - Tilim.yup()(r.ind, mesh->yend + 1, jz) = Tilim(r.ind, mesh->yend, jz); + Ti.yup()(r.ind, mesh->yend + 1, jz) = Ti(r.ind, mesh->yend, jz); Pe.yup()(r.ind, mesh->yend + 1, jz) = Ne.yup()(r.ind, mesh->yend + 1, jz) * Te.yup()(r.ind, mesh->yend + 1, jz); - Pelim.yup()(r.ind, mesh->yend + 1, jz) = - Nelim.yup()(r.ind, mesh->yend + 1, jz) * Telim.yup()(r.ind, mesh->yend + 1, jz); + Pe.yup()(r.ind, mesh->yend + 1, jz) = + Ne.yup()(r.ind, mesh->yend + 1, jz) * Te.yup()(r.ind, mesh->yend + 1, jz); Pi.yup()(r.ind, mesh->yend + 1, jz) = Ne.yup()(r.ind, mesh->yend + 1, jz) * Ti.yup()(r.ind, mesh->yend + 1, jz); - Pilim.yup()(r.ind, mesh->yend + 1, jz) = - Nelim.yup()(r.ind, mesh->yend + 1, jz) * Tilim.yup()(r.ind, mesh->yend + 1, jz); + Pi.yup()(r.ind, mesh->yend + 1, jz) = + Ne.yup()(r.ind, mesh->yend + 1, jz) * Ti.yup()(r.ind, mesh->yend + 1, jz); // No flows Vi.yup()(r.ind, mesh->yend + 1, jz) = -Vi(r.ind, mesh->yend, jz); @@ -1967,229 +1807,106 @@ int Hermes::rhs(BoutReal t) { if (parallel_sheaths){ switch (par_sheath_model) { case 0 :{ - for (const auto &bndry_par : mesh->getBoundariesPar()) { + for (const auto &bndry_par : mesh->getBoundariesPar(BoundaryParType::xout)) { for (bndry_par->first(); !bndry_par->isDone(); bndry_par->next()) { int x = bndry_par->x; int y = bndry_par->y; int z = bndry_par->z; // output.write("x: {},y: {},z: {}\n", x,y,z); - if (bndry_par->dir == 1) { - // Zero-gradient density - BoutReal nesheath = floor(Ne(x, y, z), 0.0); - - // Temperature at the sheath entrance - BoutReal tesheath = floor(Te(x, y, z), 0.0); - BoutReal tisheath = floor(Ti(x, y, z), 0.0); - - // Zero-gradient potential - BoutReal phisheath = phi(x, y, z); - BoutReal visheath = sqrt(tisheath + tesheath); + // Zero-gradient density + BoutReal nesheath = floor(Ne(x, y, z), 0.0); - if (sheath_allow_supersonic && (Vi(x, y, z) > visheath)) { - // If plasma is faster, go to plasma velocity - visheath = Vi(x, y, z); - } + // Temperature at the sheath entrance + BoutReal tesheath = floor(Te(x, y, z), 0.0); + BoutReal tisheath = floor(Ti(x, y, z), 0.0); - // Sheath current - // Note that phi/Te >= 0.0 since for phi < 0 - // vesheath is the electron saturation current - BoutReal phi_te = - floor(phisheath / tesheath, 0.0); - - BoutReal vesheath = - sqrt(tesheath) * (sqrt(mi_me) / (2. * sqrt(PI))) * exp(-phi_te); - - // J = n*(Vi - Ve) - BoutReal jsheath = nesheath * (visheath - vesheath); - if (nesheath < 1e-10) { - vesheath = visheath; - jsheath = 0.0; - } - - // Neumann conditions - Ne.yup()(x, y+bndry_par->dir, z) = nesheath; - phi.yup()(x, y+bndry_par->dir, z) = phisheath; - Vort.yup()(x, y+bndry_par->dir, z) = Vort(x, y, z); - - // Here zero-gradient Te, heat flux applied later - Te.yup()(x, y+bndry_par->dir, z) = Te(x, y, z); - Ti.yup()(x, y+bndry_par->dir, z) = Ti(x, y, z); - - Pe.yup()(x, y+bndry_par->dir, z) = Pe(x, y, z); - Pelim.yup()(x, y+bndry_par->dir, z) = Pelim(x, y, z); - Pi.yup()(x, y+bndry_par->dir, z) = Pi(x, y, z); - Pilim.yup()(x, y+bndry_par->dir, z) = Pilim(x, y, z); - - // // Dirichlet conditions - Vi.yup()(x, y+bndry_par->dir, z) = 2. * visheath - Vi(x, y, z); - if (par_sheath_ve){ - Ve.yup()(x, y+bndry_par->dir, z) = 2. * vesheath - Ve(x, y, z); - } - Jpar.yup()(x, y+bndry_par->dir, z) = - 2. * jsheath - Jpar(x, y, z); - NVi.yup()(x, y+bndry_par->dir, z) = - 2. * nesheath * visheath - NVi(x, y, z); - }else{ //backwards - // Zero-gradient density - BoutReal nesheath = floor(Ne(x, y, z), 0.0); - - // Temperature at the sheath entrance - BoutReal tesheath = floor(Te(x, y, z), 0.0); - BoutReal tisheath = floor(Ti(x, y, z), 0.0); - - // Zero-gradient potential - BoutReal phisheath = phi(x, y, z); - BoutReal visheath = -sqrt(tisheath + tesheath); + // Zero-gradient potential + BoutReal phisheath = phi(x, y, z); + BoutReal visheath = bndry_par->dir * sqrt(tisheath + tesheath); - if (sheath_allow_supersonic && (Vi(x, y, z) < visheath)) { - // If plasma is faster, go to plasma velocity - visheath = Vi(x, y, z); + if (sheath_allow_supersonic){ + if (bndry_par->dir == 1){ + if (Vi(x, y, z) > visheath){ + // If plasma is faster, go to plasma velocity + visheath = Vi(x, y, z); + } + } else { + if (Vi(x, y, z) < visheath){ + visheath = Vi(x, y, z); + } } + } - // Sheath current - // Note that phi/Te >= 0.0 since for phi < 0 - // vesheath is the electron saturation current - BoutReal phi_te = - floor(phisheath / tesheath, 0.0); + // Sheath current + // Note that phi/Te >= 0.0 since for phi < 0 + // vesheath is the electron saturation current + BoutReal phi_te = + floor(phisheath / tesheath, 0.0); - BoutReal vesheath = - -sqrt(tesheath) * (sqrt(mi_me) / (2. * sqrt(PI))) * exp(-phi_te); - - // J = n*(Vi - Ve) - BoutReal jsheath = nesheath * (visheath - vesheath); - if (nesheath < 1e-8) { - vesheath = visheath; - jsheath = 0.0; - } + BoutReal vesheath = + bndry_par->dir * sqrt(tesheath) * (sqrt(mi_me) / (2. * sqrt(PI))) * exp(-phi_te); - // Neumann conditions - Ne.ydown()(x, y+bndry_par->dir, z) = nesheath; - phi.ydown()(x, y+bndry_par->dir, z) = phisheath; - Vort.ydown()(x, y+bndry_par->dir, z) = Vort(x, y, z); - - // Here zero-gradient Te, heat flux applied later - Te.ydown()(x, y+bndry_par->dir, z) = Te(x, y, z); - Ti.ydown()(x, y+bndry_par->dir, z) = Ti(x, y, z); - - Pe.ydown()(x, y+bndry_par->dir, z) = Pe(x, y, z); - Pelim.ydown()(x, y+bndry_par->dir, z) = Pelim(x, y, z); - Pi.ydown()(x, y+bndry_par->dir, z) = Pi(x, y, z); - Pilim.ydown()(x, y+bndry_par->dir, z) = Pilim(x, y, z); - - // // Dirichlet conditions - Vi.ydown()(x, y+bndry_par->dir, z) = 2. * visheath - Vi(x, y, z); - if (par_sheath_ve){ - Ve.ydown()(x, y+bndry_par->dir, z) = 2. * vesheath - Ve(x, y, z); - } - Jpar.ydown()(x, y+bndry_par->dir, z) = - 2. * jsheath - Jpar(x, y, z); - NVi.ydown()(x, y+bndry_par->dir, z) = - 2. * nesheath * visheath - NVi(x, y, z); + // J = n*(Vi - Ve) + BoutReal jsheath = nesheath * (visheath - vesheath); + if (nesheath < 1e-10) { + vesheath = visheath; + jsheath = 0.0; } - // // Zero-gradient density - // BoutReal nesheath = floor(Ne(x, y, z), 0.0); + + // Neumann conditions + Ne.ynext(bndry_par->dir)(x, y+bndry_par->dir, z) = nesheath; + phi.ynext(bndry_par->dir)(x, y+bndry_par->dir, z) = phisheath; + Vort.ynext(bndry_par->dir)(x, y+bndry_par->dir, z) = Vort(x, y, z); - // // Temperature at the sheath entrance - // BoutReal tesheath = floor(Te(x, y, z), 0.0); - // BoutReal tisheath = floor(Ti(x, y, z), 0.0); - - // // Zero-gradient potential - // BoutReal phisheath = phi(x, y, z); - - // // Ion velocity goes to the sound speed. Note negative since out of the domain - // BoutReal visheath = -sqrt(tesheath + tisheath); - - // if (sheath_allow_supersonic && (Vi(x, y, z) < visheath)) { - // // If plasma is faster, go to plasma velocity - // visheath = Vi(x, y, z); - // } - - // if (bndry_par->dir == 1){ - // visheath = sqrt(tesheath + tisheath); - - // if (sheath_allow_supersonic && (Vi(x, y, z) > visheath)) { - // // If plasma is faster, go to plasma velocity - // visheath = Vi(x, y, z); - // } - - // } - - - // // Sheath current - // // Note that phi/Te >= 0.0 since for phi < 0 - // // vesheath is the electron saturation current - // BoutReal phi_te = - // floor(phisheath / Telim(x, y, z), 0.0); - - // BoutReal vesheath = - // -sqrt(tesheath) * (sqrt(mi_me) / (2. * sqrt(PI))) * exp(-phi_te); - - // if (bndry_par->dir == 1){ - // vesheath *= -1.; - // } - - // // J = n*(Vi - Ve) - // BoutReal jsheath = nesheath * (visheath - vesheath); - // if (nesheath < 1e-10) { - // vesheath = visheath; - // jsheath = 0.0; - // } - - // if (bndry_par->dir == 1){ - // // Apply boundary condition half-way between cells - // // Neumann conditions - // Ne.yup()(x, y+bndry_par->dir, z) = nesheath; - // phi.yup()(x, y+bndry_par->dir, z) = phisheath; - // Vort.yup()(x, y+bndry_par->dir, z) = Vort(x, y, z); - - // // Here zero-gradient Te, heat flux applied later - // // Te.yup()(x, y+bndry_par->dir, z) = Te(x, y, z); - // // Ti.yup()(x, y+bndry_par->dir, z) = Ti(x, y, z); - - // // Pe.yup()(x, y+bndry_par->dir, z) = Pe(x, y, z); - // // Pelim.yup()(x, y+bndry_par->dir, z) = Pelim(x, y, z); - // // Pi.yup()(x, y+bndry_par->dir, z) = Pi(x, y, z); - // // Pilim.yup()(x, y+bndry_par->dir, z) = Pilim(x, y, z); - - // // // Dirichlet conditions - // // Vi.yup()(x, y+bndry_par->dir, z) = 2. * visheath - Vi(x, y, z); - // if (par_sheath_ve){ - // Ve.yup()(x, y+bndry_par->dir, z) = 2. * vesheath - Ve(x, y, z); - // } - // Jpar.yup()(x, y+bndry_par->dir, z) = - // 2. * jsheath - Jpar(x, y, z); - // NVi.yup()(x, y+bndry_par->dir, z) = - // 2. * nesheath * visheath - NVi(x, y, z); - // } else if (bndry_par->dir == -1) { - // // Apply boundary condition half-way between cells - // // Neumann conditions - // Ne.ydown()(x, y+bndry_par->dir, z) = nesheath; - // phi.ydown()(x, y+bndry_par->dir, z) = phisheath; - // Vort.ydown()(x, y+bndry_par->dir, z) = Vort(x, y, z); - - // // Here zero-gradient Te, heat flux applied later - // // Te.ydown()(x, y+bndry_par->dir, z) = Te(x, y, z); - // // Ti.ydown()(x, y+bndry_par->dir, z) = Ti(x, y, z); - - // // Pe.ydown()(x, y+bndry_par->dir, z) = Pe(x, y, z); - // // Pelim.ydown()(x, y+bndry_par->dir, z) = Pelim(x, y, z); - // // Pi.ydown()(x, y+bndry_par->dir, z) = Pi(x, y, z); - // // Pilim.ydown()(x, y+bndry_par->dir, z) = Pilim(x, y, z); + // Here zero-gradient Te, heat flux applied later + Te.ynext(bndry_par->dir)(x, y+bndry_par->dir, z) = Te(x, y, z); + Ti.ynext(bndry_par->dir)(x, y+bndry_par->dir, z) = Ti(x, y, z); + + Pe.ynext(bndry_par->dir)(x, y+bndry_par->dir, z) = Pe(x, y, z); + Pi.ynext(bndry_par->dir)(x, y+bndry_par->dir, z) = Pi(x, y, z); - // // Dirichlet conditions - // // Vi.ydown()(x, y+bndry_par->dir, z) = 2. * visheath - Vi(x, y, z); - // if (par_sheath_ve) { - // Ve.ydown()(x, y+bndry_par->dir, z) = 2. * vesheath - Ve(x, y, z); - // } - // Jpar.ydown()(x, y+bndry_par->dir, z) = - // 2. * jsheath - Jpar(x, y, z); - // NVi.ydown()(x, y+bndry_par->dir, z) = - // 2. * nesheath * visheath - NVi(x, y, z); - // } + // // Dirichlet conditions + Vi.ynext(bndry_par->dir)(x, y+bndry_par->dir, z) = 2. * visheath - Vi(x, y, z); + if (par_sheath_ve){ + Ve.ynext(bndry_par->dir)(x, y+bndry_par->dir, z) = 2. * vesheath - Ve(x, y, z); + } + Jpar.ynext(bndry_par->dir)(x, y+bndry_par->dir, z) = + 2. * jsheath - Jpar(x, y, z); + NVi.ynext(bndry_par->dir)(x, y+bndry_par->dir, z) = + 2. * nesheath * visheath - NVi(x, y, z); } } + + const int ny = mesh->LocalNy; + const int nz = mesh->LocalNz; + for (const auto &bndry_par : + mesh->getBoundariesPar(BoundaryParType::xin)) { + for (bndry_par->first(); !bndry_par->isDone(); bndry_par->next()) { + const Ind3D i{(bndry_par->x * ny + bndry_par->y) * nz + bndry_par->z, + ny, nz}; + const auto inext = i.yp(bndry_par->dir); + + // Neumann conditions + Ne.ynext(bndry_par->dir)[inext] = Ne[i]; + phi.ynext(bndry_par->dir)[inext] = phi[i]; + Vort.ynext(bndry_par->dir)[inext] = Vort[i]; + + // Here zero-gradient Te, heat flux applied later + Te.ynext(bndry_par->dir)[inext] = Te[i]; + Ti.ynext(bndry_par->dir)[inext] = Ti[i]; + + Pe.ynext(bndry_par->dir)[inext] = Pe[i]; + Pi.ynext(bndry_par->dir)[inext] = Pi[i]; + + // // Dirichlet conditions + Vi.ynext(bndry_par->dir)[inext] = -Vi[i]; + Ve.ynext(bndry_par->dir)[inext] = -Ve[i]; + Jpar.ynext(bndry_par->dir)[inext] = -Jpar[i]; + NVi.ynext(bndry_par->dir)[inext] = -NVi[i]; + } + } break; } case 1:{ //insulating boundary - for (const auto &bndry_par : mesh->getBoundariesPar()) { + for (const auto &bndry_par : mesh->getBoundariesPar(BoundaryParType::xout)) { for (bndry_par->first(); !bndry_par->isDone(); bndry_par->next()) { int x = bndry_par->x; int y = bndry_par->y; int z = bndry_par->z; // Zero-gradient density @@ -2225,7 +1942,7 @@ int Hermes::rhs(BoutReal t) { // Note that phi/Te >= 0.0 since for phi < 0 // vesheath is the electron saturation current BoutReal phi_te = - floor(phisheath / Telim(x, y, z), 0.0); + floor(phisheath / Te(x, y, z), 0.0); BoutReal vesheath = visheath; // J = n*(Vi - Ve) @@ -2243,9 +1960,9 @@ int Hermes::rhs(BoutReal t) { Ti.yup()(x, y+bndry_par->dir, z) = Ti(x, y, z); Pe.yup()(x, y+bndry_par->dir, z) = Pe(x, y, z); - Pelim.yup()(x, y+bndry_par->dir, z) = Pelim(x, y, z); + Pe.yup()(x, y+bndry_par->dir, z) = Pe(x, y, z); + Pi.yup()(x, y+bndry_par->dir, z) = Pi(x, y, z); Pi.yup()(x, y+bndry_par->dir, z) = Pi(x, y, z); - Pilim.yup()(x, y+bndry_par->dir, z) = Pilim(x, y, z); // Dirichlet conditions Vi.yup()(x, y+bndry_par->dir, z) = 2. * visheath - Vi(x, y, z); @@ -2266,9 +1983,9 @@ int Hermes::rhs(BoutReal t) { Ti.ydown()(x, y+bndry_par->dir, z) = Ti(x, y, z); Pe.ydown()(x, y+bndry_par->dir, z) = Pe(x, y, z); - Pelim.ydown()(x, y+bndry_par->dir, z) = Pelim(x, y, z); + Pe.ydown()(x, y+bndry_par->dir, z) = Pe(x, y, z); + Pi.ydown()(x, y+bndry_par->dir, z) = Pi(x, y, z); Pi.ydown()(x, y+bndry_par->dir, z) = Pi(x, y, z); - Pilim.ydown()(x, y+bndry_par->dir, z) = Pilim(x, y, z); // Dirichlet conditions Vi.ydown()(x, y+bndry_par->dir, z) = 2. * visheath - Vi(x, y, z); @@ -2283,9 +2000,10 @@ int Hermes::rhs(BoutReal t) { break; } } + } else { + throw BoutException("parallel BCs are needed!"); } - - + if (!currents) { // No currents, so reset Ve to be equal to Vi // VePsi also reset, so saved in restart file correctly @@ -2293,13 +2011,6 @@ int Hermes::rhs(BoutReal t) { VePsi = Ve; } - // Ensure that Nelim, Telim and Pelim are calculated in guard cells - Nelim = floor(Ne, 1e-5); - Telim = floor(Te, 0.001 / Tnorm); - Pelim = Telim * Nelim; - Tilim = floor(Ti, 0.001 / Tnorm); - Pilim = Tilim * Nelim; - ////////////////////////////////////////////////////////////// // Plasma quantities calculated. // At this point we have calculated all boundary conditions, @@ -2319,21 +2030,47 @@ int Hermes::rhs(BoutReal t) { ////////////////////////////////////////////////////////////// // Collisions and stress tensor TRACE("Collisions"); - - // Normalised electron collision time - tau_e = (Cs0 / rho_s0) * tau_e0 * pow(Telim, 1.5) / Nelim; - - // Normalised ion-ion collision time - tau_i = (Cs0 / rho_s0) * tau_i0 * pow(Tilim, 1.5) / Nelim; + const BoutReal tau_e1 = (Cs0 / rho_s0) * tau_e0; + const BoutReal tau_i1 = (Cs0 / rho_s0) * tau_i0; + Field3D neutral_rate; if (ion_neutral && ( neutrals || (ion_neutral_rate > 0.0))) { // Include ion-neutral collisions in collision time - - Field3D neutral_rate = neutrals ? neutrals->Fperp : ion_neutral_rate; - - // Add collision frequencies (1/tau_i + neutral rate) - tau_i /= 1.0 + tau_i * neutral_rate; - } + neutral_rate = neutrals ? neutrals->Fperp : ion_neutral_rate; + } + alloc_all(tau_e); + alloc_all(tau_i); + BOUT_FOR(i, Te.getRegion("RGN_ALL")) { + // Normalised electron collision time + // tau_e[i] = mul_all((Cs0 / rho_s0) * tau_e0, div_all(mul_all(Te, + // sqrt_all(Te)), Ne)); + tau_e[i] = tau_e1 * (Te[i] * sqrt(Te[i]) / Ne[i]); + tau_e.yup()[i] = tau_e1 * (Te.yup()[i] * sqrt(Te.yup()[i]) / Ne.yup()[i]); + tau_e.ydown()[i] = + tau_e1 * (Te.ydown()[i] * sqrt(Te.ydown()[i]) / Ne.ydown()[i]); + + // Normalised ion-ion collision time + tau_i[i] = tau_i1 * (Ti[i] * sqrt(Ti[i])) / Ne[i]; + tau_i.yup()[i] = tau_i1 * (Ti.yup()[i] * sqrt(Ti.yup()[i])) / Ne.yup()[i]; + tau_i.ydown()[i] = + tau_i1 * (Ti.ydown()[i] * sqrt(Ti.ydown()[i])) / Ne.ydown()[i]; + + if (ion_neutral && (neutrals || (ion_neutral_rate > 0.0))) { + // Include ion-neutral collisions in collision time + // Add collision frequencies (1/tau_i + neutral rate) + tau_i[i] = tau_i[i] / (1 + (tau_i[i] * neutral_rate[i])); + tau_i.yup()[i] = + tau_i.yup()[i] / (1 + (tau_i.yup()[i] * neutral_rate.yup()[i])); + tau_i.ydown()[i] = + tau_i.ydown()[i] / (1 + (tau_i.ydown()[i] * neutral_rate.ydown()[i])); + } + } + // tau_e = mul_all((Cs0 / rho_s0) * tau_e0, div_all(mul_all(Te, sqrt_all(Te)), + // Ne)); tau_i = mul_all((Cs0 / rho_s0) * tau_i0, div_all(mul_all(Ti, + // sqrt_all(Ti)), Ne)); if (ion_neutral && (neutrals || (ion_neutral_rate > + // 0.0))) { + // tau_i = div_all(tau_i, add_all(1, mul_all(tau_i, neutral_rate))); + // } // Collisional damping (normalised) if (resistivity || (!electromagnetic && !FiniteElMass)) { @@ -2351,7 +2088,7 @@ int Hermes::rhs(BoutReal t) { BoutReal a0 = PI*SQ(5.29e-11); // Cross-section [m^2] // Electron thermal speed - Field3D vth_e = sqrt(mi_me*Telim); + Field3D vth_e = sqrt(mi_me*Te); // Electron-neutral collision rate Field3D nu_ne = vth_e * Nnorm * neutrals->getDensity() * a0 * rho_s0; @@ -2408,7 +2145,7 @@ int Hermes::rhs(BoutReal t) { if (thermal_conduction || sinks) { // Braginskii expression for electron parallel conduction // kappa ~ n * v_th^2 * tau - kappa_epar = 3.16 * mi_me * Telim * Nelim * tau_e; + kappa_epar = mul_all(mul_all(mul_all(mul_all(3.16, mi_me), Te), Ne), tau_e); if (kappa_limit_alpha > 0.0) { /* @@ -2421,174 +2158,30 @@ int Hermes::rhs(BoutReal t) { * R.Schneider et al. Contrib. Plasma Phys. 46, No. 1-2, 3 – 191 (2006) * DOI 10.1002/ctpp.200610001 */ - + ASSERT0(not fci_transform); // Spitzer-Harm heat flux Field3D q_SH = kappa_epar * Grad_par(Te); - Field3D q_fl = kappa_limit_alpha * Nelim * Telim * sqrt(mi_me * Telim); + Field3D q_fl = kappa_limit_alpha * Ne * Te * sqrt(mi_me * Te); kappa_epar = kappa_epar / (1. + abs(q_SH / q_fl)); // Values of kappa on cell boundaries are needed for fluxes - mesh->communicate(kappa_epar); kappa_epar.applyBoundary("dirichlet"); } // Ion parallel heat conduction - kappa_ipar = 3.9 * Tilim * Nelim * tau_i; + kappa_ipar = mul_all(mul_all(mul_all(3.9, Ti), Ne), tau_i); - // Boundary conditions on heat conduction coefficients - for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - kappa_epar(r.ind, mesh->ystart - 1, jz) = kappa_epar(r.ind, mesh->ystart, jz); - kappa_ipar(r.ind, mesh->ystart - 1, jz) = kappa_ipar(r.ind, mesh->ystart, jz); - } - } - - for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - kappa_epar(r.ind, mesh->yend + 1, jz) = kappa_epar(r.ind, mesh->yend, jz); - kappa_ipar(r.ind, mesh->yend + 1, jz) = kappa_ipar(r.ind, mesh->yend, jz); - } - } + apply_flux_bcs(kappa_epar, kappa_ipar); } - if(currents){ nu.applyBoundary(t); } + if (currents) { + nu.applyBoundary(t); + nu.applyParallelBoundary("parallel_neumann"); + } if (ion_viscosity) { - /////////////////////////////////////////////////////////// - // Ion stress tensor. Split into - // Pi_ci = Pi_ciperp + Pi_cipar - // - // In the parallel ion momentum equation the Pi_cipar term - // is solved as a parallel diffusion, so is treated separately - // All other terms are added to Pi_ciperp, even if they are - // not really parallel parts - Field3D Tifree = copy(Ti); - Tifree.applyBoundary("free_o2"); - if(fci_transform){mesh->communicate(Tifree,Pilim,Tilim,coord->Bxy);} - - // For nonlinear terms, need to evaluate qipar and qi squared - Field3D qipar = -kappa_ipar * Grad_par(Tifree); - - // Limit the maximum value of tau_i - tau_i = ceil(tau_i, 1e4); - - // Square of total heat flux, parallel and perpendicular - // The first Pi term cancels the parallel part of the second term - // Doesn't include perpendicular collisional transport - - Field3D qisq = - (SQ(kappa_ipar) - SQ((5. / 2) * Pilim)) * SQ(Grad_par(Tifree)) + - SQ((5. / 2) * Pilim) * SQ(Grad(Tifree)); // This term includes a - // parallel component which is - // cancelled in first term - - Field3D phi161Ti = phi + 1.61 * Ti; - mesh->communicate(phi161Ti,Pi); - // Perpendicular part from curvature - Pi_ciperp = - -0.5 * 0.96 * Pi * tau_i * - (fci_curvature(phi161Ti) - - fci_curvature(Pi) / Nelim); // q perpendicular - // q parallel - - if(fci_transform){ - Coordinates::FieldMetric B32kappa_ipar = B32 * kappa_ipar; - mesh->communicate(B32kappa_ipar, Tifree); - Pi_ciperp += - 0.96 * tau_i * (1.42 / B32) * - Div_par_K_Grad_par(B32kappa_ipar, Tifree); - }else{ - Pi_ciperp += - 0.96 * tau_i * (1.42 / B32) * - FV::Div_par_K_Grad_par(B32 * kappa_ipar, Tifree); - } - - Field3D logTilim = log(Tilim); - Field3D logPilim = log(Pilim); - mesh->communicate(logTilim, logPilim); - Pi_ciperp -= - 0.49 * (qipar / Pilim) * - (2.27 * Grad_par(logTilim) - Grad_par(logPilim)) + - 0.75 * (0.2 * SQ(qipar) - 0.085 * qisq) / (Pilim * Tilim); - - Field3D logB = log(coord->Bxy); - mesh->communicate(Vi, logB); - // Parallel part - Pi_cipar = -0.96 * Pi * tau_i * - (2. * Grad_par(Vi) + Vi * Grad_par(logB)); - - // Could also be written as: - // Pi_cipar = - - // 0.96*Pi*tau_i*2.*Grad_par(sqrt(coord->Bxy)*Vi)/sqrt(coord->Bxy); - - if (mesh->firstX()) { - // First cells in X subject to boundary effects. - for(int i=mesh->xstart; (i <= mesh->xend) && (i < 4); i++) { - for (int j = mesh->ystart; j <= mesh->yend; j++) { - for (int k = 0; k < mesh->LocalNz; k++) { - Pi_ciperp(i, j, k) *= float(i - mesh->xstart)/4.0; - Pi_cipar(i, j, k) *= float(i - mesh->xstart)/4.0; - } - } - } - } - - // Limit size of stress tensor components - // If the off-diagonal components of the pressure tensor are large compared - // to the scalar pressure, then the model is probably breaking down. - // This can happen in very low density regions - BOUT_FOR(i, Pi_ciperp.getRegion("RGN_ALL")) { - if (fabs(Pi_ciperp[i]) > Pi[i]) { - Pi_ciperp[i] = SIGN(Pi_ciperp[i]) * Pi[i]; - } - if (fabs(Pi_cipar[i]) > Pi[i]) { - Pi_cipar[i] = SIGN(Pi_cipar[i]) * Pi[i]; - } - } - - mesh->communicateXZ(Pi_ciperp, Pi_cipar); - if(!fci_transform){ - Pi_ciperp.clearParallelSlices(); - Pi_cipar.clearParallelSlices(); - Pi_ciperp = toFieldAligned(Pi_ciperp); - Pi_cipar = toFieldAligned(Pi_cipar); - } - - { - int jy = 1; - for (RangeIterator r = mesh->iterateBndryLowerY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - Pi_ciperp(r.ind, jy, jz) = Pi_ciperp(r.ind, jy + 1, jz); - Pi_cipar(r.ind, jy, jz) = Pi_cipar(r.ind, jy + 1, jz); - } - } - jy = mesh->yend + 1; - for (RangeIterator r = mesh->iterateBndryUpperY(); !r.isDone(); r++) { - for (int jz = 0; jz < mesh->LocalNz; jz++) { - Pi_ciperp(r.ind, jy, jz) = Pi_ciperp(r.ind, jy - 1, jz); - Pi_cipar(r.ind, jy, jz) = Pi_cipar(r.ind, jy - 1, jz); - } - } - } - - // Smooth Pi_ciperp along Y - //Pi_ciperp.applyBoundary("neumann_o2"); - Field3D Pi_ciperp_orig = copy(Pi_ciperp); - for (auto &i : Pi_ciperp.getRegion("RGN_NOBNDRY")) { - Pi_ciperp[i] = 0.5*Pi_ciperp_orig[i] + 0.25*(Pi_ciperp_orig[i.ym()] + Pi_ciperp_orig[i.yp()]); - } - if (!fci_transform){ - Pi_ciperp = fromFieldAligned(Pi_ciperp); - Pi_cipar = fromFieldAligned(Pi_cipar); - } - mesh->communicateXZ(Pi_ciperp); - - // Apply boundary conditions - Pi_ciperp.applyBoundary("neumann_o2"); - Pi_cipar.applyBoundary("neumann_o2"); - - Pi_ci = Pi_cipar + Pi_ciperp; + ASSERT0(false); // not implemented - ask Brendan } /////////////////////////////////////////////////////////// @@ -2607,63 +2200,77 @@ int Hermes::rhs(BoutReal t) { // Parallel flow if (parallel_flow) { - if (!fci_transform){ - if (currents) { - // Parallel wave speed increased to electron sound speed - // since electrostatic & electromagnetic waves are supported - ddt(Ne) -= FV::Div_par(Ne, Ve, sqrt(mi_me) * sound_speed); - }else { - // Parallel wave speed is ion sound speed - ddt(Ne) -= FV::Div_par(Ne, Ve, sound_speed); - } - }else{ - Field3D neve = mul_all(Ne,Ve); - mesh->communicate(neve); - ddt(Ne) -= Div_parP(neve); - // b = Div_parP(neve); - // Skew-symmetric form - // Field3D gparne = Grad_par(Ne); - // Field3D dparve = Div_parP(Ve); - // mesh->communicate(gparne,dparve); - // ddt(Ne) -= 0.5 * (Div_par(neve) + mul_all(Ve,gparne) + mul_all(Ne, dparve)); - // //ddt(Ne) -= 0.5 * (Div_par(neve) + Ve * Grad_par(Ne) + Ne * Div_par(Ve)); - // a = -0.5 * (Div_parP(neve) + Ve * Grad_par(Ne) + Ne * Div_par(Ve)); - // auto* coords = mesh->getCoordinates(); - // for(auto &i : Ne.getRegion("RGN_ALL")) { - // ddt(Ne)[i] -= (Ne.yup()[i.yp()] * Ve.yup()[i.yp()] / coords->Bxy[i.yp()] - // - Ne.ydown()[i.ym()] * Ve.ydown()[i.ym()] / coords->Bxy[i.ym()]) - // * coords->Bxy[i] / (coords->dy[i] * sqrt(coords->g_22[i])); - // } - - } + // if (!fci_transform){ + // if (currents) { + // // Parallel wave speed increased to electron sound speed + // // since electrostatic & electromagnetic waves are supported + // ddt(Ne) -= FV::Div_par(Ne, Ve, sqrt(mi_me) * sound_speed); + // }else { + // // Parallel wave speed is ion sound speed + // ddt(Ne) -= FV::Div_par(Ne, Ve, sound_speed); + // } + // }else{ + check_all(Ne); + check_all(Ve); + Field3D neve = mul_all(Ne, Ve); + check_all(neve); + ddt(Ne) -= Div_parP(neve); + // b = Div_parP(neve); + // Skew-symmetric form + // Field3D gparne = Grad_par(Ne); + // Field3D dparve = Div_parP(Ve); + // mesh->communicate(gparne,dparve); + // ddt(Ne) -= 0.5 * (Div_par(neve) + mul_all(Ve,gparne) + mul_all(Ne, + // dparve)); + // //ddt(Ne) -= 0.5 * (Div_par(neve) + Ve * Grad_par(Ne) + Ne * + // Div_par(Ve)); a = -0.5 * (Div_parP(neve) + Ve * Grad_par(Ne) + Ne * + // Div_par(Ve)); auto* coords = mesh->getCoordinates(); for(auto &i : + // Ne.getRegion("RGN_ALL")) { + // ddt(Ne)[i] -= (Ne.yup()[i.yp()] * Ve.yup()[i.yp()] / + // coords->Bxy[i.yp()] + // - Ne.ydown()[i.ym()] * Ve.ydown()[i.ym()] / + // coords->Bxy[i.ym()]) + // * coords->Bxy[i] / (coords->dy[i] * sqrt(coords->g_22[i])); + // } } if (j_diamag) { // Diamagnetic drift, formulated as a magnetic drift // i.e Grad-B + curvature drift - if (!fci_transform){ - ddt(Ne) -= FV::Div_f_v(Ne, -Telim * Curlb_B, ne_bndry_flux); - } else { - ddt(Ne) -= fci_curvature(Pelim); - } + ddt(Ne) -= fci_curvature(Pe); } if (ramp_mesh && (t < ramp_timescale)) { ddt(Ne) += NeTarget / ramp_timescale; } + Field3D TiTediff, tauemimeSQB; if (classical_diffusion) { // Classical perpendicular diffusion // The only term here comes from the resistive drift - Dn = (Telim + Tilim) / (tau_e * mi_me * SQ(coord->Bxy)); - if(fci_transform){ - mesh->communicate(Dn); - Field3D Ne_tauB2 = Ne / (tau_e * mi_me * SQ(coord->Bxy)); - Field3D TiTediff = Ti - 0.5 * Te; - mesh->communicate(Ne_tauB2, TiTediff); - ddt(Ne) += FV::Div_a_Laplace_perp(Dn, Ne); - ddt(Ne) += FV::Div_a_Laplace_perp(Ne_tauB2, TiTediff); - } + Field3D Ne_tauB2; + alloc_all(TiTediff); + alloc_all(tauemimeSQB); + alloc_all(Ne_tauB2); + alloc_all(Dn); + BOUT_FOR(i, Ne.getRegion("RGN_ALL")) { + tauemimeSQB[i] = tau_e[i] * mi_me * B42[i]; + Dn[i] = (Te[i] + Ti[i]) / tauemimeSQB[i]; + Ne_tauB2[i] = Ne[i] / tauemimeSQB[i]; + TiTediff[i] = Ti[i] - (0.5 * Te[i]); + + tauemimeSQB.yup()[i] = tau_e.yup()[i] * mi_me * B42.yup()[i]; + Dn.yup()[i] = (Te.yup()[i] + Ti.yup()[i]) / tauemimeSQB.yup()[i]; + Ne_tauB2.yup()[i] = Ne.yup()[i] / tauemimeSQB.yup()[i]; + TiTediff.yup()[i] = Ti.yup()[i] - (0.5 * Te.yup()[i]); + + tauemimeSQB.ydown()[i] = tau_e.ydown()[i] * mi_me * B42.ydown()[i]; + Dn.ydown()[i] = (Te.ydown()[i] + Ti.ydown()[i]) / tauemimeSQB.ydown()[i]; + Ne_tauB2.ydown()[i] = Ne.ydown()[i] / tauemimeSQB.ydown()[i]; + TiTediff.ydown()[i] = Ti.ydown()[i] - (0.5 * Te.ydown()[i]); + } + ddt(Ne) += FV::Div_a_Laplace_perp(Dn, Ne); + ddt(Ne) += FV::Div_a_Laplace_perp(Ne_tauB2, TiTediff); } if (anomalous_D > 0.0) { ddt(Ne) += FV::Div_a_Laplace_perp(a_d3d, Ne); @@ -2716,15 +2323,10 @@ int Hermes::rhs(BoutReal t) { if (low_n_diffuse) { // Diffusion which kicks in at very low density, in order to // help prevent negative density regions - if(fci_transform){ - mesh->communicate(Ne); - ddt(Ne) += Div_par_K_Grad_par(SQ(coord->dy) * coord->g_22 * 1e-4 / Nelim, Ne); - }else{ - ddt(Ne) += FV::Div_par_K_Grad_par(SQ(coord->dy) * coord->g_22 * 1e-4 / Nelim, Ne); - } + ddt(Ne) += Div_par_K_Grad_par(SQ(coord->dy) * coord->g_22 * 1e-4 / Ne, Ne); } if (low_n_diffuse_perp) { - ddt(Ne) += Div_Perp_Lap_FV_Index(1e-4 / Nelim, Ne, ne_bndry_flux); + ddt(Ne) += Div_Perp_Lap_FV_Index(1e-4 / Ne, Ne, ne_bndry_flux); } if (ne_hyper_z > 0.) { @@ -2741,7 +2343,7 @@ int Hermes::rhs(BoutReal t) { } if (numdiff > 0.0) { - for(auto &i : Ne.getRegion("RGN_NOBNDRY")) { + BOUT_FOR(i, Ne.getRegion("RGN_NOBNDRY")) { ddt(Ne)[i] += numdiff*(Ne.ydown()[i.ym()] - 2.*Ne[i] + Ne.yup()[i.yp()]); } } @@ -2760,12 +2362,10 @@ int Hermes::rhs(BoutReal t) { if (j_par) { TRACE("Vort:j_par"); // Parallel current - + // This term is central differencing so that it balances the parallel gradient // of the potential in Ohm's law - if(fci_transform){mesh->communicate(Jpar);} ddt(Vort) += Div_parP(Jpar); - a = Div_parP(Jpar); } if (j_diamag) { @@ -2773,23 +2373,14 @@ int Hermes::rhs(BoutReal t) { // Note: This term is central differencing so that it balances // the corresponding compression term in the pressure equation - if(!fci_transform){ - ddt(Vort) += Div((Pi + Pe) * Curlb_B); - }else{ - ddt(Vort) += fci_curvature(Pi + Pe); - b = fci_curvature(Pi + Pe); - } + ddt(Vort) += fci_curvature(Pi + Pe); } // Advection of vorticity by ExB if (boussinesq) { TRACE("Vort:boussinesq"); // Using the Boussinesq approximation - if(!fci_transform){ - ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(0.5 * Vort, phi, vort_bndry_flux, - poloidal_flows, false); - }else{//fci used - if (j_pol_pi){ + if (j_pol_pi){ ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(0.5 * Vort, phi, vort_bndry_flux, poloidal_flows, false) * bracket_factor; // V_ExB dot Grad(Pi) @@ -2799,27 +2390,14 @@ int Hermes::rhs(BoutReal t) { Field3D DelpPhi_2B2 = 0.5 * Delp2(phi) / SQ(Bxyz); DelpPhi_2B2.applyBoundary("free_o2"); - mesh->communicate(vEdotGradPi, DelpPhi_2B2); - if(!fci_transform){ - ddt(Vort) -= FV::Div_a_Laplace_perp(0.5 / SQ(coord->Bxy), vEdotGradPi); - }else{ - Field3D inv_2sqb = 0.5 / SQ(Bxyz); - mesh->communicate(inv_2sqb); - ddt(Vort) -= FV::Div_a_Laplace_perp(inv_2sqb, vEdotGradPi); - } + Field3D inv_2sqb = 0.5 / SQ(Bxyz); + ddt(Vort) -= FV::Div_a_Laplace_perp(inv_2sqb, vEdotGradPi); // delp2 phi v_ExB term ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(DelpPhi_2B2, phi + Pi, vort_bndry_flux, poloidal_flows) * bracket_factor; - }else if (j_pol_simplified) { - // use simplified polarization term from i.e. GBS - ddt(Vort) -= Div_n_bxGrad_f_B_XPPM(Vort, phi, vort_bndry_flux, - poloidal_flows, false) * bracket_factor; - c = Div_n_bxGrad_f_B_XPPM(Vort, phi, vort_bndry_flux, - poloidal_flows, false) * bracket_factor; - } } @@ -2835,13 +2413,11 @@ int Hermes::rhs(BoutReal t) { if (classical_diffusion) { TRACE("Vort:classical_diffusion"); // Perpendicular viscosity - Field3D tilim_3 = 0.3*Tilim; + Field3D tilim_3 = 0.3*Ti; Field3D tauisqB = tau_i * SQ(coord->Bxy); - mesh->communicate(tilim_3,tauisqB); + Field3D mu = div_all(tilim_3 , tauisqB); - if (fci_transform) {mesh->communicate(mu);} ddt(Vort) += FV::Div_a_Laplace_perp(mu, Vort); - d = FV::Div_a_Laplace_perp(mu, Vort); } if (ion_viscosity) { @@ -2882,34 +2458,12 @@ int Hermes::rhs(BoutReal t) { } if (vort_dissipation) { - if(!fci_transform){ - // Adds dissipation term like in other equations - // Maximum speed either electron sound speed or Alfven speed - Field3D max_speed = Bnorm * coord->Bxy / - sqrt(SI::mu0 * AA * SI::Mp * Nnorm * Nelim) / - Cs0; // Alfven speed (normalised by Cs0) - Field3D elec_sound = sqrt(mi_me) * sound_speed; // Electron sound speed - for (auto& i : max_speed.getRegion("RGN_ALL")) { - if (elec_sound[i] > max_speed[i]) { - max_speed[i] = elec_sound[i]; - } - - // Limit to 100x reference sound speed or light speed - BoutReal lim = BOUTMIN(100., 3e8/Cs0); - if (max_speed[i] > lim) { - max_speed[i] = lim; - } - } - ddt(Vort) -= FV::Div_par(Vort, 0.0, max_speed); - }else{ - ddt(Vort) += SQ(coord->dy)*D2DY2(Vort); - } + ddt(Vort) += SQ(coord->dy)*D2DY2(Vort); } if (phi_dissipation) { // Adds dissipation term like in other equations, but depending on gradient of potential // Note: Doesn't seem to need faster than sound speed ddt(Vort) -= SQ(coord->dy)*D2DY2(phi); - f = SQ(coord->dy)*D2DY2(phi); } } @@ -2920,69 +2474,57 @@ int Hermes::rhs(BoutReal t) { ddt(VePsi) = 0.0; - Field3D NelimVe = Nelim; + Field3D NeVe = Ne; +#warning NeVe is Ne? if (currents && (electromagnetic || FiniteElMass)) { // Evolve VePsi except for electrostatic and zero electron mass case if (resistivity) { - ddt(VePsi) -= mi_me * nu * sub_all(Ve , Vi); + ddt(VePsi) -= mi_me * nu * (Ve - Vi); // External electric field - // ddt(VePsi) += mi_me*nu*(Jpar - Jpar0)/NelimVe; + // ddt(VePsi) += mi_me*nu*(Jpar - Jpar0)/NeVe; } // Parallel electric field if (j_par) { - if(fci_transform){mesh->communicate(phi);} ddt(VePsi) += mi_me * Grad_parP(phi); } // Parallel electron pressure if (pe_par) { - if(fci_transform){mesh->communicate(Pelim);} - ddt(VePsi) -= mi_me * Grad_parP(Pelim) / NelimVe; + ddt(VePsi) -= mi_me * Grad_parP(Pe) / NeVe; } if (thermal_force) { - if(fci_transform){mesh->communicate(Te);} ddt(VePsi) -= mi_me * 0.71 * Grad_parP(Te); } if (electron_viscosity) { // Electron parallel viscosity (Braginskii) - Field3D ve_eta = 0.973 * mi_me * tau_e * Telim; + Field3D ve_eta = 0.973 * mi_me * tau_e * Te; if (eta_limit_alpha > 0.) { // SOLPS-style flux limiter // Values of alpha ~ 0.5 typically - if(fci_transform){mesh->communicate(Ve);} Field3D q_cl = ve_eta * Grad_par(Ve); // Collisional value - Field3D q_fl = eta_limit_alpha * Pelim * mi_me; // Flux limit + Field3D q_fl = eta_limit_alpha * Pe * mi_me; // Flux limit ve_eta = ve_eta / (1. + abs(q_cl / q_fl)); + if (fci_transform) { + ASSERT0("not implemented: FCI comm" == nullptr); + } mesh->communicate(ve_eta); ve_eta.applyBoundary("neumann"); } - if(fci_transform){ - mesh->communicate(ve_eta,Ve); - ddt(VePsi) += Div_par_K_Grad_par(ve_eta, Ve); - }else{ - ddt(VePsi) += FV::Div_par_K_Grad_par(ve_eta, Ve); - } + ddt(VePsi) += Div_par_K_Grad_par(ve_eta, Ve); } if (FiniteElMass) { // Finite Electron Mass. Small correction needed to conserve energy - if(!fci_transform){ - ddt(VePsi) -= Vi * Grad_par(Ve - Vi); // Parallel advection - }else{ - Field3D vdiff = sub_all(Ve,Vi); - mesh->communicate(vdiff); - ddt(VePsi) -= Vi * Grad_par(vdiff); // Parallel advection - } Field3D vdiff = sub_all(Ve,Vi); - // mesh->communicate(vdiff); + ddt(VePsi) -= Vi * Grad_par(vdiff); // Parallel advection ddt(VePsi) -= bracket(phi, vdiff, BRACKET_ARAKAWA)*bracket_factor; // ExB advection // Should also have ion polarisation advection here } @@ -2996,7 +2538,7 @@ int Hermes::rhs(BoutReal t) { } if (hyper > 0.0) { - ddt(VePsi) -= hyper * mi_me * nu * Delp2(Jpar) / Nelim; + ddt(VePsi) -= hyper * mi_me * nu * Delp2(Jpar) / Ne; } if (hyperpar > 0.0) { @@ -3014,7 +2556,7 @@ int Hermes::rhs(BoutReal t) { if (vepsi_dissipation) { // Adds dissipation term like in other equations // Maximum speed either electron sound speed or Alfven speed - Field3D max_speed = Bnorm * coord->Bxy / sqrt(SI::mu0 * AA * SI::Mp * Nnorm * Nelim) + Field3D max_speed = Bnorm * coord->Bxy / sqrt(SI::mu0 * AA * SI::Mp * Nnorm * Ne) / Cs0; // Alfven speed (normalised by Cs0) Field3D elec_sound = sqrt(mi_me) * sound_speed; // Electron sound speed for (auto& i : max_speed.getRegion("RGN_NOBNDRY")) { @@ -3029,11 +2571,7 @@ int Hermes::rhs(BoutReal t) { max_speed[i] = lim; } } - if (!fci_transform) { - ddt(VePsi) -= FV::Div_par(Ve - Vi, 0.0, max_speed); - } else { - ddt(VePsi) += SQ(coord->dy) * (D2DY2(Ve) - D2DY2(Vi)); - } + ddt(VePsi) += SQ(coord->dy) * (D2DY2(Ve) - D2DY2(Vi)); } } @@ -3043,72 +2581,42 @@ int Hermes::rhs(BoutReal t) { TRACE("Ion velocity"); if (currents) { - if(fci_transform){ // ddt(NVi) = bracket(NVi, phi, BRACKET_ARAKAWA) * bracket_factor; // ExB drift, only if electric field calculated ddt(NVi) = -Div_n_bxGrad_f_B_XPPM(NVi, phi, ne_bndry_flux, poloidal_flows) * bracket_factor; // ExB drift - - }else{ - // ExB drift, only if electric field calculated - ddt(NVi) = -Div_n_bxGrad_f_B_XPPM(NVi, phi, ne_bndry_flux, - poloidal_flows); // ExB drift - } } else { ddt(NVi) = 0.0; } if (j_diamag) { // Magnetic drift - if(!fci_transform){ - ddt(NVi) -= FV::Div_f_v(NVi, Tilim * Curlb_B, - ne_bndry_flux); // Grad-B, curvature drift - }else{ - ddt(NVi) -= fci_curvature(NVi * Tilim); - } - } - if(!fci_transform){ - ddt(NVi) -= FV::Div_par_fvv(Ne, Vi, sound_speed, false); - }else{ - Field3D nvivi = mul_all(NVi, Vi); - mesh->communicate(nvivi); - - ddt(NVi) -= Div_parP(nvivi); - // Skew-symmetric form - // ddt(NVi) -= 0.5 * (Div_par(mul_all(NVi, Vi)) + Vi * Grad_par(NVi) + NVi * Div_par(Vi)); + ddt(NVi) -= fci_curvature(NVi * Ti); } + Field3D nvivi = mul_all(NVi, Vi); + ddt(NVi) -= Div_parP(nvivi); + // Skew-symmetric form + // ddt(NVi) -= 0.5 * (Div_par(mul_all(NVi, Vi)) + Vi * Grad_par(NVi) + NVi * Div_par(Vi)); // Ignoring polarisation drift for now if (pe_par) { - if(!fci_transform){ - ddt(NVi) -= Grad_parP(Pe + Pi); - }else{ - Field3D peppi = add_all(Pe,Pi); - mesh->communicate(peppi); - ddt(NVi) -= Grad_parP(peppi); - } + Field3D peppi = add_all(Pe, Pi); + ddt(NVi) -= Grad_parP(peppi); } if (ion_viscosity) { TRACE("NVi:ion viscosity"); // Poloidal flow damping - if(fci_transform){ - // The parallel part is solved as a diffusion term - Coordinates::FieldMetric sqrtBVi = sqrtB * Vi; - Coordinates::FieldMetric Pitau_i_B = Pi * tau_i / (coord->Bxy); - mesh->communicate(sqrtBVi,Pitau_i_B); - ddt(NVi) += 1.28 * sqrtB * - Div_par_K_Grad_par(Pitau_i_B, sqrtBVi); - }else{ - ddt(NVi) += 1.28 * sqrtB * - FV::Div_par_K_Grad_par(Pi * tau_i / (coord->Bxy), sqrtB * Vi); - } - + // The parallel part is solved as a diffusion term + Coordinates::FieldMetric sqrtBVi = mul_all(B12, Vi); + Coordinates::FieldMetric Pitau_i_B = + div_all(mul_all(Pi, tau_i), (coord->Bxy)); + ddt(NVi) += 1.28 * B12 * Div_par_K_Grad_par(Pitau_i_B, sqrtBVi); + if (currents) { // Perpendicular part. B32 = B^{3/2} // This is only included if ExB flow is included - Field3D Piciperp_B32 = Pi_ciperp/B32; - mesh->communicate(Piciperp_B32); + Field3D Piciperp_B32 = div_all(Pi_ciperp, B32); ddt(NVi) -= (2. / 3) * B32 * Grad_par(Piciperp_B32); } } @@ -3119,7 +2627,7 @@ int Hermes::rhs(BoutReal t) { } if (numdiff > 0.0) { - for(auto &i : NVi.getRegion("RGN_ALL")) { + for(auto &i : NVi.getRegion("RGN_NOY")) { ddt(NVi)[i] += numdiff*(Vi.ydown()[i.ym()] - 2.*Vi[i] + Vi.yup()[i.yp()]); } // ddt(NVi) += numdiff * Div_par_diffusion_index(NVi); @@ -3137,11 +2645,8 @@ int Hermes::rhs(BoutReal t) { if (classical_diffusion) { // Using same cross-field drift as in density equation Field3D ViDn = mul_all(Vi,Dn); - mesh->communicate(ViDn); ddt(NVi) += FV::Div_a_Laplace_perp(ViDn, Ne); - Field3D NVi_tauB2 = NVi / (tau_e * mi_me * SQ(coord->Bxy)); - Field3D TiTediff = Ti - 0.5 * Te; - mesh->communicate(NVi_tauB2, TiTediff); + Field3D NVi_tauB2 = div_all(NVi, tauemimeSQB); ddt(NVi) += FV::Div_a_Laplace_perp(NVi_tauB2, TiTediff); } @@ -3160,20 +2665,16 @@ int Hermes::rhs(BoutReal t) { if (low_n_diffuse) { // Diffusion which kicks in at very low density, in order to // help prevent negative density regions - if(fci_transform){ - mesh->communicate(NVi); - ddt(NVi) += Div_par_K_Grad_par(SQ(coord->dy) * coord->g_22 * 1e-4 / Nelim, NVi); - }else{ - ddt(NVi) += FV::Div_par_K_Grad_par(SQ(coord->dy) * coord->g_22 * 1e-4 / Nelim, NVi); - } + ASSERT0(false); + ddt(NVi) += Div_par_K_Grad_par(SQ(coord->dy) * coord->g_22 * 1e-4 / Ne, NVi); } if (low_n_diffuse_perp) { - ddt(NVi) += Div_Perp_Lap_FV_Index(1e-4 / Nelim, NVi, ne_bndry_flux); + ddt(NVi) += Div_Perp_Lap_FV_Index(1e-4 / Ne, NVi, ne_bndry_flux); } if (vi_num_diff > 0.0) { // Numerical perpendicular diffusion - ddt(NVi) += Div_Perp_Lap_FV_Index(vi_num_diff * Nelim, Vi, ne_bndry_flux); + ddt(NVi) += Div_Perp_Lap_FV_Index(vi_num_diff * Ne, Vi, ne_bndry_flux); } } @@ -3184,39 +2685,22 @@ int Hermes::rhs(BoutReal t) { if (evolve_te) { if (currents) { - if(fci_transform){ - // ddt(Pe) = bracket(Pe, phi, BRACKET_ARAKAWA) * bracket_factor; - ddt(Pe) = -Div_n_bxGrad_f_B_XPPM(Pe, phi, pe_bndry_flux, poloidal_flows, true) * bracket_factor; - }else{ // Divergence of heat flux due to ExB advection - ddt(Pe) = -Div_n_bxGrad_f_B_XPPM(Pe, phi, pe_bndry_flux, poloidal_flows, true); - } + // ddt(Pe) = bracket(Pe, phi, BRACKET_ARAKAWA) * bracket_factor; + ddt(Pe) = -Div_n_bxGrad_f_B_XPPM(Pe, phi, pe_bndry_flux, poloidal_flows, true) * bracket_factor; } else { ddt(Pe) = 0.0; } if (parallel_flow_p_term) { // Parallel flow - if (fci_transform){ - Field3D peve = mul_all(Pe,Ve); - mesh->communicate(peve); - ddt(Pe) -= Div_parP(peve); - } else { - if (currents) { - ddt(Pe) -= FV::Div_par(Pe, Ve, sqrt(mi_me) * sound_speed); - } else { - ddt(Pe) -= FV::Div_par(Pe, Ve, sound_speed); - } - } + Field3D peve = mul_all(Pe,Ve); + ddt(Pe) -= Div_parP(peve); } if (j_diamag) { // Diamagnetic flow // Magnetic drift (curvature) divergence. - if (fci_transform) { - ddt(Pe) -= (5. / 3) * fci_curvature(Pe * Telim); - } else { - ddt(Pe) -= (5. / 3) * FV::Div_f_v(Pe, -Telim * Curlb_B, pe_bndry_flux); - } + ddt(Pe) -= (5. / 3) * fci_curvature(Pe * Te); // This term energetically balances diamagnetic term // in the vorticity equation @@ -3226,28 +2710,19 @@ int Hermes::rhs(BoutReal t) { // Parallel heat conduction if (thermal_conduction) { - if (fci_transform) { - mesh->communicate(kappa_epar, Te); - ddt(Pe) += (2. / 3) * Div_par_K_Grad_par(kappa_epar, Te); - } else { - ddt(Pe) += (2. / 3) * FV::Div_par_K_Grad_par(kappa_epar, Te); - } + check_all(kappa_epar); + ddt(Pe) += (2. / 3) * Div_par_K_Grad_par(kappa_epar, Te); } if (thermal_flux) { // Parallel heat convection - if (fci_transform) { - Field3D tejpar = mul_all(Te,Jpar); - mesh->communicate(tejpar); - ddt(Pe) += (2. / 3) * 0.71 * Div_parP(tejpar); - } else { - ddt(Pe) += (2. / 3) * 0.71 * Div_par(Te * Jpar); - } + Field3D tejpar = mul_all(Te,Jpar); + ddt(Pe) += (2. / 3) * 0.71 * Div_parP(tejpar); } if (currents && resistivity) { // Ohmic heating - ddt(Pe) += nu * Jpar * (Jpar - Jpar0) / Nelim; + ddt(Pe) += nu * Jpar * (Jpar - Jpar0) / Ne; } if (pe_hyper_z > 0.0) { @@ -3355,76 +2830,41 @@ int Hermes::rhs(BoutReal t) { if (parallel_sheaths){ sheath_dpe = 0.; - for (const auto &bndry_par : mesh->getBoundariesPar()) { + for (const auto &bndry_par : mesh->getBoundariesPar(BoundaryParType::xout)) { for (bndry_par->first(); !bndry_par->isDone(); bndry_par->next()) { int x = bndry_par->x; int y = bndry_par->y; int z = bndry_par->z; - if (bndry_par->dir == 1){ // forwards - // Temperature and density at the sheath entrance - BoutReal tesheath = floor( - 0.5 * (Te(x, y, z) + Te.yup()(x, y + bndry_par->dir, z)), - 0.0); - BoutReal nesheath = floor( - 0.5 * (Ne(x, y, z) + Ne.yup()(x, y + bndry_par->dir, z)), - 0.0); - BoutReal vesheath = floor( - 0.5 * (Ve(x, y, z) + Ve.yup()(x, y + bndry_par->dir, z)), - 0.0); - // BoutReal tisheath = floor( - // 0.5 * (Ti(x, y, z) + Ti.yup()(x, y + bndry_par->dir, z)), - // 0.0); - - // Sound speed (normalised units) - // BoutReal Cs = sqrt(tesheath + tisheath); - - // Heat flux - BoutReal q = (sheath_gamma_e - 1.5) * tesheath * nesheath * vesheath; - // Multiply by cell area to get power - BoutReal flux = - q - * (coord->J(x, y, z) + coord->J.yup()(x, y + bndry_par->dir, z)) - / (sqrt(coord->g_22(x, y, z)) - + sqrt(coord->g_22.yup()(x, y + bndry_par->dir, z))); - - // Divide by volume of cell, and 2/3 to get pressure - BoutReal power = - flux - / (coord->dy(x, y, z) * coord->J(x, y, z)); - // ddt(Pe)(x, y, z) -= (2. / 3) * power; - sheath_dpe(x, y, z) -= (2. / 3) * power; - } else { // backwards - // Temperature and density at the sheath entrance - BoutReal tesheath = floor( - 0.5 * (Te(x, y, z) + Te.ydown()(x, y + bndry_par->dir, z)), - 0.0); - BoutReal nesheath = floor( - 0.5 * (Ne(x, y, z) + Ne.ydown()(x, y + bndry_par->dir, z)), - 0.0); - BoutReal vesheath = floor( - 0.5 * (Ve(x, y, z) + Ve.ydown()(x, y + bndry_par->dir, z)), - 0.0); - // BoutReal tisheath = floor( - // 0.5 * (Ti(x, y, z) + Ti.ydown()(x, y + bndry_par->dir, z)), - // 0.0); - - // Sound speed (normalised units) - // BoutReal Cs = -sqrt(tesheath + tisheath); - - // Heat flux - BoutReal q = (sheath_gamma_e - 1.5) * tesheath * nesheath * vesheath; - - // Multiply by cell area to get power - BoutReal flux = - q - * (coord->J(x, y, z) + coord->J.ydown()(x, y + bndry_par->dir, z)) - / (sqrt(coord->g_22(x, y, z)) - + sqrt(coord->g_22.ydown()(x, y + bndry_par->dir, z))); - - // Divide by volume of cell, and 2/3 to get pressure - BoutReal power = - flux - / (coord->dy(x, y, z) * coord->J(x, y, z)); - sheath_dpe(x, y, z) -= (2. / 3) * power; - } + // Temperature and density at the sheath entrance + BoutReal tesheath = floor( + 0.5 * (Te(x, y, z) + Te.ynext(bndry_par->dir)(x, y + bndry_par->dir, z)), + 0.0); + BoutReal nesheath = floor( + 0.5 * (Ne(x, y, z) + Ne.ynext(bndry_par->dir)(x, y + bndry_par->dir, z)), + 0.0); + BoutReal vesheath = floor( + 0.5 * (Ve(x, y, z) + Ve.ynext(bndry_par->dir)(x, y + bndry_par->dir, z)), + 0.0); + // BoutReal tisheath = floor( + // 0.5 * (Ti(x, y, z) + Ti.ynext(bndry_par->dir)(x, y + bndry_par->dir, z)), + // 0.0); + + // Sound speed (normalised units) + // BoutReal Cs =bndry_par->dir* sqrt(tesheath + tisheath); + + // Heat flux + BoutReal q = (sheath_gamma_e - 1.5) * tesheath * nesheath * vesheath; + // Multiply by cell area to get power + BoutReal flux = + q + * (coord->J(x, y, z) + coord->J.ynext(bndry_par->dir)(x, y + bndry_par->dir, z)) + / (sqrt(coord->g_22(x, y, z)) + + sqrt(coord->g_22.ynext(bndry_par->dir)(x, y + bndry_par->dir, z))); + + // Divide by volume of cell, and 2/3 to get pressure + BoutReal power = + flux + / (coord->dy(x, y, z) * coord->J(x, y, z)); + // ddt(Pe)(x, y, z) -= (2. / 3) * power; + sheath_dpe(x, y, z) -= (2. / 3) * power; } } ddt(Pe) += sheath_dpe; @@ -3433,19 +2873,13 @@ int Hermes::rhs(BoutReal t) { // Transfer and source terms if (thermal_force) { - if (fci_transform) { - mesh->communicate(Te); - } ddt(Pe) -= (2. / 3) * 0.71 * Jpar * Grad_parP(Te); } if (pe_par_p_term) { // This term balances energetically the pressure term // in Ohm's law - if (fci_transform) { - mesh->communicate(Ve); - } - ddt(Pe) -= (2. / 3) * Pelim * Div_parP(Ve); + ddt(Pe) -= (2. / 3) * Pe * Div_parP(Ve); } if (ramp_mesh && (t < ramp_timescale)) { ddt(Pe) += PeTarget / ramp_timescale; @@ -3458,11 +2892,9 @@ int Hermes::rhs(BoutReal t) { // Combined resistive drift and cross-field heat diffusion // nu_rho2 = nu_ei * rho_e^2 in normalised units - Field3D nu_rho2 = Telim / (tau_e * mi_me * SQ(coord->Bxy)); + Field3D nu_rho2 = div_all(Te, mul_all(mul_all(tau_e, mi_me), B42)); Field3D PePi = add_all(Pe, Pi); - mesh->communicate(nu_rho2); Field3D nu_rho2Ne = mul_all(nu_rho2, Ne); - mesh->communicate(nu_rho2Ne,Te); ddt(Pe) += (2. / 3) * (FV::Div_a_Laplace_perp(nu_rho2, PePi) + (11. / 12) * FV::Div_a_Laplace_perp(nu_rho2Ne, Te)); @@ -3509,7 +2941,7 @@ int Hermes::rhs(BoutReal t) { if (energy_source) { // Add the same amount of energy to each particle - PeSource = Spe * Nelim / DC(Nelim); + PeSource = Spe * Ne / DC(Ne); } else { PeSource = Spe; } @@ -3520,7 +2952,7 @@ int Hermes::rhs(BoutReal t) { if (energy_source) { // Add the same amount of energy to each particle - PeSource = Spe * Nelim / DC(Nelim); + PeSource = Spe * Ne / DC(Ne); } else { PeSource = Spe * where(Spe, PeTarget, Pe); } @@ -3535,7 +2967,7 @@ int Hermes::rhs(BoutReal t) { if (energy_source) { // Add the same amount of energy to each particle - PeSource = Spe * Nelim / DC(Nelim); + PeSource = Spe * Ne / DC(Ne); if (source_vary_g11) { PeSource *= g11norm; @@ -3560,35 +2992,22 @@ int Hermes::rhs(BoutReal t) { if (evolve_ti) { if (currents) { - if(fci_transform){ - // ddt(Pi) = bracket(Pi, phi, BRACKET_ARAKAWA) * bracket_factor; - ddt(Pi) = -Div_n_bxGrad_f_B_XPPM(Pi, phi, pe_bndry_flux, poloidal_flows, true) * bracket_factor; - }else{ - // Divergence of heat flux due to ExB advection - ddt(Pi) = -Div_n_bxGrad_f_B_XPPM(Pi, phi, pe_bndry_flux, poloidal_flows, true); - } + // Divergence of heat flux due to ExB advection + // ddt(Pi) = bracket(Pi, phi, BRACKET_ARAKAWA) * bracket_factor; + ddt(Pi) = -Div_n_bxGrad_f_B_XPPM(Pi, phi, pe_bndry_flux, poloidal_flows, true) * bracket_factor; } else { ddt(Pi) = 0.0; } // Parallel flow if (parallel_flow_p_term) { - if (fci_transform) { - Field3D pivi = mul_all(Pi,Vi); - mesh->communicate(pivi); - ddt(Pi) -= Div_parP(pivi); - } else { - ddt(Pi) -= FV::Div_par(Pi, Vi, sound_speed); - } + Field3D pivi = mul_all(Pi,Vi); + ddt(Pi) -= Div_parP(pivi); } if (j_diamag) { // Diamagnetic flow // Magnetic drift (curvature) divergence - if (fci_transform) { - ddt(Pi) -= (5. / 3) * fci_curvature(Pi * Tilim); - } else { - ddt(Pi) -= (5. / 3) * FV::Div_f_v(Pi, Tilim * Curlb_B, pe_bndry_flux); - } + ddt(Pi) -= (5. / 3) * fci_curvature(Pi * Ti); // Compression of ExB flow // These terms energetically balances diamagnetic term @@ -3596,49 +3015,32 @@ int Hermes::rhs(BoutReal t) { // ddt(Pi) -= (2. / 3) * Pi * (Curlb_B * Grad(phi)); ddt(Pi) -= (2. / 3) * Pi * fci_curvature(phi); - if (fci_transform) { - // Vector3D pipecb = (Pi + Pe); - // mesh->communicate(pipe); - ddt(Pi) += Pi * fci_curvature(Pi + Pe); - } else { - ddt(Pi) += Pi * Div((Pe + Pi) * Curlb_B); - } + ddt(Pi) += Pi * fci_curvature(Pi + Pe); } if (j_par) { - if (fci_transform) { - mesh->communicate(Pi); - } if (boussinesq) { ddt(Pi) -= (2. / 3) * Jpar * Grad_parP(Pi); } else { - ddt(Pi) -= (2. / 3) * Jpar * Grad_parP(Pi) / Nelim; + ddt(Pi) -= (2. / 3) * Jpar * Grad_parP(Pi) / Ne; } } // Parallel heat conduction if (thermal_conduction) { - if (fci_transform) { - mesh->communicate(kappa_ipar, Ti); - ddt(Pi) += (2. / 3) * Div_par_K_Grad_par(kappa_ipar, Ti); - } else { - ddt(Pi) += (2. / 3) * FV::Div_par_K_Grad_par(kappa_ipar, Ti); - } + ddt(Pi) += (2. / 3) * Div_par_K_Grad_par(kappa_ipar, Ti); } // Parallel pressure gradients (sound waves) if (pe_par_p_term) { // This term balances energetically the pressure term // in the parallel momentum equation - if (fci_transform) { - mesh->communicate(Vi); - } - ddt(Pi) -= (2. / 3) * Pilim * Div_parP(Vi); + ddt(Pi) -= (2. / 3) * Pi * Div_parP(Vi); } if (electron_ion_transfer) { // Electron-ion heat transfer - Wi = (3. / mi_me) * Nelim * (Te - Ti) / tau_e; + Wi = (3. / mi_me) * Ne * (Te - Ti) / tau_e; ddt(Pi) += (2. / 3) * Wi; ddt(Pe) -= (2. / 3) * Wi; } @@ -3647,50 +3049,66 @@ int Hermes::rhs(BoutReal t) { // Classical diffusion if (classical_diffusion) { - // Cross-field heat conduction - // kappa_perp = 2 * n * nu_ii * rho_i^2 - Field3D Pi_B2tau = 2. * Pilim / (SQ(coord->Bxy) * tau_i); - mesh->communicate(Pi_B2tau); + Field3D Pi_B2tau, PePi, nu_rho2, nu_rho2Ne; + alloc_all(Pi_B2tau); + alloc_all(PePi); + alloc_all(nu_rho2); + alloc_all(nu_rho2Ne); + BOUT_FOR(i, Pi.getRegion("RGN_ALL")) { + // Cross-field heat conduction + // kappa_perp = 2 * n * nu_ii * rho_i^2 + Pi_B2tau[i] = (2. * Pi[i]) / (B42[i] * tau_i[i]); + nu_rho2[i] = Te[i] / (tau_e[i] * mi_me * B42[i]); + PePi[i] = Pe[i] + Pi[i]; + nu_rho2Ne[i] = nu_rho2[i] * Ne[i]; + + Pi_B2tau.yup()[i] = + (2. * Pi.yup()[i]) / (B42.yup()[i] * tau_i.yup()[i]); + nu_rho2.yup()[i] = + Te.yup()[i] / (tau_e.yup()[i] * mi_me * B42.yup()[i]); + PePi.yup()[i] = Pe.yup()[i] + Pi.yup()[i]; + nu_rho2Ne.yup()[i] = nu_rho2.yup()[i] * Ne.yup()[i]; + + Pi_B2tau.ydown()[i] = + (2. * Pi.ydown()[i]) / (B42.ydown()[i] * tau_i.ydown()[i]); + nu_rho2.ydown()[i] = + Te.ydown()[i] / (tau_e.ydown()[i] * mi_me * B42.ydown()[i]); + PePi.ydown()[i] = Pe.ydown()[i] + Pi.ydown()[i]; + nu_rho2Ne.ydown()[i] = nu_rho2.ydown()[i] * Ne.ydown()[i]; + } + + // BOUT_FOR(i, Pi.getRegion("RGN_NOBNDRY")) { ddt(Pi) += (2. / 3) * FV::Div_a_Laplace_perp(Pi_B2tau, Ti); // Resistive drift terms - // nu_rho2 = (Ti/Te) * nu_ei * rho_e^2 in normalised units - Field3D nu_rho2 = Tilim / (tau_e * mi_me * SQ(coord->Bxy)); - Field3D PePi = add_all(Pe , Pi); - mesh->communicate(nu_rho2,PePi); - Field3D nu_rho2Ne = mul_all(nu_rho2, Ne); - mesh->communicate(nu_rho2Ne,Te); - ddt(Pi) += (5. / 3) - * (FV::Div_a_Laplace_perp(nu_rho2, PePi) - - (3. / 2) * FV::Div_a_Laplace_perp(nu_rho2Ne, Te)); + // mesh->communicate(nu_rho2Ne,Te); + ddt(Pi) += (5. / 3) * (FV::Div_a_Laplace_perp(nu_rho2, PePi) - + (1.5) * FV::Div_a_Laplace_perp(nu_rho2Ne, Te)); // Collisional heating from perpendicular viscosity // in the vorticity equation if (currents) { Vector3D Grad_perp_vort = Grad(Vort); - Field3D phiPi = phi+Pi; - mesh->communicate(phiPi); + Field3D phiPi = add_all(phi, Pi); Grad_perp_vort.y = 0.0; // Zero parallel component - ddt(Pi) -= (2. / 3) * (3. / 10) * Tilim / (SQ(coord->Bxy) * tau_i) + ddt(Pi) -= (2. / 3) * (3. / 10) * Ti / (SQ(coord->Bxy) * tau_i) * (Grad_perp_vort * Grad(phiPi)); } } if (ion_viscosity) { // Collisional heating due to parallel viscosity - Field3D sqrtBVi = mul_all(sqrtB,Vi); - mesh->communicate(sqrtBVi,Vi); - ddt(Pi) += - (2. / 3) * 1.28 * (Pi * tau_i / sqrtB) * Grad_par(sqrtBVi) * Div_parP(Vi); + Field3D sqrtBVi = mul_all(B12, Vi); + ddt(Pi) += (2. / 3) * 1.28 * (Pi * tau_i / B12) * Grad_par(sqrtBVi) * + Div_parP(Vi); if (currents) { ddt(Pi) -= (4. / 9) * Pi_ciperp * Div_parP(Vi); //(4. / 9) * Vi * B32 * Grad_par(Pi_ciperp / B32); - Field3D phiPi = phi + Pi; - mesh->communicate(phiPi); + Field3D phiPi = add_all(phi, Pi); ddt(Pi) -= (2. / 6) * Pi_ci * fci_curvature(phiPi);//Curlb_B * Grad(phiPi); ddt(Pi) += (2. / 9) * bracket(Pi_ci, phi + Pi, BRACKET_ARAKAWA) * bracket_factor; } @@ -3807,81 +3225,45 @@ int Hermes::rhs(BoutReal t) { if (parallel_sheaths){ sheath_dpi = 0.0; - for (const auto &bndry_par : mesh->getBoundariesPar()) { + for (const auto &bndry_par : mesh->getBoundariesPar(BoundaryParType::xout)) { for (bndry_par->first(); !bndry_par->isDone(); bndry_par->next()) { int x = bndry_par->x; int y = bndry_par->y; int z = bndry_par->z; - if (bndry_par->dir == 1){ // forwards - // Temperature and density at the sheath entrance - BoutReal tisheath = floor( - 0.5 * (Ti(x, y, z) + Ti.yup()(x, y + bndry_par->dir, z)), - 0.0); - BoutReal tesheath = floor( - 0.5 * (Te(x, y, z) + Te.yup()(x, y + bndry_par->dir, z)), - 0.0); - BoutReal nesheath = floor( - 0.5 * (Ne(x, y, z) + Ne.yup()(x, y + bndry_par->dir, z)), - 0.0); - BoutReal visheath = floor( - 0.5 * (Vi(x, y, z) + Vi.yup()(x, y + bndry_par->dir, z)), - 0.0); - - // BoutReal tesheath = floor(Te(x,y,z),0.); - // BoutReal tisheath = floor(Te(x,y,z),0.); - // BoutReal nesheath = floor(Ne(x,y,z), 0.); - // Sound speed (normalised units) - BoutReal Cs = sqrt(tesheath + tisheath); - - // Heat flux - BoutReal q = (sheath_gamma_i - 1.5) * tisheath * nesheath * visheath; - - // Multiply by cell area to get power - BoutReal flux = - q - * (coord->J(x, y, z) + coord->J.yup()(x, y + bndry_par->dir, z)) - / (sqrt(coord->g_22(x, y, z)) - + sqrt(coord->g_22.yup()(x, y + bndry_par->dir, z))); - - // Divide by volume of cell, and 2/3 to get pressure - BoutReal power = - flux - / (coord->dy(x, y, z) * coord->J(x, y, z)); - sheath_dpi(x, y, z) -= (3. / 2) * power; - - } else { // backwards - // // Temperature and density at the sheath entrance - BoutReal tisheath = floor( - 0.5 * (Ti(x, y, z) + Ti.ydown()(x, y + bndry_par->dir, z)), - 0.0); - BoutReal tesheath = floor( - 0.5 * (Te(x, y, z) + Te.ydown()(x, y + bndry_par->dir, z)), - 0.0); - BoutReal nesheath = floor( - 0.5 * (Ne(x, y, z) + Ne.ydown()(x, y + bndry_par->dir, z)), - 0.0); - BoutReal visheath = floor( - 0.5 * (Vi(x, y, z) + Vi.ydown()(x, y + bndry_par->dir, z)), - 0.0); - // BoutReal tesheath = floor(Te(x,y,z),0.); - // BoutReal tisheath = floor(Te(x,y,z),0.); - // BoutReal nesheath = floor(Ne(x,y,z), 0.); - - // Sound speed (normalised units) - BoutReal Cs = -sqrt(tesheath + tisheath); - - // Heat flux - BoutReal q = (sheath_gamma_i - 1.5) * tisheath * nesheath * visheath; - - // Multiply by cell area to get power - BoutReal flux = q * (coord->J(x, y, z) + coord->J.ydown()(x, y + bndry_par->dir, z)) - / (sqrt(coord->g_22(x, y, z)) + sqrt(coord->g_22.ydown()(x, y + bndry_par->dir, z))); - - // Divide by volume of cell, and 2/3 to get pressure - BoutReal power = - flux - / (coord->dy(x, y, z) * coord->J(x, y, z)); - sheath_dpi(x, y, z) -= (3. / 2) * power; - - } + // Temperature and density at the sheath entrance + BoutReal tisheath = floor( + 0.5 * (Ti(x, y, z) + Ti.ynext(bndry_par->dir)(x, y + bndry_par->dir, z)), + 0.0); + BoutReal tesheath = floor( + 0.5 * (Te(x, y, z) + Te.ynext(bndry_par->dir)(x, y + bndry_par->dir, z)), + 0.0); + BoutReal nesheath = floor( + 0.5 * (Ne(x, y, z) + Ne.ynext(bndry_par->dir)(x, y + bndry_par->dir, z)), + 0.0); + BoutReal visheath = floor( + 0.5 * (Vi(x, y, z) + Vi.ynext(bndry_par->dir)(x, y + bndry_par->dir, z)), + 0.0); + + // BoutReal tesheath = floor(Te(x,y,z),0.); + // BoutReal tisheath = floor(Te(x,y,z),0.); + // BoutReal nesheath = floor(Ne(x,y,z), 0.); + // Sound speed (normalised units) + BoutReal Cs = bndry_par->dir * sqrt(tesheath + tisheath); + + // Heat flux + BoutReal q = (sheath_gamma_i - 1.5) * tisheath * nesheath * visheath; + + // Multiply by cell area to get power + BoutReal flux = + q + * (coord->J(x, y, z) + coord->J.ynext(bndry_par->dir)(x, y + bndry_par->dir, z)) + / (sqrt(coord->g_22(x, y, z)) + + sqrt(coord->g_22.ynext(bndry_par->dir)(x, y + bndry_par->dir, z))); + + // Divide by volume of cell, and 2/3 to get pressure + BoutReal power = + flux + / (coord->dy(x, y, z) * coord->J(x, y, z)); + sheath_dpi(x, y, z) -= (3. / 2) * power; + } } ddt(Pi) += sheath_dpi; @@ -3918,7 +3300,7 @@ int Hermes::rhs(BoutReal t) { if (energy_source) { // Add the same amount of energy to each particle - PiSource = Spi * Nelim / DC(Nelim); + PiSource = Spi * Ne / DC(Ne); } else { PiSource = Spi; } @@ -3929,7 +3311,7 @@ int Hermes::rhs(BoutReal t) { if (energy_source) { // Add the same amount of energy to each particle - PiSource = Spi * Nelim / DC(Nelim); + PiSource = Spi * Ne / DC(Ne); } else { PiSource = Spi * where(Spi, PiTarget, Pi); } @@ -3944,7 +3326,7 @@ int Hermes::rhs(BoutReal t) { if (energy_source) { // Add the same amount of energy to each particle - PiSource = Spi * Nelim / DC(Nelim); + PiSource = Spi * Ne / DC(Ne); if (source_vary_g11) { PiSource *= g11norm; @@ -4102,8 +3484,10 @@ int Hermes::rhs(BoutReal t) { if (neutral_friction) { // Vorticity if (boussinesq) { - ddt(Vort) -= FV::Div_a_Laplace_perp( - neutrals->Fperp / (Nelim * SQ(coord->Bxy)), phi); + Field3D tmp = neutrals->Fperp / (Ne * SQ(coord->Bxy)); + mesh->communicate(tmp); + tmp.applyParallelBoundary("parallel_neumann"); + ddt(Vort) -= FV::Div_a_Laplace_perp(tmp, phi); } else { ddt(Vort) -= FV::Div_a_Laplace_perp(neutrals->Fperp / SQ(coord->Bxy), phi); @@ -4258,13 +3642,14 @@ int Hermes::rhs(BoutReal t) { if (sinks) { // Sink terms for 2D simulations - // Field3D nsink = 0.5*Ne*sqrt(Telim)*sink_invlpar; // n C_s/ (2L) // + // Field3D nsink = 0.5*Ne*sqrt(Te)*sink_invlpar; // n C_s/ (2L) // // Sound speed flow to targets Field3D nsink = 0.5 * sqrt(Ti) * Ne * sink_invlpar; nsink = floor(nsink, 0.0); ddt(Ne) -= nsink; + ASSERT0(not fci_transform); Field3D conduct = (2. / 3) * kappa_epar * Te * SQ(sink_invlpar); conduct = floor(conduct, 0.0); ddt(Pe) -= conduct // Heat conduction @@ -4275,8 +3660,8 @@ int Hermes::rhs(BoutReal t) { /////////////////////////// // Sheath dissipation closure - Field3D phi_te = floor(phi / Telim, 0.0); - Field3D jsheath = Nelim * sqrt(Telim) * + Field3D phi_te = floor(phi / Te, 0.0); + Field3D jsheath = Ne * sqrt(Te) * (1 - (sqrt(mi_me) / (2. * sqrt(PI))) * exp(-phi_te)); ddt(Vort) += jsheath * sink_invlpar; @@ -4293,8 +3678,8 @@ int Hermes::rhs(BoutReal t) { // Include a drift-wave closure in the core region // as in Hasegawa-Wakatani and SOLT models - Field2D Tedc = DC(Telim); // Zonal average - Field3D Q = phi - Telim * log(Nelim); + Field2D Tedc = DC(Te); // Zonal average + Field3D Q = phi - Te * log(Ne); // Drift-wave operator Field3D Adw = alpha_dw * pow(Tedc, 1.5) * (Q - DC(Q)); @@ -4331,7 +3716,7 @@ int Hermes::rhs(BoutReal t) { } return 0; -} +} // rhs? /*! * Preconditioner. Solves the heat conduction @@ -4377,11 +3762,15 @@ const Field3D Hermes::Div_parP(const Field3D &f) { auto* coords = mesh->getCoordinates(); Field3D result; result.allocate(); - for(auto &i : f.getRegion("RGN_NOBNDRY")) { + const auto fup = f.yup(); + const auto fdown = f.ydown(); + BOUT_FOR(i, f.getRegion("RGN_NOBNDRY")) { + // for(auto &i : f.getRegion("RGN_NOBNDRY")) { auto yp = i.yp(); auto ym = i.ym(); - result[i] = (f.yup()[yp] / coords->Bxy.yup()[yp] - f.ydown()[ym] / coords->Bxy.ydown()[ym]) - * coords->Bxy[i] / (coords->dy[i] * sqrt(coords->g_22[i])); + result[i] = (fup[yp] / coords->Bxy.yup()[yp] - + fdown[ym] / coords->Bxy.ydown()[ym]) * + coords->Bxy[i] / (coords->dy[i] * sqrt(coords->g_22[i])); } return result; //return Div_par(f) + 0.5*beta_e*coord->Bxy*bracket(psi, f/coord->Bxy, BRACKET_ARAKAWA); diff --git a/hermes-2.hxx b/hermes-2.hxx index 6fe15da..4b93b1c 100644 --- a/hermes-2.hxx +++ b/hermes-2.hxx @@ -276,14 +276,12 @@ private: // Mesh quantities - Coordinates::FieldMetric B32, sqrtB; + Coordinates::FieldMetric B12, B32, B42; bool fci_transform; Field3D Bxyz, logB; Field3D bracket_factor; const Field3D fci_curvature(const Field3D &f); - - Field3D a,b,c,d,f; //Debugging variables }; /// Fundamental constants diff --git a/makefile b/makefile index 361d7a1..69b794c 100644 --- a/makefile +++ b/makefile @@ -1,10 +1,39 @@ -BOUT_TOP = ../.. +BOUT_TOP ?= ../.. TARGET = hermes-2 DIRS = atomicpp -SOURCEC = hermes-2.cxx div_ops.cxx loadmetric.cxx radiation.cxx neutral-model.cxx diffusion2d.cxx recycling.cxx full-velocity.cxx mixed.cxx +SOURCEC = hermes-2.cxx div_ops.cxx loadmetric.cxx radiation.cxx neutral-model.cxx \ + diffusion2d.cxx recycling.cxx full-velocity.cxx mixed.cxx \ + atomicpp/ImpuritySpecies.cxx atomicpp/Prad.cxx \ + atomicpp/RateCoefficient.cxx atomicpp/sharedFunctions.cxx + +OBJ = $(SOURCEC:%.cxx=%.o) +TARGET ?= $(SOURCEC:%.cxx=%) + +ifeq (, $(shell which bout-config)) +$(error No bout-config in $$PATH ($(PATH))) +endif + +# Use the bout-config script to get compiler, flags +CXX:=$(shell bout-config --cxx) + +CFLAGS:=$(shell bout-config --cflags) +LD:=$(shell bout-config --ld) +LDFLAGS:=$(shell bout-config --libs) + +$(TARGET): makefile $(OBJ) + @echo " Linking" $(TARGET) + $(LD) -o $(TARGET) $(OBJ) $(LDFLAGS) -lbout++ # -lfmt + +%.o: %.cxx + @echo " Compiling " $< + @echo "$<" + @$(CXX) $(CFLAGS) -c $< -o $@ + +.PHONY: clean +clean: + rm -f $(OBJ) $(TARGET) -include $(BOUT_TOP)/make.config diff --git a/neutral-model.cxx b/neutral-model.cxx index 94d04f7..d04090e 100644 --- a/neutral-model.cxx +++ b/neutral-model.cxx @@ -16,6 +16,7 @@ NeutralModel *NeutralModel::create(Solver *solver, Mesh *mesh, Options &options) { // Decide which neutral model to use std::string type = options["type"].withDefault("none"); + options["bulk"].setConditionallyUsed(); if (type == "none") { // Neutral model which does nothing diff --git a/radiation.cxx b/radiation.cxx index acc0077..85921ac 100644 --- a/radiation.cxx +++ b/radiation.cxx @@ -1,4 +1,3 @@ - #include #include // for mesh object #include @@ -211,17 +210,21 @@ BoutReal HydrogenRadiatedPower::excitation(BoutReal Te) { // Collision rate coefficient [m3/s] BoutReal UpdatedRadiatedPower::ionisation(BoutReal T) { double fION; // Rate coefficient - double TT; - TT = T; + const double TT = T; - double ioncoeffs[9] = {-3.271397E1, 1.353656E1, -5.739329, - 1.563155, -2.877056E-1, 3.482560e-2, - -2.631976E-3, 1.119544E-4, -2.039150E-6}; + const double ioncoeffs[9] = {-3.271397E1, 1.353656E1, -5.739329, + 1.563155, -2.877056E-1, 3.482560e-2, + -2.631976E-3, 1.119544E-4, -2.039150E-6}; + + const double logT = log(TT); double lograte = 0.0; - for (int i = 0; i <= 8; i++) { - lograte = lograte + ioncoeffs[i] * pow(log(TT), i); + double powT = 1; + for (int i = 0; i < 9; i++) { + // lograte += ioncoeffs[i]*pow(logT,i); + lograte += ioncoeffs[i] * powT; + powT *= logT; } fION = exp(lograte) * 1.0E-6; @@ -231,95 +234,109 @@ BoutReal UpdatedRadiatedPower::ionisation(BoutReal T) { // [m3/s] BoutReal UpdatedRadiatedPower::recombination(BoutReal n, BoutReal T) { - double TT, RDNE, RTE, DNE, E, RN, RT, RNJ, RTI, suma, fHAV, fRAD; - int i, j, i1, j1; + // double TT,RDNE,RTE,DNE,E,RN,RT,RNJ,RTI,suma,fHAV,fRAD; + // int i,j,i1,j1; if (n < 1e3) // Log(n) used, so prevent NaNs return 0.0; - - double MATA[9][9] = { - { - -2.855728479302E+01, -7.664042607917E-01, -4.930424003280E-03, - -5.386830982777E-03, -1.626039237665E-04, 6.080907650243E-06, - 2.101102051942E-05, -2.770717597683E-06, 1.038235939800E-07, - }, - { - 3.488563234375E-02, -3.583233366133E-03, -3.620245352252E-03, - -9.532840484460E-04, 1.888048628708E-04, -1.014890683861E-05, - 2.245676563601E-05, -4.695982369246E-06, 2.523166611507E-07, - }, - { - -2.799644392058E-02, -7.452514292790E-03, 6.958711963182E-03, - 4.631753807534E-04, 1.288577690147E-04, -1.145028889459E-04, - -2.245624273814E-06, 3.250878872873E-06, -2.145390398476E-07, - }, - { - 1.209545317879E-02, 2.709299760454E-03, -2.139257298118E-03, - -5.371179699661E-04, -1.634580516353E-05, 5.942193980802E-05, - -2.944873763540E-06, -9.387290785993E-07, 7.381435237585E-08, - }, - { - -2.436630799820E-03, -7.745129766167E-04, 4.603883706734E-04, - 1.543350502150E-04, -9.601036952725E-06, -1.211851723717E-05, - 1.002105099354E-06, 1.392391630459E-07, -1.299713684966E-08, - }, - { - 2.837893719800E-04, 1.142444698207E-04, -5.991636837395E-05, - -2.257565836876E-05, 3.425262385387E-06, 1.118965496365E-06, - -1.291320799814E-07, -1.139093288575E-08, 1.265189576423E-09, - }, - { - -1.886511169084E-05, -9.382783518064E-06, 4.729262545726E-06, - 1.730782954588E-06, -4.077019941998E-07, -4.275321573501E-08, - 7.786155463269E-09, 5.178505597480E-10, -6.854203970018E-11, - }, - { - 6.752155602894E-07, 3.902800099653E-07, -1.993485395689E-07, - -6.618240780594E-08, 2.042041097083E-08, 3.708616111085E-10, - -2.441127783437E-10, -9.452402157390E-12, 1.836615031798E-12, - }, - { - -1.005893858779E-08, -6.387411585521E-09, 3.352589865190E-09, - 1.013364275013E-09, -3.707977721109E-10, 7.068450112690E-12, - 3.773208484020E-12, -4.672724022059E-14, -1.640492364811E-14, - }, - }; - - RDNE = n; - RTE = T; - - DNE = RDNE; - TT = RTE; - E = DNE * 1.0E-14; - - RN = log(E); - RT = log(TT); - - suma = 0.0; - for (i = 1; i <= 9; i++) { - i1 = i - 1; - for (j = 1; j <= 9; j++) { - j1 = j - 1; - RNJ = pow(RN, j1); - if ((RN == 0.0) && (j1 == 0)) - RNJ = 1.0; - RTI = pow(RT, i1); - if ((RT == 0.0) && (i1 == 0)) - RTI = 1.0; - suma = suma + MATA[j - 1][i - 1] * RNJ * RTI; + // use transpose, its ~ 10% faster + /*const double MATA[9][9]= { + {-2.855728479302E+01, -7.664042607917E-01, -4.930424003280E-03, + -5.386830982777E-03, -1.626039237665E-04, 6.080907650243E-06, + 2.101102051942E-05, -2.770717597683E-06, 1.038235939800E-07,}, + { 3.488563234375E-02, -3.583233366133E-03, -3.620245352252E-03, + -9.532840484460E-04, 1.888048628708E-04, -1.014890683861E-05, + 2.245676563601E-05, -4.695982369246E-06, 2.523166611507E-07,}, + {-2.799644392058E-02, -7.452514292790E-03, 6.958711963182E-03, + 4.631753807534E-04, 1.288577690147E-04, -1.145028889459E-04, + -2.245624273814E-06, 3.250878872873E-06, -2.145390398476E-07,}, + { 1.209545317879E-02, 2.709299760454E-03, -2.139257298118E-03, + -5.371179699661E-04, -1.634580516353E-05, 5.942193980802E-05, + -2.944873763540E-06, -9.387290785993E-07, 7.381435237585E-08,}, + {-2.436630799820E-03, -7.745129766167E-04, 4.603883706734E-04, + 1.543350502150E-04, -9.601036952725E-06, -1.211851723717E-05, + 1.002105099354E-06, 1.392391630459E-07, -1.299713684966E-08,}, + { 2.837893719800E-04, 1.142444698207E-04, -5.991636837395E-05, + -2.257565836876E-05, 3.425262385387E-06, 1.118965496365E-06, + -1.291320799814E-07, -1.139093288575E-08, 1.265189576423E-09,}, + {-1.886511169084E-05, -9.382783518064E-06, 4.729262545726E-06, + 1.730782954588E-06, -4.077019941998E-07, -4.275321573501E-08, + 7.786155463269E-09, 5.178505597480E-10, -6.854203970018E-11,}, + { 6.752155602894E-07, 3.902800099653E-07, -1.993485395689E-07, + -6.618240780594E-08, 2.042041097083E-08, 3.708616111085E-10, + -2.441127783437E-10, -9.452402157390E-12, 1.836615031798E-12,}, + {-1.005893858779E-08, -6.387411585521E-09, 3.352589865190E-09, + 1.013364275013E-09, -3.707977721109E-10, 7.068450112690E-12, + 3.773208484020E-12, -4.672724022059E-14, -1.640492364811E-14,}, + };*/ + const double MATA[9][9] = { + {-28.55728479302, 0.03488563234375, -0.02799644392058, 0.01209545317879, + -0.00243663079982, 0.00028378937198, -1.886511169084E-05, 6.752155602894E-07, + -1.005893858779E-08}, + {-0.7664042607917, -0.00358323336613, -0.00745251429279, 0.00270929976045, + -0.00077451297662, 0.00011424446982, -9.382783518064E-06, 3.902800099653E-07, + -6.387411585521E-09}, + {-0.00493042400328, -0.00362024535225, 0.00695871196318, -0.00213925729812, + 0.00046038837067, -5.991636837395E-05, 4.729262545726E-06, -1.993485395689E-07, + 3.35258986519E-09}, + {-0.00538683098278, -0.00095328404845, 0.00046317538075, -0.00053711796997, + 0.00015433505022, -2.257565836876E-05, 1.730782954588E-06, -6.618240780594E-08, + 1.013364275013E-09}, + {-0.00016260392377, 0.00018880486287, 0.00012885776901, -1.634580516353E-05, + -9.601036952725E-06, 3.425262385387E-06, -4.077019941998E-07, 2.042041097083E-08, + -3.707977721109E-10}, + {6.080907650243E-06, -1.014890683861E-05, -0.00011450288895, 5.942193980802E-05, + -1.211851723717E-05, 1.118965496365E-06, -4.275321573501E-08, 3.708616111085E-10, + 7.06845011269E-12}, + {2.101102051942E-05, 2.245676563601E-05, -2.245624273814E-06, -2.94487376354E-06, + 1.002105099354E-06, -1.291320799814E-07, 7.786155463269E-09, -2.441127783437E-10, + 3.77320848402E-12}, + {-2.770717597683E-06, -4.695982369246E-06, 3.250878872873E-06, -9.387290785993E-07, + 1.392391630459E-07, -1.139093288575E-08, 5.17850559748E-10, -9.45240215739E-12, + -4.672724022059E-14}, + {1.0382359398E-07, 2.523166611507E-07, -2.145390398476E-07, 7.381435237585E-08, + -1.299713684966E-08, 1.265189576423E-09, -6.854203970018E-11, 1.836615031798E-12, + -1.640492364811E-14}}; + + const double RN = log(n * 1.0E-14); + const double RT = log(T); + + double suma = 0.0; + // double RNJ[9]; + // for (int j=0;j<9;j++){ + // //if ((RN==0.0) && (j==0)) + // //RNJ[j]=1.0; + // } + { + double RTI = 1; + for (int i = 0; i < 9; i++) { + double RNJ = 1; + double sumb = 0; + for (int j = 0; j < 9; j++) { + sumb += MATA[i][j] * RNJ; + RNJ *= RN; + } + suma = suma + (sumb * RTI); + RTI *= RT; } } - fHAV = exp(suma) * 1.0E-6 / (1.0 + 0.125 * TT); + double fHAV = exp(suma) * 1.0E-6 / (1.0 + 0.125 * T); - double A = 3.92E-20; - double B = 3.0E-124 * pow(1.6E-19, -4.5); - double Ry = 13.60569; - double chi = 0.35; + //const double A = 3.92E-20; + // const double B = 3.0E-124 * pow(1.6E-19,-4.5); + const double B = 1.144409179687500325061682019168e-39; + const double Ry = 13.60569; + // const double Ry15 = pow(Ry,1.5); + //const double Ry15 = 50.18580066442199694165537948720; + const double chi = 0.35; - fRAD = A * pow(Ry, 1.5) * 1.0 / (sqrt(TT) * (Ry + chi * TT)); + const double ARy15 = 1.9672833860453424198e-18; - double fREC = fHAV + fRAD + B * DNE * pow(TT, -5.0); + double fRAD = ARy15 / (sqrt(T) * (Ry + chi * T)); + + // double fREC = fHAV + fRAD + B*n*pow(T,-5.0); + double fREC = fHAV + fRAD + B * n / (T * T * T * T * T); return fREC; } @@ -327,69 +344,72 @@ BoutReal UpdatedRadiatedPower::recombination(BoutReal n, BoutReal T) { // [m3/s] BoutReal UpdatedRadiatedPower::chargeExchange(BoutReal T) { double fCX; - double TT; - - TT = T; - - double E = 10.0; - - double cxcoeffs[9][9] = { - { - -1.829079582E1, 1.640252721E-1, 3.364564509E-2, 9.530225559E-3, - -8.519413900E-4, -1.247583861E-3, 3.014307546E-4, -2.499323170E-5, - 6.932627238E-7, - }, - { - 2.169137616E-1, -1.106722014E-1, -1.382158680E-3, 7.348786287E-3, - -6.343059502E-4, -1.919569450E-4, 4.075019352E-5, -2.850044983E-6, - 6.966822400E-8, - }, - { - 4.307131244E-2, 8.948693625E-3, -1.209480567E-2, -3.675019470E-4, - 1.039643391E-3, -1.553840718E-4, 2.670827249E-6, 7.695300598E-7, - -3.783302282E-8, - }, - { - -5.754895093E-4, 6.062141761E-3, 1.075907882E-3, -8.119301728E-4, - 8.911036876E-6, 3.175388950E-5, -4.515123642E-6, 2.187439284E-7, - -2.911233952E-9, - }, - { - -1.552077120E-3, -1.210431588E-3, 8.297212634E-4, 1.361661817E-4, - -1.008928628E-4, 1.080693990E-5, 5.106059414E-7, -1.299275586E-7, - 5.117133050E-9, - }, - { - -1.876800283E-4, -4.052878752E-5, -1.907025663E-4, 1.141663042E-5, - 1.775681984E-5, -3.149286924E-6, 3.105491555E-8, 2.274394089E-8, - -1.130988251E-9, - }, - { - 1.125490271E-4, 2.875900436E-5, 1.338839629E-5, -4.340802793E-6, - -7.003521917E-7, 2.318308730E-7, -6.030983538E-9, -1.755944926E-9, - 1.005189187E-10, - }, - { - -1.238982763E-5, -2.616998140E-6, -1.171762874E-7, 3.517971869E-7, - -4.928692833E-8, 1.756388999E-10, -1.446756796E-10, 7.143183138E-11, - -3.989884106E-12, - }, - { - 4.163596197E-7, 7.558092849E-8, -1.328404104E-8, -9.170850254E-9, - 3.208853884E-9, -3.952740759E-10, 2.739558476E-11, -1.693040209E-12, - 6.388219930E-14, - }, - }; - - double lograte = 0.0; - for (int i = 0; i <= 8; i++) { - for (int j = 0; j <= 8; j++) { - lograte = lograte + cxcoeffs[i][j] * pow(log(TT), i) * pow(log(E), j); + if (T < 0.14660) { // avoid unphysical part of function + fCX = 2.02022322361958e-14; + } else { + + const double logT = log(T); + + // double E = 10.0; // used for computing the vector from the commented matrix + + /*const double cxcoeffs[9][9] = { + {-1.829079582E1, 1.640252721E-1, 3.364564509E-2,\ + 9.530225559E-3, -8.519413900E-4, -1.247583861E-3,\ + 3.014307546E-4, -2.499323170E-5, 6.932627238E-7,}, + { 2.169137616E-1, -1.106722014E-1, -1.382158680E-3,\ + 7.348786287E-3, -6.343059502E-4, -1.919569450E-4,\ + 4.075019352E-5, -2.850044983E-6, 6.966822400E-8,}, + { 4.307131244E-2, 8.948693625E-3, -1.209480567E-2,\ + -3.675019470E-4, 1.039643391E-3, -1.553840718E-4,\ + 2.670827249E-6, 7.695300598E-7, -3.783302282E-8,}, + {-5.754895093E-4, 6.062141761E-3, 1.075907882E-3,\ + -8.119301728E-4, 8.911036876E-6, 3.175388950E-5,\ + -4.515123642E-6, 2.187439284E-7, -2.911233952E-9,}, + {-1.552077120E-3, -1.210431588E-3, 8.297212634E-4,\ + 1.361661817E-4, -1.008928628E-4, 1.080693990E-5,\ + 5.106059414E-7, -1.299275586E-7, 5.117133050E-9,}, + {-1.876800283E-4, -4.052878752E-5, -1.907025663E-4,\ + 1.141663042E-5, 1.775681984E-5, -3.149286924E-6,\ + 3.105491555E-8, 2.274394089E-8, -1.130988251E-9,}, + { 1.125490271E-4, 2.875900436E-5, 1.338839629E-5,\ + -4.340802793E-6, -7.003521917E-7, 2.318308730E-7,\ + -6.030983538E-9, -1.755944926E-9, 1.005189187E-10,}, + {-1.238982763E-5, -2.616998140E-6, -1.171762874E-7,\ + 3.517971869E-7, -4.928692833E-8, 1.756388999E-10,\ + -1.446756796E-10, 7.143183138E-11, -3.989884106E-12,}, + { 4.163596197E-7, 7.558092849E-8, -1.328404104E-8,\ + -9.170850254E-9, 3.208853884E-9, -3.952740759E-10,\ + 2.739558476E-11, -1.693040209E-12, 6.388219930E-14,}, + };*/ + const double cxcoeffs[9] = {-1.76861862430790530709e+01, 1.93633165585571261069e-02, + 1.48638830739889656052e-02, 1.08809966343519897575e-02, + -3.78840793338734954480e-04, -8.45866042005527008808e-04, + 1.90656650772697163192e-04, -1.61164747567831752339e-05, + 4.76171332457707722191e-07}; + + double lograte = 0.0; + double plogT=1; + for (int i = 0; i < 9; i++) { + /*lograte = 0.0; + for (int j=0;j<=8;j++) + { + lograte = lograte + cxcoeffs[i][j] * pow(log(E),j); + } + printf("%.20e,\t",lograte);*/ + // pow(log(TT),i); + lograte += cxcoeffs[i] * plogT; + plogT*=logT; + } + // throw BoutException("I'm done ...\n"); + // std::cerr << lograte << "\t"; + BoutReal logratelim = 650; + if (lograte > logratelim) { + lograte = logratelim; + throw BoutException("UpdatedRadiatedPower::chargeExchange limiting lograte!\n"); } + fCX = 1.0E-6 * exp(lograte); } - fCX = 1.0E-6 * exp(lograte); - return fCX; } diff --git a/rotating-ellipse-tiny.fci.nc b/rotating-ellipse-tiny.fci.nc new file mode 100644 index 0000000..fcf3663 Binary files /dev/null and b/rotating-ellipse-tiny.fci.nc differ