bspline working ok?

This commit is contained in:
Anton Ljungdahl 2024-10-12 16:29:06 +02:00
parent ed31ba62ef
commit dd89a7b6de
8 changed files with 1504 additions and 15 deletions

Binary file not shown.

View File

@ -0,0 +1,11 @@
{
"folders":
[
{
"path": "."
},
{
"path": "src"
}
]
}

View File

@ -0,0 +1,304 @@
{
"auto_complete":
{
"selected_items":
[
]
},
"buffers":
[
{
"contents": "#include <stdlib.h>\n\n// ---\n// Header includes\n#include \"base/base_inc.h\"\n#include \"os/os_inc.h\"\n\n\n// CUDA headers\n#include \"kernels.h\"\n\n// ---\n// .C includes\n#include \"base/base_inc.c\"\n#include \"os/os_inc.c\"\n#include \"os/os_entry_point.c\"\n\n\n#define grid_file_path_bin \"D:\\\\dev\\\\eigsol_gpu\\\\out\\\\grid.bin\"\n#define grid_file_path \"D:\\\\dev\\\\eigsol_gpu\\\\out\\\\grid.dat\"\n\nglobal U32 enable_logging = 1;\n\n// Complex number with double precision\ntypedef struct Z64 Z64;\nstruct Z64\n{\n F64 re;\n F64 im;\n};\n\ntypedef struct Grid Grid;\nstruct Grid \n{\n F64 start;\n F64 end;\n U32 num_steps;\n F64 *points;\n};\n\n\ntypedef struct BSplineCtx BSplineCtx;\nstruct BSplineCtx \n{\n U32 order;\n U32 num_knotpoints;\n U32 num_bsplines;\n U32 num_phys_points;\n};\n\n\n\n//~ Globals\nglobal BSplineCtx bspline_ctx = {0};\n\n//~ Functions\n\nfunction inline void\nlogger_debug(String8 msg)\n{\n OutputDebugString(msg.str);\n}\n\n#define LOG(msg) if(enable_logging) { logger_debug(msg) };\n\n\nfunction void\nwrite_array_binary_F64(String8 path_to_file, F64 *values, U32 array_size) \n{\n OS_Handle file_handle = OS_file_open(OS_AccessFlag_Write | OS_AccessFlag_CreateNew, \n path_to_file);\n { \n ArenaTemp scratch = scratch_get(0, 0);\n String8List list = {0};\n String8 temp = {0};\n temp.str = (U8*)values;\n temp.size = sizeof(F64)*array_size;\n str8_list_push(scratch.arena, &list, temp);\n OS_file_write(scratch.arena, file_handle, 0, list, 0);\n \n String8List log_list = {0};\n str8_list_push(scratch.arena, &log_list, str8_lit(\"Wrote binary array data to\"));\n str8_list_push(scratch.arena, &log_list, path_to_file);\n StringJoin join = {0};\n join.sep = str8_lit(\" \");\n join.post = str8_lit(\"\\n\");\n String8 log_msg = str8_list_join(scratch.arena, log_list, &join);\n OutputDebugString(log_msg.str);\n }\n OS_file_close(file_handle);\n}\n\nfunction void\nwrite_array_F64(String8 path_to_file, F64 *values, U32 array_size, char* fmt) \n{\n OS_Handle file_handle = OS_file_open(OS_AccessFlag_Write | OS_AccessFlag_CreateNew, \n path_to_file);\n {\n ArenaTemp scratch = scratch_get(0, 0);\n String8List list = {0};\n for(U32 i = 0; i < array_size; i++) \n {\n str8_list_pushf(scratch.arena, &list, fmt, values[i]);\n }\n OS_file_write(scratch.arena, file_handle, 0, list, 0);\n \n String8List log_list = {0};\n str8_list_push(scratch.arena, &log_list, str8_lit(\"Wrote array to\"));\n str8_list_push(scratch.arena, &log_list, path_to_file);\n StringJoin join = {0};\n join.sep = str8_lit(\" \");\n join.post = str8_lit(\"\\n\");\n String8 log_msg = str8_list_join(scratch.arena, log_list, &join);\n OutputDebugString(log_msg.str);\n }\n OS_file_close(file_handle);\n \n}\n\n\n\nfunction\nvoid EntryPoint() \n{\n OS_InitReceipt os_receipt = OS_init();\n OS_InitGfxReceipt os_gfx_receipt = OS_gfx_init(os_receipt);\n\n Arena *arena = m_make_arena();\n\n Grid grid = {0};\n grid.start = 0.0;\n grid.end = 10.0;\n grid.num_steps = 1000;\n \n grid.points = PushArrayZero(arena, F64, grid.num_steps);\n F64 step_size = (grid.end-grid.start)/(F64)grid.num_steps;\n grid.points[0] = grid.start;\n grid.points[grid.num_steps-1] = grid.end;\n for(U32 i = 1; i < grid.num_steps-1; i++)\n {\n grid.points[i] = grid.points[i-1] + step_size;\n }\n\n write_array_binary_F64(str8_lit(grid_file_path_bin), grid.points, grid.num_steps);\n write_array_F64(str8_lit(grid_file_path), grid.points, grid.num_steps, \"%13.6e\\n\");\n \n // Create bsplines\n U32 k = 4;\n U32 N = 14;\n bspline_ctx.order = k; \n bspline_ctx.num_knotpoints = N; \n U32 num_bsplines = N-k;\n U32 num_phys_points = N-2*k+2; // Remove k points at each end, and then add back the first and last points.\n\n\n\n}\n\n\n",
"file": "src/main.c",
"file_size": 3834,
"file_write_time": 133732027810885494,
"settings":
{
"buffer_size": 3679,
"encoding": "UTF-8",
"line_ending": "Windows"
}
}
],
"build_system": "Packages/C++/C Single File.sublime-build",
"build_system_choices":
[
[
[
[
"Packages/C++/C Single File.sublime-build",
""
],
[
"Packages/C++/C Single File.sublime-build",
"Run"
]
],
[
"Packages/C++/C Single File.sublime-build",
""
]
]
],
"build_varint": "",
"command_palette":
{
"height": 0.0,
"last_filter": "",
"selected_items":
[
[
"View",
"View Package File"
],
[
"Color",
"Colorsublime: Install Theme"
],
[
"Pack",
"Package Control: Install Package"
],
[
"Install Pack",
"Package Control: Install Package"
],
[
"Package Insta",
"Package Control: Install Package"
],
[
"Color Sch",
"UI: Select Color Scheme"
],
[
"Install Pac",
"Package Control: Install Package"
],
[
"Package",
"Install Package Control"
],
[
"View ",
"View Package File"
]
],
"width": 0.0
},
"console":
{
"height": 0.0,
"history":
[
]
},
"distraction_free":
{
"menu_visible": true,
"show_minimap": false,
"show_open_files": false,
"show_tabs": false,
"side_bar_visible": false,
"status_bar_visible": false
},
"expanded_folders":
[
"/W/eigsol_gpu",
"/W/eigsol_gpu/src"
],
"file_history":
[
"/C/Users/anton/AppData/Roaming/Sublime Text 3/Packages/Default/Default (Windows).sublime-keymap",
"/C/sbs/sb1/java/guidesign/src/com/comsol/guidesign/views/%USER%",
"/C/Users/antonlj/AppData/Roaming/Sublime Text 3/Packages/Colorsublime/Colorsublime.sublime-settings",
"/C/Users/antonlj/AppData/Roaming/Sublime Text 3/Packages/User/Colorsublime.sublime-settings",
"/C/sbs/sb1/java/guidesign/src/com/comsol/guidesign/actions/DesignCreatorActionFactory.java",
"/C/sbs/sb1/java/design/src/com/comsol/design/operations/DesignOperation.java",
"/C/Users/antonlj/AppData/Roaming/Sublime Text 3/Packages/Default/symbol.py",
"/C/Program Files/Sublime Text 3/Packages/Default.sublime-package",
"/C/Users/antonlj/AppData/Roaming/Sublime Text 3/Packages/User/Preferences.sublime-settings",
"/C/sbs/sb1/java/testgui/src/com/comsol/testgui/builder/TBuilderFeatures.java",
"/C/Users/antonlj/AppData/Roaming/Sublime Text 3/Packages/Default/Preferences.sublime-settings"
],
"find":
{
"height": 40.0
},
"find_in_files":
{
"height": 103.333333333,
"where_history":
[
]
},
"find_state":
{
"case_sensitive": false,
"find_history":
[
"main.c"
],
"highlight": true,
"in_selection": false,
"preserve_case": false,
"regex": false,
"replace_history":
[
],
"reverse": false,
"show_context": true,
"use_buffer2": true,
"whole_word": false,
"wrap": true
},
"groups":
[
{
"selected": 0,
"sheets":
[
{
"buffer": 0,
"file": "src/main.c",
"semi_transient": false,
"settings":
{
"buffer_size": 3679,
"regions":
{
},
"selection":
[
[
2699,
2699
]
],
"settings":
{
"syntax": "Packages/C++/C.sublime-syntax",
"tab_size": 2,
"translate_tabs_to_spaces": true
},
"translation.x": 0.0,
"translation.y": 1552.0,
"zoom_level": 1.0
},
"stack_index": 0,
"type": "text"
}
]
}
],
"incremental_find":
{
"height": 25.0
},
"input":
{
"height": 0.0
},
"layout":
{
"cells":
[
[
0,
0,
1,
1
]
],
"cols":
[
0.0,
1.0
],
"rows":
[
0.0,
1.0
]
},
"menu_visible": true,
"output.exec":
{
"height": 267.0
},
"output.find_results":
{
"height": 0.0
},
"pinned_build_system": "",
"project": "eigsol_gpu.sublime-project",
"replace":
{
"height": 46.0
},
"save_all_on_build": true,
"select_file":
{
"height": 0.0,
"last_filter": "",
"selected_items":
[
[
"main",
"src\\main.c"
],
[
"BuilderInf",
"design\\src\\com\\comsol\\design\\util\\link\\BuilderInfoVisitor.java"
],
[
"DesignPhysicsM",
"design\\src\\com\\comsol\\design\\DesignPhysicsMaker.java"
],
[
"TBuilderFeat",
"testgui\\src\\com\\comsol\\testgui\\builder\\TBuilderFeatures.java"
]
],
"width": 0.0
},
"select_project":
{
"height": 0.0,
"last_filter": "",
"selected_items":
[
],
"width": 0.0
},
"select_symbol":
{
"height": 392.0,
"last_filter": "DesignPhysicsMak",
"selected_items":
[
[
"DesignPhysicsMak",
"DesignPhysicsMaker"
],
[
"DeviceModelFeatu",
"DeviceModelFeatureOperation"
]
],
"width": 644.0
},
"selected_group": 0,
"settings":
{
},
"show_minimap": false,
"show_open_files": false,
"show_tabs": true,
"side_bar_visible": false,
"side_bar_width": 145.0,
"status_bar_visible": true,
"template_settings":
{
}
}

1000
out/Bspline3.dat Normal file

File diff suppressed because it is too large Load Diff

35
out/knotpoints.dat Normal file
View File

@ -0,0 +1,35 @@
0.000000e+00
0.000000e+00
0.000000e+00
0.000000e+00
3.571429e-01
7.142857e-01
1.071429e+00
1.428571e+00
1.785714e+00
2.142857e+00
2.500000e+00
2.857143e+00
3.214286e+00
3.571429e+00
3.928571e+00
4.285714e+00
4.642857e+00
5.000000e+00
5.357143e+00
5.714286e+00
6.071429e+00
6.428571e+00
6.785714e+00
7.142857e+00
7.500000e+00
7.857143e+00
8.214286e+00
8.571429e+00
8.928571e+00
9.285714e+00
9.642857e+00
1.000000e+01
1.000000e+01
1.000000e+01
1.000000e+01

File diff suppressed because one or more lines are too long

View File

@ -18,22 +18,57 @@
#define grid_file_path_bin "D:\\dev\\eigsol_gpu\\out\\grid.bin"
#define grid_file_path "D:\\dev\\eigsol_gpu\\out\\grid.dat"
#define knotpoints_file_path "D:\\dev\\eigsol_gpu\\out\\knotpoints.dat"
global U32 enable_logging = 1;
// Complex number with double precision
typedef struct Z64 Z64;
struct Z64
{
F64 re;
F64 im;
};
typedef struct Grid Grid;
struct Grid {
struct Grid
{
F64 start;
F64 end;
U32 num_steps;
F64 *points;
};
typedef struct BSplineCtx BSplineCtx;
struct BSplineCtx {
struct BSplineCtx
{
Arena* arena;
U32 order;
F64 *knotpoints;
U32 num_knotpoints;
U32 num_bsplines;
U32 num_phys_points;
};
void write_array_binary_F64(String8 path_to_file, F64 *values, U32 array_size)
//~ Globals
global BSplineCtx g_bspline_ctx = {0};
//~ Functions
function inline void
logger_debug(String8 msg)
{
OutputDebugString(msg.str);
}
#define LOG(msg) if(enable_logging) { logger_debug(msg) };
function void
write_array_binary_F64(String8 path_to_file, F64 *values, U32 array_size)
{
OS_Handle file_handle = OS_file_open(OS_AccessFlag_Write | OS_AccessFlag_CreateNew,
path_to_file);
@ -45,12 +80,21 @@ void write_array_binary_F64(String8 path_to_file, F64 *values, U32 array_size)
temp.size = sizeof(F64)*array_size;
str8_list_push(scratch.arena, &list, temp);
OS_file_write(scratch.arena, file_handle, 0, list, 0);
String8List log_list = {0};
str8_list_push(scratch.arena, &log_list, str8_lit("Wrote binary array data to"));
str8_list_push(scratch.arena, &log_list, path_to_file);
StringJoin join = {0};
join.sep = str8_lit(" ");
join.post = str8_lit("\n");
String8 log_msg = str8_list_join(scratch.arena, log_list, &join);
OutputDebugString(log_msg.str);
}
OS_file_close(file_handle);
}
void write_array_F64(String8 path_to_file, F64 *values, U32 array_size, char* fmt)
function void
write_array_F64(String8 path_to_file, F64 *values, U32 array_size, char* fmt)
{
OS_Handle file_handle = OS_file_open(OS_AccessFlag_Write | OS_AccessFlag_CreateNew,
path_to_file);
@ -62,13 +106,56 @@ void write_array_F64(String8 path_to_file, F64 *values, U32 array_size, char* fm
str8_list_pushf(scratch.arena, &list, fmt, values[i]);
}
OS_file_write(scratch.arena, file_handle, 0, list, 0);
String8List log_list = {0};
str8_list_push(scratch.arena, &log_list, str8_lit("Wrote array to"));
str8_list_push(scratch.arena, &log_list, path_to_file);
StringJoin join = {0};
join.sep = str8_lit(" ");
join.post = str8_lit("\n");
String8 log_msg = str8_list_join(scratch.arena, log_list, &join);
OutputDebugString(log_msg.str);
}
OS_file_close(file_handle);
}
function F64 bspline_recursion(F64 x, U32 k, U32 i)
{
F64 *t = g_bspline_ctx.knotpoints;
if(k == 1)
{
if(i < g_bspline_ctx.num_bsplines && x >= t[i] && x < t[i+1])
{
return 1.0;
} else {
return 0.0;
}
}
else
{
F64 term1 = (x - t[i])/(t[i+k-1]-t[i])*bspline_recursion(x, k-1, i);
F64 term2 = (t[i+k]-x)/(t[i+k]-t[i+1])*bspline_recursion(x, k-1, i+1);
return term1 + term2;
}
}
void EntryPoint() {
function F64 get_bspline_F64(F64 x_coord, U32 index)
{
U32 k = g_bspline_ctx.order;
F64 out = bspline_recursion(x_coord, k, index);
return out;
}
function
void EntryPoint()
{
OS_InitReceipt os_receipt = OS_init();
OS_InitGfxReceipt os_gfx_receipt = OS_gfx_init(os_receipt);
@ -78,20 +165,68 @@ void EntryPoint() {
grid.start = 0.0;
grid.end = 10.0;
grid.num_steps = 1000;
grid.points = PushArrayZero(arena, F64, grid.num_steps);
F64 step_size = (grid.end-grid.start)/(F64)grid.num_steps;
grid.points[0] = grid.start;
grid.points[grid.num_steps-1] = grid.end;
for(U32 i = 1; i < grid.num_steps-1; i++) {
for(U32 i = 1; i < grid.num_steps-1; i++)
{
grid.points[i] = grid.points[i-1] + step_size;
}
write_array_binary_F64(str8_lit(grid_file_path_bin), grid.points, grid.num_steps);
write_array_F64(str8_lit(grid_file_path), grid.points, grid.num_steps, "%13.6e\n");
// Create knotpoint sequence.
U32 k = 4;
U32 N = 35;
g_bspline_ctx.order = k;
g_bspline_ctx.num_knotpoints = N;
g_bspline_ctx.num_bsplines = N-k;
g_bspline_ctx.num_phys_points = N-(2*k)+2; // Remove k points at each end, and then add back the first and last points.
g_bspline_ctx.arena = arena;
g_bspline_ctx.knotpoints = (F64*)PushArrayZero(arena, F64, g_bspline_ctx.num_knotpoints);
// Set up physical points;
F64 delta = (grid.end-grid.start)/(g_bspline_ctx.num_phys_points-1);
// Set ghost points including first physical
U32 phys_point_index = g_bspline_ctx.num_phys_points + k-1;
for(U32 i = 0; i < k; i++)
{
g_bspline_ctx.knotpoints[i] = grid.start;
}
for(U32 i = k; i < phys_point_index; i++)
{
g_bspline_ctx.knotpoints[i] = g_bspline_ctx.knotpoints[i-1] + delta;
}
// Set the last points
F64 last_physical = grid.end;
for(U32 i = phys_point_index; i < g_bspline_ctx.num_knotpoints; i++)
{
g_bspline_ctx.knotpoints[i] = last_physical;
}
write_array_F64(str8_lit(knotpoints_file_path), g_bspline_ctx.knotpoints, g_bspline_ctx.num_knotpoints, "%13.6e\n");
String8 Bspline3_fn = str8_lit("D:\\dev\\eigsol_gpu\\out\\Bspline3.dat");
F64 *Bspline3 = (F64 *)PushArrayZero(arena, F64, grid.num_steps);
for(U32 i = 0; i < grid.num_steps; i++)
{
Bspline3[i] = get_bspline_F64(grid.points[i], 3);
}
write_array_F64(Bspline3_fn, Bspline3, grid.num_steps, "%13.6e\n");
}

Binary file not shown.