.. _subroutines-examples:

********
Examples
********

.. _use-pressure:

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:

.. code-block:: sierrainput

   #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 :numref:`use-topoq`) on groups of elements, though this would increase subroutine complexity.

.. code-block:: fortran

   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

.. _use-error:

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:

.. code-block:: sierrainput

   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:

.. code-block:: sierrainput

   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:

.. code-block:: fortran

         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

.. _use-transf:

Transform Output Stresses to a Cylindrical Coordinate System
============================================================

The following code is a user subroutine to transform element stresses in global :math:`x`, :math:`y`, and :math:`z` coordinates to a global cylindrical coordinate system. This subroutine could be used to transform the relatively meaningless shell stress in :math:`x`, :math:`y`, and :math:`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.

.. code-block:: sierrainput

   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:

.. code-block:: fortran

         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
