Hi everyone,
I’ll get straight to the point with my problem and example code, then add a bit more information at the bottom about why I’m doing this for those interested in the motivation behind it.
I’m solving a moving boundary problem iteratively that I’m re-meshing for on each iteration using mshr’s Polygon mesh function and the new position of my moving boundary. The re-meshing every iteration is proving to be very time consuming, especially due to the high resolution I need on my moving boundary. Is there a way to optimise this?
Ideally I’d like my code to run in parallel, however it takes even longer re-meshing in parallel. I’ve also tried using gmsh, but that is also slower (though I may not be doing it the best way). Any ideas?
My example code:
from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np
import time
rank = MPI.comm_world.rank
mesh_file = File('MeshFolder/moving_mesh.pvd')
#Set inital position of moving boundary with high resolution
points_y = np.linspace(0.1, -0.1, num=int(5001), endpoint=True).reshape(-2,1)
points_x = np.zeros(len(points_y)).reshape(-2,1)
points = np.concatenate((points_x,points_y),axis=1)
#First and last points have alpha=0 so they hold the corner positions
alpha = np.ones(len(points))
alpha[[0,len(points)-1]] = 0
#Create initial mesh
domain_vertices = [Point(0,-5),Point(5,-5),Point(5,5),Point(0,5)]
for idx in range(len(points)):
domain_vertices.append(Point(points[idx]))
polygon_domain = Polygon(domain_vertices)
mesh = generate_mesh(polygon_domain,100)
mesh_file << (mesh,0)
#Iterate
for count in range(5):
t0 = time.time()
#Calculate new boundary position
points[:,0] -= 0.1*alpha
#Create new mesh
domain_vertices = [Point(0,-5),Point(5,-5),Point(5,5),Point(0,5)]
for idx in range(len(points)):
domain_vertices.append(Point(points[idx]))
polygon_domain = Polygon(domain_vertices)
mesh = generate_mesh(polygon_domain,100)
mesh_file << (mesh,count+1)
#Display time taken
if rank == 0:
print(count, time.time()-t0)
Motivation:
I’m developing a front-tracking numerical method for simulating a fluid flowing through a pore network. I start with a random pore network (Fig1 - red is solid, green is empty space) with a liquid above it. I have some equations I solve in the liquid which affect the movement of the liquid. I start by meshing the liquid phase then iteratively solve the equations in the liquid, then move the parts of the boundary touching the green regions. As this results in my mesh distorting and I need the moving boundary to be able to split and move in different directions and re-join as it moves through the pore network, I then create a new mesh and interpolate the solution onto it. In order to accurately model the movement through the pore network I need a very high resolution in this region (Fig2).
