usnistgov / fipy

FiPy is a Finite Volume PDE solver written in Python
http://pages.nist.gov/fipy/en/latest
Other
489 stars 148 forks source link

Is there a way to make mesh generation faster? #1047

Open CalebDmArcher opened 1 month ago

CalebDmArcher commented 1 month ago

Currently, I am doing simulation with a shape that has a large top and bottom surface area but it is thin. The top and bottom surfaces are 2cm * 5cm, but the thicknesses of the layers are 5e-4m and 5e-5m.

I have tried to make the mesh size 5e-4 which has quite a fast mesh generation with Gmsh but the result is very rough (compare to the theoretical answer). I've tried to decrease the mesh size a little bit and the result does get closer to the theoretical value but not enough and I cannot decrease the mesh size further due to the heavy load of mesh generation.

For this special type of shape, do you have any suggestions for me?

Below is the code for the mesh part:

from fipy import Gmsh3D
import matplotlib
matplotlib.use('TkAgg')

t_cu = 0.02
t_fr4 = 0.02

# Gmsh script as a string
gmsh_script = """
SetFactory("OpenCASCADE");

// Define layer dimensions
Lx = 0.05;
Ly = 0.02;
t_cu = 5.e-5;
t_fr4 = 5.e-4;

// Mesh sizes for each layer
dx_cu = 5.e-4;
dx_fr4 = 5.e-4;

// Create Layer 1
Box(1) = {0, 0, 0, Lx, Ly, t_cu};

// Extrude Layer 2 from Layer 1
Extrude {0, 0, t_fr4} {
  Surface{6};
}

// Extrude Layer 3 from Layer 2
Extrude {0, 0, t_fr4} {
  Surface{11};
}

Coherence;

// Meshing properties for each layer
Field[1] = Box;
Field[1].VIn = dx_cu;
Field[1].VOut = dx_cu;
Field[1].XMin = 0; Field[1].XMax = Lx;
Field[1].YMin = 0; Field[1].YMax = Ly;
Field[1].ZMin = 0; Field[1].ZMax = t_cu;

Field[2] = Box;
Field[2].VIn = dx_fr4;
Field[2].VOut = dx_fr4;
Field[2].XMin = 0; Field[2].XMax = Lx;
Field[2].YMin = 0; Field[2].YMax = Ly;
Field[2].ZMin = t_cu; Field[2].ZMax = t_cu + t_fr4;

Field[3] = Box;
Field[3].VIn = dx_cu;
Field[3].VOut = dx_cu;
Field[3].XMin = 0; Field[1].XMax = Lx;
Field[3].YMin = 0; Field[1].YMax = Ly;
Field[3].ZMin = t_cu + t_fr4; Field[1].ZMax = t_cu + t_fr4 + t_cu;

Field[4] = Min;
Field[4].FieldsList = {1,2,3};

Background Field = 4;

// Generate 3D mesh
Mesh 4;
"""

# Use Gmsh2D to create a mesh from the Gmsh script
mesh = Gmsh3D(gmsh_script)

# Set spatial variables for calculation and plot
xc, yc, zc = mesh.cellCenters
x, y, z = mesh.faceCenters

# top_faces = mesh.facesTop & (z > (t_cu + t_fr4)) # Identify the top faces
top_faces = mesh.facesFront
# top_faces = mesh.facesTop
cell_areas = mesh._faceAreas[top_faces]  # Calculate the areas of the top faces
A = sum(cell_areas)  # Calculate the uniform heat flux q/A
print(f"Top surface area is {A} m^2")

# print(top_faces)

exterior_faces = mesh.exteriorFaces
exterior_areas = mesh._faceAreas[exterior_faces]
total_exterior_area = sum(exterior_areas)
print(f"Total exterior surface area is {total_exterior_area} m^2")

bottom_faces = mesh.facesBack
# bottom_faces = mesh.facesBottom # Identify the bottom faces
cell_areas_bottom = mesh._faceAreas[bottom_faces]  # Calculate the areas of the top faces
Ab = sum(cell_areas_bottom)  # Calculate the uniform heat flux q/A
print(f"Bottom surface area is {Ab} m^2")

bottom_faces = mesh.facesBottom

# import tempfile
# # Write the Gmsh script to a temporary file
# with tempfile.NamedTemporaryFile(suffix='.geo', delete=False) as temp_geo_file:
#     temp_geo_file.write(gmsh_script.encode('utf-8'))
#     temp_geo_file_path = temp_geo_file.name
#
# import gmsh
# # draw the mesh
# # Initialize the Gmsh Python API
# gmsh.initialize()
# # gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 2)
# # Merge the .geo file
# # gmsh.merge('cube.geo')
# gmsh.merge(temp_geo_file_path)
# # Generate the 3D mesh
# gmsh.model.mesh.generate(3)
# # Launch the Gmsh GUI to visualize the mesh (optional)
# gmsh.fltk.run()
# # Finalize the Gmsh Python API
# gmsh.finalize()
CalebDmArcher commented 1 month ago

A Gmsh developer suggested using gmsh your_script.geo -nt 8 -algo hxt - to mesh with a new 3D algo with 8 threads. So the code below has been tried which has some problems. The biggest problem is that Gmsh3D cannot use the '.msh' file generated from Gmsh. Seems like Gmsh3D has its own way of using the '.geo' file to generate something that Fipy could use. So I guess I need to find a way to transfer the '.msh' file to something that is just like what Gmsh3D generates so Fipy can work on it.

Do you have any suggestions? Or are there any mistakes?

import os
from fipy import Gmsh3D
import matplotlib
import subprocess
import tempfile
import gmsh

matplotlib.use('TkAgg')

# Gmsh script as a string
gmsh_script = """
SetFactory("OpenCASCADE");

// Define layer dimensions
Lx = 0.05;
Ly = 0.02;
t_cu = 5.e-5;
t_fr4 = 5.e-4;

// Mesh sizes for each layer
dx_cu = 5.e-4;
dx_fr4 = 5.e-4;

// Create Layer 1
Box(1) = {0, 0, 0, Lx, Ly, t_cu};

// Extrude Layer 2 from Layer 1
Extrude {0, 0, t_fr4} {
  Surface{6};
}

// Extrude Layer 3 from Layer 2
Extrude {0, 0, t_cu} {
  Surface{11};
}

Coherence;

// Meshing properties for each layer
Field[1] = Box;
Field[1].VIn = dx_cu;
Field[1].VOut = dx_cu;
Field[1].XMin = 0; Field[1].XMax = Lx;
Field[1].YMin = 0; Field[1].YMax = Ly;
Field[1].ZMin = 0; Field[1].ZMax = t_cu;

Field[2] = Box;
Field[2].VIn = dx_fr4;
Field[2].VOut = dx_fr4;
Field[2].XMin = 0; Field[2].XMax = Lx;
Field[2].YMin = 0; Field[2].YMax = Ly;
Field[2].ZMin = t_cu; Field[2].ZMax = t_cu + t_fr4;

Field[3] = Box;
Field[3].VIn = dx_cu;
Field[3].VOut = dx_cu;
Field[3].XMin = 0; Field[1].XMax = Lx;
Field[3].YMin = 0; Field[1].YMax = Ly;
Field[3].ZMin = t_cu + t_fr4; Field[1].ZMax = t_cu + t_fr4 + t_cu;

Field[4] = Min;
Field[4].FieldsList = {1,2,3};

Background Field = 4;

// Generate 3D mesh
Mesh 4;
"""

# Write the Gmsh script to a temporary file
with tempfile.NamedTemporaryFile(suffix='.geo', delete=False) as temp_geo_file:
    temp_geo_file.write(gmsh_script.encode('utf-8'))
    temp_geo_file_path = temp_geo_file.name

# Define the path for the .msh file
msh_file_path = temp_geo_file_path.replace('.geo', '.msh')

# Run Gmsh as a subprocess with the specified options
subprocess.run(['gmsh', temp_geo_file_path, '-3', '-nt', '8', '-algo', 'hxt', '-o', msh_file_path])

# Check if the .msh file was created
if not os.path.exists(msh_file_path):
    raise FileNotFoundError(f"Mesh file {msh_file_path} not found. Gmsh might have failed to generate it.")

# Initialize the Gmsh API
gmsh.initialize()
gmsh.open(msh_file_path)

# Print some basic information about the mesh
model = gmsh.model
print(f"Number of nodes: {model.mesh.getNodes()[0].size}")
print(f"Number of elements: {model.mesh.getElements()[1][0].size}")

# Optionally, visualize the mesh
gmsh.fltk.run()

# Finalize the Gmsh API
gmsh.finalize()

# Load the mesh generated by Gmsh
mesh = Gmsh3D(msh_file_path)

# Set spatial variables for calculation and plot
xc, yc, zc = mesh.cellCenters
x, y, z = mesh.faceCenters

top_faces = mesh.facesFront
cell_areas = mesh._faceAreas[top_faces]  # Calculate the areas of the top faces
A = sum(cell_areas)  # Calculate the uniform heat flux q/A
print(f"Top surface area is {A} m^2")

exterior_faces = mesh.exteriorFaces
exterior_areas = mesh._faceAreas[exterior_faces]
total_exterior_area = sum(exterior_areas)
print(f"Total exterior surface area is {total_exterior_area} m^2")

bottom_faces = mesh.facesBack
cell_areas_bottom = mesh._faceAreas[bottom_faces]  # Calculate the areas of the top faces
Ab = sum(cell_areas_bottom)  # Calculate the uniform heat flux q/A
print(f"Bottom surface area is {Ab} m^2")
CalebDmArcher commented 1 month ago

May I close this thread?

guyer commented 1 month ago

If your question has been addressed, sure. I have been away, so haven't had a chance to look at your question, yet.