How to use the underworld.systems function in underworld

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) ] )