diff --git a/.github/c-h-after-uncrustify.sh b/.github/c-h-after-uncrustify.sh new file mode 100755 index 0000000000..93f750c3af --- /dev/null +++ b/.github/c-h-after-uncrustify.sh @@ -0,0 +1,50 @@ +#!/bin/bash +# ----------------------------------------------------------------------------- +# Post-process C/H files after uncrustify formatting +# ----------------------------------------------------------------------------- +# Fixes: Function return type pointer style +# FROM: struct foo * +# function_name(...) +# TO: struct foo* +# function_name(...) +# +# Usage: +# .github/c-after-uncrustify.sh [file2.h ...] +# .github/c-after-uncrustify.sh core/zero/*.h +# ----------------------------------------------------------------------------- + +if [ $# -eq 0 ]; then + echo "Usage: $0 [file2.h ...]" + echo "Post-processes C/H files after uncrustify formatting." + exit 1 +fi + +for file in "$@"; do + if [ ! -f "$file" ]; then + echo "Warning: File not found: $file" + continue + fi + + echo "Processing: $file" + + # Fix function return type pointer style: + # When a line ends with " *" and the next line is a function name, + # change " *" to "*" (remove the space before the asterisk) + # + # Pattern: "word *\n" followed by a line starting with a function name + # This handles: + # struct gkyl_dg_array_mask * + # gkyl_dg_array_mask_acquire(...) + # Becomes: + # struct gkyl_dg_array_mask* + # gkyl_dg_array_mask_acquire(...) + + perl -i -0pe ' + # Match: (type) space asterisk newline (function_name) + # Replace with: (type) asterisk newline (function_name) + s/(\w) \*\n(\w+\s*\()/$1*\n$2/g; + ' "$file" + +done + +echo "Done." diff --git a/.github/cuda-after-uncrustify.sh b/.github/cuda-after-uncrustify.sh new file mode 100755 index 0000000000..c67d2e867b --- /dev/null +++ b/.github/cuda-after-uncrustify.sh @@ -0,0 +1,71 @@ +#!/bin/bash +# ----------------------------------------------------------------------------- +# Post-process CUDA files after uncrustify formatting +# ----------------------------------------------------------------------------- +# Fixes issues caused by uncrustify not understanding CUDA syntax: +# 1. Removes spaces in/around kernel launch operators: << < -> <<<, >> > -> >>> +# 2. Fixes indentation of arguments after kernel launch calls +# +# Usage: +# .github/cuda-after-uncrustify.sh [file2.cu ...] +# .github/cuda-after-uncrustify.sh core/zero/*.cu +# ----------------------------------------------------------------------------- + +if [ $# -eq 0 ]; then + echo "Usage: $0 [file2.cu ...]" + echo "Post-processes CUDA files after uncrustify formatting." + exit 1 +fi + +for file in "$@"; do + if [ ! -f "$file" ]; then + echo "Warning: File not found: $file" + continue + fi + + echo "Processing: $file" + + # Fix 1: Remove spaces in and around <<< and >>> operators + # Handles: << < -> <<<, >> > -> >>>, and removes surrounding spaces + sed -i 's/<< > >/>>>/g' "$file" + sed -i 's/ <<<\s*/<<>> />>>/' "$file" + + # Fix 2: Fix indentation of continuation lines after kernel launches + # Find lines with <<<...>>>( and fix indentation of following lines + # until we hit the closing ); + perl -i -pe ' + BEGIN { $in_kernel_call = 0; $base_indent = ""; } + + # Detect kernel launch line ending with ( + if (/^(\s*)(\S+)\s*<<<.*>>>\s*\(\s*$/) { + $in_kernel_call = 1; + $base_indent = $1 . " "; # 2 spaces from function start + next; + } + + # If we are in a kernel call, fix indentation + if ($in_kernel_call) { + # Check if this line ends the call + if (/\);\s*$/) { + s/^(\s*)/$base_indent/; + $in_kernel_call = 0; + } else { + # Fix indentation of continuation lines + s/^(\s*)/$base_indent/; + } + } + ' "$file" + + # Fix 3: Function return type pointer style: + # When a line ends with " *" and the next line is a function name, + # change " *" to "*" (remove the space before the asterisk) + perl -i -0pe ' + # Match: (type) space asterisk newline (function_name) + # Replace with: (type) asterisk newline (function_name) + s/(\w) \*\n(\w+\s*\()/$1*\n$2/g; + ' "$file" + +done + +echo "Done." diff --git a/.github/uncrustify-wrapper.sh b/.github/uncrustify-wrapper.sh new file mode 100755 index 0000000000..daa686f654 --- /dev/null +++ b/.github/uncrustify-wrapper.sh @@ -0,0 +1,90 @@ +#!/bin/bash +# ----------------------------------------------------------------------------- +# Uncrustify wrapper for VS Code +# ----------------------------------------------------------------------------- +# This script wraps uncrustify to add post-processing for CUDA files. +# Set this as the "Uncrustify › Executable Path" in VS Code settings. +# +# It passes all arguments to uncrustify, then runs cuda-after-uncrustify.sh +# if the file being formatted is a .cu file. +# +# Works with: +# - Format Document +# - Format Selection (uncrustify handles this via --frag flag) +# +# How to set up in VS Code: +# +# Navigate to .vscode settings in this repository and open .vscode/settings.json. +# Then, set the "uncrustify.executablePath" setting to point to this script +# using an ABSOLUTE path, e.g. +# +# "uncrustify.configPath.linux": "/path/to/gkeyll/.github/uncrustify.cfg", +# "uncrustify.executablePath.linux": "/path/to/gkeyll/.github/uncrustify-wrapper.sh", +# +# To enable formatting of .h and .cu files, add these lines to .vscode/settings.json: +# "files.associations": { +# "*.h": "cpp", +# "*.cu": "cpp", +# ----------------------------------------------------------------------------- + +SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)" +CUDA_POST_SCRIPT="$SCRIPT_DIR/cuda-after-uncrustify.sh" +C_H_POST_SCRIPT="$SCRIPT_DIR/c-h-after-uncrustify.sh" + +# Find the actual uncrustify binary; fall back to known locations if PATH is stale +# (VS Code may have a cached PATH from before uncrustify was installed). +UNCRUSTIFY_BIN=$(which uncrustify 2>/dev/null) +if [ -z "$UNCRUSTIFY_BIN" ]; then + for candidate in /usr/bin/uncrustify /usr/local/bin/uncrustify; do + if [ -x "$candidate" ]; then + UNCRUSTIFY_BIN="$candidate" + break + fi + done +fi +if [ -z "$UNCRUSTIFY_BIN" ]; then + echo "Error: uncrustify not found in PATH or known locations" >&2 + exit 1 +fi + +# Run uncrustify with all passed arguments +"$UNCRUSTIFY_BIN" "$@" +UNCRUSTIFY_EXIT=$? + +# If uncrustify failed, exit with its error code +if [ $UNCRUSTIFY_EXIT -ne 0 ]; then + exit $UNCRUSTIFY_EXIT +fi + +# Check if any argument is a .cu, .c, or .h file and if --replace or -o was used +# (meaning the file was modified in place or output was written) +CUDA_FILE="" +C_H_FILE="" +HAS_REPLACE=0 + +for arg in "$@"; do + case "$arg" in + --replace|--no-backup) + HAS_REPLACE=1 + ;; + *.cu) + CUDA_FILE="$arg" + ;; + *.c|*.h) + C_H_FILE="$arg" + ;; + esac +done + +# Postprocessing scripts +# .c and .h files +if [ -n "$C_H_FILE" ] && [ $HAS_REPLACE -eq 1 ] && [ -x "$C_H_POST_SCRIPT" ]; then + "$C_H_POST_SCRIPT" "$C_H_FILE" >/dev/null 2>&1 +fi + +# CUDA files +if [ -n "$CUDA_FILE" ] && [ $HAS_REPLACE -eq 1 ] && [ -x "$CUDA_POST_SCRIPT" ]; then + "$CUDA_POST_SCRIPT" "$CUDA_FILE" >/dev/null 2>&1 +fi + +exit 0 diff --git a/.github/uncrustify.cfg b/.github/uncrustify.cfg new file mode 100644 index 0000000000..d45cc4deb3 --- /dev/null +++ b/.github/uncrustify.cfg @@ -0,0 +1,230 @@ +# Uncrustify Configuration for Gkeyll +# ----------------------------------------------------------------------------- +# This configuration file defines the coding style for the Gkeyll codebase. +# Based on K&R style with project-specific customizations. +# +# Usage: +# uncrustify -c .github/uncrustify.cfg --replace +# uncrustify -c .github/uncrustify.cfg --check +# +# VScode: install the uncrustify extension and set the path to .github/uncrustify.cfg +# +# NOTE: when using uncrustify on .cu files, it may incorrectly add spaces in CUDA kernel +# syntax (<<<...>>>) by adding spaces. Postprocess the .cu files manually. +# ----------------------------------------------------------------------------- + +# General options +newlines = lf # Unix line endings +input_tab_size = 2 # Tab size for input +output_tab_size = 2 # Tab size for output +indent_columns = 2 # Number of columns for indentation +indent_with_tabs = 0 # Use spaces, not tabs + +# File filtering - only C and header files (NOT .cu - breaks CUDA syntax) +file_ext = c,h # Only process these extensions + +# ----------------------------------------------------------------------------- +# Brace style: K&R (opening brace on same line as control statement) +# ----------------------------------------------------------------------------- +nl_if_brace = remove # if () { on same line +nl_else_brace = remove # else { on same line +nl_elseif_brace = remove # else if () { on same line +nl_for_brace = remove # for () { on same line +nl_while_brace = remove # while () { on same line +nl_do_brace = remove # do { on same line +nl_switch_brace = remove # switch () { on same line +nl_brace_else = force # } else on separate lines +nl_brace_while = remove # } while on same line (for do-while) + +# Function braces: opening brace on new line (K&R for C functions) +nl_fdef_brace = force # Function definition: brace on new line +nl_func_paren = remove # Keep function name and ( on same line +nl_func_def_paren = remove # Keep function def name and ( on same line +nl_func_decl_start = ignore # Don't force newline after return type +nl_func_decl_end = ignore # Don't force newline before closing ) +nl_func_def_end = ignore # Don't force newline before ) in definitions + +# Struct/enum/union braces +nl_struct_brace = remove # struct { on same line +nl_union_brace = remove # union { on same line +nl_enum_brace = remove # enum { on same line + +# ----------------------------------------------------------------------------- +# Spacing around operators and keywords +# ----------------------------------------------------------------------------- +sp_arith = force # Space around arithmetic operators +sp_assign = force # Space around assignment +sp_compare = force # Space around comparison operators +sp_bool = force # Space around boolean operators +sp_after_comma = force # Space after comma +sp_before_comma = remove # No space before comma + +# Exception: allow no space around + in certain contexts (like poly_order+1) +sp_arith_additive = ignore # Allow no space around + and - + +# Pointer and reference spacing: +# - Variable declarations: * attached to name (e.g., struct foo *bar) +# - Function return types: * attached to type (e.g., struct foo*) +# - Function parameters: * attached to name (e.g., struct foo *bar) +# Note: sp_before_ptr_star_func and sp_after_ptr_star_func control function return types +sp_before_ptr_star = force # Space before * in variable declarations +sp_after_ptr_star = remove # No space after * (attached to name) +sp_between_ptr_star = remove # No space between multiple * +sp_before_ptr_star_func = remove # No space before * in function return types (attached to type) +sp_after_ptr_star_func = ignore # Don't force space after * in return types (newline follows) +sp_ptr_star_paren = remove # No space between * and ( in casts +sp_after_ptr_star_qualifier = force # Space after * before qualifiers like GKYL_RESTRICT (produces: double* GKYL_RESTRICT) + +# Function call/definition spacing +sp_func_call_paren = remove # No space between function name and ( +sp_func_def_paren = remove # No space between function name and ( +sp_func_proto_paren = remove # No space in prototypes +sp_func_class_paren = remove # No space for class methods +sp_inside_paren = remove # No space inside parentheses +sp_inside_fparen = remove # No space inside function parentheses + +# Control statement spacing +sp_before_sparen = force # Space before ( in if/for/while +sp_inside_sparen = remove # No space inside control statement parens +sp_after_sparen = force # Space after ) in control statements +sp_sparen_brace = force # Space between ) and { in control statements + +# Brace spacing +sp_inside_braces = force # Space inside { } for arrays and initializers +sp_inside_braces_struct = remove # No space inside struct { } (keep struct fields tight) +sp_inside_braces_enum = remove # No space inside enum { } +sp_brace_else = force # Space between } and else +sp_else_brace = force # Space between else and { + +# Struct initialization and compound literals +nl_brace_struct_var = remove # Keep struct var = { on same line +pos_bool = trail # Position of closing brace for compound literals + +# Semicolon spacing +sp_before_semi = remove # No space before semicolon +sp_after_semi = force # Space after semicolon in for loops +sp_after_semi_for_empty = remove # No space in empty for clauses + +# Cast spacing +sp_after_cast = remove # No space after cast (void*)ptr +sp_inside_paren_cast = remove # No space inside cast parens + +# Square bracket spacing +sp_before_square = remove # No space before [ +sp_inside_square = remove # No space inside [ ] + +# ----------------------------------------------------------------------------- +# Newlines and blank lines +# ----------------------------------------------------------------------------- +nl_after_brace_open = false # No newline after { in single-line +nl_after_func_body = 2 # Blank line after function body +nl_after_func_proto = 1 # Single newline after prototype +nl_max = 2 # Maximum consecutive blank lines +nl_before_block_comment = 2 # Blank line before block comments +nl_after_struct = 1 # Newline after struct definition + +# Collapse empty blocks +nl_collapse_empty_body = true # Collapse empty function bodies + +# ----------------------------------------------------------------------------- +# Alignment (minimal - Gkeyll uses simple alignment) +# ----------------------------------------------------------------------------- +align_var_def_span = 0 # Don't align variable definitions +align_var_def_star_style = 1 # Align * with variable name +align_assign_span = 0 # Don't align assignments +align_enum_equ_span = 0 # Don't align enum assignments +align_struct_init_span = 0 # Don't align struct initializers +align_right_cmt_span = 0 # Don't align right-side comments +align_pp_define_span = 0 # Don't align preprocessor defines + +# Function call parentheses and struct literals +indent_paren_nl = true # Indent when parentheses are split +indent_paren_close = 0 # Don't add extra indent to closing paren +indent_paren_after_func_def = false # Don't extra-indent function def params +pos_arith = trail # Position arithmetic operators at end of line +nl_func_call_paren = remove # Keep opening paren on same line +nl_func_call_paren_empty = remove # Keep empty parens together + +# Compound literal / brace initializer indentation +indent_brace_parent = false # Don't indent braces relative to parent +indent_continue = 2 # Continue line indent +indent_shift = 0 # Additional shift for continuation +nl_ds_struct_enum_cmt = false # Don't add newline after struct/enum comment +eat_blanks_before_close_brace = true # Remove blank lines before } +nl_func_call_end = remove # Keep closing ) on same line when possible +nl_func_call_end_multi_line = false # Don't force newline before ) in multi-line calls + +# ----------------------------------------------------------------------------- +# Comments +# ----------------------------------------------------------------------------- +cmt_indent_multi = true # Indent multi-line comments +cmt_star_cont = true # Add * at start of continued comments +cmt_sp_before_star_cont = 0 # Space before continued comment * +cmt_sp_after_star_cont = 1 # Space after continued comment * +sp_cmt_cpp_start = force # Space after // in comments + +# ----------------------------------------------------------------------------- +# Preprocessor +# ----------------------------------------------------------------------------- +pp_indent = remove # Don't indent preprocessor directives +pp_space_after = remove # No space after # in directives +pp_indent_if = 0 # Don't indent #if blocks +pp_if_indent_code = false # Don't indent code in #if blocks + +# ----------------------------------------------------------------------------- +# Line length and wrapping +# ----------------------------------------------------------------------------- +code_width = 100 # Maximum line width. 100 is best for Vim +ls_for_split_full = true # Split for loops at semicolons +ls_func_split_full = false # Don't force split on long function calls + +# Control where line breaks can occur +ls_code_width = false # Disable automatic line splitting +pos_arith = trail # Keep operators at end of line +pos_assign = trail # Keep assignment operators at end of line +pos_bool = trail # Keep boolean operators at end of line + +# ----------------------------------------------------------------------------- +# Single-line statements (require braces per coding guidelines) +# ----------------------------------------------------------------------------- +mod_full_brace_if = add # Always add braces to if +mod_full_brace_for = add # Always add braces to for +mod_full_brace_while = add # Always add braces to while +mod_full_brace_do = add # Always add braces to do +mod_full_brace_if_chain = 2 # Add braces to entire if-else chain (2 = add) + +# ----------------------------------------------------------------------------- +# Parentheses +# ----------------------------------------------------------------------------- +mod_paren_on_return = remove # No parens on return: return x; not return (x); + +# ----------------------------------------------------------------------------- +# Sorting (disabled - maintain original order) +# ----------------------------------------------------------------------------- +mod_sort_include = false # Don't sort includes +mod_sort_using = false # Don't sort using statements + +# ----------------------------------------------------------------------------- +# Type handling +# ----------------------------------------------------------------------------- +indent_class = true # Indent class bodies +indent_func_def_param = true # Indent continued function params +indent_func_call_param = true # Indent continued call params +indent_switch_case = 2 # Indent case labels inside switch statements + +# Long parameter lists: continuation indent +indent_param = 2 # Continuation indent for parameters +indent_func_def_param_paren_pos_threshold = 0 + +# Function parameter wrapping +ls_func_split_full = false # Don't force function parameter splitting +ls_code_width = false # Disable automatic line splitting based on code_width + +# ----------------------------------------------------------------------------- +# Special Gkeyll patterns +# ----------------------------------------------------------------------------- +# Handle common macros +set MACRO GKYL_MAX2 GKYL_MIN2 NELM NCOM NSIZE + +# Function attributes - tell uncrustify these are function attributes +set ATTRIBUTE GKYL_CU_DH GKYL_CU_D GKYL_CU_H GKYL_RESTRICT \ No newline at end of file diff --git a/gyrokinetic/unit/ctest_position_map.c b/gyrokinetic/unit/ctest_position_map.c index 7a12623471..8909cc4b13 100644 --- a/gyrokinetic/unit/ctest_position_map.c +++ b/gyrokinetic/unit/ctest_position_map.c @@ -14,7 +14,8 @@ #include void -test_nonuniform_position_map(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx) +test_nonuniform_position_map(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, + void *ctx) { double poly_order = 2; double z = xn[0]; @@ -23,13 +24,14 @@ test_nonuniform_position_map(double t, const double *GKYL_RESTRICT xn, double *G if (z < -left) fout[0] = z; else if (z < right) - fout[0] = - pow(z - right, poly_order)/fabs(pow(left-right, poly_order-1)) + right; + fout[0] = -pow(z - right, poly_order) / fabs(pow(left - right, poly_order - 1)) + right; else fout[0] = z; } void -test_nonuniform_position_map_slope(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx) +test_nonuniform_position_map_slope(double t, const double *GKYL_RESTRICT xn, + double *GKYL_RESTRICT fout, void *ctx) { double poly_order = 2; double z = xn[0]; @@ -38,53 +40,54 @@ test_nonuniform_position_map_slope(double t, const double *GKYL_RESTRICT xn, dou if (z < -left) fout[0] = 1.0; else if (z < right) - fout[0] = - poly_order * pow(z - right, poly_order-1)/fabs(pow(left-right, poly_order-1)); + fout[0] = -poly_order* pow(z - right, poly_order - 1) / fabs(pow(left - right, poly_order - 1)); else fout[0] = 1.0; } void -test_identity_position_map(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx) +test_identity_position_map(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, + void *ctx) { fout[0] = xn[0]; } void -test_nonuniform_position_map_3x(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx) +test_nonuniform_position_map_3x(double t, const double *GKYL_RESTRICT xn, + double *GKYL_RESTRICT fout, void *ctx) { double poly_order = 2; double left = 0.25; double right = 0.75; - for (int i = 0; i<3; i++) - { + for (int i = 0; i < 3; i++) { double z = xn[i]; if (z < -left) fout[i] = z; else if (z < right) - fout[i] = - pow(z - right, poly_order)/fabs(pow(left-right, poly_order-1)) + right; + fout[i] = -pow(z - right, poly_order) / fabs(pow(left - right, poly_order - 1)) + right; else fout[i] = z; } } -void +void bmag_func(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx) { double x = xn[0]; double s = 0.6; double c = 0.; // double B = (4*pow(s*(x-c),2) - 0.3*pow(s*(x-c),4) + 1)*exp(-pow(s*(x-c),2)); - double B = 1/(1+100*pow(x-M_PI/2,2)) + 1/(1+100*pow(x+M_PI/2,2)); + double B = 1 / (1 + 100 * pow(x - M_PI / 2, 2)) + 1 / (1 + 100 * pow(x + M_PI / 2, 2)); fout[0] = B; } void test_position_map_init_1x() { - int cells[] = {32}; + int cells[] = { 32 }; int poly_order = 1; - double lower[] = {0.0}, upper[] = {1.0}; - int dim = sizeof(lower)/sizeof(lower[0]); + double lower[] = { 0.0 }, upper[] = { 1.0 }; + int dim = sizeof(lower) / sizeof(lower[0]); // Grids. struct gkyl_rect_grid grid; gkyl_rect_grid_init(&grid, dim, lower, upper, cells); @@ -92,14 +95,15 @@ test_position_map_init_1x() int ghost[] = { 1, 1 }; struct gkyl_range localRange, localRange_ext; // local, local-ext ranges. gkyl_create_grid_ranges(&grid, ghost, &localRange_ext, &localRange); - + // Basis functions. struct gkyl_basis basis; gkyl_cart_modal_serendip(&basis, dim, poly_order); struct gkyl_position_map_inp pos_map_inp = { - .maps = {test_nonuniform_position_map, test_nonuniform_position_map, test_nonuniform_position_map}, - .ctxs = {NULL, NULL, NULL}, + .maps = { test_nonuniform_position_map, test_nonuniform_position_map, + test_nonuniform_position_map }, + .ctxs = { NULL, NULL, NULL }, }; struct gkyl_position_map *pos_map = gkyl_position_map_new(pos_map_inp, @@ -115,14 +119,13 @@ test_position_map_init_1x() gkyl_position_map_release(pos_map); } - void test_position_map_init_1x_null() { - int cells[] = {8}; + int cells[] = { 8 }; int poly_order = 1; - double lower[] = {0.0}, upper[] = {1.0}; - int dim = sizeof(lower)/sizeof(lower[0]); + double lower[] = { 0.0 }, upper[] = { 1.0 }; + int dim = sizeof(lower) / sizeof(lower[0]); // Grids. struct gkyl_rect_grid grid; gkyl_rect_grid_init(&grid, dim, lower, upper, cells); @@ -130,19 +133,19 @@ test_position_map_init_1x_null() int ghost[] = { 1, 1 }; struct gkyl_range localRange, localRange_ext; // local, local-ext ranges. gkyl_create_grid_ranges(&grid, ghost, &localRange_ext, &localRange); - + // Basis functions. struct gkyl_basis basis; gkyl_cart_modal_serendip(&basis, dim, poly_order); struct gkyl_position_map_inp pos_map_inp = { }; - + struct gkyl_position_map *pos_map = gkyl_position_map_new(pos_map_inp, \ grid, localRange, localRange_ext, localRange, localRange_ext, basis); TEST_CHECK(pos_map->id == GKYL_PMAP_USER_INPUT); - for (double i = 0; i < 1; i = i+0.1){ - double x[1] = {i}; + for (double i = 0; i < 1; i = i + 0.1) { + double x[1] = { i }; double y[1]; pos_map->maps[0](0.0, x, y, pos_map->ctxs[0]); TEST_CHECK(y[0] == x[0]); @@ -166,10 +169,10 @@ test_position_map_init_1x_null() void test_position_map_init_2x() { - int cells[] = {8,8}; + int cells[] = { 8, 8 }; int poly_order = 1; - double lower[] = {0.0, 0.0}, upper[] = {1.0, 1.0}; - int dim = sizeof(lower)/sizeof(lower[0]); + double lower[] = { 0.0, 0.0 }, upper[] = { 1.0, 1.0 }; + int dim = sizeof(lower) / sizeof(lower[0]); // Grids. struct gkyl_rect_grid grid; gkyl_rect_grid_init(&grid, dim, lower, upper, cells); @@ -177,14 +180,15 @@ test_position_map_init_2x() int ghost[] = { 1, 1 }; struct gkyl_range localRange, localRange_ext; // local, local-ext ranges. gkyl_create_grid_ranges(&grid, ghost, &localRange_ext, &localRange); - + // Basis functions. struct gkyl_basis basis; gkyl_cart_modal_serendip(&basis, dim, poly_order); struct gkyl_position_map_inp pos_map_inp = { - .maps = {test_nonuniform_position_map, test_nonuniform_position_map, test_nonuniform_position_map}, - .ctxs = {0, 0, 0}, + .maps = { test_nonuniform_position_map, test_nonuniform_position_map, + test_nonuniform_position_map }, + .ctxs = { 0, 0, 0 }, }; struct gkyl_position_map *pos_map = gkyl_position_map_new(pos_map_inp, \ @@ -203,25 +207,26 @@ test_position_map_init_2x() void test_position_map_init_3x() { - int cells[] = {8, 8, 8}; + int cells[] = { 8, 8, 8 }; int poly_order = 1; - double lower[] = {0.0, 0.0, 0.0}, upper[] = {1.0, 1.0, 1.0}; - int dim = sizeof(lower)/sizeof(lower[0]); + double lower[] = { 0.0, 0.0, 0.0 }, upper[] = { 1.0, 1.0, 1.0 }; + int dim = sizeof(lower) / sizeof(lower[0]); // Grids. struct gkyl_rect_grid grid; gkyl_rect_grid_init(&grid, dim, lower, upper, cells); // Ranges - int ghost[] = { 1, 1, 1}; + int ghost[] = { 1, 1, 1 }; struct gkyl_range localRange, localRange_ext; // local, local-ext ranges. gkyl_create_grid_ranges(&grid, ghost, &localRange_ext, &localRange); - + // Basis functions. struct gkyl_basis basis; gkyl_cart_modal_serendip(&basis, dim, poly_order); struct gkyl_position_map_inp pos_map_inp = { - .maps = {test_nonuniform_position_map, test_nonuniform_position_map, test_nonuniform_position_map}, - .ctxs = {0, 0, 0}, + .maps = { test_nonuniform_position_map, test_nonuniform_position_map, + test_nonuniform_position_map }, + .ctxs = { 0, 0, 0 }, }; struct gkyl_position_map *pos_map = gkyl_position_map_new(pos_map_inp, \ @@ -240,94 +245,97 @@ test_position_map_init_3x() void test_position_map_set() { - int cells[] = {8, 8, 8}; + int cells[] = { 8, 8, 8 }; int poly_order = 1; - double lower[] = {0.0, 0.0, 0.0}, upper[] = {1.0, 1.0, 1.0}; - int dim = sizeof(lower)/sizeof(lower[0]); + double lower[] = { 0.0, 0.0, 0.0 }, upper[] = { 1.0, 1.0, 1.0 }; + int dim = sizeof(lower) / sizeof(lower[0]); // Grids. struct gkyl_rect_grid grid; gkyl_rect_grid_init(&grid, dim, lower, upper, cells); // Ranges - int ghost[] = { 1, 1, 1}; + int ghost[] = { 1, 1, 1 }; struct gkyl_range localRange, localRange_ext; // local, local-ext ranges. gkyl_create_grid_ranges(&grid, ghost, &localRange_ext, &localRange); - + // Basis functions. struct gkyl_basis basis; gkyl_cart_modal_serendip(&basis, dim, poly_order); - + struct gkyl_position_map_inp pos_map_inp = { - .maps = {test_nonuniform_position_map, test_nonuniform_position_map, test_nonuniform_position_map}, - .ctxs = {0, 0, 0}, + .maps = { test_nonuniform_position_map, test_nonuniform_position_map, + test_nonuniform_position_map }, + .ctxs = { 0, 0, 0 }, }; struct gkyl_position_map *pos_map = gkyl_position_map_new(pos_map_inp, \ grid, localRange, localRange_ext, localRange, localRange_ext, basis); struct gkyl_array *pmap_arr_set = gkyl_array_new(GKYL_DOUBLE, \ - 3*pos_map->basis.num_basis, pos_map->local_ext.volume); + 3 * pos_map->basis.num_basis, pos_map->local_ext.volume); gkyl_array_clear(pmap_arr_set, 1.0); gkyl_position_map_set_mc2nu(pos_map, pmap_arr_set); - double *pos_map_i = pos_map->mc2nu->data; - for (unsigned i=0; imc2nu->size; ++i) - TEST_CHECK( gkyl_compare(pos_map_i[i], 1.0, 1e-14) ); + double *pos_map_i = pos_map->mc2nu->data; + for (unsigned i = 0; i < pos_map->mc2nu->size; ++i) { + TEST_CHECK(gkyl_compare(pos_map_i[i], 1.0, 1e-14) ); + } gkyl_array_release(pmap_arr_set); gkyl_position_map_release(pos_map); } - void test_gkyl_position_map_eval_mc2nu() { - int cells[] = {8, 8, 8}; + int cells[] = { 8, 8, 8 }; int poly_order = 2; - double lower[] = {0.0, 0.0, 0.0}, upper[] = {1.0, 1.0, 1.0}; - int dim = sizeof(lower)/sizeof(lower[0]); + double lower[] = { 0.0, 0.0, 0.0 }, upper[] = { 1.0, 1.0, 1.0 }; + int dim = sizeof(lower) / sizeof(lower[0]); // Grids. struct gkyl_rect_grid grid; gkyl_rect_grid_init(&grid, dim, lower, upper, cells); // Ranges - int ghost[] = { 1, 1, 1}; + int ghost[] = { 1, 1, 1 }; struct gkyl_range localRange, localRange_ext; // local, local-ext ranges. gkyl_create_grid_ranges(&grid, ghost, &localRange_ext, &localRange); - + // Basis functions. struct gkyl_basis basis; gkyl_cart_modal_serendip(&basis, dim, poly_order); - + struct gkyl_position_map_inp pos_map_inp = { - .maps = {test_nonuniform_position_map, test_nonuniform_position_map, test_nonuniform_position_map}, - .ctxs = {0, 0, 0}, + .maps = { test_nonuniform_position_map, test_nonuniform_position_map, + test_nonuniform_position_map }, + .ctxs = { 0, 0, 0 }, }; struct gkyl_position_map *pos_map = gkyl_position_map_new(pos_map_inp, \ grid, localRange, localRange_ext, localRange, localRange_ext, basis); struct gkyl_array *pmap_arr_set = gkyl_array_new(GKYL_DOUBLE, \ - 3*pos_map->basis.num_basis, pos_map->local_ext.volume); + 3 * pos_map->basis.num_basis, pos_map->local_ext.volume); gkyl_proj_on_basis *projDistf = gkyl_proj_on_basis_new(&grid, &basis, - poly_order+1, 3, test_nonuniform_position_map_3x, 0); + poly_order + 1, 3, test_nonuniform_position_map_3x, 0); gkyl_proj_on_basis_advance(projDistf, 0.0, &localRange, pmap_arr_set); gkyl_proj_on_basis_release(projDistf); gkyl_position_map_set_mc2nu(pos_map, pmap_arr_set); - for (int i=0; i<3; i++) { - for (int j=0; j<3; j++) { - for (int k=0; k<5; k++) { - double x[3] = {i/10.0, j/10.0, k/10.0}; + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + for (int k = 0; k < 5; k++) { + double x[3] = { i / 10.0, j / 10.0, k / 10.0 }; double x_fa[3]; gkyl_position_map_eval_mc2nu(pos_map, x, x_fa); double x_analytic[3]; test_nonuniform_position_map(0.0, &x[0], &x_analytic[0], 0); test_nonuniform_position_map(0.0, &x[1], &x_analytic[1], 0); test_nonuniform_position_map(0.0, &x[2], &x_analytic[2], 0); - for (int d=0; d<3; ++d) - TEST_CHECK( gkyl_compare(x_fa[d], x_analytic[d], 1e-12) ); + for (int d = 0; d < 3; ++d) { + TEST_CHECK(gkyl_compare(x_fa[d], x_analytic[d], 1e-12) ); + } } } } @@ -336,48 +344,48 @@ test_gkyl_position_map_eval_mc2nu() gkyl_position_map_release(pos_map); } - void test_gkyl_position_map_slope() { - int cells[] = {8, 8, 8}; + int cells[] = { 8, 8, 8 }; int poly_order = 2; - double lower[] = {0.0, 0.0, 0.0}, upper[] = {1.0, 1.0, 1.0}; - int dim = sizeof(lower)/sizeof(lower[0]); + double lower[] = { 0.0, 0.0, 0.0 }, upper[] = { 1.0, 1.0, 1.0 }; + int dim = sizeof(lower) / sizeof(lower[0]); // Grids. struct gkyl_rect_grid grid; gkyl_rect_grid_init(&grid, dim, lower, upper, cells); // Ranges - int ghost[] = { 1, 1, 1}; + int ghost[] = { 1, 1, 1 }; struct gkyl_range localRange, localRange_ext; // local, local-ext ranges. gkyl_create_grid_ranges(&grid, ghost, &localRange_ext, &localRange); - + // Basis functions. struct gkyl_basis basis; gkyl_cart_modal_serendip(&basis, dim, poly_order); - + struct gkyl_position_map_inp pos_map_inp = { - .maps = {test_nonuniform_position_map, test_nonuniform_position_map, test_nonuniform_position_map}, - .ctxs = {0, 0, 0}, + .maps = { test_nonuniform_position_map, test_nonuniform_position_map, + test_nonuniform_position_map }, + .ctxs = { 0, 0, 0 }, }; struct gkyl_position_map *pos_map = gkyl_position_map_new(pos_map_inp, \ grid, localRange, localRange_ext, localRange, localRange_ext, basis); struct gkyl_array *pmap_arr_set = gkyl_array_new(GKYL_DOUBLE, \ - 3*pos_map->basis.num_basis, pos_map->local_ext.volume); + 3 * pos_map->basis.num_basis, pos_map->local_ext.volume); gkyl_proj_on_basis *projDistf = gkyl_proj_on_basis_new(&grid, &basis, - poly_order+1, 3, test_nonuniform_position_map_3x, 0); + poly_order + 1, 3, test_nonuniform_position_map_3x, 0); gkyl_proj_on_basis_advance(projDistf, 0.0, &localRange, pmap_arr_set); gkyl_proj_on_basis_release(projDistf); gkyl_position_map_set_mc2nu(pos_map, pmap_arr_set); - for (int i=0; i<8; i++) { - for (int j=0; j<8; j++) { - for (int k=0; k<8; k++) { - double x[3] = {i/8.0, j/8.0, k/8.0}; + for (int i = 0; i < 8; i++) { + for (int j = 0; j < 8; j++) { + for (int k = 0; k < 8; k++) { + double x[3] = { i / 8.0, j / 8.0, k / 8.0 }; if (x[0] == 0.25 || x[0] == 0.75) continue; if (x[1] == 0.25 || x[1] == 0.75) @@ -392,8 +400,9 @@ test_gkyl_position_map_slope() slope[0] = gkyl_position_map_slope(pos_map, 0, x[0], 1e-6, i, &localRange); slope[1] = gkyl_position_map_slope(pos_map, 1, x[1], 1e-6, j, &localRange); slope[2] = gkyl_position_map_slope(pos_map, 2, x[2], 1e-6, k, &localRange); - for (int d=0; d<3; ++d) - TEST_CHECK( gkyl_compare(slope[d], x_analytic[d], 1e-6) ); + for (int d = 0; d < 3; ++d) { + TEST_CHECK(gkyl_compare(slope[d], x_analytic[d], 1e-6) ); + } } } } @@ -404,10 +413,10 @@ test_gkyl_position_map_slope() void test_position_polynomial_map_optimize_1x() { - int cells[] = {64}; + int cells[] = { 64 }; int poly_order = 1; - double lower[] = {-M_PI+1e-2}, upper[] = {M_PI-1e-2}; - int dim = sizeof(lower)/sizeof(lower[0]); + double lower[] = { -M_PI + 1e-2 }, upper[] = { M_PI - 1e-2 }; + int dim = sizeof(lower) / sizeof(lower[0]); // Grids. struct gkyl_rect_grid grid; gkyl_rect_grid_init(&grid, dim, lower, upper, cells); @@ -415,7 +424,7 @@ test_position_polynomial_map_optimize_1x() int ghost[] = { 1, 1 }; struct gkyl_range localRange, localRange_ext; // local, local-ext ranges. gkyl_create_grid_ranges(&grid, ghost, &localRange_ext, &localRange); - + // Basis functions. struct gkyl_basis basis; gkyl_cart_modal_serendip(&basis, dim, poly_order); @@ -425,20 +434,22 @@ test_position_polynomial_map_optimize_1x() .map_strength = 1.0, }; - struct gkyl_position_map *pos_map = gkyl_position_map_new(pos_map_inp,\ + struct gkyl_position_map *pos_map = gkyl_position_map_new(pos_map_inp, \ grid, localRange, localRange_ext, localRange, localRange_ext, basis); // Project bmag_func onto bmag_global - struct gkyl_array *bmag_global = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, localRange_ext.volume); - gkyl_proj_on_basis *projB = gkyl_proj_on_basis_new(&grid, &basis, poly_order+1, 1, bmag_func, 0); + struct gkyl_array *bmag_global = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, + localRange_ext.volume); + gkyl_proj_on_basis *projB = gkyl_proj_on_basis_new(&grid, &basis, poly_order + 1, 1, bmag_func, + 0); gkyl_proj_on_basis_advance(projB, 0.0, &localRange, bmag_global); gkyl_proj_on_basis_release(projB); - + struct gkyl_rect_grid grid3D; - double lower3D[] = {0.4, -0.1, lower[0]}, upper3D[] = {0.6, 0.1, upper[0]}; - int cells3D[] = { 1, 1, cells[0]}; + double lower3D[] = { 0.4, -0.1, lower[0] }, upper3D[] = { 0.6, 0.1, upper[0] }; + int cells3D[] = { 1, 1, cells[0] }; gkyl_rect_grid_init(&grid3D, 3, lower3D, upper3D, cells3D); - int ghost3D[] = { 1, 1 , 1}; + int ghost3D[] = { 1, 1, 1 }; struct gkyl_range localRange3D, localRange3D_ext; // local, local-ext ranges. gkyl_create_grid_ranges(&grid3D, ghost3D, &localRange3D_ext, &localRange3D); @@ -447,14 +458,14 @@ test_position_polynomial_map_optimize_1x() gkyl_position_map_optimize(pos_map, grid3D, localRange3D); TEST_CHECK(pos_map->to_optimize == true); - TEST_CHECK( gkyl_compare(pos_map->constB_ctx->theta_throat, 1.565796, 1e-6) ); - TEST_CHECK( gkyl_compare(pos_map->constB_ctx->Bmag_throat, 1.093613, 1e-6) ); - TEST_CHECK( gkyl_compare(pos_map->constB_ctx->psi, 0.5, 1e-6) ); - TEST_CHECK( gkyl_compare(pos_map->constB_ctx->alpha, 0.0, 1e-6) ); - TEST_CHECK( gkyl_compare(pos_map->constB_ctx->map_strength, 1.0, 1e-6) ); - TEST_CHECK( pos_map->constB_ctx->map_order_center == 2 ); - TEST_CHECK( pos_map->constB_ctx->map_order_expander == 3 ); - TEST_CHECK( pos_map->constB_ctx->N_theta_boundaries == 65 ); + TEST_CHECK(gkyl_compare(pos_map->constB_ctx->theta_throat, 1.565796, 1e-6) ); + TEST_CHECK(gkyl_compare(pos_map->constB_ctx->Bmag_throat, 1.093613, 1e-6) ); + TEST_CHECK(gkyl_compare(pos_map->constB_ctx->psi, 0.5, 1e-6) ); + TEST_CHECK(gkyl_compare(pos_map->constB_ctx->alpha, 0.0, 1e-6) ); + TEST_CHECK(gkyl_compare(pos_map->constB_ctx->map_strength, 1.0, 1e-6) ); + TEST_CHECK(pos_map->constB_ctx->map_order_center == 2); + TEST_CHECK(pos_map->constB_ctx->map_order_expander == 3); + TEST_CHECK(pos_map->constB_ctx->N_theta_boundaries == 65); gkyl_position_map_release(pos_map); gkyl_array_release(bmag_global); @@ -463,10 +474,10 @@ test_position_polynomial_map_optimize_1x() void test_position_map_numeric_optimize_1x() { - int cells[] = {64}; + int cells[] = { 64 }; int poly_order = 1; - double lower[] = {-M_PI+1e-2}, upper[] = {M_PI-1e-2}; - int dim = sizeof(lower)/sizeof(lower[0]); + double lower[] = { -M_PI + 1e-2 }, upper[] = { M_PI - 1e-2 }; + int dim = sizeof(lower) / sizeof(lower[0]); // Grids. struct gkyl_rect_grid grid; gkyl_rect_grid_init(&grid, dim, lower, upper, cells); @@ -474,7 +485,7 @@ test_position_map_numeric_optimize_1x() int ghost[] = { 1, 1 }; struct gkyl_range localRange, localRange_ext; // local, local-ext ranges. gkyl_create_grid_ranges(&grid, ghost, &localRange_ext, &localRange); - + // Basis functions. struct gkyl_basis basis; gkyl_cart_modal_serendip(&basis, dim, poly_order); @@ -484,20 +495,22 @@ test_position_map_numeric_optimize_1x() .map_strength = 1.0, }; - struct gkyl_position_map *pos_map = gkyl_position_map_new(pos_map_inp,\ + struct gkyl_position_map *pos_map = gkyl_position_map_new(pos_map_inp, \ grid, localRange, localRange_ext, localRange, localRange_ext, basis); // Project bmag_func onto bmag_global - struct gkyl_array *bmag_global = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, localRange_ext.volume); - gkyl_proj_on_basis *projB = gkyl_proj_on_basis_new(&grid, &basis, poly_order+1, 1, bmag_func, 0); + struct gkyl_array *bmag_global = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, + localRange_ext.volume); + gkyl_proj_on_basis *projB = gkyl_proj_on_basis_new(&grid, &basis, poly_order + 1, 1, bmag_func, + 0); gkyl_proj_on_basis_advance(projB, 0.0, &localRange, bmag_global); gkyl_proj_on_basis_release(projB); struct gkyl_rect_grid grid3D; - double lower3D[] = {0.4, -0.1, lower[0]}, upper3D[] = {0.6, 0.1, upper[0]}; - int cells3D[] = {1, 1, cells[0]}; + double lower3D[] = { 0.4, -0.1, lower[0] }, upper3D[] = { 0.6, 0.1, upper[0] }; + int cells3D[] = { 1, 1, cells[0] }; gkyl_rect_grid_init(&grid3D, 3, lower3D, upper3D, cells3D); - int ghost3D[] = { 1, 1 , 1}; + int ghost3D[] = { 1, 1, 1 }; struct gkyl_range localRange3D, localRange3D_ext; // local, local-ext ranges. gkyl_create_grid_ranges(&grid3D, ghost3D, &localRange3D_ext, &localRange3D); @@ -505,27 +518,31 @@ test_position_map_numeric_optimize_1x() gkyl_position_map_set_bmag(pos_map, NULL, bmag_global); gkyl_position_map_optimize(pos_map, grid3D, localRange3D); - double theta_extrema_analytic[5] = {lower[0], lower[0]/2, 0.0, upper[0]/2, upper[0]}; + double theta_extrema_analytic[5] = { lower[0], lower[0] / 2, 0.0, upper[0] / 2, upper[0] }; - TEST_CHECK( pos_map->constB_ctx->num_extrema == 5 ); - TEST_CHECK( gkyl_compare(pos_map->constB_ctx->theta_extrema[0], theta_extrema_analytic[0], 1e-15) ); - TEST_CHECK( gkyl_compare(pos_map->constB_ctx->theta_extrema[1], theta_extrema_analytic[1], 1e-15) ); - TEST_CHECK( gkyl_compare(pos_map->constB_ctx->theta_extrema[2], theta_extrema_analytic[2], 1e-15) ); - TEST_CHECK( gkyl_compare(pos_map->constB_ctx->theta_extrema[3], theta_extrema_analytic[3], 1e-15) ); - TEST_CHECK( gkyl_compare(pos_map->constB_ctx->theta_extrema[4], theta_extrema_analytic[4], 1e-15) ); + TEST_CHECK(pos_map->constB_ctx->num_extrema == 5); + TEST_CHECK(gkyl_compare(pos_map->constB_ctx->theta_extrema[0], theta_extrema_analytic[0], + 1e-15) ); + TEST_CHECK(gkyl_compare(pos_map->constB_ctx->theta_extrema[1], theta_extrema_analytic[1], + 1e-15) ); + TEST_CHECK(gkyl_compare(pos_map->constB_ctx->theta_extrema[2], theta_extrema_analytic[2], + 1e-15) ); + TEST_CHECK(gkyl_compare(pos_map->constB_ctx->theta_extrema[3], theta_extrema_analytic[3], + 1e-15) ); + TEST_CHECK(gkyl_compare(pos_map->constB_ctx->theta_extrema[4], theta_extrema_analytic[4], + 1e-15) ); gkyl_position_map_release(pos_map); gkyl_array_release(bmag_global); } - void test_position_map_numeric_calculate_1x() { - int cells[] = {64}; + int cells[] = { 64 }; int poly_order = 1; - double lower[] = {-M_PI+1e-2}, upper[] = {M_PI-1e-2}; - int dim = sizeof(lower)/sizeof(lower[0]); + double lower[] = { -M_PI + 1e-2 }, upper[] = { M_PI - 1e-2 }; + int dim = sizeof(lower) / sizeof(lower[0]); // Grids. struct gkyl_rect_grid grid; gkyl_rect_grid_init(&grid, dim, lower, upper, cells); @@ -533,7 +550,7 @@ test_position_map_numeric_calculate_1x() int ghost[] = { 1, 1 }; struct gkyl_range localRange, localRange_ext; // local, local-ext ranges. gkyl_create_grid_ranges(&grid, ghost, &localRange_ext, &localRange); - + // Basis functions. struct gkyl_basis basis; gkyl_cart_modal_serendip(&basis, dim, poly_order); @@ -543,20 +560,22 @@ test_position_map_numeric_calculate_1x() .map_strength = 1.0, }; - struct gkyl_position_map *pos_map = gkyl_position_map_new(pos_map_inp,\ + struct gkyl_position_map *pos_map = gkyl_position_map_new(pos_map_inp, \ grid, localRange, localRange_ext, localRange, localRange_ext, basis); // Project bmag_func onto bmag_global - struct gkyl_array *bmag_global = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, localRange_ext.volume); - gkyl_proj_on_basis *projB = gkyl_proj_on_basis_new(&grid, &basis, poly_order+1, 1, bmag_func, 0); + struct gkyl_array *bmag_global = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, + localRange_ext.volume); + gkyl_proj_on_basis *projB = gkyl_proj_on_basis_new(&grid, &basis, poly_order + 1, 1, bmag_func, + 0); gkyl_proj_on_basis_advance(projB, 0.0, &localRange, bmag_global); gkyl_proj_on_basis_release(projB); struct gkyl_rect_grid grid3D; - double lower3D[] = {0.4, -0.1, lower[0]}, upper3D[] = {0.6, 0.1, upper[0]}; - int cells3D[] = {1, 1, cells[0]}; + double lower3D[] = { 0.4, -0.1, lower[0] }, upper3D[] = { 0.6, 0.1, upper[0] }; + int cells3D[] = { 1, 1, cells[0] }; gkyl_rect_grid_init(&grid3D, 3, lower3D, upper3D, cells3D); - int ghost3D[] = { 1, 1 , 1}; + int ghost3D[] = { 1, 1, 1 }; struct gkyl_range localRange3D, localRange3D_ext; // local, local-ext ranges. gkyl_create_grid_ranges(&grid3D, ghost3D, &localRange3D_ext, &localRange3D); @@ -566,7 +585,7 @@ test_position_map_numeric_calculate_1x() double theta_map = 1.0; pos_map->maps[2](0.0, &theta_map, &theta_map, pos_map->ctxs[2]); - TEST_CHECK( gkyl_compare(theta_map, 1.505924, 1e-5) ); + TEST_CHECK(gkyl_compare(theta_map, 1.505924, 1e-5) ); gkyl_position_map_release(pos_map); gkyl_array_release(bmag_global); @@ -578,7 +597,7 @@ TEST_LIST = { { "test_position_map_init_2x", test_position_map_init_2x }, { "test_position_map_init_3x", test_position_map_init_3x }, { "test_position_map_set", test_position_map_set }, - { "test_gkyl_position_map_eval_mc2nu", test_gkyl_position_map_eval_mc2nu }, + { "test_gkyl_position_map_eval_mc2nu", test_gkyl_position_map_eval_mc2nu }, { "test_gkyl_position_map_slope", test_gkyl_position_map_slope }, { "test_position_polynomial_map_optimize_1x", test_position_polynomial_map_optimize_1x }, { "test_position_map_numeric_optimize_1x", test_position_map_numeric_optimize_1x }, diff --git a/gyrokinetic/zero/gkyl_position_map.h b/gyrokinetic/zero/gkyl_position_map.h index 8d88742e85..dc46563c41 100644 --- a/gyrokinetic/zero/gkyl_position_map.h +++ b/gyrokinetic/zero/gkyl_position_map.h @@ -15,13 +15,14 @@ enum gkyl_position_map_id { GKYL_PMAP_XPT_COMPRESSION, // Compresses cells near X-point (For use in MB Tokamaks) }; -typedef void (*mc2nu_t)(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, void *ctx); +typedef void (*mc2nu_t)(double t, const double *GKYL_RESTRICT xn, double *GKYL_RESTRICT fout, + void *ctx); struct gkyl_position_map_inp { enum gkyl_position_map_id id; - mc2nu_t maps[3]; // Position mapping in each position direction. This is defined in full 3x, + mc2nu_t maps[3]; // Position mapping in each position direction. This is defined in full 3x, // not in deflated coordinates. - mc2nu_t map_derivs[3]; // Derivative of mapping in each position direction. This is defined in full 3x, + mc2nu_t map_derivs[3]; // Derivative of mapping in each position direction. This is defined in full 3x, // not in deflated coordinates. void *ctxs[3]; // Context for each position mapping function. double map_strength; // Zero is uniform mapping, one is fully nonuniform mapping. How strong the nonuniformity is @@ -62,14 +63,14 @@ struct gkyl_position_map { // Stuff for constant B mapping struct gkyl_bmag_ctx *bmag_ctx; // Context for magnetic field calculation struct gkyl_position_map_const_B_ctx *constB_ctx; // Context for constant B mapping - struct gkyl_position_map_xpt_ctx *xpt_ctx; // Context for X-point compression mapping + struct gkyl_position_map_xpt_ctx *xpt_ctx; // Context for X-point compression mapping }; struct gkyl_position_map_const_B_ctx { mc2nu_t maps_backup[3]; // Backup of the position mapping functions. void *ctxs_backup[3]; // Backup of the context for each position mapping function. - - double psi, alpha; // The psi and alpha values for the middle flux surface to identify the 1D line we are optimizing + + double psi, alpha; // The psi and alpha values for the middle flux surface to identify the 1D line we are optimizing double psi_min, psi_max; // The max and min psi values for the simulation double alpha_min, alpha_max; // The max and min alpha values for the simulation double theta_min, theta_max; // The max and min theta values for the simulation @@ -107,8 +108,8 @@ struct gkyl_position_map_xpt_ctx { }; /** - * Create a new position map object. A position map is a function that maps - * uniform computational coordinates to non-uniform coordinates in the + * Create a new position map object. A position map is a function that maps + * uniform computational coordinates to non-uniform coordinates in the * same coordinate space as computational coordinates. (e.g. uniform field * aligned -> non-uniform field aligned). * @@ -121,7 +122,7 @@ struct gkyl_position_map_xpt_ctx { * @return New position map object. */ struct gkyl_position_map* gkyl_position_map_new(struct gkyl_position_map_inp pmap_info, - struct gkyl_rect_grid grid, struct gkyl_range local, struct gkyl_range local_ext, + struct gkyl_rect_grid grid, struct gkyl_range local, struct gkyl_range local_ext, struct gkyl_range global, struct gkyl_range global_ext, struct gkyl_basis basis); /** @@ -132,7 +133,6 @@ struct gkyl_position_map* gkyl_position_map_new(struct gkyl_position_map_inp pma struct gkyl_position_map* gkyl_position_map_inew(struct gkyl_position_map_inew_inp inp); - /** Create a new null position map object. This is a position map that does nothing. * All maps are identity maps. * @return New null position map object. @@ -142,13 +142,13 @@ gkyl_position_map_null_new(); /** * Set the position map object. Copy the non-uniform map array to the position map object. - * + * * @param gpm Position map object. * @param mc2nu Position map array. - * + * * @note This function is used to set the position map array in the position map object. */ -void gkyl_position_map_set_mc2nu(struct gkyl_position_map* gpm, struct gkyl_array* mc2nu); +void gkyl_position_map_set_mc2nu(struct gkyl_position_map *gpm, struct gkyl_array *mc2nu); /** * Set the magnetic field array in the position map object. This is used to set the magnetic field @@ -159,12 +159,12 @@ void gkyl_position_map_set_mc2nu(struct gkyl_position_map* gpm, struct gkyl_arra * @param bmag Magnetic field array. */ void -gkyl_position_map_set_bmag(struct gkyl_position_map* gpm, struct gkyl_comm* comm, - struct gkyl_array* bmag); +gkyl_position_map_set_bmag(struct gkyl_position_map *gpm, struct gkyl_comm *comm, + struct gkyl_array *bmag); /** * Set the function paramters for the map object. - * + * * @param gpm Position map object. * @param zcut half wavelength of sinusoidal mapping. * @param zcenter location of largest cells. @@ -172,8 +172,8 @@ gkyl_position_map_set_bmag(struct gkyl_position_map* gpm, struct gkyl_comm* comm * @param psisep separatrix psi value. */ void -gkyl_position_map_set_compression(struct gkyl_position_map* gpm, double zcut, - double zcenter, double w, double psisep); +gkyl_position_map_set_compression(struct gkyl_position_map *gpm, double zcut, + double zcenter, double w, double psisep); /** * Evaluate the position mapping at a specific computational (position) coordinate. @@ -184,11 +184,11 @@ gkyl_position_map_set_compression(struct gkyl_position_map* gpm, double zcut, * @param xnu Resulting non-uniform position coordinates. */ void -gkyl_position_map_eval_mc2nu(const struct gkyl_position_map* gpm, const double *xc, double *xnu); +gkyl_position_map_eval_mc2nu(const struct gkyl_position_map *gpm, const double *xc, double *xnu); /** * Evaluate the slope of the position mapping at a specific computational (position) coordinate. - * + * * @param gpm Gkyl position map object. * @param ix_map Index of the map to evaluate. Calls gpm->maps[index]. * @param x Computational position coordinates. @@ -198,7 +198,7 @@ gkyl_position_map_eval_mc2nu(const struct gkyl_position_map* gpm, const double * * @return Slope of the position mapping. */ double -gkyl_position_map_slope(const struct gkyl_position_map* gpm, int ix_map, +gkyl_position_map_slope(const struct gkyl_position_map *gpm, int ix_map, double x, double dx, int ix_comp, const struct gkyl_range *nrange); /** @@ -207,22 +207,21 @@ gkyl_position_map_slope(const struct gkyl_position_map* gpm, int ix_map, * * @param gpm Position map object. */ -struct gkyl_position_map* gkyl_position_map_acquire(const struct gkyl_position_map* gpm); +struct gkyl_position_map* gkyl_position_map_acquire(const struct gkyl_position_map *gpm); /** * Optimize the position map object for constant B mapping. - * + * * @param gpm Position map object. * @param grid 3D Position space grid. * @param global 3D Global position range. */ -void gkyl_position_map_optimize(struct gkyl_position_map* gpm, struct gkyl_rect_grid grid, +void gkyl_position_map_optimize(struct gkyl_position_map *gpm, struct gkyl_rect_grid grid, struct gkyl_range global); - /** * Release pointer to (and eventually memory associated with) - * the position map object. + * the position map object. * * @param Position map object. */ diff --git a/gyrokinetic/zero/gkyl_position_map_priv.h b/gyrokinetic/zero/gkyl_position_map_priv.h index ab3b3ae3c5..dd73dbe3a5 100644 --- a/gyrokinetic/zero/gkyl_position_map_priv.h +++ b/gyrokinetic/zero/gkyl_position_map_priv.h @@ -2,8 +2,7 @@ #include // Context for numeric root finding B mapping -struct opt_Theta_ctx -{ +struct opt_Theta_ctx { struct gkyl_position_map *gpm; struct gkyl_bmag_ctx *bmag_ctx; double dB_target; // How much B should change in 1 cell @@ -24,19 +23,20 @@ static void gkyl_position_map_free(const struct gkyl_ref_count *ref); /** * Calculates the location of the throat of the mirror in the theta direction * and the value of Bmag at that location. - * + * * @param constB_ctx Context for the constant B mapping * @param bmag_ctx Context for the magnetic field calculation */ static void -calculate_mirror_throat_location_polynomial(struct gkyl_position_map_const_B_ctx *constB_ctx, struct gkyl_bmag_ctx *bmag_ctx) +calculate_mirror_throat_location_polynomial(struct gkyl_position_map_const_B_ctx *constB_ctx, + struct gkyl_bmag_ctx *bmag_ctx) { // Parameters to use for the midpoint rule root finding algorithm to find the throat of the mirror int itterations = 10; int points_per_level = 20; // Assumes symmetry along theta, centered at 0, and two local maxima in Bmag that are symmetric - enum { X_IDX, Y_IDX, Z_IDX }; // arrangement of cartesian coordinates + enum {X_IDX, Y_IDX, Z_IDX}; // arrangement of cartesian coordinates double psi = constB_ctx->psi; double alpha = constB_ctx->alpha; double xp[3]; @@ -48,19 +48,16 @@ calculate_mirror_throat_location_polynomial(struct gkyl_position_map_const_B_ctx double maximum_Bmag = 0.0; double maximum_Bmag_location = 0.0; double fout[3]; - for (int j = 0; j < itterations; j++) - { + for (int j = 0; j < itterations; j++) { double dz = (interval_right - interval_left) / points_per_level; maximum_Bmag = 0.0; maximum_Bmag_location = 0.0; - for (int i = 0; i < points_per_level; i++) - { + for (int i = 0; i < points_per_level; i++) { double z = interval_left + i * dz; xp[Z_IDX] = z; gkyl_calc_bmag_global(0.0, xp, fout, bmag_ctx); double Bmag = fout[0]; - if (Bmag > maximum_Bmag) - { + if (Bmag > maximum_Bmag) { maximum_Bmag = Bmag; maximum_Bmag_location = z; } @@ -76,7 +73,7 @@ calculate_mirror_throat_location_polynomial(struct gkyl_position_map_const_B_ctx * Converts our uniform coordinate along field line length to a non-uniform coordinate * according to the polynomial mapping of arbitrary order. * Notation: We switch from theta to z here. Both are the third computational coordinate. - * + * * @param t Time * @param xn Uniform coordinate * @param fout Non-uniform coordinate @@ -95,36 +92,31 @@ position_map_constB_z_polynomial(double t, const double *xn, double *fout, void double uniform_coordinate = xn[0]; double nonuniform_coordinate, left, right; int n; - if (uniform_coordinate >= z_min && uniform_coordinate <= z_max) - { - if (uniform_coordinate <= -z_m) - { + if (uniform_coordinate >= z_min && uniform_coordinate <= z_max) { + if (uniform_coordinate <= -z_m) { left = -z_m; right = z_min; n = n_ex; } - else if (uniform_coordinate <= 0.0) - { + else if (uniform_coordinate <= 0.0) { left = -z_m; right = 0.0; n = n_ct; } - else if (uniform_coordinate <= z_m) - { + else if (uniform_coordinate <= z_m) { left = z_m; right = 0.0; n = n_ct; } - else - { + else { left = z_m; right = z_max; n = n_ex; } - nonuniform_coordinate = (pow(right - left, 1 - n) * pow(uniform_coordinate - left, n) + left) * frac + uniform_coordinate * (1 - frac); + nonuniform_coordinate = (pow(right - left, 1 - n) * pow(uniform_coordinate - left, + n) + left) * frac + uniform_coordinate * (1 - frac); } - else - { + else { nonuniform_coordinate = uniform_coordinate; } fout[0] = nonuniform_coordinate; @@ -133,16 +125,17 @@ position_map_constB_z_polynomial(double t, const double *xn, double *fout, void /** * Calculates the optimal orders for the polynomail mapping of the B field. * The optimal orders are the orders that minimize the maximum dB/dTheta in that region. - * + * * @param constB_ctx Context for the constant B mapping * @param bmag_ctx Context for the magnetic field calculation */ static void -calculate_optimal_mapping_polynomial(struct gkyl_position_map_const_B_ctx *constB_ctx, struct gkyl_bmag_ctx *bmag_ctx) +calculate_optimal_mapping_polynomial(struct gkyl_position_map_const_B_ctx *constB_ctx, + struct gkyl_bmag_ctx *bmag_ctx) { // Could be refined further by doing midpoint root finding for maximum dB/dz // Expander region - enum { X_IDX, Y_IDX, Z_IDX }; // arrangement of cartesian coordinates + enum {X_IDX, Y_IDX, Z_IDX}; // arrangement of cartesian coordinates double psi = constB_ctx->psi; double alpha = constB_ctx->alpha; double xp[3]; @@ -158,12 +151,10 @@ calculate_optimal_mapping_polynomial(struct gkyl_position_map_const_B_ctx *const double max_dB_dCell_prior = 99999999.99; double max_dB_dCell; double max_dB_dCell_order1 = 0.0; - while (1) - { + while (1) { max_dB_dCell = 0.0; constB_ctx->map_order_expander = expander_order; - for (int iz = 0; iz < scan_cells; iz++) - { + for (int iz = 0; iz < scan_cells; iz++) { double left_xi = scan_left + iz * scan_dxi; double right_xi = scan_left + (iz + 1) * scan_dxi; double psi = constB_ctx->psi; @@ -180,42 +171,35 @@ calculate_optimal_mapping_polynomial(struct gkyl_position_map_const_B_ctx *const gkyl_calc_bmag_global(0.0, xp, fout, bmag_ctx); double Bmag_right = fout[0]; double dB_dCell = (Bmag_right - Bmag_left); - if (fabs(dB_dCell) > max_dB_dCell) - { + if (fabs(dB_dCell) > max_dB_dCell) { max_dB_dCell = fabs(dB_dCell); } } double improvement = max_dB_dCell_prior - max_dB_dCell; - if (improvement > 1e-3) - { + if (improvement > 1e-3) { expander_order++; max_dB_dCell_prior = max_dB_dCell; } - else if (improvement < 0) - { + else if (improvement < 0) { expander_order--; constB_ctx->map_order_expander = expander_order; break; } - else - { + else { break; } - } double max_dB_dCell_expander = max_dB_dCell; - //Center region + // Center region scan_left = 0.0; scan_right = constB_ctx->theta_throat; scan_dxi = (scan_right - scan_left) / scan_cells; int center_order = 1; max_dB_dCell_prior = 99999999.99; - while (1) - { + while (1) { max_dB_dCell = 0.0; constB_ctx->map_order_center = center_order; - for (int iz = 0; iz < scan_cells; iz++) - { + for (int iz = 0; iz < scan_cells; iz++) { double left_xi = scan_left + iz * scan_dxi; double right_xi = scan_left + (iz + 1) * scan_dxi; @@ -231,37 +215,32 @@ calculate_optimal_mapping_polynomial(struct gkyl_position_map_const_B_ctx *const double Bmag_right = fout[0]; double dB_dCell = (Bmag_right - Bmag_left); - if (fabs(dB_dCell) > max_dB_dCell) - { + if (fabs(dB_dCell) > max_dB_dCell) { max_dB_dCell = fabs(dB_dCell); } } double improvement = max_dB_dCell_prior - max_dB_dCell; - if (improvement > 1e-3) - { + if (improvement > 1e-3) { center_order++; max_dB_dCell_prior = max_dB_dCell; } - else if (improvement < 0) - { + else if (improvement < 0) { center_order--; constB_ctx->map_order_center = center_order; break; } - else - { + else { break; } } } - // Utility functions for numeric root finding B mapping /** * Calculates dB/dTheta numerically at a given xn value. Calculates * the derivative to the left of the point. - * + * * @param theta The theta value to calculate the derivative at * @param ctx The context for the position map */ @@ -270,7 +249,8 @@ calc_bmag_global_derivative(double theta, void *ctx) { struct gkyl_position_map *gpm = ctx; struct gkyl_bmag_ctx *bmag_ctx = gpm->bmag_ctx; - double dtheta_cell = (gpm->constB_ctx->theta_max - gpm->constB_ctx->theta_min)/gpm->constB_ctx->N_theta_boundaries; + double dtheta_cell = (gpm->constB_ctx->theta_max - gpm->constB_ctx->theta_min) / + gpm->constB_ctx->N_theta_boundaries; double h = 1e-2 * dtheta_cell; double xh[3]; double fout[3]; @@ -279,7 +259,7 @@ calc_bmag_global_derivative(double theta, void *ctx) xh[2] = theta - h; gkyl_calc_bmag_global(0.0, xh, fout, bmag_ctx); double Bmag_plus = fout[0]; - xh[2] = theta - 2*h; + xh[2] = theta - 2 * h; gkyl_calc_bmag_global(0.0, xh, fout, bmag_ctx); double Bmag_minus = fout[0]; return (Bmag_plus - Bmag_minus) / (h); @@ -288,7 +268,7 @@ calc_bmag_global_derivative(double theta, void *ctx) /** * Finds the local min and max of the B field along the field line * specified by the input psi and alpha values. - * + * * @param gpm The position map object */ static void @@ -297,7 +277,7 @@ find_B_field_extrema(struct gkyl_position_map *gpm) // Assumes we are P1 in z, which means maxima and minima can only be in the center or edge of cells struct gkyl_position_map_const_B_ctx *constB_ctx = gpm->constB_ctx; struct gkyl_bmag_ctx *bmag_ctx = gpm->bmag_ctx; - enum { X_IDX, Y_IDX, Z_IDX }; // arrangement of cartesian coordinates + enum {X_IDX, Y_IDX, Z_IDX}; // arrangement of cartesian coordinates double psi = constB_ctx->psi; double alpha = constB_ctx->alpha; double xp[3]; @@ -314,46 +294,44 @@ find_B_field_extrema(struct gkyl_position_map *gpm) double *theta_extrema = gkyl_malloc(sizeof(double) * (npts + 1)); double *bmag_extrema = gkyl_malloc(sizeof(double) * (npts + 1)); - for (int i = 0; i <= npts; i++){ + for (int i = 0; i <= npts; i++) { double theta = theta_lo + i * theta_dxi; xp[Z_IDX] = theta; gkyl_calc_bmag_global(0.0, xp, &bmag_vals[i], bmag_ctx); dbmag_vals[i] = calc_bmag_global_derivative(theta, gpm); - if (i==0) continue; + if (i == 0) { + continue; + } // Minima - if (dbmag_vals[i] > 0 && dbmag_vals[i-1] < 0){ - if (bmag_vals[i] < bmag_vals[i-1]) - { + if (dbmag_vals[i] > 0 && dbmag_vals[i - 1] < 0) { + if (bmag_vals[i] < bmag_vals[i - 1]) { theta_extrema[extrema] = theta; bmag_extrema[extrema] = bmag_vals[i]; extrema++; } - else - { + else { theta_extrema[extrema] = theta - theta_dxi; - bmag_extrema[extrema] = bmag_vals[i-1]; + bmag_extrema[extrema] = bmag_vals[i - 1]; extrema++; } } // Maxima - if (dbmag_vals[i] < 0 && dbmag_vals[i-1] > 0){ - if (bmag_vals[i] > bmag_vals[i-1]) - { + if (dbmag_vals[i] < 0 && dbmag_vals[i - 1] > 0) { + if (bmag_vals[i] > bmag_vals[i - 1]) { theta_extrema[extrema] = theta; bmag_extrema[extrema] = bmag_vals[i]; extrema++; } - else - { + else { theta_extrema[extrema] = theta - theta_dxi; - bmag_extrema[extrema] = bmag_vals[i-1]; + bmag_extrema[extrema] = bmag_vals[i - 1]; extrema++; } } } - + // Set final extrema after the loop. MR April 22 2025 theta_extrema[0] = constB_ctx->theta_min; xp[Z_IDX] = constB_ctx->theta_min; @@ -365,8 +343,7 @@ find_B_field_extrema(struct gkyl_position_map *gpm) extrema++; gpm->constB_ctx->num_extrema = extrema; - for (int i = 0; i < extrema; i++) - { + for (int i = 0; i < extrema; i++) { gpm->constB_ctx->theta_extrema[i] = theta_extrema[i]; gpm->constB_ctx->bmag_extrema[i] = bmag_extrema[i]; } @@ -374,31 +351,39 @@ find_B_field_extrema(struct gkyl_position_map *gpm) // Identify 1 for maxima, 0 for minima // Left edge - if (bmag_extrema[0] > bmag_extrema[1]) - { gpm->constB_ctx->min_or_max[0] = 1; } // Maximum - else if (bmag_extrema[0] < bmag_extrema[1]) - { gpm->constB_ctx->min_or_max[0] = 0; } // Minimum - else - { printf("Error: Extrema is not an extrema. Position_map optimization failed\n"); } + if (bmag_extrema[0] > bmag_extrema[1]) { + gpm->constB_ctx->min_or_max[0] = 1; + } // Maximum + else if (bmag_extrema[0] < bmag_extrema[1]) { + gpm->constB_ctx->min_or_max[0] = 0; + } // Minimum + else { + printf("Error: Extrema is not an extrema. Position_map optimization failed\n"); + } // Middle points - for (int i = 1; i < extrema - 1; i++) - { - if (bmag_extrema[i] > bmag_extrema[i-1] && bmag_extrema[i] > bmag_extrema[i+1]) - { gpm->constB_ctx->min_or_max[i] = 1; } // Maximum - else if (bmag_extrema[i] < bmag_extrema[i-1] && bmag_extrema[i] < bmag_extrema[i+1]) - { gpm->constB_ctx->min_or_max[i] = 0; } // Minimum - else - { printf("Error: Extrema is not an extrema. Position_map optimization failed\n"); } + for (int i = 1; i < extrema - 1; i++) { + if (bmag_extrema[i] > bmag_extrema[i - 1] && bmag_extrema[i] > bmag_extrema[i + 1]) { + gpm->constB_ctx->min_or_max[i] = 1; + } // Maximum + else if (bmag_extrema[i] < bmag_extrema[i - 1] && bmag_extrema[i] < bmag_extrema[i + 1]) { + gpm->constB_ctx->min_or_max[i] = 0; + } // Minimum + else { + printf("Error: Extrema is not an extrema. Position_map optimization failed\n"); + } } // Right edge - if (bmag_extrema[extrema-1] > bmag_extrema[extrema-2]) - { gpm->constB_ctx->min_or_max[extrema-1] = 1; } // Maximum - else if (bmag_extrema[extrema-1] < bmag_extrema[extrema-2]) - { gpm->constB_ctx->min_or_max[extrema-1] = 0; } // Minimum - else - { printf("Error: Extrema is not an extrema. Position_map optimization failed\n"); } + if (bmag_extrema[extrema - 1] > bmag_extrema[extrema - 2]) { + gpm->constB_ctx->min_or_max[extrema - 1] = 1; + } // Maximum + else if (bmag_extrema[extrema - 1] < bmag_extrema[extrema - 2]) { + gpm->constB_ctx->min_or_max[extrema - 1] = 0; + } // Minimum + else { + printf("Error: Extrema is not an extrema. Position_map optimization failed\n"); + } // Free mallocs gkyl_free(bmag_vals); @@ -410,7 +395,7 @@ find_B_field_extrema(struct gkyl_position_map *gpm) /** * Refines the extrema found in the B field along the field line * specified by the input psi and alpha values. - * + * * @param gpm The position map object */ static void @@ -421,7 +406,7 @@ refine_B_field_extrema(struct gkyl_position_map *gpm) struct gkyl_position_map_const_B_ctx *constB_ctx = gpm->constB_ctx; struct gkyl_bmag_ctx *bmag_ctx = gpm->bmag_ctx; - enum { X_IDX, Y_IDX, Z_IDX }; // arrangement of cartesian coordinates + enum {X_IDX, Y_IDX, Z_IDX}; // arrangement of cartesian coordinates double psi = constB_ctx->psi; double alpha = constB_ctx->alpha; double xp[3]; @@ -432,8 +417,7 @@ refine_B_field_extrema(struct gkyl_position_map *gpm) double theta_hi = constB_ctx->theta_max; double theta_dxi = (theta_hi - theta_lo) / npts; - for (int i = 1; i < gpm->constB_ctx->num_extrema - 1; i++) - { + for (int i = 1; i < gpm->constB_ctx->num_extrema - 1; i++) { double theta = gpm->constB_ctx->theta_extrema[i]; xp[Z_IDX] = theta; double bmag_cent, bmag_left, bmag_right; @@ -449,38 +433,38 @@ refine_B_field_extrema(struct gkyl_position_map *gpm) double extrema_Bmag_location; double bmag_out; bool is_maximum; - if (bmag_cent > bmag_left && bmag_cent > bmag_right) - { is_maximum = true; } // Local maxima - else if (bmag_cent < bmag_left && bmag_cent < bmag_right) - { is_maximum = false; } // Local minima - else - { printf("Error: Extrema is not an extrema. Position_map optimization failed\n"); + if (bmag_cent > bmag_left && bmag_cent > bmag_right) { + is_maximum = true; + } // Local maxima + else if (bmag_cent < bmag_left && bmag_cent < bmag_right) { + is_maximum = false; + } // Local minima + else { + printf("Error: Extrema is not an extrema. Position_map optimization failed\n"); break; } // Midpoint rule refinement - for (int j = 0; j < num_iterations; j++) - { + for (int j = 0; j < num_iterations; j++) { double dz = (interval_right - interval_left) / num_points_per_level; - if (is_maximum) - { extrema_Bmag = 0.0; } - else - { extrema_Bmag = 99999999999999999.; } + if (is_maximum) { + extrema_Bmag = 0.0; + } + else { + extrema_Bmag = 99999999999999999.; + } extrema_Bmag_location = 0.0; - for (int k = 0; k <= num_points_per_level; k++) - { + for (int k = 0; k <= num_points_per_level; k++) { double z = interval_left + k * dz; xp[Z_IDX] = z; gkyl_calc_bmag_global(0.0, xp, &bmag_out, bmag_ctx); - if (is_maximum && bmag_out > extrema_Bmag) - { + if (is_maximum && bmag_out > extrema_Bmag) { extrema_Bmag = bmag_out; extrema_Bmag_location = z; } - else if (!is_maximum && bmag_out < extrema_Bmag) - { + else if (!is_maximum && bmag_out < extrema_Bmag) { extrema_Bmag = bmag_out; extrema_Bmag_location = z; } @@ -494,9 +478,8 @@ refine_B_field_extrema(struct gkyl_position_map *gpm) // Find the change in B over each cell double B_total_change = 0.0; // Total change in magnetic field - for (int i = 1; i < gpm->constB_ctx->num_extrema; i++) - { - B_total_change += fabs(gpm->constB_ctx->bmag_extrema[i] - gpm->constB_ctx->bmag_extrema[i-1]); + for (int i = 1; i < gpm->constB_ctx->num_extrema; i++) { + B_total_change += fabs(gpm->constB_ctx->bmag_extrema[i] - gpm->constB_ctx->bmag_extrema[i - 1]); } gpm->constB_ctx->dB_cell = B_total_change / (gpm->constB_ctx->N_theta_boundaries); } @@ -504,7 +487,7 @@ refine_B_field_extrema(struct gkyl_position_map *gpm) /** * Function used for root finding to determine the optimal theta value * for the numeric constant dB mapping. - * + * * @param theta The theta value to evaluate * @param ctx The context for the root finder. Type opt_Theta_ctx */ @@ -535,7 +518,7 @@ position_map_numeric_optimization_function(double theta, void *ctx) /** * Maps the uniform computational coordinate to a non-uniform coordinate * according to the numeric constant B mapping. - * + * * @param t Time * @param xn Uniform coordinate * @param fout Non-uniform coordinate @@ -560,8 +543,7 @@ position_map_constB_z_numeric(double t, const double *xn, double *fout, void *ct // Set strict floor and ceiling limits for theta // This is to prevent the root finding algorithm from going out of bounds // Not fout[0] = theta because of the finite differences and can lead to jumps - if (it <= 0 || it >= num_boundaries) - { + if (it <= 0 || it >= num_boundaries) { fout[0] = (it <= 0) ? theta_lo : theta_hi; return; } @@ -571,14 +553,11 @@ position_map_constB_z_numeric(double t, const double *xn, double *fout, void *ct // Initial guess is not accurate because the theta_extrema are not Theta_extrema // We use itteration to further refine this, but it's a good initial guess int region = 0; - for (int i = 1; i <= num_extrema-2; i++) - { - if (theta >= theta_extrema[i]) - { + for (int i = 1; i <= num_extrema - 2; i++) { + if (theta >= theta_extrema[i]) { region = i; } - else - { + else { break; } } @@ -592,12 +571,11 @@ position_map_constB_z_numeric(double t, const double *xn, double *fout, void *ct dB_target = dB_cell * it; bool outside_region = true; // Asuume that we identified the region incorrectly - while (outside_region) - { + while (outside_region) { dB_global_lower = 0.0; - for (int i = 0; i < region; i++) - { - dB_global_lower += fabs(gpm->constB_ctx->bmag_extrema[i+1] - gpm->constB_ctx->bmag_extrema[i]); + for (int i = 0; i < region; i++) { + dB_global_lower += fabs(gpm->constB_ctx->bmag_extrema[i + 1] - + gpm->constB_ctx->bmag_extrema[i]); } B_lower_region = gpm->constB_ctx->bmag_extrema[region]; @@ -606,7 +584,7 @@ position_map_constB_z_numeric(double t, const double *xn, double *fout, void *ct ridders_ctx.B_lower_region = B_lower_region; interval_lower = theta_extrema[region]; - interval_upper = theta_extrema[region+1]; + interval_upper = theta_extrema[region + 1]; interval_lower_eval = position_map_numeric_optimization_function(interval_lower, &ridders_ctx); interval_upper_eval = position_map_numeric_optimization_function(interval_upper, &ridders_ctx); @@ -619,7 +597,8 @@ position_map_constB_z_numeric(double t, const double *xn, double *fout, void *ct // Just use the corresponding endpoint if (fabs(interval_lower_eval) < fabs(interval_upper_eval)) { fout[0] = interval_lower; - } else { + } + else { fout[0] = interval_upper; } return; @@ -638,7 +617,7 @@ position_map_constB_z_numeric(double t, const double *xn, double *fout, void *ct else if (interval_lower_eval < 0.0 && interval_upper_eval < 0.0) { // If the bounds on the interval are both negative, we should move up a region to make it pass through zero region++; - if (region > num_extrema-2) { + if (region > num_extrema - 2) { // If we can't move up any regions and leave the simulation domain, we are likely on the upper limit of the domain and should just return the input theta fout[0] = theta_hi; return; @@ -655,7 +634,8 @@ position_map_constB_z_numeric(double t, const double *xn, double *fout, void *ct return; } else { - fprintf(stderr, "Warning: Unexpected interval evaluation state in position_map_constB_z_numeric. Using theta directly.\n"); + fprintf(stderr, + "Warning: Unexpected interval evaluation state in position_map_constB_z_numeric. Using theta directly.\n"); fout[0] = theta; return; } @@ -665,28 +645,25 @@ position_map_constB_z_numeric(double t, const double *xn, double *fout, void *ct struct gkyl_qr_res res = gkyl_ridders(position_map_numeric_optimization_function, &ridders_ctx, interval_lower, interval_upper, interval_lower_eval, interval_upper_eval, 10, 1e-6); double Theta = res.res; - fout[0] = Theta*gpm->constB_ctx->map_strength + theta*(1-gpm->constB_ctx->map_strength); + fout[0] = Theta * gpm->constB_ctx->map_strength + theta * (1 - gpm->constB_ctx->map_strength); bool enable_limits_min_B = gpm->constB_ctx->enable_maximum_slope_limits_at_min_B; bool enable_limits_max_B = gpm->constB_ctx->enable_maximum_slope_limits_at_max_B; - if (enable_limits_min_B || enable_limits_max_B) - { + if (enable_limits_min_B || enable_limits_max_B) { // Set a minimum cell size on the edges // Assume that at inflection points, Theta = theta. This should be true - double Theta_left = interval_lower; + double Theta_left = interval_lower; double Theta_right = interval_upper; double theta_middle = 0.5 * (interval_lower + interval_upper); bool left_is_maximum = gpm->constB_ctx->min_or_max[region]; - bool right_is_maximum = gpm->constB_ctx->min_or_max[region+1]; + bool right_is_maximum = gpm->constB_ctx->min_or_max[region + 1]; - if (theta > theta_middle && left_is_maximum) - { + if (theta > theta_middle && left_is_maximum) { enable_limits_max_B = false; } - if (theta < theta_middle && right_is_maximum) - { + if (theta < theta_middle && right_is_maximum) { enable_limits_max_B = false; } @@ -694,39 +671,36 @@ position_map_constB_z_numeric(double t, const double *xn, double *fout, void *ct double max_slope_max_B = gpm->constB_ctx->maximum_slope_at_max_B; double right_straight_line_value, left_straight_line_value; - if (left_is_maximum){ - left_straight_line_value = max_slope_max_B * theta + (1-max_slope_max_B) * Theta_left; + if (left_is_maximum) { + left_straight_line_value = max_slope_max_B * theta + (1 - max_slope_max_B) * Theta_left; } else { - left_straight_line_value = max_slope_min_B * theta + (1-max_slope_min_B) * Theta_left; + left_straight_line_value = max_slope_min_B * theta + (1 - max_slope_min_B) * Theta_left; } - if (right_is_maximum){ - right_straight_line_value = max_slope_max_B * theta + (1-max_slope_max_B) * Theta_right; + if (right_is_maximum) { + right_straight_line_value = max_slope_max_B * theta + (1 - max_slope_max_B) * Theta_right; } else { - right_straight_line_value = max_slope_min_B * theta + (1-max_slope_min_B) * Theta_right; + right_straight_line_value = max_slope_min_B * theta + (1 - max_slope_min_B) * Theta_right; } - if ( fout[0] < right_straight_line_value && - ((right_is_maximum && enable_limits_max_B) || - ((!right_is_maximum) && enable_limits_min_B))) - { + if (fout[0] < right_straight_line_value && + ((right_is_maximum && enable_limits_max_B) || + ((!right_is_maximum) && enable_limits_min_B))) { fout[0] = right_straight_line_value; } - if (fout[0] > left_straight_line_value && + if (fout[0] > left_straight_line_value && ((left_is_maximum && enable_limits_max_B) || - ((!left_is_maximum) && enable_limits_min_B))) - { + ((!left_is_maximum) && enable_limits_min_B))) { fout[0] = left_straight_line_value; } } } // Context for Gaussian-weighted integration -struct gaussian_weight_ctx -{ +struct gaussian_weight_ctx { struct gkyl_position_map *gpm; double theta_c; // Center point for Gaussian double wd2; // Half-width of averaging window @@ -739,11 +713,11 @@ position_map_constB_z_numeric_dbl_exp_wrapper(double z, void *ctx) struct gaussian_weight_ctx *gw_ctx = ctx; double fout[3]; position_map_constB_z_numeric(0.0, &z, fout, gw_ctx->gpm); - + // Apply Gaussian weight: exp(-(z-theta_c)^2 / (2*sigma^2)) double dz = z - gw_ctx->theta_c; double weight = exp(-dz * dz / (2.0 * gw_ctx->sigma * gw_ctx->sigma)); - + return fout[0] * weight; } @@ -751,18 +725,18 @@ double gaussian_norm_wrapper(double z, void *ctx) { struct gaussian_weight_ctx *gw_ctx = ctx; - + // Return just the Gaussian weight for normalization double dz = z - gw_ctx->theta_c; double weight = exp(-dz * dz / (2.0 * gw_ctx->sigma * gw_ctx->sigma)); - + return weight; } /** * Maps the uniform computational coordinate to a non-uniform coordinate * according to the numeric constant B mapping. - * + * * @param t Time * @param xn Uniform coordinate * @param fout Non-uniform coordinate @@ -772,8 +746,7 @@ static void position_map_constB_z_numeric_moving_average(double t, const double *xn, double *fout, void *ctx) { struct gkyl_position_map *gpm = ctx; - if (gpm->constB_ctx->gaussian_std == 0.0) - { + if (gpm->constB_ctx->gaussian_std == 0.0) { position_map_constB_z_numeric(t, xn, fout, ctx); return; } @@ -783,20 +756,19 @@ position_map_constB_z_numeric_moving_average(double t, const double *xn, double const double max_width = gpm->constB_ctx->gaussian_max_integration_width; const double tmin = gpm->constB_ctx->theta_min; const double tmax = gpm->constB_ctx->theta_max; - + // Shrink the half-width symmetrically to stay within bounds // This ensures the integration window is always centered at theta_c double dist_to_min = theta_c - tmin; double dist_to_max = tmax - theta_c; - double wd2 = fmin(fmin(dist_to_min, dist_to_max) * 0.99, max_width/2); - + double wd2 = fmin(fmin(dist_to_min, dist_to_max) * 0.99, max_width / 2); + // If the symmetric window is too small, fall back to unsmoothed - if (wd2 < 1e-6) - { + if (wd2 < 1e-6) { position_map_constB_z_numeric(t, xn, fout, ctx); return; } - + double rng_lo = theta_c - wd2; double rng_up = theta_c + wd2; @@ -825,7 +797,7 @@ position_map_constB_z_numeric_moving_average(double t, const double *xn, double /** * Converts our uniform coordinate along field line length to a non-uniform coordinate * according to a sinusoidal mapping with a specified compression factor at the ends - * + * * @param t Time * @param xn Uniform coordinate * @param fout Non-uniform coordinate @@ -835,19 +807,20 @@ static void position_map_xpt_compression(double t, const double *xn, double *fout, void *ctx) { struct gkyl_position_map_xpt_ctx *app = ctx; - double uniform_coordinate = xn[0]; + double uniform_coordinate = xn[0]; double F = 1.0 / (1.0 - app->compression_factor); - double A = 1.0/F; + double A = 1.0 / F; double zcut = app->zcut; double zshift = uniform_coordinate - app->zcenter; - double nonuniform_coordinate = A * (sin(M_PI*zshift/zcut)*zcut/M_PI + F*zshift) + app->zcenter; + double nonuniform_coordinate = A * (sin(M_PI * zshift / zcut) * zcut / M_PI + F * zshift) + + app->zcenter; fout[0] = nonuniform_coordinate; } /** * Converts our uniform coordinate psi to a non-uniform coordinate * according to a sinusoidal mapping with a specified compression factor at the ends - * + * * @param t Time * @param xn Uniform coordinate * @param fout Non-uniform coordinate @@ -857,19 +830,20 @@ static void position_map_sep_compression(double t, const double *xn, double *fout, void *ctx) { struct gkyl_position_map_xpt_ctx *app = ctx; - double uniform_coordinate = xn[0]; + double uniform_coordinate = xn[0]; double F = 1.0 / (1.0 - app->compression_factor); - double A = 1.0/F; + double A = 1.0 / F; double w = app->w; double xshift = uniform_coordinate - app->psisep; - double nonuniform_coordinate = A * (-sin(M_PI*xshift/w)*w/M_PI + F*xshift) + app->psisep; + double nonuniform_coordinate = A * (-sin(M_PI * xshift / w) * w / M_PI + F * xshift) + + app->psisep; fout[0] = nonuniform_coordinate; } /** - * Evaluates the derivative of the nonuniform coordinate wrt the uniform coordinate + * Evaluates the derivative of the nonuniform coordinate wrt the uniform coordinate * according to a sinusoidal mapping with a specified compression factor at the ends - * + * * @param t Time * @param xn Uniform coordinate * @param fout Non-uniform coordinate @@ -879,19 +853,19 @@ static void position_map_deriv_xpt_compression(double t, const double *xn, double *fout, void *ctx) { struct gkyl_position_map_xpt_ctx *app = ctx; - double uniform_coordinate = xn[0]; + double uniform_coordinate = xn[0]; double F = 1.0 / (1.0 - app->compression_factor); - double A = 1.0/F; + double A = 1.0 / F; double zcut = app->zcut; double zshift = uniform_coordinate - app->zcenter; - double deriv = A * (cos(M_PI*zshift/zcut) + F); + double deriv = A * (cos(M_PI * zshift / zcut) + F); fout[0] = deriv; } /** - * Evaluates the derivative of the nonuniform coordinate wrt the uniform coordinate + * Evaluates the derivative of the nonuniform coordinate wrt the uniform coordinate * according to a sinusoidal mapping with a specified compression factor at the ends - * + * * @param t Time * @param xn Uniform coordinate * @param fout Non-uniform coordinate @@ -901,12 +875,11 @@ static void position_map_deriv_sep_compression(double t, const double *xn, double *fout, void *ctx) { struct gkyl_position_map_xpt_ctx *app = ctx; - double uniform_coordinate = xn[0]; + double uniform_coordinate = xn[0]; double F = 1.0 / (1.0 - app->compression_factor); - double A = 1.0/F; + double A = 1.0 / F; double w = app->w; double xshift = uniform_coordinate - app->psisep; - double deriv = A * (-cos(M_PI*xshift/w) + F); + double deriv = A * (-cos(M_PI * xshift / w) + F); fout[0] = deriv; } - diff --git a/gyrokinetic/zero/position_map.c b/gyrokinetic/zero/position_map.c index 47b0e95230..e814ea7882 100644 --- a/gyrokinetic/zero/position_map.c +++ b/gyrokinetic/zero/position_map.c @@ -25,7 +25,6 @@ gkyl_position_map_identity_slope(double t, const double *xn, double *fout, void fout[0] = 1.0; } - struct gkyl_position_map* gkyl_position_map_null_new() { @@ -40,8 +39,8 @@ gkyl_position_map_null_new() gpm->bmag_ctx = gkyl_malloc(sizeof(struct gkyl_bmag_ctx)); gpm->bmag_ctx->bmag = gkyl_array_new(GKYL_DOUBLE, 1, 1); gpm->ref_count = gkyl_ref_count_init(gkyl_position_map_free); - - for (int i = 0; i < 3; i++){ + + for (int i = 0; i < 3; i++) { gpm->maps[i] = gkyl_position_map_identity; gpm->map_derivs[i] = gkyl_position_map_identity_slope; gpm->ctxs[i] = 0; @@ -62,7 +61,8 @@ gkyl_position_map_inew(struct gkyl_position_map_inew_inp inp) struct gkyl_position_map* gkyl_position_map_new(struct gkyl_position_map_inp pmap_info, struct gkyl_rect_grid grid, - struct gkyl_range local, struct gkyl_range local_ext, struct gkyl_range global, struct gkyl_range global_ext, + struct gkyl_range local, struct gkyl_range local_ext, struct gkyl_range global, + struct gkyl_range global_ext, struct gkyl_basis basis) { struct gkyl_position_map *gpm = gkyl_malloc(sizeof(*gpm)); @@ -71,13 +71,13 @@ gkyl_position_map_new(struct gkyl_position_map_inp pmap_info, struct gkyl_rect_g gpm->bmag_ctx = gkyl_malloc(sizeof(struct gkyl_bmag_ctx)); gpm->bmag_ctx->bmag = gkyl_array_new(GKYL_DOUBLE, basis.num_basis, global_ext.volume); gpm->to_optimize = false; - gpm->use_map_derivs = (pmap_info.id == GKYL_PMAP_XPT_COMPRESSION || pmap_info.id == GKYL_PMAP_USER_INPUT_W_DERIVATIVE) ? true : false; - + gpm->use_map_derivs = (pmap_info.id == GKYL_PMAP_XPT_COMPRESSION || + pmap_info.id == GKYL_PMAP_USER_INPUT_W_DERIVATIVE) ? true : false; gpm->constB_ctx = gkyl_malloc(sizeof(struct gkyl_position_map_const_B_ctx)); gpm->xpt_ctx = gkyl_malloc(sizeof(struct gkyl_position_map_xpt_ctx)); - for (int i = 0; i < 3; i++){ + for (int i = 0; i < 3; i++) { gpm->maps[i] = gkyl_position_map_identity; gpm->map_derivs[i] = gkyl_position_map_identity_slope; gpm->ctxs[i] = 0; @@ -87,20 +87,19 @@ gkyl_position_map_new(struct gkyl_position_map_inp pmap_info, struct gkyl_rect_g gpm->xpt_ctx->ctxs_backup[i] = 0; } - switch (pmap_info.id) - { + switch (pmap_info.id) { case GKYL_PMAP_USER_INPUT: - for (int i = 0; i < 3; i++){ - if (pmap_info.maps[i] != 0) - { gpm->maps[i] = pmap_info.maps[i]; + for (int i = 0; i < 3; i++) { + if (pmap_info.maps[i] != 0) { + gpm->maps[i] = pmap_info.maps[i]; gpm->ctxs[i] = pmap_info.ctxs[i]; } } case GKYL_PMAP_USER_INPUT_W_DERIVATIVE: - for (int i = 0; i < 3; i++){ - if (pmap_info.maps[i] != 0) - { gpm->maps[i] = pmap_info.maps[i]; + for (int i = 0; i < 3; i++) { + if (pmap_info.maps[i] != 0) { + gpm->maps[i] = pmap_info.maps[i]; gpm->map_derivs[i] = pmap_info.map_derivs[i]; gpm->ctxs[i] = pmap_info.ctxs[i]; } @@ -108,9 +107,9 @@ gkyl_position_map_new(struct gkyl_position_map_inp pmap_info, struct gkyl_rect_g case GKYL_PMAP_CONSTANT_DB_POLYNOMIAL: - for (int i = 0; i < 2; i++){ - if (pmap_info.maps[i] != 0) - { gpm->constB_ctx->maps_backup[i] = pmap_info.maps[i]; + for (int i = 0; i < 2; i++) { + if (pmap_info.maps[i] != 0) { + gpm->constB_ctx->maps_backup[i] = pmap_info.maps[i]; gpm->constB_ctx->ctxs_backup[i] = pmap_info.ctxs[i]; } } @@ -118,33 +117,37 @@ gkyl_position_map_new(struct gkyl_position_map_inp pmap_info, struct gkyl_rect_g case GKYL_PMAP_CONSTANT_DB_NUMERIC: - for (int i = 0; i < 2; i++){ - if (pmap_info.maps[i] != 0) - { gpm->constB_ctx->maps_backup[i] = pmap_info.maps[i]; + for (int i = 0; i < 2; i++) { + if (pmap_info.maps[i] != 0) { + gpm->constB_ctx->maps_backup[i] = pmap_info.maps[i]; gpm->constB_ctx->ctxs_backup[i] = pmap_info.ctxs[i]; } } gpm->constB_ctx->map_strength = pmap_info.map_strength; - if (pmap_info.maximum_slope_at_min_B == 0.) - { gpm->constB_ctx->enable_maximum_slope_limits_at_min_B = false; } - else - { gpm->constB_ctx->enable_maximum_slope_limits_at_min_B = true; } + if (pmap_info.maximum_slope_at_min_B == 0.) { + gpm->constB_ctx->enable_maximum_slope_limits_at_min_B = false; + } + else { + gpm->constB_ctx->enable_maximum_slope_limits_at_min_B = true; + } gpm->constB_ctx->maximum_slope_at_min_B = pmap_info.maximum_slope_at_min_B; - if (pmap_info.maximum_slope_at_max_B == 0.) - { gpm->constB_ctx->enable_maximum_slope_limits_at_max_B = false; } - else - { gpm->constB_ctx->enable_maximum_slope_limits_at_max_B = true; } + if (pmap_info.maximum_slope_at_max_B == 0.) { + gpm->constB_ctx->enable_maximum_slope_limits_at_max_B = false; + } + else { + gpm->constB_ctx->enable_maximum_slope_limits_at_max_B = true; + } gpm->constB_ctx->maximum_slope_at_max_B = pmap_info.maximum_slope_at_max_B; gpm->constB_ctx->gaussian_std = pmap_info.gaussian_std; gpm->constB_ctx->gaussian_max_integration_width = pmap_info.gaussian_max_integration_width; case GKYL_PMAP_XPT_COMPRESSION: - for (int i = 0; i < 2; i++){ - if (pmap_info.maps[i] != 0) - { gpm->xpt_ctx->maps_backup[i] = pmap_info.maps[i]; + for (int i = 0; i < 2; i++) { + if (pmap_info.maps[i] != 0) { + gpm->xpt_ctx->maps_backup[i] = pmap_info.maps[i]; gpm->xpt_ctx->ctxs_backup[i] = pmap_info.ctxs[i]; } } @@ -159,8 +162,8 @@ gkyl_position_map_new(struct gkyl_position_map_inp pmap_info, struct gkyl_rect_g gpm->global = global; gpm->global_ext = global_ext; gpm->basis = basis; - gpm->cdim = grid.ndim; - gpm->mc2nu = gkyl_array_new(GKYL_DOUBLE, 3*gpm->basis.num_basis, gpm->local_ext.volume); + gpm->cdim = grid.ndim; + gpm->mc2nu = gkyl_array_new(GKYL_DOUBLE, 3 * gpm->basis.num_basis, gpm->local_ext.volume); gpm->ref_count = gkyl_ref_count_init(gkyl_position_map_free); struct gkyl_position_map *gpm_out = gpm; @@ -168,14 +171,14 @@ gkyl_position_map_new(struct gkyl_position_map_inp pmap_info, struct gkyl_rect_g } void -gkyl_position_map_set_mc2nu(struct gkyl_position_map* gpm, struct gkyl_array* mc2nu) +gkyl_position_map_set_mc2nu(struct gkyl_position_map *gpm, struct gkyl_array *mc2nu) { gkyl_array_copy(gpm->mc2nu, mc2nu); } void -gkyl_position_map_set_bmag(struct gkyl_position_map* gpm, struct gkyl_comm* comm, - struct gkyl_array* bmag) +gkyl_position_map_set_bmag(struct gkyl_position_map *gpm, struct gkyl_comm *comm, + struct gkyl_array *bmag) { gpm->to_optimize = true; int N_boundaries = gpm->constB_ctx->N_theta_boundaries; @@ -189,19 +192,20 @@ gkyl_position_map_set_bmag(struct gkyl_position_map* gpm, struct gkyl_comm* comm } else { gkyl_comm_array_allgather_host(comm, &gpm->local, \ - &gpm->global, bmag, (struct gkyl_array*) gpm->bmag_ctx->bmag); + &gpm->global, bmag, (struct gkyl_array *)gpm->bmag_ctx->bmag); } } void -gkyl_position_map_set_compression(struct gkyl_position_map* gpm, double zcut, double zcenter, double w, double psisep) +gkyl_position_map_set_compression(struct gkyl_position_map *gpm, double zcut, double zcenter, + double w, double psisep) { gpm->xpt_ctx->zcut = zcut; gpm->xpt_ctx->zcenter = zcenter; gpm->xpt_ctx->w = w; gpm->xpt_ctx->psisep = psisep; - if (gpm->xpt_ctx->radial_compression_factor!=0.0) { + if (gpm->xpt_ctx->radial_compression_factor != 0.0) { gpm->maps[0] = position_map_sep_compression; gpm->map_derivs[0] = position_map_deriv_sep_compression; gpm->ctxs[0] = gpm->xpt_ctx; @@ -214,7 +218,7 @@ gkyl_position_map_set_compression(struct gkyl_position_map* gpm, double zcut, do gpm->maps[1] = gpm->xpt_ctx->maps_backup[1]; gpm->ctxs[1] = gpm->xpt_ctx->ctxs_backup[1]; - if (gpm->xpt_ctx->compression_factor!=0.0) { + if (gpm->xpt_ctx->compression_factor != 0.0) { gpm->maps[2] = position_map_xpt_compression; gpm->map_derivs[2] = position_map_deriv_xpt_compression; gpm->ctxs[2] = gpm->xpt_ctx; @@ -225,12 +229,14 @@ gkyl_position_map_set_compression(struct gkyl_position_map* gpm, double zcut, do } } -void -gkyl_position_map_eval_mc2nu(const struct gkyl_position_map* gpm, const double *x_comp, double *x_fa) +void +gkyl_position_map_eval_mc2nu(const struct gkyl_position_map *gpm, const double *x_comp, + double *x_fa) { int cidx[GKYL_MAX_CDIM]; - for(int i = 0; i < gpm->grid.ndim; i++){ - int idxtemp = gpm->global.lower[i] + (int) floor((x_comp[i] - (gpm->grid.lower[i]) )/gpm->grid.dx[i]); + for (int i = 0; i < gpm->grid.ndim; i++) { + int idxtemp = gpm->global.lower[i] + + (int)floor((x_comp[i] - (gpm->grid.lower[i]) ) / gpm->grid.dx[i]); idxtemp = GKYL_MAX2(GKYL_MIN2(idxtemp, gpm->local.upper[i]), gpm->local.lower[i]); cidx[i] = idxtemp; } @@ -239,34 +245,33 @@ gkyl_position_map_eval_mc2nu(const struct gkyl_position_map* gpm, const double * double cxc[gpm->grid.ndim]; double x_log[gpm->grid.ndim]; gkyl_rect_grid_cell_center(&gpm->grid, cidx, cxc); - for(int i = 0; i < gpm->grid.ndim; i++){ - x_log[i] = (x_comp[i]-cxc[i])/(gpm->grid.dx[i]*0.5); + for (int i = 0; i < gpm->grid.ndim; i++) { + x_log[i] = (x_comp[i] - cxc[i]) / (gpm->grid.dx[i] * 0.5); } double xyz_fa[3]; - for(int i = 0; i < 3; i++){ - xyz_fa[i] = gpm->basis.eval_expand(x_log, &pmap_coeffs[i*gpm->basis.num_basis]); + for (int i = 0; i < 3; i++) { + xyz_fa[i] = gpm->basis.eval_expand(x_log, &pmap_coeffs[i * gpm->basis.num_basis]); } - for (int i=0; igrid.ndim; i++) { + for (int i = 0; i < gpm->grid.ndim; i++) { x_fa[i] = xyz_fa[i]; } - x_fa[gpm->grid.ndim-1] = xyz_fa[2]; + x_fa[gpm->grid.ndim - 1] = xyz_fa[2]; } void -gkyl_position_map_optimize(struct gkyl_position_map* gpm, struct gkyl_rect_grid grid, +gkyl_position_map_optimize(struct gkyl_position_map *gpm, struct gkyl_rect_grid grid, struct gkyl_range global) { - enum { PSI_IDX, AL_IDX, TH_IDX }; // arrangement of computational coordinates - gpm->constB_ctx->psi_max = grid.upper[PSI_IDX]; - gpm->constB_ctx->psi_min = grid.lower[PSI_IDX]; + enum {PSI_IDX, AL_IDX, TH_IDX}; // arrangement of computational coordinates + gpm->constB_ctx->psi_max = grid.upper[PSI_IDX]; + gpm->constB_ctx->psi_min = grid.lower[PSI_IDX]; gpm->constB_ctx->alpha_max = grid.upper[AL_IDX]; gpm->constB_ctx->alpha_min = grid.lower[AL_IDX]; gpm->constB_ctx->theta_max = grid.upper[TH_IDX]; gpm->constB_ctx->theta_min = grid.lower[TH_IDX]; gpm->constB_ctx->N_theta_boundaries = global.upper[TH_IDX] - global.lower[TH_IDX] + 2; - if (gpm->id == GKYL_PMAP_CONSTANT_DB_POLYNOMIAL && gpm->to_optimize == true) - { + if (gpm->id == GKYL_PMAP_CONSTANT_DB_POLYNOMIAL && gpm->to_optimize == true) { double psi_center = 0.5 * (gpm->constB_ctx->psi_min + gpm->constB_ctx->psi_max); double alpha_center = 0.5 * (gpm->constB_ctx->alpha_min + gpm->constB_ctx->alpha_max); @@ -281,15 +286,15 @@ gkyl_position_map_optimize(struct gkyl_position_map* gpm, struct gkyl_rect_grid gpm->bmag_ctx->cbasis = &gpm->basis; gpm->bmag_ctx->cgrid = &gpm->grid; - gpm->constB_ctx->psi = psi_center; - gpm->constB_ctx->alpha = alpha_center; + gpm->constB_ctx->psi = psi_center; + gpm->constB_ctx->alpha = alpha_center; calculate_mirror_throat_location_polynomial(gpm->constB_ctx, gpm->bmag_ctx); calculate_optimal_mapping_polynomial(gpm->constB_ctx, gpm->bmag_ctx); } - else if (gpm->id == GKYL_PMAP_CONSTANT_DB_NUMERIC && gpm->to_optimize == true) - { - double psi_center = pow(0.5 * (sqrt(gpm->constB_ctx->psi_min) + sqrt(gpm->constB_ctx->psi_max)), 2.0); + else if (gpm->id == GKYL_PMAP_CONSTANT_DB_NUMERIC && gpm->to_optimize == true) { + double psi_center = pow(0.5 * (sqrt(gpm->constB_ctx->psi_min) + sqrt(gpm->constB_ctx->psi_max)), + 2.0); double alpha_center = 0.5 * (gpm->constB_ctx->alpha_min + gpm->constB_ctx->alpha_max); gpm->maps[0] = gpm->constB_ctx->maps_backup[0]; @@ -300,11 +305,11 @@ gkyl_position_map_optimize(struct gkyl_position_map* gpm, struct gkyl_rect_grid gpm->ctxs[2] = gpm; gpm->bmag_ctx->crange_global = &gpm->global; - gpm->bmag_ctx->cbasis = &gpm->basis; - gpm->bmag_ctx->cgrid = &gpm->grid; + gpm->bmag_ctx->cbasis = &gpm->basis; + gpm->bmag_ctx->cgrid = &gpm->grid; - gpm->constB_ctx->psi = psi_center; - gpm->constB_ctx->alpha = alpha_center; + gpm->constB_ctx->psi = psi_center; + gpm->constB_ctx->alpha = alpha_center; find_B_field_extrema(gpm); refine_B_field_extrema(gpm); @@ -312,11 +317,10 @@ gkyl_position_map_optimize(struct gkyl_position_map* gpm, struct gkyl_rect_grid } double -gkyl_position_map_slope(const struct gkyl_position_map* gpm, int ix_map, +gkyl_position_map_slope(const struct gkyl_position_map *gpm, int ix_map, double x, double dx, int ix_comp, const struct gkyl_range *nrange) { - if (gpm->use_map_derivs) - { + if (gpm->use_map_derivs) { double slope; gpm->map_derivs[ix_map](0.0, &x, &slope, gpm->ctxs[ix_map]); return slope; @@ -328,28 +332,25 @@ gkyl_position_map_slope(const struct gkyl_position_map* gpm, int ix_map, gpm->maps[ix_map](0.0, &x_left, &f_left, gpm->ctxs[ix_map]); gpm->maps[ix_map](0.0, &x_right, &f_right, gpm->ctxs[ix_map]); double slope; - if (ix_comp == nrange->lower[ix_map]) - { + if (ix_comp == nrange->lower[ix_map]) { gpm->maps[ix_map](0.0, &x, &f, gpm->ctxs[ix_map]); slope = (f_right - f) / dx; } - else if (ix_comp == nrange->upper[ix_map]) - { + else if (ix_comp == nrange->upper[ix_map]) { gpm->maps[ix_map](0.0, &x, &f, gpm->ctxs[ix_map]); slope = (f - f_left) / dx; } - else - { + else { slope = (f_right - f_left) / (2.0 * dx); } return slope; } struct gkyl_position_map* -gkyl_position_map_acquire(const struct gkyl_position_map* gpm) +gkyl_position_map_acquire(const struct gkyl_position_map *gpm) { gkyl_ref_count_inc(&gpm->ref_count); - return (struct gkyl_position_map*) gpm; + return (struct gkyl_position_map *)gpm; } void @@ -364,8 +365,7 @@ gkyl_position_map_free(const struct gkyl_ref_count *ref) struct gkyl_position_map *gpm = container_of(ref, struct gkyl_position_map, ref_count); gkyl_array_release(gpm->mc2nu); gkyl_array_release(gpm->bmag_ctx->bmag); - if (gpm->to_optimize == true) - { + if (gpm->to_optimize == true) { gkyl_free(gpm->constB_ctx->theta_extrema); gkyl_free(gpm->constB_ctx->bmag_extrema); gkyl_free(gpm->constB_ctx->min_or_max);