Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def _calibrate_pressureField(self):
surfaceArea = uw.utils.Integral(fn=1.0, mesh=self.mesh,
integrationType='surface',
surfaceIndexSet=self.topWall)
surfacepressureFieldIntegral = uw.utils.Integral(
fn=self.pressureField,
mesh=self.mesh,
integrationType='surface',
surfaceIndexSet=self.topWall
)
(area,) = surfaceArea.evaluate()
(p0,) = surfacepressureFieldIntegral.evaluate()
offset = p0/area
self.pressureField.data[:] -= offset
self.pressSmoother.smooth()
for material in self.materials:
if material.viscosity:
# **Nusselt number**
#
# The Nusselt number is the ratio between convective and conductive heat transfer
#
# \\[
# Nu = -h \frac{ \int_0^l \partial_z T (x, z=h) dx}{ \int_0^l T (x, z=0) dx}
# \\]
#
#
#
#
# In[15]:
nuTop = uw.utils.Integral( fn=temperatureField.fn_gradient[1],
mesh=mesh, integrationType='Surface',
surfaceIndexSet=mesh.specialSets["MaxJ_VertexSet"])
nuBottom = uw.utils.Integral( fn=temperatureField,
mesh=mesh, integrationType='Surface',
surfaceIndexSet=mesh.specialSets["MinJ_VertexSet"])
# In[16]:
nu = - nuTop.evaluate()[0]/nuBottom.evaluate()[0]
print('Nusselt number = {0:.6f}'.format(nu))
# **RMS velocity**
#
advector = uw.systems.SwarmAdvector( swarm=swarm, velocityField=velocityField, order=2 )
# Analysis tools
# -----
#
# **RMS velocity**
#
# Set up integrals used to calculate the RMS velocity.
# In[15]:
vdotv = fn.math.dot( velocityField, velocityField )
v2sum_integral = uw.utils.Integral( mesh=mesh, fn=vdotv )
volume_integral = uw.utils.Integral( mesh=mesh, fn=1. )
# Main simulation loop
# -----
#
# The main time stepping loop begins here. Before this the time and timestep are initialised to zero and the output statistics arrays are set up. Also the frequency of outputting basic statistics to the screen is set in steps_output.
#
# Note that there are two ``advector.integrate`` steps, one for each swarm, that need to be done each time step.
# In[16]:
# Stepping. Initialise time and timestep.
time = 0.
step = 0
nsteps = 10
if(uw.rank()==0):
--------
>>> mesh = uw.mesh.FeMesh_Cartesian(minCoord=(0.0,0.0), maxCoord=(1.0,2.0))
>>> fn_1 = uw.function.misc.constant(2.0)
>>> np.allclose( mesh.integrate( fn_1 )[0], 4 )
True
>>> fn_2 = uw.function.misc.constant(2.0) * (0.5, 1.0)
>>> np.allclose( mesh.integrate( fn_2 ), [2,4] )
True
"""
# julian, i've dumbed this down for now as i'm not sure how to handle the
# circular dependency. i think the performance hit will be generally
# negligible
_volume_integral = uw.utils.Integral(mesh=self, fn=fn)
return _volume_integral.evaluate()
# open hdf5 file
h5f = h5py.File(name=filename, mode="r", driver='mpio', comm=MPI.COMM_WORLD)
dset = h5f.get('data')
if dset == None:
raise RuntimeError("Can't find 'data' in file '{0}'.\n".format(filename))
if dset.shape[1] != self.particleCoordinates.data.shape[1]:
raise RuntimeError("Cannot load file data on current swarm. Data in file '{0}', " \
"has {1} components -the particlesCoords has {2} components".format(filename, dset.shape[1], self.particleCoordinates.data.shape[1]))
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
nProcs = comm.Get_size()
if rank == 0 and verbose:
bar = uw.utils._ProgressBar( start=0, end=dset.shape[0]-1, title="loading "+filename)
# try and read the procCount attribute & assume that if nProcs in .h5 file
# is equal to the current no. procs then the particles will be distributed the
# same across the processors. (Danger if different discretisations are used... i think)
# else try and load the whole .h5 file.
# we set the 'offset' & 'size' variables to achieve the above
offset = 0
totalsize = size = dset.shape[0] # number of particles in h5 file
if try_optimise:
procCount = h5f.attrs.get('proc_offset')
if procCount is not None and nProcs == len(procCount):
for p_i in xrange(rank):
offset += procCount[p_i]
size = procCount[rank]
def solve_velocityField(self):
fnVel = (-self._pressureField.fn_gradient + self.fn_bodyforce) * self._kMatTerm.fn
if self._swarmVarVelocity:
self._swarmVarVelocity.data[:] = fnVel.evaluate(self._swarm)
if self._velocityField:
velproj = uw.utils.MeshVariable_Projection(self._velocityField,fnVel,self._swarm)
velproj.solve()
if not self.velocityField:
raise ValueError("Please link a Velocity Field to the Isostasy solver")
if not self.materialIndexField:
raise ValueError("Please link a Material Field to the Isostasy solver")
if not self.materialIndexField:
raise ValueError("Please link a Material Field to the Isostasy solver")
if not self.densityFn:
raise ValueError("Please link a Density Field to the Isostasy solver")
if not self.reference_mat:
raise ValueError("Please define a reference Material for the Isostasy solver")
# Create utilities
self.MaterialIndexFieldFloat = self.swarm.add_variable(dataType="double", count=1 )
self.DensityVar = uw.mesh.MeshVariable(self.mesh, nodeDofCount=1)
self.MaterialVar = uw.mesh.MeshVariable(self.mesh, nodeDofCount=1)
self.projectorDensity = uw.utils.MeshVariable_Projection(self.DensityVar,
self.densityFn,
type=0 )
self.projectorMaterial = uw.utils.MeshVariable_Projection(self.MaterialVar,
self.MaterialIndexFieldFloat, type=0 )
if not self.mesh._cself.isRegular:
raise TypeError("You are using an irregular mesh: \
isostasy module only works with regular meshes")
if self.surface is not None and not isinstance(self.surface, uw.swarm._swarm.Swarm):
raise TypeError("'surface' must be a tuple'")
self.initialized = True
"""
# use barrier as there are some file open operations below
# and we need to ensure that all procs have finished writing
# before we try and open any files.
uw.mpi.barrier()
if uw.mpi.rank == 0:
if not isinstance(varname, str):
raise ValueError("'varname' must be of type str")
if not isinstance(swarmname, str):
raise ValueError("'swarmname' must be of type str")
if not isinstance(filename, str):
raise ValueError("'filename' must be of type str")
if not isinstance(swarmSavedData, uw.utils.SavedFileData ):
raise ValueError("'swarmSavedData' must be of type SavedFileData")
if not isinstance(varSavedData, uw.utils.SavedFileData ):
raise ValueError("'varSavedData' must be of type SavedFileData")
if not isinstance(modeltime, (int,float)):
raise ValueError("'modeltime' must be of type int or float")
modeltime = float(modeltime) # make modeltime a float
# get the elementMesh - if self is a subMeshed variable get the parent
if self.swarm != swarmSavedData.pyobj:
raise RuntimeError("'swarmSavedData file doesn't correspond to the object's swarm")
if not filename.lower().endswith('.xdmf'):
filename += '.xdmf'
# the xmf file is stored in 'string'
# 1st write header
string = uw.utils._xdmfheader()
"""
swarmVgrad = velocityField.fn_gradient.evaluate(swarm)
stretching.data[:,0] += dt * (swarmVgrad[:,0] * stretching.data[:,0] + swarmVgrad[:,1] * stretching.data[:,1])
stretching.data[:,1] += dt * (swarmVgrad[:,2] * stretching.data[:,0] + swarmVgrad[:,3] * stretching.data[:,1])
# update plastic strain on those swarm particles that yielded
swarmYield = viscosityFn.evaluate(swarm) < etaV
swarmStrainRateInv = strainRate_2ndInvariantFn.evaluate(swarm)
particleIsViscPlastic = materialVariable.evaluate(swarm) == materialV
plasticStrainIncrement = dt * np.where(swarmYield, np.where(particleIsViscPlastic, swarmStrainRateInv, 0.0) , 0.0 )
plasticStrain.data[:] += plasticStrainIncrement
if (step%5 ==0):
# figVelocityPressure.save_image( outputPath + "figVP-" + str(step).zfill(4))
projectorStress = uw.utils.MeshVariable_Projection( meshDevStress,
fn.tensor.second_invariant(devStressFn), type=0 )
projectorStress.solve()
# figStrain.save_image( outputPath + "figStrain-" + str(step).zfill(4))
if uw.rank()==0:
print('step = {0:3d}; time = {1:.3e}; xmax = {2:.3f}; pmax = {3:.4f}; cpu time = {4:.2e}'
.format(step, time, xmax, plasticStrain.evaluate(swarm).max(), cpuTime.clock()-startTime))
time += dt
step += 1
# Post simulation images
# -----
#
def vertical_gradient_value(self, fn=None):
"""Returns vertical gradient within the shell of scalar uw function"""
## Need to check here that the supplied function is
## valid and is a scalar. Mask function is to eliminate the
## effect of the values in the core.
rGrad = uw.utils.Integral(fn * (self.unit_heightFn-self._c0)*self.maskFn, self).evaluate()[0] / self._c1
return rGrad