Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
solA = fn.analytic.SolCx()
stokesSystem = uw.systems.Stokes(velocityField,pressureField,solA.fn_viscosity,solA.fn_bodyforce,conditions=[freeslip,], rtolerance=1.e-5)
# To set up Uzawa solver with mumps
#from libUnderworld import petsc
#petsc.OptionsInsertString("-Uzawa_velSolver_pc_factor_mat_solver_package mumps -Uzawa_velSolver_ksp_type preonly -Uzawa_velSolver_pc_type lu -Uzawa_velSolver_ksp_converged_reason -Uzawa_velSolver_ksp_view")
#Running Uzawa solve (current default)
#stokesSystem.solve()
# In[3]:
#Run the BSSCR Solver
# can optionally set penalty this way
solver=uw.systems.Solver(stokesSystem)
solver.options.A11.ksp_rtol=1e-4
solver.options.scr.ksp_rtol=1e-3
solver.options.main.Q22_pc_type='uwscale'
#solver.options.mg.active=False
#solver.options.A11.set_direct()
#solver.options.A11.list()
#solve
solver.solve()
velocityField = uw.mesh.MeshVariable(linearMesh,2)
velocityField.data[:] = (0.,0.)
pressureField = uw.mesh.MeshVariable(constantMesh,1)
pressureField.data[:] = 0.
# freeslip
IWalls = linearMesh.specialSets["MinI_VertexSet"] + linearMesh.specialSets["MaxI_VertexSet"]
JWalls = linearMesh.specialSets["MinJ_VertexSet"] + linearMesh.specialSets["MaxJ_VertexSet"]
freeslip = uw.conditions.DirichletCondition(velocityField, (IWalls, JWalls))
solA = fn.analytic.SolCx()
stokesSystem = uw.systems.Stokes(velocityField,pressureField,solA.fn_viscosity,solA.fn_bodyforce,conditions=[freeslip,], rtolerance=1.e-5)
# To set up Uzawa solver with mumps
#from libUnderworld import petsc
#petsc.OptionsInsertString("-Uzawa_velSolver_pc_factor_mat_solver_package mumps -Uzawa_velSolver_ksp_type preonly -Uzawa_velSolver_pc_type lu -Uzawa_velSolver_ksp_converged_reason -Uzawa_velSolver_ksp_view")
#Running Uzawa solve (current default)
#stokesSystem.solve()
# In[3]:
#Run the BSSCR Solver
# can optionally set penalty this way
solver=uw.systems.Solver(stokesSystem)
solver.options.A11.ksp_rtol=1e-4
solver.options.scr.ksp_rtol=1e-3
solver.options.main.Q22_pc_type='uwscale'
mswarm_home_pts.data[:] = mswarm.particleCoordinates.data[:]
# if np.any(localID == -1):
# print("{} - particles missing: {}".format(uw.mpi.rank, np.where(localID == -1).shape[0]), flush=True )
# Make these variables accessible
self._mswarm = mswarm
self._mswarm_global_particles = mswarm.particleGlobalCount
self._mswarm_map = mswarm_map
self._mswarm_home_pts = mswarm_home_pts
self._mswarm_phiStar = mswarm_phiStar
if self._mswarm_advector == None:
madvector = uw.systems.SwarmAdvector(velocityField=self.vField, swarm=self._mswarm)
self._mswarm_advector = madvector
return
# In[17]:
stokes = uw.systems.Stokes( velocityField = velocityField,
pressureField = pressureField,
conditions = [freeslipBC,],
fn_viscosity = viscosity,
fn_bodyforce = buoyancyFn )
# Implicit Stokes solver
solver = uw.systems.Solver( stokes )
solver.set_inner_method("mumps")
solver.set_penalty(1.0e7)
# advDiff uses an explicit / timestepping approach
advDiff = uw.systems.AdvectionDiffusion( phiField = temperatureField,
phiDotField = temperatureDotField,
velocityField = velocityField,
fn_diffusivity = 1.0,
conditions = [tempBC,] )
# System setup
# -----
#
# **Setup a Stokes system**
#
# In[13]:
stokesPIC = uw.systems.Stokes( velocityField = velocityField,
pressureField = pressureField,
conditions = [freeslipBC,],
fn_viscosity = viscosity,
fn_bodyforce = buoyancyFn )
# get the default stokes equation solver
solver = uw.systems.Solver( stokesPIC )
# **Create an advection diffusion system**
#
# In[14]:
advDiff = uw.systems.AdvectionDiffusion( phiField = temperatureField,
phiDotField = temperatureDotField,
velocityField = velocityField,
fn_diffusivity = 1.0,
conditions = [tempBC,] )
# Analysis tools
# -----
def init_stokes_system(self):
self.buoyancyFn = self.densityFn * self.gravity
if any([material.viscosity for material in self.materials]):
stokes = uw.systems.Stokes(velocityField=self.velocityField,
pressureField=self.pressureField,
conditions=self._velocityBCs,
fn_viscosity=self.viscosityFn,
fn_bodyforce=self.buoyancyFn,
fn_one_on_lambda = None)
self.solver = uw.systems.Solver(stokes)
fn_pos = fn.coord() - centre
fn_conds.append( (fn.math.dot( fn_pos, fn_pos ) < sphereRadius**2, 1.5) )
fn_conds.append( (True, 1.0) )
fn_density = fn.branching.conditional( fn_conds )
fn_gravity = fn.misc.constant(9.8) * (0., 0., -1.)
fn_buoyancy = fn_gravity * fn_density
v_const = fn.misc.constant([4.0,2.,0.])
f_const = fn.misc.constant([1.0,2.,0.])
e_const = fn.misc.constant(1.0)
model = uw.systems.pl_StokesModel(filename=None)
model.SetViscosity(f_const[0])
model.SetRHS(fn_buoyancy)
model.Solve()
self.empty = False
# Should do some checking first
self.mesh = mesh
self.velocity = velocityField
self.thickness = fthickness
self.ID = fID
self.insidePt = insidePt
self.director = None
# Set up the swarm and variables on all procs
self.swarm = uw.swarm.Swarm( mesh=self.mesh, particleEscape=True )
self.director = self.swarm.add_variable( dataType="double", count=2)
self._swarm_advector = uw.systems.SwarmAdvector( swarm=self.swarm,
velocityField=self.velocity, order=2 )
self.swarm.add_particles_with_coordinates(np.stack((pointsX, pointsY)).T)
self.director.data[...] = 0.0
self._update_kdtree()
self._update_surface_normals()
return
import underworld as uw
from underworld import libUnderworld
from libUnderworld import petsc_layer as pl
from underworld import function as fn
import numpy as np
model = uw.systems.pl_PoissonModel()
# these following 2 fn are made available by the model
temperature = model.fn_temperature
'''
PI = fn.misc.constant(3.14159265359)
my_k = fn.math.sin( PI * ( fn_r / 5. - 1. ) )
'''
someFn = fn.misc.constant(-1)
fn_r = fn.math.sqrt( fn.math.dot( fn.coord(), fn.coord() ) )
fn_r2 = (10.*temperature+1.)*fn_r
benFn = fn.branching.conditional( [ (fn.coord()[2] > 0.0, fn_r),
(True , fn_r2) ] )