10.2.4. Examples
10.2.4.1. Pressure as a Function of Space and Time
(The following example provides functionality – a blast load on a surface – more applicable to Presto than Adagio. It is included in both manuals as it is instructive in the general use of a user-defined subroutine.)
The following code is an example of a user subroutine to compute blast pressures on a group of faces. The blast pressure simulates a blast occurring at a specified position and time. The blast wave radiates out from the center of the blast and dissipates over time. This subroutine is included in the input file as follows:
#In the SIERRA scope:
user subroutine file = blast_load.F
#In the region scope:
begin pressure
surface = surface_1
surface subroutine = blast_load
subroutine real parameter: pos_x = 5.0
subroutine real parameter: pos_y = 5.0
subroutine real parameter: pos_z = 1.6
subroutine real parameter: wave_speed = 1.5e+02
subroutine real parameter: blast_time = 0.0
subroutine real parameter: blast_energy = 1.0e+09
subroutine real parameter: blast_wave_width = 0.75
end pressure
The FORTRAN 77 subroutine listing follows. It would be possible to increase the speed of this subroutine by calling the topology functions (see Section 10.2.1.2.6) on groups of elements, though this would increase subroutine complexity.
c
c Subroutine to simulate a blast load on a surface
c
subroutine blast_load(num_faces, num_vals,
& eval_time, faceID, pressure, flags, err_code)
implicit none
c
c Subroutine input arguments
c
integer num_faces
double precision eval_time
integer faceID(num_faces)
integer num_vals
c
c Subroutine output arguments
c
double precision pressure(num_vals, num_faces)
integer flags(num_faces)
integer err_code
c
c Variables to hold the subroutine parameters
c
double precision pos_x, pos_y, pos_z, wave_speed,
& blast_time, blast_energy,
& blast_wave_width
c
c Local variables
c
integer iface, inode
integer cur_face_id, face_topo, num_nodes
integer num_comp_check
double precision dist, blast_o_rad, blast_i_rad
double precision blast_volume, blast_pressure
integer query_error
double precision face_center(3)
c
c Create some static variables to hold queried
c information. Assume no face has more than 10
c nodes
c
double precision face_nodes(10)
double precision face_coords(3, 10)
c
c Extract the subroutine parameters
c
call aupst_get_real_param("pos_x",pos_x,query_error)
call aupst_get_real_param("pos_y",pos_y,query_error)
call aupst_get_real_param("pos_z",pos_z,query_error)
call aupst_get_real_param("wave_speed",wave_speed,
& query_error)
call aupst_get_real_param("blast_energy",
& blast_energy,query_error)
call aupst_get_real_param("blast_time",
& blast_time,query_error)
call aupst_get_real_param("blast_wave_width",
& blast_wave_width, query_error)
c
c Determine the outer radius of the blast wave
c
blast_o_rad = (eval_time - blast_time) * wave_speed
if(blast_o_rad .le. 0.0) return;
c
c Determine the inner radius of the blast wave
c
blast_i_rad = blast_o_rad - blast_wave_width
if(blast_i_rad .le. 0.0) blast_i_rad = 0.0
c
c Determine the total volume the blast wave occupies
c
blast_volume = 3.1415 * (4.0/3.0) *
& (blast_o_rad**2 - blast_i_rad**2)
c
c Determine the total pressure on faces inside the
c blast wave
c
blast_pressure = blast_energy / blast_volume
c
c Loop over all faces in the set
c
do iface = 1, num_faces
c
c Extract the topology of the current face
c
cur_face_id = faceID(iface)
call aupst_get_face_topology(1, cur_face_id,
& face_topo, query_error)
c
c Determine the number of nodes of the current face
c
num_nodes = mod(face_topo,100)
c
c Extract the node ids for nodes contained in the current
c face
c
call aupst_get_face_nodes(1, cur_face_id,
& face_nodes, query_error)
c
c Extract the nodal coordinates of the face nodes
c
call aupst_get_node_var(num_nodes, 3, face_nodes,
& face_coords, "coordinates", query_error)
c
c Compute the centroid of the face
c
face_center(1) = 0.0
face_center(2) = 0.0
face_center(3) = 0.0
do inode = 1, num_nodes
face_center(1) = face_center(1) +
& face_coords(1,inode)
face_center(2) = face_center(2) +
& face_coords(2,inode)
face_center(3) = face_center(3) +
& face_coords(3,inode)
enddo
face_center(1) = face_center(1)/num_nodes
face_center(2) = face_center(2)/num_nodes
face_center(3) = face_center(3)/num_nodes
c
c Determine the distance from the current face
c to the blast center
c
dist = sqrt((face_center(1) - pos_x)**2 +
& (face_center(2) - pos_y)**2 +
& (face_center(3) - pos_z)**2)
c
c Apply pressure to the current face if it falls within
c the blast wave
c
if(dist .ge. blast_i_rad .and.
& dist .le. blast_o_rad) then
pressure(1,iface) = blast_pressure
else
pressure(1,iface) = 0.0
endif
enddo
err_code = 0
end
10.2.4.2. Error Between a Computed and an Analytic Solution
The following code is a user subroutine to compute the error between Sierra/SM-computed results and results from an analytic manufactured solution. This subroutine is called by a USER OUTPUT command block immediately before producing an output Exodus file. The error for the mesh is computed by taking the squared difference between the computed and analytic displacements at every node. Finally, a global sum of the error is produced along with the square root norm of the error.
This user subroutine requires a user variable, which is defined in the Sierra/SM input file. The command block for the user variable specified in this user subroutine is as follows:
begin user variable conv_error
type = global real length = 1
global operator = sum
initial value = 0.0
end user variable conv_error
The subroutine is called in the Sierra/SM input file as follows:
begin user output
node set = nodelist_10
node set subroutine = conv0_error
subroutine real parameter: char_length = 1.0
subroutine real parameter: char_time = 1.0e-3
subroutine real parameter: x_offset = 0.0
subroutine real parameter: y_offset = 0.0
subroutine real parameter: z_offset = 0.0
subroutine real parameter: t_offset = 0.0
subroutine real parameter: u0 = 0.01
subroutine real parameter: v0 = 0.02
subroutine real parameter: w0 = 0.03
subroutine real parameter: alpha = 1.0
subroutine real parameter: youngs_modulus = 10.0e6
subroutine real parameter: poissons_ratio = 0.3
subroutine real parameter: density = 0.0002588
subroutine real parameter: num_nodes = 125.0
end user output
The FORTRAN listing for the subroutine is as follows:
subroutine conv0_error(num_nodes, num_vals,
& eval_time, nodeID, values, flags, ierror)
implicit none
integer num_nodes
integer num_vals
double precision eval_time
integer nodeID(num_nodes)
double precision values(1)
integer flags(1)
integer ierror
c
c Local vars
c
integer inode
integer error_code
double precision clength, ctime, xoff, yoff, zoff, toff
double precision zero, one, two, three, four, nine
double precision mod_coords(3,3000)
double precision cdispl(3,3000)
integer num_comp_check
double precision expat
double precision x, y, z, t
double precision u0, v0, w0, alpha
double precision pi
double precision half
double precision mdisplx, mdisply, mdisplz
double precision xdiff, ydiff, zdiff
double precision conv_error
double precision numnod
pi = 3.141592654
half = 0.5
zero = 0.0
one = 1.0
two = 2.0
three = 3.0
four = 4.0
nine = 9.0
c
c Check that the nodal coordinates will fit into the
c statically allocated array
c
if(num_nodes .gt. 3000) then
write(6,*) ERROR in sphere disp, cdispl,
& num_nodes exceeds static array size
ierror = 1
return
endif
c
c Extract the model coordinates for all nodes
c
call aupst_check_node_var(num_nodes, num_comp_check,
& nodeID, "model_coordinates",
& ierror)
if(ierror .ne. 0) return
if(num_comp_check .ne. 3) return
call aupst_get_node_var(num_nodes, num_comp_check,
& nodeID, mod_coords, "model_coordinates",
& ierror)
c
c Extract the computed displacements for all nodes
c
call aupst_check_node_var(num_nodes, num_comp_check,
& nodeID, "displacement",
& ierror)
if(ierror .ne. 0) return
if(num_comp_check .ne. 3) return
call aupst_get_node_var(num_nodes, num_comp_check,
& nodeID, cdispl, "displacement",
& ierror)
c
c Extract the subroutine parameters.
c
call aupst_get_real_param("char_length",
& clength,error_code)
call aupst_get_real_param("char_time",
& ctime,error_code)
call aupst_get_real_param("x_offset",xoff,error_code)
call aupst_get_real_param("y_offset",yoff,error_code)
call aupst_get_real_param("z_offset",zoff,error_code)
call aupst_get_real_param("t_offset",toff,error_code)
call aupst_get_real_param("u0",u0,error_code)
call aupst_get_real_param("v0",v0,error_code)
call aupst_get_real_param("w0",w0,error_code)
call aupst_get_real_param("alpha",alpha,error_code)
call aupst_get_real_param("num_nodes",
& numnod,error_code)
c
c Calculate a solution scaling factor
c
expat = half * ( one - cos( pi * eval_time / ctime ) )
c
c Compute the expected solution at each node and do a
c sum of the differences from the analytic solution
c
conv_error = zero
do inode = 1, num_nodes
c
c Set the displacement value from the manufactured solution
c
x = ( mod_coords(1,inode) - xoff ) / clength
y = ( mod_coords(2,inode) - yoff ) / clength
z = ( mod_coords(3,inode) - zoff ) / clength
c
mdisplx = u0 * sin(x) * cos(two*y) * cos(three*z)
* * expat
mdisply = v0 * cos(three*x) * sin(y) * cos(two*z)
* * expat
mdisplz = w0 * cos(two*x) * cos(three*y) * sin(z)
* * expat
c
xdiff = mdisplx - cdispl(1,inode)
ydiff = mdisply - cdispl(2,inode)
zdiff = mdisplz - cdispl(3,inode)
conv_error = conv_error + xdiff*xdiff
* + ydiff*ydiff
* + zdiff*zdiff
c
enddo
c
ierror = 0
c
c Do a parallel sum of the squared errors and extract
c the total summed value on all processors
c
call aupst_put_global_var(1,conv_error,
& "conv_error","sum",ierror)
call aupst_get_global_var(1,conv_error,
& "conv_error",ierror)
c
c Take the square root of the errors and store that as
c the net error norm
c
conv_error = sqrt(conv_error) / sqrt(numnod)
call aupst_put_global_var(1,conv_error,
& "conv_error","none",ierror)
c
return
end
10.2.4.3. Transform Output Stresses to a Cylindrical Coordinate System
The following code is a user subroutine to transform element stresses in global \(x\), \(y\), and \(z\) coordinates to a global cylindrical coordinate system. This subroutine could be used to transform the relatively meaningless shell stress in \(x\), \(y\), and \(z\) coordinates to more meaningful tangential, hoop, and radial stresses. The subroutine is called from a USER OUTPUT command block. It reads in the old stresses, transforms them, and writes them back out to a user-created scratch variable, defined via a USER VARIABLE command block, for output.
begin user variable cyl_stress
type = element sym_tensor length = 1
initial value = 0.0 0.0 0.0 0.0 0.0 0.0
end user variable
begin user output
block = block_1
element block subroutine = aupst_cyl_transform
subroutine string parameter: origin_point = Point_O
subroutine string parameter: z_point = Point_Z
subroutine string parameter: xz_point = Point_XZ
subroutine string parameter: input_stress = memb_stress
subroutine string parameter: output_stress = cyl_stress
end user output
The FORTRAN listing for the subroutine is as follows:
subroutine aupst_cyl_transform(num_elems, num_vals,
* eval_time, elemID, values, flags, ierror)
implicit none
#include<framewk/Fmwk_type_sizes_decl.par>
#include<framewk/Fmwk_type_sizes.par>
c
c Subroutine Arguments
c
c num_elems: Input: Number of elements to calculate on
c num_vals : Input: Ignored
c eval_time: Input: Time at which to evaluate the stress.
c elemID : Input: Global sierra IDs of the input elements
c values : I/O : Ignored, stress will be stored manually
c flags : I/O : Ignored
c ierror :Output: Returns non-zero if an error occurs
c
integer num_elems
integer num_vals
double precision eval_time
integer elemID(num_elems)
double precision values(1)
integer flags(1)
integer ierror
c
c Fortran cannot dynamically allocate memory, thus worksets
c will be iterated over by chucks each of size chunk_size.
c
integer chunk_size
parameter (chunk_size = 100)
integer chunk_ids(chunk_size)
c
c Subroutine parameter data
c
character*80 origin_point_name
double precision origin_point(3)
character*80 z_point_name
double precision z_point(3)
character*80 xz_point_name
double precision xz_point(3)
character*80 input_stress_name
character*80 output_stress_name
c
c Local element data for centroids and rotation vectors
c
double precision cent(3)
double precision centerline_pos(3)
double precision dot_prod
double precision z_vec(3)
double precision r_vec(3)
double precision theta_vec(3)
double precision rotation_tensor(9)
c
c Chunk data storage
c
double precision elem_centroid(3, chunk_size)
double precision input_stress_val(6, chunk_size)
double precision output_stress_val(6, chunk_size)
c
c Simple iteration variables
c
integer error_code
integer ichunk, ielem
integer zero_elem, nel
c
c Extract the current subroutine parameters. origin_point
c is the origin of the coordinate system
c z_point is a point on the z axis of the coordinate system
c xz_point is a point on the xz plane
c
call aupst_get_string_param("origin_point",
& origin_point_name,
& error_code)
call aupst_get_string_param("z_point",
& z_point_name,
& error_code)
call aupst_get_string_param("xz_point",
& xz_point_name,
& error_code)
call aupst_get_string_param("input_stress",
& input_stress_name,
& error_code)
call aupst_get_string_param("output_stress",
& output_stress_name,
& error_code)
c
c Use the point names to look up the coordinates of each
c relevant point
c
call aupst_get_point(origin_point_name, origin_point,
& error_code)
call aupst_get_point(z_point_name, z_point,
& error_code)
call aupst_get_point(xz_point_name, xz_point,
& error_code)
c
c Compute the z axis vector
c
z_vec(1) = z_point(1) - origin_point(1)
z_vec(2) = z_point(2) - origin_point(2)
z_vec(3) = z_point(3) - origin_point(3)
c
c Transform z_vec into a unit vector, abort if it is invalid
c
call aupst_unitize_vector(z_vec, ierror)
if(ierror .ne. 0) return
c
c Loop over chunks of the data arrays
c
do ichunk = 1, (num_elems/chunk_size + 1)
c
c Determine the first and last element number for the
c current chunk of elements
c
zero_elem = (ichunk-1) * chunk_size
if((zero_elem + chunk_size) .gt. num_elems) then
nel = num_elems - zero_elem
else
nel = chunk_size
endif
c
c Copy the elemIDs for all elems in the current chunk to a
c temporary array
c
do ielem = 1, nel
chunk_ids(ielem) = elemID(zero_elem + ielem)
enddo
c
c Extract the element centroids and stresses
c
call aupst_get_elem_centroid(nel, chunk_ids,
& elem_centroid,
& ierror)
call aupst_get_elem_var(nel, 6, chunk_ids,
& input_stress_val,
& input_stress_name, ierror)
c
c Loop over each element in the current chunk
c
do ielem = 1, nel
c
c Find the closest point on the cylinder centerline axis
c to the element centroid
c
cent(1) = elem_centroid(1, ielem) - origin_point(1)
cent(2) = elem_centroid(2, ielem) - origin_point(2)
cent(3) = elem_centroid(3, ielem) - origin_point(3)
dot_prod = cent(1) * z_vec(1) +
& cent(2) * z_vec(2) +
& cent(3) * z_vec(3)
centerline_pos(1) = z_vec(1) * dot_prod
centerline_pos(2) = z_vec(2) * dot_prod
centerline_pos(3) = z_vec(3) * dot_prod
c
c Compute the current normal radial vector
c
r_vec(1) = cent(1) - centerline_pos(1)
r_vec(2) = cent(2) - centerline_pos(2)
r_vec(3) = cent(3) - centerline_pos(3)
call aupst_unitize_vector(r_vec, ierror)
if(ierror .ne. 0) return
c
c Compute the current hoop vector
c
theta_vec(1) = z_vec(2)*r_vec(3) - r_vec(2)*z_vec(3)
theta_vec(2) = z_vec(3)*r_vec(1) - r_vec(3)*z_vec(1)
theta_vec(3) = z_vec(1)*r_vec(2) - r_vec(1)*z_vec(2)
c
c The r, theta, and z vectors describe the new stress
c coordinate system, Transform the input stress tensor
c in x,y,z coords to the output stress tensor in r, theta,
c and z coords use the unit vectors to create a rotation
c tensor
c
rotation_tensor(k_f36xx) = r_vec(1)
rotation_tensor(k_f36yx) = r_vec(2)
rotation_tensor(k_f36zx) = r_vec(3)
rotation_tensor(k_f36xy) = theta_vec(1)
rotation_tensor(k_f36yy) = theta_vec(2)
rotation_tensor(k_f36zy) = theta_vec(3)
rotation_tensor(k_f36xz) = z_vec(1)
rotation_tensor(k_f36yz) = z_vec(2)
rotation_tensor(k_f36zz) = z_vec(3)
c
c Rotate the current stress tensor to the new configuration
c
call fmth_rotate_symten33(1, 1, 0, rotation_tensor,
& input_stress_val(1,ielem),
& output_stress_val(1,ielem))
enddo
c
c Store the new stress
c
call aupst_put_elem_var(nel, 6, chunk_ids,
& output_stress_val,
& output_stress_name, ierror)
enddo
ierror = 0
end