Skip to content

Commit

Permalink
Spherical harmonics nn (#40)
Browse files Browse the repository at this point in the history 
* added new ansatz

* debug new ansatz

* branch should be abandoned here

* debug new ansatz

* fix scaling error new ansatz

* fix config errors

* temporarily fixed dynamic ansatz rotation bug

* fix inheritance error for new ansatz

* mesh acceleration

* add structured hohlraum

* added spherical harmonics basis generator

* parallelize datagen

* added flag for dynamic ansatz

* added rotated monomial basis

* fixed compilation errors

* fixed rotated basis linesource

* added rotated hohlraum testcase

* code cleanup

* restructure old files

* added regularized monomial m1 networks

* added spherical harmonics inference and started tests on rotation

* added sh models

* added rotated sh networks

* beautifier

* added rotated sh 2d and 3d models. EXPERIMENTALLY

* re-add mirroring

* did dynamic ansatz after rotation

* added symmetry correction

* fixed 2d sh projection

* added 2d sh to datagenerator

* remove rotation from data-gen

* update 2d mn solver

* delete redundant logs

* add rotation to 2d data-gen

* removed debugging output

* parallelize sampler again

* added 2d sh rotator up do degree 4

* fixed rotation matrix

* test run for sh2d rotated m1 closure

* fix compile error

* attempt to fix vector size bug

* add option to turn of rotation

* change m1_2d model

* change m1_2d models

* delete tf history files

* bugfix: u_2 was manually set to 0 in 2d sh network evaluation - removed line

* bugfix: dynamic closure was missing in non-rotated 2d monomial and sh NN closures. Removed rotation correction for m1 monomial NN closure

* bugfix: false naming of 2d m1 sh model

* bugfix: vector length for rotated sh neural network closure

* renamed some vectors

* fix: formatting, BUGFIX: rotation angle

* bugfix: made several function arguments in optimizer const for type safety

* enhance: made rotated datagen mirrored

* ongoing: fix rotated mn bug for spherical harmonics

* bufgfix: Rotated m1 works now

* enhance: rotated m1 models added

* bugfix: removed column deletion for higher order rotated models

* formatting

* modelchange: m1_2d_rot

* change: dockerfile

* enhance: added sh m2 models. change: newton optimizer stopping criterion to u recons. change: mn_normalized_ml solver

* change: default values for newton

* change: changed error measure of neural cells to relative norm error

* enhance: add M3 spherical harmonics models

* added measure for inaccuraate nn cells

---------

Co-authored-by: Steffen <[email protected]>
  • Loading branch information
@ScSteffen
ScSteffen and Steffen authored Nov 15, 2023
1 parent 22a5c21 commit 2cd1be9
Show file tree
Hide file tree
Showing 524 changed files with 134,292 additions and 62,973 deletions.
147 changes: 50 additions & 97 deletions examples/meshes/create_structured_hohlraum_2d_mesh.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,131 +5,84 @@
from optparse import OptionParser


def add_block(x0, y0, x_len, y_len, char_length, geom):
def add_block(x0,y0,lengthX, lengthY,char_length,geom):
coords = np.array([
[x0, y0, 0.0],
[x0 + x_len, y0, 0.0],
[x0 + x_len, y0 + y_len, 0.0],
[x0, y0 + y_len, 0.0]
[x0+lengthX, y0, 0.0],
[x0+lengthX, y0+lengthY, 0.0],
[x0, y0+lengthY, 0.0]
])
return geom.add_polygon(coords, char_length)



def main():
print("---------- Start Creating the mesh ------------")
print("Parsing options")
# --- parse options ---
parser = OptionParser()
parser.add_option("-o", "--output_name", dest="output_name", default="hohlraum")
parser.add_option("-c", "--char_length", dest="char_length", default=0.025)
parser.add_option("-o", "--output_name", dest="output_name", default="struct_2dmesh")
parser.add_option("-c", "--char_length", dest="char_length", default=0.01)
parser.add_option("-s", "--start_pt", dest="start_pt", nargs=2, default=(0,0))
parser.add_option("-l", "--length", dest="length", nargs=2, default=(1,1))
parser.add_option("-b", "--boundary", dest="b_type", default="void")
(options, args) = parser.parse_args()

options.output_name = str(options.output_name)
options.char_length = float(options.char_length)
options.b_type = str(options.b_type)
char_length = options.char_length
lengthX = float(options.length[0])
lengthY = float(options.length[1])
x0 = float(options.start_pt[0])
y0 = float(options.start_pt[1])

geom = pg.opencascade.Geometry()
domain = add_block(-0.1, -0.05, 1.45, 1.4, char_length, geom)
domain = add_block(x0,y0,lengthX, lengthY,char_length,geom)

boxes = [domain]
lines = []
if options.b_type == "void":
geom.add_physical(domain.lines, label="void")

# add interior polygon
coords = np.array([
[0.85, 0.25, 0.0],
[0.45, 0.25, 0.0],
[0.45, 1.05, 0.0],
[0.85, 1.05, 0.0],
[0.85, 1.0, 0.0],
[0.5, 1.0, 0.0],
[0.5, 0.3, 0.0],
[0.85, 0.3, 0.0],
])
boxes.append(geom.add_polygon(coords, char_length))
elif options.b_type == "voidAdiabatic":
adiabatic = list()
void = list()
adiabatic.append(domain.lines[0])
adiabatic.append(domain.lines[2])

# add outer polygon
coords = np.array([
[0.0, 0.0, 0.0],
[1.3, 0.0, 0.0],
[1.3, 1.3, 0.0],
[0.0, 1.3, 0.0],
[0.0, 1.25, 0.0],
[1.25, 1.25, 0.0],
[1.25, 0.05, 0.0],
[0.0, 0.05, 0.0],
])
boxes.append(geom.add_polygon(coords, char_length))

"""
points = []
points.append(geom.add_point([1.0, 1.05, 0.0], lcar=char_length))
points.append(geom.add_point([1.05, 1.05, 0.0], lcar=char_length))
points.append(geom.add_point([1.05, 1.0, 0.0], lcar=char_length))
points.append(geom.add_point([1.0, 1.0, 0.0], lcar=char_length))
lines = []
lines.append(geom.add_line(points[0], points[1]))
lines.append(geom.add_line(points[1], points[2]))
lines.append(geom.add_line(points[2], points[3]))
lines.append(geom.add_line(points[3], points[0]))
lineloop = geom.add_line_loop(lines)
surf = geom.add_surface(lineloop)
geom.add_physical(surf)
"""
# points = []
# points.append(geom.add_point([0.0, 0.0, 0.0], lcar=char_length))
# points.append(geom.add_point([1.3, 0.0, 0.0], lcar=char_length))
# points.append(geom.add_point([1.3, 1.3, 0.0], lcar=char_length))
# points.append(geom.add_point([0.0, 1.3, 0.0], lcar=char_length))
# points.append(geom.add_point([0.0, 1.25, 0.0], lcar=char_length))
# points.append(geom.add_point([0.05, 1.25, 0.0], lcar=char_length))
# points.append(geom.add_point([1.25, 1.25, 0.0], lcar=char_length))
# points.append(geom.add_point([1.25, 1.25, 0.0], lcar=char_length))
# points.append(geom.add_point([1.25, 0.05, 0.0], lcar=char_length))
# points.append(geom.add_point([0.05, 0.05, 0.0], lcar=char_length))
# points.append(geom.add_point([0.0, 0.05, 0.0], lcar=char_length))
#
# for i in range(len(points) - 1):
# lines.append(geom.add_line(points[i], points[i + 1]))
# lines.append(geom.add_line(points[len(points) - 1], points[0]))
#
# lineloop = geom.add_line_loop(lines)
# surf = geom.add_surface(lineloop)
# geom.add_physical(surf)

# add left side absorption region
boxes.append(add_block(0, 0.25, 0.05, 0.8, char_length, geom))
boxes.append(add_block(-0.05, 0.05, 0.05, 0.2, char_length, geom))
boxes.append(add_block(-0.05, 1.05, 0.05, 0.2, char_length, geom))

# add line that denotes the inflow ghost cells
# coords = np.array([
# [-0.01, -0.01, 0.0],
# [-0.01, 1.31, 0.0]])

# point1 = geom.add_point((0.0, 0.05, 0.0), lcar=char_length)
# point2 = geom.add_point((0.0, 0.25, 0.0), lcar=char_length)
# t = geom.add_line(point1, point2)
# geom.add_physical_line([t], label="void")

# lines.append(t)

geom.boolean_fragments(boxes, [])
geom.add_physical(domain.lines, label="void")

# domain.lines.append(t)
# geom.add_raw_code("Transfinite Surface{:};")
# geom.add_raw_code('Mesh.Algorithm= 3;')
# geom.add_raw_code("Recombine Surface {:} = 0;")
void.append(domain.lines[1])
void.append(domain.lines[3])

geom.add_physical(adiabatic, label="adiabatic")
geom.add_physical(void, label="void")

elif options.b_type == "dirichletNeumann":
dirichlet = list()
wall_low = list()
wall_up = list()
wall_low.append(domain.lines[0])
wall_up.append(domain.lines[2])

dirichlet.append(domain.lines[1])
dirichlet.append(domain.lines[3])

geom.add_physical(dirichlet, label="dirichlet")
geom.add_physical(wall_low, label="wall_low")
geom.add_physical(wall_up, label="wall_up")
print("| Neumann marker has tag wall_low and wall_up")
print("| Dirichlet marker has tag dirichlet")

else:
print("Boundary condition not yet implemented")

geom.add_raw_code("Transfinite Surface{:};")
geom.add_raw_code('Mesh.Algorithm= 3;')
geom.add_raw_code("Recombine Surface {:} = 0;")

mesh_code = geom.get_code()
with open(options.output_name + ".geo", "w") as mesh_file:
mesh_file.write(mesh_code)
os.system('gmsh ' + options.output_name + '.geo -2 -format su2 -save_all')
# os.system('rm ' + options.output_name + '.geo')
os.system('rm ' + options.output_name + '.geo')
print("---------- Successfully created the mesh ------------")


Expand Down
Loading

0 comments on commit 2cd1be9

Please sign in to comment.