Skip to content

Commit

Permalink
Merge branch 'development' of https://github.com/underworldcode/under…
Browse files Browse the repository at this point in the history 
…world3 into development
  • Loading branch information
@bknight1
bknight1 committed Nov 9, 2023
2 parents 004f236 + 592298c commit da1d18d
Show file tree
Hide file tree
Showing 3 changed files with 72 additions and 212 deletions.
235 changes: 58 additions & 177 deletions Jupyterbook/Notebooks/Examples-NavierStokes/Ex_Navier_Stokes_Benchmarks_NS_DFG_2d_SLCN.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,9 @@
import underworld3 as uw
from underworld3 import timing

import nest_asyncio
nest_asyncio.apply()

import numpy as np
import sympy

Expand All @@ -53,8 +56,8 @@
# These can be set when launching the script as
# mpirun python3 scriptname -uw_resolution=0.1 etc

resolution = uw.options.getReal("model_resolution", default=10)
refinement = uw.options.getInt("model_refinement", default=2)
resolution = uw.options.getReal("model_resolution", default=15)
refinement = uw.options.getInt("model_refinement", default=0)
model = uw.options.getInt("model_number", default=3)
maxsteps = uw.options.getInt("max_steps", default=1000)
restart_step = uw.options.getInt("restart_step", default=-1)
Expand Down Expand Up @@ -197,50 +200,6 @@ def pipemesh_return_coords_to_bounds(coords):
Vb = (4.0 * U0 * y * (0.41 - y)) / 0.41**2


# +
# check the mesh if in a notebook / serial

if uw.mpi.size == 1:
import numpy as np
import pyvista as pv
import vtk

pv.global_theme.background = "white"
pv.global_theme.window_size = [1250, 1250]
pv.global_theme.anti_aliasing = "msaa"
pv.global_theme.jupyter_backend = "panel"
pv.global_theme.smooth_shading = True
pv.global_theme.camera["viewup"] = [0.0, 1.0, 0.0]
pv.global_theme.camera["position"] = [0.0, 0.0, 1.0]

pipemesh.vtk("tmpmesh.vtk")
pvmesh = pv.read(f"tmpmesh.vtk")

# point sources at cell centres

points = np.zeros((pipemesh._centroids.shape[0], 3))
points[:, 0] = pipemesh._centroids[:, 0]
points[:, 1] = pipemesh._centroids[:, 1]
point_cloud = pv.PolyData(points)

pl = pv.Plotter()

pl.add_mesh(
pvmesh,
cmap="coolwarm",
edge_color="Black",
show_edges=True,
use_transparency=False,
opacity=1.0,
)

# pl.add_mesh(pvmesh,'Black', 'wireframe', opacity=0.75)
# pl.add_mesh(pvstream)

# pl.remove_scalar_bar("mag")

pl.show()

# +
v_soln = uw.discretisation.MeshVariable("U", pipemesh, pipemesh.dim, degree=2)
vs_soln = uw.discretisation.MeshVariable("Us", pipemesh, pipemesh.dim, degree=2)
Expand Down Expand Up @@ -299,7 +258,7 @@ def pipemesh_return_coords_to_bounds(coords):
order=2,
)

navier_stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel(v_soln)
navier_stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel

# Constant visc

Expand Down Expand Up @@ -384,43 +343,20 @@ def pipemesh_return_coords_to_bounds(coords):
# check the mesh if in a notebook / serial

if uw.mpi.size == 1:
import numpy as np
import pyvista as pv
import vtk
import underworld3 as uw
import underworld3.visualisation

pv.global_theme.background = "white"
pv.global_theme.window_size = [1250, 1250]
pv.global_theme.anti_aliasing = "msaa"
pv.global_theme.jupyter_backend = "panel"
pv.global_theme.smooth_shading = True
pv.global_theme.camera["viewup"] = [0.0, 1.0, 0.0]
pv.global_theme.camera["position"] = [0.0, 0.0, 1.0]
pvmesh = uw.visualisation.mesh_to_pv_mesh(pipemesh)
pvmesh.point_data["V"] = uw.visualisation.vector_fn_to_pv_points(pvmesh, v_soln.sym)
pvmesh.point_data["Vmag"] = uw.visualisation.scalar_fn_to_pv_points(pvmesh, v_soln.sym.dot(v_soln.sym))

pipemesh.vtk("tmpmesh.vtk")
pvmesh = pv.read(f"tmpmesh.vtk")
velocity_points = underworld3.visualisation.meshVariable_to_pv_cloud(v_soln)
velocity_points.point_data["V"] = uw.visualisation.vector_fn_to_pv_points(velocity_points, v_soln.sym)

with pipemesh.access():
# usol = navier_stokes._u_star_projector.u.data.copy()
usol = v_soln.data.copy()
pl = pv.Plotter(window_size=(1000, 750))

with pipemesh.access():
pvmesh.point_data["Vmag"] = uw.function.evalf(
sympy.sqrt(v_soln.sym.dot(v_soln.sym)), pipemesh.data
)
pvmesh.point_data["P"] = uw.function.evalf(p_soln.fn, pipemesh.data)

v_vectors = np.zeros((pipemesh.data.shape[0], 3))
v_vectors[:, 0] = uw.function.evalf(v_soln[0].sym, pipemesh.data)
v_vectors[:, 1] = uw.function.evalf(v_soln[1].sym, pipemesh.data)
pvmesh.point_data["V"] = v_vectors

arrow_loc = np.zeros((v_soln.coords.shape[0], 3))
arrow_loc[:, 0:2] = v_soln.coords[...]

arrow_length = np.zeros((v_soln.coords.shape[0], 3))
arrow_length[:, 0:2] = usol[...]

# point sources at cell centres
# point sources at cell centres for streamlines

points = np.zeros((pipemesh._centroids.shape[0], 3))
points[:, 0] = pipemesh._centroids[:, 0]
Expand All @@ -431,14 +367,7 @@ def pipemesh_return_coords_to_bounds(coords):
point_cloud, vectors="V", integration_direction="forward", max_steps=10
)

pl = pv.Plotter()

pl.add_arrows(arrow_loc, arrow_length, mag=0.025 / U0, opacity=0.75)

# pl.add_points(spoint_cloud,color="Black",
# render_points_as_spheres=False,
# point_size=5, opacity=0.66
# )
pl.add_arrows(velocity_points.points, velocity_points.point_data["V"], mag=0.025 / U0, opacity=0.75)

pl.add_mesh(
pvmesh,
Expand All @@ -456,98 +385,42 @@ def pipemesh_return_coords_to_bounds(coords):
# pl.remove_scalar_bar("mag")

pl.show()
# -


# +
if uw.mpi.size == 1:
pv.global_theme.background = "white"
pv.global_theme.window_size = [1250, 1000]
pv.global_theme.anti_aliasing = "msaa"
pv.global_theme.jupyter_backend = "panel"
pv.global_theme.smooth_shading = True
# pv.global_theme.camera['viewup'] = [0.0, 1.0, 0.0]
# pv.global_theme.camera['position'] = [0.0, 0.0, 2.0]

pl = pv.Plotter()


def plot_V_mesh(filename):
import underworld3 as uw


if uw.mpi.size == 1:
import numpy as np

import pyvista as pv
import vtk

## Plotting into existing pl (memory leak in pyvista)
pl.clear()

pvmesh = pv.read(f".meshes/ns_pipe_flow_{resolution}.msh")

with passive_swarm.access():
points = np.zeros((passive_swarm.data.shape[0], 3))
points[:, 0] = passive_swarm.data[:, 0]
points[:, 1] = passive_swarm.data[:, 1]

point_cloud = pv.PolyData(points)

import underworld3.visualisation
pvmesh = uw.visualisation.mesh_to_pv_mesh(pipemesh)
pvmesh.point_data["V"] = uw.visualisation.vector_fn_to_pv_points(pvmesh, v_soln.sym)
pvmesh.point_data["Omega"] = uw.visualisation.scalar_fn_to_pv_points(pvmesh, vorticity.sym)
pvmesh.point_data["Vmag"] = uw.visualisation.scalar_fn_to_pv_points(pvmesh, v_soln.sym.dot(v_soln.sym))
velocity_points = underworld3.visualisation.meshVariable_to_pv_cloud(v_soln)
velocity_points.point_data["V"] = uw.visualisation.vector_fn_to_pv_points(velocity_points, v_soln.sym)

pl = pv.Plotter(window_size=(1000, 750))
# point sources at cell centres for streamlines
points = np.zeros((pipemesh._centroids.shape[0], 3))
points[:, 0] = pipemesh._centroids[:, 0]
points[:, 1] = pipemesh._centroids[:, 1]

c_point_cloud = pv.PolyData(points)

# with swarm.access():
# spoints = np.zeros((swarm.particle_coordinates.data.shape[0], 3))
# spoints[:, 0] = swarm.particle_coordinates.data[:, 0]
# spoints[:, 1] = swarm.particle_coordinates.data[:, 1]
# spoint_cloud = pv.PolyData(spoints)

with pipemesh.access():
pvmesh.point_data["P"] = uw.function.evalf(p_soln.sym[0], pipemesh.data)
pvmesh.point_data["Omega"] = uw.function.evalf(
vorticity.sym[0], pipemesh.data
)

with pipemesh.access():
usol = v_soln.data.copy()

v_vectors = np.zeros((pipemesh.data.shape[0], 3))
v_vectors[:, 0] = uw.function.evalf(v_soln[0].sym, pipemesh.data)
v_vectors[:, 1] = uw.function.evalf(v_soln[1].sym, pipemesh.data)
pvmesh.point_data["V"] = v_vectors

arrow_loc = np.zeros((v_soln.coords.shape[0], 3))
arrow_loc[:, 0:2] = v_soln.coords[...]

arrow_length = np.zeros((v_soln.coords.shape[0], 3))
arrow_length[:, 0:2] = usol[...]

pl.add_arrows(arrow_loc, arrow_length, mag=0.033 / U0, opacity=0.5)

point_cloud = pv.PolyData(points)

pvstream = pvmesh.streamlines_from_source(
c_point_cloud, vectors="V", integration_direction="both", max_time=0.25
)

# pl.add_mesh(pvmesh,'Black', 'wireframe', opacity=0.75)

# pl.add_points(
# spoint_cloud,
# color="Grey",
# render_points_as_spheres=True,
# point_size=3,
# opacity=0.5,
# )

pl.add_mesh(
pvmesh,
cmap="coolwarm",
edge_color="Black",
show_edges=False,
scalars="Omega",
use_transparency=False,
opacity=0.5,
point_cloud, vectors="V", integration_direction="forward", max_steps=10
)


pl.add_arrows(velocity_points.points, velocity_points.point_data["V"], mag=0.025 / U0, opacity=0.75)
pl.add_mesh(pvstream)

pl.add_points(
Expand All @@ -562,7 +435,7 @@ def plot_V_mesh(filename):
pl.camera.SetFocalPoint(0.75, 0.2, 0.0)
pl.camera.SetClippingRange(1.0, 8.0)

pl.remove_scalar_bar("Omega")
# pl.remove_scalar_bar("Omega")
pl.remove_scalar_bar("mag")
pl.remove_scalar_bar("V")

Expand All @@ -572,8 +445,6 @@ def plot_V_mesh(filename):
window_size=(2560, 1280),
return_img=False,
)


# -

ts = 0
Expand All @@ -585,10 +456,6 @@ def plot_V_mesh(filename):
delta_t_cfl = navier_stokes.estimate_dt()
delta_t = dt_ns # min(delta_t_swarm, dt_ns)

print(
f"Memory usage [1] = {python_process.memory_info().rss//1000000} Mb", flush=True
)

navier_stokes.solve(timestep=dt_ns, zero_init_guess=False)

# update passive swarm
Expand Down Expand Up @@ -633,14 +500,26 @@ def plot_V_mesh(filename):
# check the mesh if in a notebook / serial

if uw.mpi.size == 1:
import numpy as np

import pyvista as pv
import vtk
import underworld3 as uw
import underworld3.visualisation

pvmesh = uw.visualisation.mesh_to_pv_mesh(pipemesh)
pvmesh.point_data["V"] = uw.visualisation.vector_fn_to_pv_points(pvmesh, v_soln.sym)
pvmesh.point_data["Vmag"] = uw.visualisation.scalar_fn_to_pv_points(pvmesh, v_soln.sym.dot(v_soln.sym))
pvmesh.point_data["Omega"] = uw.visualisation.scalar_fn_to_pv_points(pvmesh, vorticity.sym)

velocity_points = underworld3.visualisation.meshVariable_to_pv_cloud(v_soln)
velocity_points.point_data["V"] = uw.visualisation.vector_fn_to_pv_points(velocity_points, v_soln.sym)

pl = pv.Plotter(window_size=(1000, 750))


pv.global_theme.background = "white"
pv.global_theme.window_size = [1250, 1250]
pv.global_theme.anti_aliasing = "msaa"
pv.global_theme.jupyter_backend = "panel"
pv.global_theme.jupyter_backend = "client"
pv.global_theme.smooth_shading = True
pv.global_theme.camera["viewup"] = [0.0, 1.0, 0.0]
pv.global_theme.camera["position"] = [0.0, 0.0, 1.0]
Expand All @@ -650,7 +529,7 @@ def plot_V_mesh(filename):
with pipemesh.access():
usol = v_soln.data.copy()

ustar = navier_stokes.DuDt.psi_star[0].sym
ustar = navier_stokes.Unknowns.DuDt.psi_star[0].sym

with pipemesh.access():
pvmesh.point_data["Vmag"] = uw.function.evalf(
Expand Down Expand Up @@ -723,3 +602,5 @@ def plot_V_mesh(filename):
print(navier_stokes.DFDt.psi_star[0].data.max())

navier_stokes._u_f1


4 changes: 2 additions & 2 deletions underworld3/cython/petsc_generic_snes_solvers.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1542,7 +1542,7 @@ class SNES_Stokes_SaddlePt(Solver):
self.is_setup = False

def _setup_history_terms(self):
self.Unknowns.DuDt = uw.swarm.SemiLagrange_Updater(
self.Unknowns.DuDt = uw.swarm.SemiLagrange_D_Dt(
self.mesh,
self.u.sym,
self.u.sym,
Expand All @@ -1556,7 +1556,7 @@ class SNES_Stokes_SaddlePt(Solver):
smoothing=0.0,
)

self.Unknowns.DFDt = uw.swarm.SemiLagrange_Updater(
self.Unknowns.DFDt = uw.swarm.SemiLagrange_D_Dt(
self.mesh,
sympy.Matrix.zeros(self.mesh.dim, self.mesh.dim),
self.u.sym,
Expand Down
Loading

0 comments on commit da1d18d

Please sign in to comment.