From 8781b98732a88239ef4ea68cc0104d305219cedb Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 11 Apr 2024 00:26:01 +0000 Subject: [PATCH] deploy: 0d2ed110d63eeff96a5d1c79bb7c6e8054ec52ab --- ...x_AdvectionDiffusionSLCN_RotationTest.html | 2 - .../Ex_NavierStokesRotationTest.html | 155 +++++++++++------- ...vier_Stokes_Benchmarks_NS_DFG_2d_SLCN.html | 16 +- .../Ex_Navier_Stokes_Visualise_Disc.html | 97 ++++++----- .../Ex_Sheared_Layer_Elastic.html | 7 + .../Ex_Sheared_Layer_Test.html | 46 +++++- .../Ex_Stokes_Cartesian_SolC.html | 9 +- .../Ex_AdvectionDiffusionSLCN_RotationTest.py | 7 +- .../Ex_Convection_Disc_InternalHeat.py | 2 + .../Ex_NavierStokesRotationTest.py | 135 +++++++++------ ...Navier_Stokes_Benchmarks_NS_DFG_2d_SLCN.py | 8 +- .../Ex_Navier_Stokes_Visualise_Disc.py | 98 +++++------ .../Ex_Sheared_Layer_Elastic.py | 2 + .../Ex_Sheared_Layer_Test.py | 37 ++++- .../Ex_Stokes_Cartesian_SolC.py | 4 +- main/objects.inv | Bin 3428 -> 3435 bytes main/searchindex.js | 2 +- main/slideshows/index.html | 2 +- 18 files changed, 400 insertions(+), 229 deletions(-) diff --git a/main/Notebooks/Examples-Convection/Ex_AdvectionDiffusionSLCN_RotationTest.html b/main/Notebooks/Examples-Convection/Ex_AdvectionDiffusionSLCN_RotationTest.html index 7763046..4eedd56 100644 --- a/main/Notebooks/Examples-Convection/Ex_AdvectionDiffusionSLCN_RotationTest.html +++ b/main/Notebooks/Examples-Convection/Ex_AdvectionDiffusionSLCN_RotationTest.html @@ -454,8 +454,6 @@

Field (SemiLagrange) Advection solver test
-
import petsc4py
-from petsc4py import PETSc
-
-import underworld3 as uw
+
import underworld3 as uw
 from underworld3.systems import Stokes
 from underworld3.systems import NavierStokesSLCN
-
 from underworld3 import function
 
+
 import numpy as np
 import sympy
 

@@ -496,11 +493,17 @@ 
Navier Stokes test: boundary driven ring with step change in boundary condit
 # These can be set when launching the script as
 # mpirun python3 scriptname -uw_resolution=0.1 etc
 
-resolution = uw.options.getInt("model_resolution", default=25)
+resolution = uw.options.getInt("model_resolution", default=10)
 refinement = uw.options.getInt("model_refinement", default=0)
-maxsteps = uw.options.getInt("max_steps", default=1000)
+maxsteps = uw.options.getInt("max_steps", default=25)
 restart_step = uw.options.getInt("restart_step", default=-1)
 rho = uw.options.getReal("rho", default=1000)
+
+outdir="output"
+
+if uw.mpi.rank == 0:
+    print(f"restart: {restart_step}")
+    print(f"resolution: {resolution}")
@@ -508,7 +511,7 @@

Navier Stokes test: boundary driven ring with step change in boundary condit 
meshball = uw.meshing.Annulus(
-    radiusOuter=1.0, radiusInner=0.5, cellSize=1/resolution, qdegree=3
+    radiusOuter=1.0, radiusInner=0.0, cellSize=1/resolution, qdegree=3
 )
@@ -516,6 +519,13 @@

Navier Stokes test: boundary driven ring with step change in boundary condit

+
meshball.view()
+
+
+
+
+
+
# Define some functions on the mesh
 
 import sympy
@@ -536,8 +546,8 @@ 
Navier Stokes test: boundary driven ring with step change in boundary condit
 # Rigid body rotation v_theta = constant, v_r = 0.0
 
 theta_dot = 2.0 * np.pi  # i.e one revolution in time 1.0
-v_x = -1.0 * r * theta_dot * sympy.sin(th)  # * y # to make a convergent / divergent bc
-v_y = r * theta_dot * sympy.cos(th)  # * y
+v_x = -1.0 * r * theta_dot * sympy.sin(th) * y # to make a convergent / divergent bc
+v_y = r * theta_dot * sympy.cos(th) * y
@@ -547,7 +557,7 @@

Navier Stokes test: boundary driven ring with step change in boundary condit 
v_soln = uw.discretisation.MeshVariable("U", meshball, meshball.dim, degree=2)
 p_soln = uw.discretisation.MeshVariable("P", meshball, 1, degree=1)
 vorticity = uw.discretisation.MeshVariable(
-    "\omega", meshball, 1, degree=1, continuous=False
+    "\omega", meshball, 1, degree=1, continuous=True
 )
@@ -555,7 +565,7 @@

Navier Stokes test: boundary driven ring with step change in boundary condit

-
navier_stokes = NavierStokesSLCN(
+
navier_stokes = uw.systems.NavierStokes(
     meshball,
     velocityField=v_soln,
     pressureField=p_soln,
@@ -609,16 +619,6 @@ 
Navier Stokes test: boundary driven ring with step change in boundary condit
 passive_swarm.populate(
     fill_param=3,
 )
-
-# add new points at the 12 o'clock position
-
-npoints = 100
-passive_swarm.dm.addNPoints(npoints)
-with passive_swarm.access(passive_swarm.particle_coordinates):
-    for i in range(npoints):
-        passive_swarm.particle_coordinates.data[-1 : -(npoints + 1) : -1, :] = np.array(
-            [-0.05, 0.9] + 0.1 * np.random.random((npoints, 2))
-        )
@@ -636,8 +636,32 @@

Navier Stokes test: boundary driven ring with step change in boundary condit navier_stokes.bodyforce = sympy.Matrix([0, 0]) # Velocity boundary conditions -navier_stokes.add_dirichlet_bc((v_x, v_y), "Upper", (0, 1)) -navier_stokes.add_dirichlet_bc((0.0, 0.0), "Lower", (0, 1)) +# navier_stokes.add_dirichlet_bc((v_x, v_y), "Upper") +# navier_stokes.add_dirichlet_bc((0.0, 0.0), "Lower") + +# Try this one: + + +# upper_mask = meshball.meshVariable_mask_from_label("UW_Boundaries", meshball.boundaries.Upper.value ) +# lower_mask = meshball.meshVariable_mask_from_label("UW_Boundaries", meshball.boundaries.Lower.value ) + +# vbc_xn = 1000 * ((v_soln.sym[0] - v_x) * upper_mask.sym[0] + v_soln.sym[0] * lower_mask.sym[0]) +# vbc_yn = 1000 * ((v_soln.sym[1] - v_y) * upper_mask.sym[0] + v_soln.sym[1] * lower_mask.sym[0]) + +# vbc_x = v_x * upper_mask.sym[0] + 0 * lower_mask.sym[0] +# vbc_y = v_y * upper_mask.sym[0] + 0 * lower_mask.sym[0] + +# navier_stokes.add_natural_bc( (vbc_xn, vbc_yn), "All_Boundaries") + +# navier_stokes.add_natural_bc( (1000 * (v_soln.sym[0] - v_x), 1000* ( v_soln.sym[1] - v_y)), "Upper") +# navier_stokes.add_natural_bc( (1000 * v_x, 1000*v_y), "Lower") + + +# navier_stokes.add_natural_bc( (vbc_xn, vbc_yn), "Upper") +# navier_stokes.add_natural_bc( (vbc_xn, vbc_yn), "Lower") + +navier_stokes.add_dirichlet_bc((v_x, v_y), "Upper") +# navier_stokes.add_dirichlet_bc((0.0, 0.0), "Lower") expt_name = f"Cylinder_NS_rho_{navier_stokes.rho}_{resolution}"

@@ -646,22 +670,44 @@

Navier Stokes test: boundary driven ring with step change in boundary condit

-
navier_stokes.solve(timestep=0.1)
-navier_stokes.estimate_dt()
+
navier_stokes.delta_t_physical = 0.1
+navier_stokes.solve(timestep=0.1, verbose=False, evalf=True, order=1)
+# navier_stokes.rho = rho
-
nodal_vorticity_from_v.solve()
+
navier_stokes.estimate_dt()
+

+

+
+
+
+
+
if restart_step > 0:
+    if uw.mpi.rank == 0:
+        print(f"Reading step {restart_step}", flush=True)
+
+    passive_swarm = uw.swarm.Swarm(mesh=meshball)
+    passive_swarm.read_timestep(
+        expt_name, "passive_swarm", restart_step, outputPath=outdir
+    )
+    
+    v_soln.read_timestep(expt_name, "U", restart_step, outputPath=outdir)
+    p_soln.read_timestep(expt_name, "P", restart_step, outputPath=outdir)
+    
+    # Flux history variable might be a good idea
-
# check the mesh if in a notebook / serial
+
nodal_vorticity_from_v.solve()
+
+# check the mesh if in a notebook / serial
 if uw.mpi.size == 1:
     import pyvista as pv
     import underworld3.visualisation as vis
@@ -699,7 +745,7 @@ 
Navier Stokes test: boundary driven ring with step change in boundary condit
     pl.add_arrows(
         velocity_points.points,
         velocity_points.point_data["V"],
-        mag=1.0e-2,
+        mag=2.0e-2,
         opacity=0.75,
     )
 
@@ -720,13 +766,9 @@ 
Navier Stokes test: boundary driven ring with step change in boundary condit
     pl.remove_scalar_bar("V")
 
     pl.show(jupyter_backend="client")
-

-

-
-
-
-
-
def plot_V_mesh(filename):
+
+
+def plot_V_mesh(filename):
     if uw.mpi.size == 1:
         import pyvista as pv
         import underworld3.visualisation as vis
@@ -812,7 +854,10 @@ 
Navier Stokes test: boundary driven ring with step change in boundary condit
-
ts = 0
+
if restart_step > 0:
+    ts = restart_step
+else:
+    ts = 0
@@ -820,35 +865,31 @@

Navier Stokes test: boundary driven ring with step change in boundary condit 
# Time evolution model / update in time
+navier_stokes.delta_t_physical = 0.1
+
+delta_t = 0.1 #  5.0 * navier_stokes.estimate_dt()
 
-for step in range(0, 100):  # 250
-    delta_t = 2.0 * navier_stokes.estimate_dt()
-    navier_stokes.solve(timestep=delta_t, zero_init_guess=False)
+for step in range(0, maxsteps+1):  # 250
+    
+    # if step%10 == 0:
+    #     delta_t = 5.0 * navier_stokes.estimate_dt()
 
+    navier_stokes.solve(timestep=delta_t, zero_init_guess=False, evalf=True)    
     passive_swarm.advection(v_soln.sym, delta_t, order=2, corrector=False, evalf=False)
 
     nodal_vorticity_from_v.solve()
 
-    npoints = 100
-    passive_swarm.dm.addNPoints(npoints)
-    with passive_swarm.access(passive_swarm.particle_coordinates):
-        for i in range(npoints):
-            passive_swarm.particle_coordinates.data[-1 : -(npoints + 1) : -1, :] = np.array(
-                [-0.05, 0.9] + 0.1 * np.random.random((npoints, 2))
-            )
-
-
     if uw.mpi.rank == 0:
-        print("Timestep {}, dt {}".format(ts, delta_t))
+        print("Timestep {}, dt {}".format(ts, delta_t), flush=True)
 
     if ts % 5 == 0:
-        plot_V_mesh(filename="output/{}_step_{}".format(expt_name, ts))
-
+        plot_V_mesh(filename=f"{outdir}/{expt_name}_step_{ts}")
+        
         meshball.write_timestep(
             expt_name,
             meshUpdates=True,
-            meshVars=[p_soln, v_soln, vorticity, St],
-            outputPath="output",
+            meshVars=[p_soln, v_soln, vorticity],
+            outputPath=outdir,
             index=ts,
         )
 
@@ -856,7 +897,7 @@ 
Navier Stokes test: boundary driven ring with step change in boundary condit
             expt_name,
             "passive_swarm",
             swarmVars=None,
-            outputPath="output",
+            outputPath=outdir,
             index=ts,
             force_sequential=True,
         )
@@ -866,14 +907,14 @@ 
Navier Stokes test: boundary driven ring with step change in boundary condit
+

! open .

-
pvmesh.point_data["V"].min()
+
navier_stokes._u_f0
-

! open .