Merge pull request #38 from prkkumar/exchange_bc

new Laplacian for exchange BC

Examples/Tests/Macroscopic_Maxwell/PSSW_halfCos_exci/analysis_getline.bash

@@ -0,0 +1,11 @@
+#!/bin/bash
+FILE="$1"
+LINE_NO=$2
+i=0
+while read line; do
+  i=$(( i + 1 ))
+  case $i in $LINE_NO) echo "$line"; break;; esac
+done <"$FILE"
+
+### how to run this script: 
+### $ bash squashTime.bash FILE LINE_NO

Examples/Tests/Macroscopic_Maxwell/PSSW_halfCos_exci/analysis_readfield.bash

@@ -0,0 +1,6 @@
+#!/bin/bash
+for timestep in $(seq 0 500 5000000 ); do
+    newnumber='000000000'${timestep}      # get number, pack with zeros
+    newnumber=${newnumber:(-9)}       # the last seven characters
+    ./WritePlotfileToASCII3d.gnu.ex infile=diags/plt${newnumber} | tee raw_data/${newnumber}.txt
+done

Examples/Tests/Macroscopic_Maxwell/PSSW_halfCos_exci/analysis_squashTime.bash

@@ -0,0 +1,6 @@
+#!/bin/bash
+for timestep in $(seq 0 500 5000000); do
+    newnumber='000000000'${timestep}      # get number, pack with zeros
+    newnumber=${newnumber:(-9)}       # the last seven characters
+    bash getLine.bash raw_data/${newnumber}.txt 128 | tee -a Mx_xface_center_128.txt
+done

Examples/Tests/Macroscopic_Maxwell/PSSW_halfCos_exci/inputs_LLG_PSSW_asymm_exci

@@ -0,0 +1,129 @@
+# set up to duplicate perpendicular standing spin wave for exchange boundary condition demo
+# max step
+max_step = 5000000
+
+# number of grid points
+amr.n_cell = n_cellx n_celly n_cellz
+
+# Maximum allowable size of each subdomain
+amr.max_grid_size = grid_number
+amr.blocking_factor = 8
+
+amr.max_level = 0
+
+# Geometry
+geometry.coord_sys   = 0                     # 0: Cartesian
+geometry.prob_lo     =  -width/2 -length/2 -thickness
+geometry.prob_hi     =   width/2  length/2  thickness
+
+# Boundary condition
+boundary.field_lo = periodic periodic pec
+boundary.field_hi = periodic periodic pec
+warpx.serialize_ics = 1
+
+# Verbosity
+warpx.verbose = 1
+warpx.use_filter = 0
+warpx.cfl = 1000.0 # We can use a large cfl number here because the Maxwell equations are not included.
+warpx.mag_time_scheme_order = 2 # default 1
+warpx.mag_M_normalization = 1
+warpx.mag_LLG_coupling = 0
+warpx.mag_LLG_exchange_coupling = 1
+warpx.mag_LLG_anisotropy_coupling = 0
+
+# Algorithms
+algo.current_deposition = esirkepov
+algo.em_solver_medium = macroscopic # vacuum/macroscopic
+
+my_constants.thickness = 345.0e-9
+my_constants.width = 6900.0e-9
+my_constants.length = 6900.0e-9
+my_constants.n_cellx = 8
+my_constants.n_celly = 8
+my_constants.n_cellz = 256
+my_constants.grid_number = 256
+
+my_constants.epsilon_r = 13.0
+my_constants.pi = 3.14159265359
+my_constants.mu0 = 1.25663706212e-6
+my_constants.epsilon0 = 8.84e-12
+my_constants.Ms_ga = 1750.
+
+my_constants.flag_none = 0
+my_constants.flag_hs = 1
+my_constants.flag_ss = 2
+
+algo.macroscopic_sigma_method = laxwendroff # laxwendroff or backwardeuler
+macroscopic.sigma_function(x,y,z) = "0.0"
+macroscopic.epsilon_function(x,y,z) = "epsilon0 * epsilon_r"
+macroscopic.mu_function(x,y,z) = "mu0"
+
+#unit conversion: 1 Gauss = (1000/4pi) A/m
+macroscopic.mag_Ms_init_style = "parse_mag_Ms_function" # parse or "constant"
+macroscopic.mag_Ms_function(x,y,z) = "(z < thickness/2 + 1e-12) * (z > -thickness/2 - 1e-12) * Ms_ga*1000/4/pi" # in unit A/m
+
+macroscopic.mag_alpha_init_style = "parse_mag_alpha_function" # parse or "constant"
+macroscopic.mag_alpha_function(x,y,z) = "(z < thickness/2 + 1e-12) * (z > -thickness/2 - 1e-12) * 0.005" # alpha is unitless, typical values range from 1e-3 ~ 1e-5
+
+macroscopic.mag_gamma_init_style = "parse_mag_gamma_function" # parse or "constant"
+macroscopic.mag_gamma_function(x,y,z) = " -1.759e11 " # gyromagnetic ratio is constant for electrons in all materials
+
+macroscopic.mag_exchange_init_style = "parse_mag_exchange_function" # parse or "constant"
+macroscopic.mag_exchange_function(x,y,z) = "(z < thickness/2 + 1e-12) * (z > -thickness/2 - 1e-12) * 3.1e-12" # YIG
+
+#macroscopic.mag_anisotropy_init_style = "parse_mag_anisotropy_function" # parse or "constant"
+#macroscopic.mag_anisotropy_function(x,y,z) = "0.0" # Ku = mu0 * Ms * H_anis / 2.0, all fields in SI units, H_anis = 20 Oe
+
+macroscopic.mag_max_iter = 100 # maximum number of M iteration in each time step
+macroscopic.mag_tol = 1.e-6 # M magnitude relative error tolerance compared to previous iteration
+macroscopic.mag_normalized_error = 0.1 # if M magnitude relatively changes more than this value, raise a red flag
+#macroscopic.mag_LLG_anisotropy_axis = 0.0 1.0 0.0
+
+#################################
+############ FIELDS #############
+#################################
+my_constants.Hy_bias_oe = 300 # in Oe
+my_constants.Hx_rf_oe = 3 # in Oe
+my_constants.c = 299792458.
+my_constants.frequency = 2.5e9 # frequency of input microwave H
+my_constants.TP = 1/frequency
+my_constants.flag_none = 0
+my_constants.flag_hs = 1
+my_constants.flag_ss = 2
+
+#warpx.E_excitation_on_grid_style = "parse_E_excitation_grid_function"
+#warpx.Ex_excitation_grid_function(x,y,z,t) = "0.0"
+#warpx.Ey_excitation_grid_function(x,y,z,t) = "0.0"
+#warpx.Ez_excitation_grid_function(x,y,z,t) = "0.0"
+#warpx.Ex_excitation_flag_function(x,y,z) = "flag_none"
+# warpx.Ey_excitation_flag_function(x,y,z) = "flag_ss * (z > -domain_size/2.0 + thickness - domain_size/grid_number/2.0) * (z <= -domain_size/2.0 + thickness + domain_size/grid_number/2.0)"
+#warpx.Ey_excitation_flag_function(x,y,z) = "flag_none"
+#warpx.Ez_excitation_flag_function(x,y,z) = "flag_none"
+
+#unit conversion: 1 Gauss = 1 Oersted = (1000/4pi) A/m
+#calculation of H_bias: H_bias (oe) = frequency / 2.8e6
+
+warpx.H_bias_excitation_on_grid_style = parse_h_bias_excitation_grid_function
+#warpx.Hx_bias_excitation_grid_function(x,y,z,t) = "Hx_rf_oe * 1000/4/pi * (exp(-(t-3*TP)**2/(2*TP**2))*cos(2*pi*frequency*t))"
+warpx.Hx_bias_excitation_grid_function(x,y,z,t) = "Hx_rf_oe * 1000/4/pi * (exp(-(t-3*TP)**2/(2*TP**2))*cos(2*pi*frequency*t)) * cos(z/thickness * pi)"
+
+#warpx.Hx_bias_excitation_grid_function(x,y,z,t) = "0.0"
+warpx.Hy_bias_excitation_grid_function(x,y,z,t) = "Hy_bias_oe * 1000/4/pi"
+warpx.Hz_bias_excitation_grid_function(x,y,z,t) = "0.0"
+
+warpx.Hx_bias_excitation_flag_function(x,y,z) = "flag_hs"
+warpx.Hy_bias_excitation_flag_function(x,y,z) = "flag_hs"
+warpx.Hz_bias_excitation_flag_function(x,y,z) = "flag_none"
+
+warpx.M_ext_grid_init_style = parse_M_ext_grid_function
+warpx.Mx_external_grid_function(x,y,z)= "0.0"
+warpx.My_external_grid_function(x,y,z)= " (z < thickness/2 + 1e-12) * (z > -thickness/2 - 1e-12) * Ms_ga*1000/4/pi "
+warpx.Mz_external_grid_function(x,y,z)= "0.0"
+
+
+# Diagnostics
+diagnostics.diags_names = plt
+plt.intervals = 500
+plt.diag_type = Full
+plt.fields_to_plot = Mx_xface Mx_yface Mx_zface My_xface My_yface My_zface Mz_xface Mz_yface Mz_zface
+plt.file_min_digits = 9

Examples/Tests/Macroscopic_Maxwell/inputs_3d_LLG_exchangeBC

@@ -0,0 +1,88 @@
+################################
+####### GENERAL PARAMETERS ######
+#################################
+max_step = 6000
+amr.n_cell = 64 64 64 # number of cells spanning the domain in each coordinate direction at level 0
+amr.max_grid_size = 32 # maximum size of each AMReX box, used to decompose the domain
+amr.blocking_factor = 16
+geometry.coord_sys = 0
+
+geometry.prob_lo = -10e-9 -10e-9 -10.0e-9
+geometry.prob_hi =  10e-9  10e-9  10.0e-9
+boundary.field_lo = periodic periodic periodic
+boundary.field_hi = periodic periodic periodic
+amr.max_level = 0
+
+#################################
+############ NUMERICS ###########
+#################################
+warpx.verbose = 1
+warpx.use_filter = 0
+warpx.cfl = 1000
+warpx.mag_time_scheme_order = 2 # default 1
+warpx.mag_M_normalization = 1 # 1 is saturated
+warpx.mag_LLG_coupling = 0
+warpx.mag_LLG_exchange_coupling = 1
+
+algo.em_solver_medium = macroscopic # vacuum/macroscopic
+algo.macroscopic_sigma_method = laxwendroff # laxwendroff or backwardeuler
+
+macroscopic.sigma_function(x,y,z) = "0.0"
+macroscopic.epsilon_function(x,y,z) = "8.8541878128e-12"
+macroscopic.mu_function(x,y,z) = "1.25663706212e-06"
+
+my_constants.mag_lo_x = -5.e-9
+my_constants.mag_hi_x = 5.e-9
+my_constants.mag_lo_y = -5.e-9
+my_constants.mag_hi_y = 5.e-9
+my_constants.mag_lo_z = -5.e-9
+my_constants.mag_hi_z = 5.e-9
+
+#unit conversion: 1 Gauss = (1000/4pi) A/m
+macroscopic.mag_Ms_init_style = "parse_mag_Ms_function" # parse or "constant"
+macroscopic.mag_Ms_function(x,y,z) = "1.4e5 * (x >= mag_lo_x - 1e-12) * (x <= mag_hi_x + 1e-12) * (y >= mag_lo_y - 1e-12) * (y <= mag_hi_y + 1e-12) * (z >= mag_lo_z - 1e-12) * (z <= mag_hi_z + 1e-12)"  # in unit A/m, equal to 1750 Gauss; Ms must be nonzero for LLG
+
+macroscopic.mag_alpha_init_style = "parse_mag_alpha_function" # parse or "constant"
+macroscopic.mag_alpha_function(x,y,z) = "0.058" # * (x >= mag_lo_x - 1e-12) * (x <= mag_hi_x + 1e-12) * (y >= mag_lo_y - 1e-12) * (y <= mag_hi_y + 1e-12) * (z >= mag_lo_z - 1e-12) * (z <= mag_hi_z + 1e-12)"  # alpha is unitless, calculated from linewidth Delta_H = 40 Oersted
+
+macroscopic.mag_gamma_init_style = "parse_mag_gamma_function" # parse or "constant"
+macroscopic.mag_gamma_function(x,y,z) = "-1.759e11" # * (x >= mag_lo_x - 1e-12) * (x <= mag_hi_x + 1e-12) * (y >= mag_lo_y - 1e-12) * (y <= mag_hi_y + 1e-12) * (z >= mag_lo_z - 1e-12) * (z <= mag_hi_z + 1e-12)" # gyromagnetic ratio is constant for electrons in all materials
+
+macroscopic.mag_exchange_init_style = "parse_mag_exchange_function" # parse or "constant"
+macroscopic.mag_exchange_function(x,y,z) = "3.76e-12" # * (x >= mag_lo_x) * (x <= mag_hi_x) * (y >= mag_lo_y - 1e-12) * (y <= mag_hi_y + 1e-12) * (z >= mag_lo_z - 1e-12) * (z <= mag_hi_z + 1e-12)" # exchange coupling constanit; Must be non-zero when exchange coupling is ON
+
+macroscopic.mag_max_iter = 100 # maximum number of M iteration in each time step
+macroscopic.mag_tol = 1.e-6 # M magnitude relative error tolerance compared to previous iteration
+macroscopic.mag_normalized_error = 0.1 # if M magnitude relatively changes more than this value, raise a red flag
+
+#################################
+############ FIELDS #############
+#################################
+warpx.E_ext_grid_init_style = parse_E_ext_grid_function
+warpx.Ex_external_grid_function(x,y,z) = 0.
+warpx.Ey_external_grid_function(x,y,z) = 0.
+warpx.Ez_external_grid_function(x,y,z) = 0.
+
+warpx.H_ext_grid_init_style = parse_H_ext_grid_function
+warpx.Hx_external_grid_function(x,y,z)= 0.
+warpx.Hy_external_grid_function(x,y,z) = 0.
+warpx.Hz_external_grid_function(x,y,z) = 0.
+
+#unit conversion: 1 Gauss = 1 Oersted = (1000/4pi) A/m
+#calculation of H_bias: H_bias (oe) = frequency / 2.8e6
+warpx.H_bias_ext_grid_init_style = parse_H_bias_ext_grid_function
+warpx.Hx_bias_external_grid_function(x,y,z)= 0.
+warpx.Hy_bias_external_grid_function(x,y,z)= "3.7e4" # in A/m, equal to 464 Oersted
+warpx.Hz_bias_external_grid_function(x,y,z)= 0.
+
+warpx.M_ext_grid_init_style = parse_M_ext_grid_function
+warpx.Mx_external_grid_function(x,y,z)= "1.4e5 * (x >= mag_lo_x - 1e-12) * (x <= mag_hi_x + 1e-12) * (y >= mag_lo_y - 1e-12) * (y < 0.) * (z >= mag_lo_z - 1e-12) * (z <= mag_hi_z + 1e-12)"
+warpx.My_external_grid_function(x,y,z) = 0.
+warpx.Mz_external_grid_function(x,y,z) = "1.4e5 * (x >= mag_lo_x - 1e-12) * (x <= mag_hi_x + 1e-12) * (y >= 0.) * (y <= mag_hi_y + 1e-12) * (z >= mag_lo_z - 1e-12) * (z <= mag_hi_z + 1e-12)"
+
+#Diagnostics
+diagnostics.diags_names = plt
+plt.intervals = 100
+plt.diag_type = Full
+plt.fields_to_plot = Ex Ey Ez Hx Hy Hz Bx By Bz Mx_xface My_xface Mz_xface Mx_yface My_yface Mz_yface Mx_zface My_zface Mz_zface
+plt.plot_raw_fields = 0

Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CartesianYeeAlgorithm.H

@@ -15,6 +15,9 @@
 #include <AMReX_Array4.H>
 #include <AMReX_Gpu.H>
 #include <AMReX_REAL.H>
+#include "Utils/WarpXUtil.H"
+#include "WarpX.H"
+#include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H"
 
 #include <array>
 #include <cmath>
@@ -239,6 +242,109 @@ struct CartesianYeeAlgorithm {
 
 #endif
 
+#ifdef WARPX_MAG_LLG
+
+    /**
+     * Perform divergence of gradient along x on M field when exchange coupling is on */
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static amrex::Real LaplacianDx_Mag (
+        amrex::Array4<amrex::Real> const& F,
+        amrex::Real const * const coefs_x, int const n_coefs_x, amrex::Real const Ms_lo_x, amrex::Real const Ms_hi_x,
+        int const i, int const j, int const k, int const ncomp=0, int const nodality=0) {
+
+        amrex::Real const inv_dx = coefs_x[0];
+        if (nodality == 0){ // at x face (normal face). dM/dx = 0
+           if (Ms_hi_x == 0.){
+               return 0.5 * inv_dx * inv_dx * (8. * F(i-1, j, k, ncomp) - F(i-2, j, k, ncomp) - 7. * F(i, j, k, ncomp));
+           } else if (Ms_lo_x == 0){
+               return 0.5 * inv_dx * inv_dx * (8. * F(i+1, j, k, ncomp) - F(i+2, j, k, ncomp) - 7. * F(i, j, k, ncomp));
+           } else {
+               return inv_dx*(UpwardDx(F, coefs_x, n_coefs_x, i, j, k, ncomp) - DownwardDx(F, coefs_x, n_coefs_x, i, j, k, ncomp));
+           }
+        } else { // at y or z faces
+           if (Ms_hi_x == 0.){
+               return inv_dx*(0. - DownwardDx(F, coefs_x, n_coefs_x, i, j, k, ncomp));
+           } else if (Ms_lo_x == 0.){
+               return inv_dx*(UpwardDx(F, coefs_x, n_coefs_x, i, j, k, ncomp) - 0.);
+           } else {
+               return inv_dx*(UpwardDx(F, coefs_x, n_coefs_x, i, j, k, ncomp) - DownwardDx(F, coefs_x, n_coefs_x, i, j, k, ncomp));
+           }
+        }
+    }
+
+    /**
+     * Perform divergence of gradient along y on M field when exchange coupling is on*/
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static amrex::Real LaplacianDy_Mag (
+        amrex::Array4<amrex::Real> const& F,
+        amrex::Real const * const coefs_y, int const n_coefs_y, amrex::Real const Ms_lo_y, amrex::Real const Ms_hi_y,
+        int const i, int const j, int const k, int const ncomp=0, int const nodality=0) {
+
+        amrex::Real const inv_dy = coefs_y[0];
+        if (nodality == 1){ // y face (normal face), dM/dy = 0
+           if (Ms_hi_y == 0.){
+               return 0.5 * inv_dy * inv_dy * (8. * F(i, j-1, k, ncomp) - F(i, j-2, k, ncomp) - 7. * F(i, j, k, ncomp));
+           } else if (Ms_lo_y == 0.) {
+               return 0.5 * inv_dy * inv_dy * (8. * F(i, j+1, k, ncomp) - F(i, j+2, k, ncomp) - 7. * F(i, j, k, ncomp));
+           } else {
+               return inv_dy*(UpwardDy(F, coefs_y, n_coefs_y, i, j, k, ncomp) - DownwardDy(F, coefs_y, n_coefs_y, i, j, k, ncomp));
+           }
+
+        } else { // x or z faces
+           if (Ms_hi_y == 0.){
+               return inv_dy*(0. - DownwardDy(F, coefs_y, n_coefs_y, i, j, k, ncomp));
+           } else if (Ms_lo_y == 0.) {
+               return inv_dy*(UpwardDy(F, coefs_y, n_coefs_y, i, j, k, ncomp) - 0.);
+           } else {
+               return inv_dy*(UpwardDy(F, coefs_y, n_coefs_y, i, j, k, ncomp) - DownwardDy(F, coefs_y, n_coefs_y, i, j, k, ncomp));
+           }
+        }
+    }
+
+     /**
+     * Perform divergence of gradient along z on M field when exchange coupling is on*/
+   AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static amrex::Real LaplacianDz_Mag (
+        amrex::Array4<amrex::Real> const& F,
+        amrex::Real const * const coefs_z, int const n_coefs_z, amrex::Real const Ms_lo_z, amrex::Real const Ms_hi_z,
+        int const i, int const j, int const k, int const ncomp=0, int const nodality=0) {
+
+        amrex::Real const inv_dz = coefs_z[0];
+        if (nodality == 2){ // z face (normal face) dM/dz = 0
+           if ( Ms_hi_z == 0.) {
+               return 0.5 * inv_dz * inv_dz * (8. * F(i, j, k-1, ncomp) - F(i, j, k-2, ncomp) - 7. * F(i, j, k, ncomp));
+           } else if ( Ms_lo_z == 0.){
+               return 0.5 * inv_dz * inv_dz * (8. * F(i, j, k+1, ncomp) - F(i, j, k+2, ncomp) - 7. * F(i, j, k, ncomp));
+           } else {
+               return inv_dz*(UpwardDz(F, coefs_z, n_coefs_z, i, j, k, ncomp) - DownwardDz(F, coefs_z, n_coefs_z, i, j, k, ncomp));
+           }
+
+        } else { // x or y face
+           if ( Ms_hi_z == 0.) {
+               return inv_dz*(0. - DownwardDz(F, coefs_z, n_coefs_z, i, j, k, ncomp));
+           } else if ( Ms_lo_z == 0.){
+               return inv_dz*(UpwardDz(F, coefs_z, n_coefs_z, i, j, k, ncomp) - 0.);
+           } else {
+               return inv_dz*(UpwardDz(F, coefs_z, n_coefs_z, i, j, k, ncomp) - DownwardDz(F, coefs_z, n_coefs_z, i, j, k, ncomp));
+           }
+        }
+    }
+
+     /**
+     * Compute the sum to get Laplacian of M field when exchange coupling is on*/
+   AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    static amrex::Real Laplacian_Mag (
+        amrex::Array4<amrex::Real> const& F,
+        amrex::Real const * const coefs_x, amrex::Real const * const coefs_y, amrex::Real const * const coefs_z,
+        int const n_coefs_x, int const n_coefs_y, int const n_coefs_z,
+        amrex::Real const Ms_lo_x, amrex::Real const Ms_hi_x, amrex::Real const Ms_lo_y, amrex::Real const Ms_hi_y, amrex::Real const Ms_lo_z, amrex::Real const Ms_hi_z,
+        int const i, int const j, int const k, int const ncomp=0, int const nodality=0) {
+
+        return LaplacianDx_Mag(F, coefs_x, n_coefs_x, Ms_lo_x, Ms_hi_x, i, j, k, ncomp, nodality) + LaplacianDy_Mag(F, coefs_y, n_coefs_y, Ms_lo_y, Ms_hi_y, i, j, k, ncomp, nodality) + LaplacianDz_Mag(F, coefs_z, n_coefs_z, Ms_lo_z, Ms_hi_z, i, j, k, ncomp, nodality);
+    }
+
+#endif
+
 };
 
 #endif // WARPX_FINITE_DIFFERENCE_ALGORITHM_CARTESIAN_YEE_H_

Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H

@@ -142,6 +142,7 @@ class FiniteDifferenceSolver
                        std::unique_ptr<MacroscopicProperties> const &macroscopic_properties);
 
         void MacroscopicEvolveHM_2nd (
+                       int lev,
                        std::array<std::unique_ptr<amrex::MultiFab>, 3> &Mfield,
                        std::array<std::unique_ptr<amrex::MultiFab>, 3> &Hfield,            // H Maxwell
                        std::array< std::unique_ptr<amrex::MultiFab>, 3>& Bfield,
@@ -327,6 +328,7 @@ class FiniteDifferenceSolver
 
         template< typename T_Algo >
         void MacroscopicEvolveHMCartesian_2nd(
+            int lev,
             std::array<std::unique_ptr<amrex::MultiFab>, 3> &Mfield,
             std::array<std::unique_ptr<amrex::MultiFab>, 3> &Hfield,            // H Maxwell
             std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield,