Playing with the Diablo mesh

Author: lmoresi

Notebooks/Examples-Meshing/Batmesh_number_one.h5

@@ -39,12 +39,14 @@ class bd(Enum):
 batmesh.dm.view()
 
 # %%
-Rayleigh_number = uw.function.expression(r"\textrm{Ra}", sympy.sympify(10)**6, "Rayleigh number")
+Rayleigh_number = uw.function.expression(r"\textrm{Ra}", sympy.sympify(10000000), "Rayleigh number")
 
 # %%
-v_soln = uw.discretisation.MeshVariable("U", batmesh, batmesh.dim, degree=2, continuous=True)
-p_soln = uw.discretisation.MeshVariable("P", batmesh, 1, degree=1, continuous=True)
-t_soln = uw.discretisation.MeshVariable("T", batmesh, 1, degree=2)
+v_soln  = uw.discretisation.MeshVariable("U", batmesh, batmesh.dim, degree=2, continuous=True)
+v_soln1  = uw.discretisation.MeshVariable("U1", batmesh, batmesh.dim, degree=2, continuous=True)
+p_soln  = uw.discretisation.MeshVariable("P", batmesh, 1, degree=1, continuous=True)
+t_soln  = uw.discretisation.MeshVariable("T", batmesh, 1, degree=2, varsymbol=r"T_0")
+t_soln1 = uw.discretisation.MeshVariable("T1", batmesh, 1, degree=2, varsymbol=r"T_1")
 cell_properties = uw.discretisation.MeshVariable("Bat", batmesh, 1, degree=0)
 
 
@@ -68,18 +70,24 @@ class bd(Enum):
 
 stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
 stokes.constitutive_model.Parameters.shear_viscosity_0 = 1
-stokes.tolerance = 1.0e-3
+stokes.tolerance = 1.0e-5
 
 unit_r_vec = batmesh.CoordinateSystem.unit_e_0
+unit_y_vec = batmesh.CoordinateSystem.unit_j 
 
 # free slip.
 # note with petsc we always need to provide a vector of correct cardinality.
 
-stokes.add_natural_bc(10000 * unit_r_vec.dot(v_soln.sym) * unit_r_vec, "Upper")
-stokes.bodyforce = Rayleigh_number * unit_r_vec * t_soln.sym[0] 
+stokes.add_essential_bc([0,0], "Upper")
+# stokes.add_natural_bc(1e6  * Rayleigh_number * unit_r_vec.dot(v_soln.sym) * unit_r_vec, "Upper")
 
-# Make the bat stagnant
-# stokes.bodyforce -= 10000 * (1-cell_properties.sym[0]) * v_soln.sym
+stokes.bodyforce = Rayleigh_number * unit_y_vec * t_soln.sym[0] + \
+                   Rayleigh_number * unit_y_vec * t_soln1.sym[0]
+
+
+stokes.petsc_options.setValue("ksp_monitor", None)
+stokes.petsc_options.setValue("snes_monitor", None)
+stokes.petsc_options.setValue("snes_min_it", 1)
 
 stokes.view()
 
@@ -92,48 +100,65 @@ class bd(Enum):
 
 adv_diff.constitutive_model = uw.constitutive_models.DiffusionModel
 adv_diff.constitutive_model.Parameters.diffusivity = 1
-adv_diff.f = 10 * (1-cell_properties.sym[0])
-
 adv_diff.add_dirichlet_bc(0.0, "Upper")
 
+
+adv_diff1 = uw.systems.AdvDiffusionSLCN(
+    batmesh,
+    u_Field=t_soln1,
+    V_fn=v_soln1,
+)
+
+adv_diff1.constitutive_model = uw.constitutive_models.DiffusionModel
+adv_diff1.constitutive_model.Parameters.diffusivity = 1
+adv_diff1.add_dirichlet_bc(0.0, "Upper")
+
 t_init = 1-cell_properties.sym[0]
 
-with batmesh.access(t_soln):
+with batmesh.access(t_soln, t_soln1):
     t_soln.data[:,0] = uw.function.evalf(t_init, t_soln.coords)
+    t_soln1.data[:,0] = uw.function.evalf(t_init, t_soln1.coords)
 
 
-# %%
-stokes.solve(zero_init_guess=True)
 
 # %%
-v_soln.view()
+stokes.solve(zero_init_guess=True, picard=0, verbose=False)
+with batmesh.access(v_soln1):
+    v_soln1.data[...] = -v_soln.data[...]
 
 # %%
+
 t_step = 0
 
 # %%
 
-for step in range(0, 50):
+for step in range(0, 250):
     stokes.solve(zero_init_guess=False)
-    delta_t = 2.0 * stokes.estimate_dt()
-    adv_diff.solve(timestep=delta_t, zero_init_guess=False)
+    with batmesh.access(v_soln1):
+        v_soln1.data[...] = -v_soln.data[...]
+
+    
+    delta_t = 5.0 * stokes.estimate_dt()
+    adv_diff.solve(timestep=delta_t, zero_init_guess=True)
+    adv_diff1.solve(timestep=delta_t, zero_init_guess=True)
 
     # stats then loop
     tstats = t_soln.stats()
 
     if uw.mpi.rank == 0:
         print("Timestep {}, dt {}".format(step, delta_t))
 
-    batmesh.write_timestep(
-        "Batmesh",
-        meshUpdates=True,
-        meshVars=[p_soln, v_soln, t_soln],
-        outputPath="output",
-        index=t_step,
-    )
+    if t_step%1==0:
+        batmesh.write_timestep(
+            "Batmesh",
+            meshUpdates=True,
+            meshVars=[p_soln, v_soln, t_soln, t_soln1],
+            outputPath="output",
+            index=t_step,
+        )
 
-    with batmesh.access(t_soln):
-        t_soln.data[:,0] = np.maximum(t_soln.data[:,0], uw.function.evalf(t_init, t_soln.coords))
+    # with batmesh.access(t_soln):
+    #     t_soln.data[:,0] = np.maximum(t_soln.data[:,0], uw.function.evalf(t_init, t_soln.coords))
 
     t_step += 1
 
@@ -146,69 +171,133 @@ class bd(Enum):
 
     pvmesh = vis.mesh_to_pv_mesh(batmesh)
     pvmesh.point_data["T"] = vis.scalar_fn_to_pv_points(pvmesh, t_soln.sym[0])
+    pvmesh.point_data["R"] = vis.scalar_fn_to_pv_points(pvmesh, batmesh.CoordinateSystem.R[0])
+    pvmesh.point_data["T1"] = vis.scalar_fn_to_pv_points(pvmesh, t_soln1.sym[0])
     pvmesh.point_data["V"] = vis.vector_fn_to_pv_points(pvmesh, v_soln.sym)
     pvmesh.point_data["L"] = vis.scalar_fn_to_pv_points(pvmesh, cell_properties.sym[0])
 
+    pvmesh_t = vis.meshVariable_to_pv_mesh_object(t_soln)
+    pvmesh_t.point_data["T"] = vis.scalar_fn_to_pv_points(pvmesh_t, t_soln.sym[0])
+    pvmesh_t.point_data["T1"] = vis.scalar_fn_to_pv_points(pvmesh_t, t_soln1.sym[0])
+
+    pvmesh_t.point_data["To"] = (pvmesh_t.point_data["T"] > 0.25).astype(float)
+    pvmesh_t.point_data["T1o"] = (pvmesh_t.point_data["T1"] > 0.25).astype(float)
+    pvmesh_t.point_data["TD"] = pvmesh_t.point_data["T"] - pvmesh_t.point_data["T1"]
+
+
     pl = pv.Plotter(window_size=(750, 750))
+    pl.enable_depth_peeling()
 
+    doughnut = pvmesh.clip_scalar(scalars="R", 
+                                  value=0.3, 
+                                  invert=False)
+    
     pvstream = pvmesh.streamlines_from_source(
-        pvmesh.cell_centers(), 
+        pv.PointSet(doughnut.cell_centers().points[::5]), 
         vectors="V", 
         integrator_type=45,
         surface_streamlines=True, 
         max_steps=1000,
         max_time=0.25,
     )
 
+
+    # pl.add_mesh(pvmesh, 
+    #             style="wireframe",
+    #             color="#FEFEF0",
+    #             opacity=0.1)
+
+    pl.add_mesh(pvstream, 
+                cmap=["orange", "white", "green"], 
+                show_scalar_bar=False,
+                opacity=0.1
+               )
+
+
     pl.add_mesh(
-                pvmesh,
+                pvmesh_t,
                 cmap="Oranges",
                 scalars="T",
+                opacity="sigmoid",
                 edge_color="Grey",
                 show_edges=False,
                 use_transparency=False,
-                opacity=1.0,
-                show_scalar_bar=True,
-                clim=[0,1]
+                show_scalar_bar=False,
                )
 
     pl.add_mesh(
-                pvmesh,
-                scalars="L",
-                cmap="Greys_r",
-                edge_color="Black",
-                opacity="L",
+                pvmesh_t,
+                copy_mesh=True,
+                cmap="Blues",
+                scalars="T1",
+                opacity="sigmoid",
+                edge_color="Grey",
                 show_edges=False,
-                use_transparency=True,
+                use_transparency=False,
                 show_scalar_bar=False,
-                
+                # clim=[0,2]
                )
-
-    # pl.add_mesh(pvstream, cmap="jet", show_scalar_bar=False)
+  
 
 
+    # pl.add_mesh(
+    #             pvmesh,
+    #             scalars="L",
+    #             cmap="Greys_r",
+    #             edge_color="Black",
+    #             opacity="L",
+    #             show_edges=False,
+    #             use_transparency=True,
+    #             show_scalar_bar=False,
+    #            )
+
 
     # pl.add_arrows(pvmesh.points, pvmesh.point_data["V"],
-    #               mag=float(20.0/Rayleigh_number.sym),
+    #               mag=float(0.2/Rayleigh_number.sym),
     #              show_scalar_bar=False)
 
 
     pl.show(jupyter_backend='html')
 
+# %%
+0/0
+
+# %%
+
 # %%
 ## Check the results saved to file
 
 expt_name = "Batmesh"
 output_path="output"
-step = 20
+step = 50
 
 v_soln.read_timestep(expt_name, "U", step, outputPath=output_path)
 p_soln.read_timestep(expt_name, "P", step, outputPath=output_path)
 t_soln.read_timestep(expt_name, "T", step, outputPath=output_path)
+t_soln1.read_timestep(expt_name, "T1", step, outputPath=output_path)
+
 
 
+# %%
+t_soln.stats()
+
+# %%
+t_soln1.stats()
 
 # %%
 adv_diff
 
 # %%
+adv_diff1
+
+# %%
+adv_diff1.DuDt.bdf(1)
+
+# %%
+adv_diff.DuDt.bdf(1)
+
+# %%
+batmesh.dm.view()
+
+# %% [markdown]
+# #