Added version documents and example boundary condition file for DOUAR-WSMP v0.1

c1427018 · Dave Whipp · 19af5802 · c1427018 · c1427018 · c1427018
Commit c1427018 authored 14 years ago by Dave Whipp
--- a/src/ChangeLog_v0.1
+++ b/src/ChangeLog_v0.1
+--------------------------------------------------------------------------------
+DOUAR-WSMP v0.1 - 23/03/2011
+--------------------------------------------------------------------------------
+Changes listed with revision number in the douar_newiso branch (r###) in the
+following format:
+--------------------------------------------------------------------------------
+r###:
+    <Description of the change or additional feature>
+    - <Files affected, file2, file3...>
+--------------------------------------------------------------------------------
+
+Major changes since 2009:
+--------------------------------------------------------------------------------
+r175-183,195-196,198,216,218-223,225,227-233,236,244,247,249-255,258-261:
+    Added new define_isostasy_bc.f90 routine to allow velocity boundary
+    conditions to adjust to the crustal subsidence/uplift.
+    - define_cloud.f90, define_isostasy_bc.f90, DOUAR.f90,
+      interpolate_ov_on_osolve.f90, isostasy.f90, Makefile,
+      module_definitions.f90, octree_interpolate_six.f90 (new@251),
+      read_input_file.f90, write_global_output.f90
+
+r234:
+    Added modifications suggested by J. Braun to be able to vertically compress
+    the FE mesh by factor vex
+    - define_cloud.f90, define_surface.f90, erosion.f90, make_matrix.f90,
+      make_pressure.f90, module_definitions.f90, move_cloud.f90,
+      move_surface.f90, read_input_file.f90, write_global_output.f90
+
+r257,262,267,271-274,292-293,296:
+    Calculate and apply isostasy after surfaces are advected
+    - build_system_wsmp.f90, DOUAR.f90, define_cloud.f90, define_ov.f90,
+      define_surface.f90, interpolate_ov_on_osolve.f90, isostasy.f90,
+      make_cut.f90, make_matrix.f90, module_definitions.f90, move_surface.f90,
+      write_global_output.f90
+
+r297-302:
+    Added ability of points to be removed when density becomes too high in a
+    given surface (surface level + 4)
+    - DOUAR.f90, refine_surface.f90
+
+r306,310-311,315 (make_cut.f90 rewrite),317 (pressure_cut.f90 rewrite),327,330:
+    Added ability to choose to use a volumetric majority or minority to define
+    elemental material properties, in addition to using divFEM
+    - make_cut.f90, module_definitions.f90, pressure_cut.f90, read_input_file.f90
+
+r331-333,336:
+    Revised vrm.f90 to have consistent formatting and 7 Drucker-Prager options.
+    Some are known to be bad, but available for comparison with older versions
+    of the code.
+    - vrm.f90
+--------------------------------------------------------------------------------
+
+Minor changes:
+--------------------------------------------------------------------------------
+r164:
+    e2d is now output as sqrt(J2d), not J2d (r164)
+    - do_leaf_measurements.f90
+
+r165-173,184-188,192,194,197,199,256,263:
+    Elemental pressure, smoothed pressure, viscosity, pre-isostasy velocities,
+    isostasy-only velocities, is_plastic and countnode now output
+    - define_ov.f90, DOUAR.f90, input_of_outputs.txt, isostasy.f90,
+      module_definitions.f90, post.f90, write_global_output.f90
+
+r189-190:
+    Gave post a Makefile like the rest of the code
+    - Makefile, VTK/Makefile
+
+r212:
+    Fixed output formatting to be more flexible
+    - post.f90, 
+
+r235,238-242,328:
+    Added surface type(s)
+    12: rectangular emboss with specified slope
+    13: Simple plateau geometry with variable width
+    - create_surfaces.f90
+
+r265:
+    Modified MPI usage to call mpi_ssend() rather than mpi_send()
+    - solve_with_pwgsmp.f90, solve_with_pwssmp.f90
+
+r270:
+    Implemented fix from Anshul to make the temperature calculation work
+    properly
+    - solve_with_pwgsmp.f90
+
+r305,309,314,321,323,325-327,329,441:
+    Added storage/output of elemental dilatation rate
+    - define_cloud.f90, define_ov.f90, do_leaf_measurements.f90, DoRuRe.f90,
+      DOUAR.f90, input_of_outputs.txt, module_definitions.f90, post.f90,
+      write_global_output.f90
+
+r307,311,316,318,326-327:
+    Added ability to use an initial strain rate to place plastic materials on
+    yield
+    - DOUAR.f90, make_matrix.f90, module_definitions.f90, read_input_file.f90,
+      vrm.f90
+
+r307, 311-312, 327:
+    Added option to specify an elemental bulk viscosity rather than penalty
+    value
+    - make_matrix.f90, make_pressure.f90, module_definitions.f90,
+      read_input_file.f90
+
+r319:
+    Revised nonlinear iteration convergence criterion to determine convergence
+    has been reached when the maximum vector difference over the maximum vector
+    length from the previous iteration is below the given tolerance value
+    - compute_convergence_criterion.f90
+
+r436:
+    Removed output of 'done visu.vtk' from post.f90 since it is not part of the
+    standard VTK file format
+    - post.f90
+
+r439, 441-443:
+    Added input option for choosing to use a 2D form of the strain rate tensor
+    invariants. Assumes plane strain in the x-z plane.
+    - do_leaf_measurements.f90, improve_osolve.f90, make_matrix.f90,
+      module_definitions.f90, module_invariants.f90, read_input_file.f90
+--------------------------------------------------------------------------------
\ No newline at end of file
--- a/src/EXAMPLE_v0.1
+++ b/src/EXAMPLE_v0.1
+--------------------------------------------------------------------------------
+DOUAR-WSMP v0.1 - 23/03/2011
+--------------------------------------------------------------------------------
+
+Running a DOUAR-WSMP example calculation
+--------------------------------------------------------------------------------
+To facilitate testing of DOUAR-WSMP, we have included a simple example input.txt
+and define_bc.f90 file (input_DOUAR_v0.1.txt and define_bc_DOUAR_v0.1.f90). You
+can use these files to run an example calculation that will ensure DOUAR is
+built properly and functional. To run the example, follow the steps below.
+
+1.  Copy the define_bc_DOUAR_v0.1.f90 file in the DOUAR-WSMP base directory to
+    'define_bc.f90'.
+2.  Build the 'douar_wsmp' executable following the directions provided in the
+    INSTALL file.
+3.  Create a new directory to store the DOUAR-WSMP executable, input file and
+    output (e.g., mkdir ~/DOUAR-WSMP-test).
+4.  Copy the DOUAR-WSMP executable (douar_wsmp) and input file
+    (input_DOUAR_v0.1.txt) to the new directory. Note that the input file should
+    be renamed 'input.txt'.
+5.  Create subdirectories needed to store the DOUAR-WSMP output in the new
+    directory:
+        OUT
+        DEBUG/mpilogs
+6.  Assuming you have a MPI-enabled system, you can run the test from the
+    command line. For example:
+        mpirun -np 8 ./douar_wsmp
+    or submit a job to the queue for systems with job management software. The
+    example calculation should less than 10 minutes on 8 modern processor cores.
+7.  This should produce output in the OUT subdirectory that can be converted to
+    VTK format using the 'post' post-processing program. VTK output can be
+    viewed in Mayavi, Paraview or other suitable viewers.
\ No newline at end of file
--- a/src/INSTALL_v0.1
+++ b/src/INSTALL_v0.1
+--------------------------------------------------------------------------------
+DOUAR-WSMP v0.1 - 23/03/2011
+--------------------------------------------------------------------------------
+
+Installation instructions
+--------------------------------------------------------------------------------
+There are several sets needed to build DOUAR-WSMP, the number of which varies
+depending on the compiler selected and libraries used.
+
+In general, the process is:
+1.  Select an appropriate Makefile.inc from the available options (or create
+    your own) and copy that file to Makefile.inc in the base DOUAR directory.
+2.  Edit the Makefile.inc to ensure the compiler and library locations are
+    correct.
+3.  Copy your desired boundary condition file to define_bc.f90 in the DOUAR base
+    directory.
+4.  Build DOUAR by running 'make all'. If successful, this will create a binary
+    called 'douar_wsmp'. If unsuccessful, consult the *.compile files in the
+    base directory and subdirectories, or the douar_wsmp.link.std* files.
+5.  Build the post-processing program by changing to the VTK directory and
+    running 'make all'. If successful, this will create a binary called 'post'.
+6.  The following environment variables should be set prior to running a
+    DOUAR-WSMP calculation:
+        MALLOC_TRIM_THRESHOLD_=100000000    (for WSMP)
+        MALLOC_MMAP_MAX_=0                  (for WSMP)
+        MKL_NUM_THREADS=1                   (if using the Intel MKL)
+        ulimit -s 131072                    (sets stack size, larger values may
+                                            resolve crashes in DOUAR-WSMP)
+--------------------------------------------------------------------------------
+
+Makefile.inc.* - Makefile include files
+                 Known working combinations of compilers and libraries
+                 (<compiler; MPI; additional libraries>)
+--------------------------------------------------------------------------------
+Makefile.inc.geodyncomp.intel.mpich2 (ifort/icc v12; MPICH2; WSMP, Intel MKL)
+Makefile.inc.geodyncomp.intel.openmpi (ifort/icc v12; Open MPI; WSMP, Intel MKL)
+Makefile.inc.gpc.intel.intelmpi (ifort/icc v11/12; IntelMPI; WSMP, Intel MKL)
+Makefile.inc.gpc.intel.mvapich2 (ifort/icc v11/12; MVAPICH2; WSMP, Intel MKL)
+Makefile.inc.gpc.intel.openmpi (ifort/icc v11; Open MPI; WSMP, Intel MKL)
+Makefile.inc.mahone.pgi.openmpi (PGI v8; Open MPI; WSMP, ScaLAPACK, BLACS)
+Makefile.inc.placentia.intel.mvapich2 (ifort/icc v12; MVAPICH2; WSMP, Intel MKL)
+Makefile.inc.placentia.intel.openmpi (ifort/icc v12; Open MPI; WSMP, Intel MKL)
+Makefile.inc.placentia.pgi.openmpi (PGI v8; Open MPI; WSMP, ScaLAPACK, BLACS)
+
+Obsolete Makefile.* from previous code versions (Available for reference)
+                 (<compiler; MPI; additional libraries>), when known
+--------------------------------------------------------------------------------
+Makefile.inc.mahone.orig (similar to Makefile.inc.mahone.pgi.openmpi)
+Makefile.inc.p690 (DOUAR-WSMP for IBM p690 running AIX)
+Makefile.cray (DOUAR-WSMP for Cray)
+Makefile.ifort (DOUAR-MUMPS for Linux; ifort/icc)
+Makefile.macosx (DOUAR-MUMPS for Mac; ifort/icc)
+Makefile.macosx.cedric (DOUAR-MUMPS for Mac)
+Makefile.pathscale (DOUAR-WSMP for Linux; Pathscale)
+--------------------------------------------------------------------------------
+
+Makefile - List of source files to compile
+--------------------------------------------------------------------------------
+Makefile (For use with DOUAR-WSMP v0.1)
+
+Obsolete Makefiles for reference
+--------------------------------------------------------------------------------
+Makefile_Cedric
+Makefile_DavidWhipp
+Makefile_Jean
+--------------------------------------------------------------------------------
\ No newline at end of file
--- a/src/README_v0.1
+++ b/src/README_v0.1
+--------------------------------------------------------------------------------
+DOUAR-WSMP v0.1 - 23/03/2011
+--------------------------------------------------------------------------------
+
+This is the README for DOUAR-WSMP version 0.1. This version of the code has been
+modified primarily by D. Whipp and D. Guptill (Dalhousie) from a version of the
+code handed down by J. Braun in early 2009. This version uses the WSMP solver,
+and has diverged from the original version of the code that used the MUMPS
+solver. Below, I list the changes that have been implemented in the code since
+we began tracking them in July of 2009. There are still a large number of
+features that need to be added, but hopefully this is helpful to new DOUAR
+users.
+--------------------------------------------------------------------------------
+
+- Changes are summarized in the ChangeLog file.
+- Detail on building the code can be found in the INSTALL file.
+- Instructions for running a simple test model can be found in the EXAMPLE file.
+--------------------------------------------------------------------------------
+
+Known bugs:
+--------------------------------------------------------------------------------
+1.  DOUAR-WSMP does not run reliably when linked with Open MPI. The cause of
+    this trouble is unclear.
+2.  The subdirectory DEBUG/mpilogs must exist when using debug level 1 or
+    higher. It seems this should only need to exist for level 2 or higher.
+--------------------------------------------------------------------------------
+
+Known limitations:
+--------------------------------------------------------------------------------
+1.  Constant temperature boundary conditions for the thermal calculation are
+    applied at the top set of nodes on the free surface. Since this surface is
+    always a cut surface, this will not represent the actual surface
+    temperature.
+2.  The largest successful simulation using the current code version was ~2.7
+    million elements on 128 cores with 8GB of RAM per core.
+3.  Some MPI implementations, including Open MPI do not support passing of
+    command line options to the code. This means that the use of the different
+    define_bc options is not functioning properly. This should be resolved in
+    version 0.2.
+4.  The code will crash if any surface is advected into the base of the model.
+    It is not clear how best to handle these surfaces.
+--------------------------------------------------------------------------------
+
+Near-term planned feature additions/modifications:
+--------------------------------------------------------------------------------
+1.  Boundary condition selection and input parameters need to be input values
+    from the input.txt file.
+2.  There is no need to refine the surfaces following isostatic adjustment, so
+    the code will be modified to turn that off.
+3.  Need to further modify the velocity boundary conditions for the HT models so
+    that they adjust more accurately to crustal subsidence. Rather than looking
+    at the isostatic subsidence array, the Moho surface position should be used
+    to modify the basal velocity B/Cs.
+4.  Strain softening in vrm.f90 should be a function so that it could be more
+    flexible and apply to either the friction angle or cohesion.
+5.  More output should be reported to the stdout to ensure proper conversion of
+    input values. This may already be done in some of the debugging files, but
+    that is not currently clear.
+6.  There are numerous unused variables that could be pruned from the code and
+    unintialized variables that need initial values.
+--------------------------------------------------------------------------------
+
+Long-term planned feature additions/modifications:
+--------------------------------------------------------------------------------
+1.  The solver should be selectable by the user so that either MUMPS or WSMP
+    could be used for people who don't have a WSMP license.
\ No newline at end of file
--- a/src/define_bc_DOUAR_v0.1.f90
+++ b/src/define_bc_DOUAR_v0.1.f90
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+!                                                                              |
+!              ||===\\                                                         | 
+!              ||    \\                                                        |
+!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
+!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
+!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
+!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
+!                                                                              |
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+!                                                                              |
+!              DEFINE_BC   Feb. 2009                                           |
+!                                                                              |
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+
+subroutine define_bc (params,osolve,vo,zi) 
+
+!------------------------------------------------------------------------------|
+!((((((((((((((((  Purpose of the routine  )))))))))))))))))))))))))))))))))))))
+!------------------------------------------------------------------------------|
+! this routine assigns the boundary condition for the Tibet experiment
+
+!------------------------------------------------------------------------------|
+!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
+!------------------------------------------------------------------------------|
+
+use definitions
+
+implicit none
+
+type (parameters) params
+type (octreesolve) osolve
+type (void) vo
+type (ziso) zi
+
+!------------------------------------------------------------------------------|
+!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
+!------------------------------------------------------------------------------|
+
+integer i,iproc,nproc,ierr
+double precision eps,lsf0,pi,lorig,h,x1,x2,phi,spu,spv,yend,cper,cscl,xstart
+double precision ystart,theta,l,vin,vzfluxscl,cntvel,dxy,xend,xsym,ymax,xwidth
+double precision ywidth
+double precision,dimension(:),allocatable :: x0,ldisp
+integer ie,ij,j,jp,nelemx,nelemz
+logical firstcall
+
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+
+INCLUDE 'mpif.h'
+
+call mpi_comm_size (mpi_comm_world,nproc,ierr)
+call mpi_comm_rank (mpi_comm_world,iproc,ierr)
+
+eps=1.d-10
+osolve%kfix=0
+osolve%kfixt=0
+pi=atan(1.d0)*4.d0
+vin=1.d0                                                                        ! Reference velocity
+
+select case(trim(params%infile))
+
+case('input.txt','input.small.txt')
+   do i=1,osolve%nnode
+     if (osolve%x(i).lt.eps) then
+       osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=1.d0
+       osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+       osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
+     endif
+     if (osolve%x(i).gt.1.d0-eps) then
+       osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
+!       osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+!       osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
+     endif
+     if (osolve%y(i).lt.eps) then
+!       osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
+       osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+!       osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
+     endif
+     if (osolve%y(i).gt.1.d0-eps) then
+!       osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
+       osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+!       osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
+     endif
+     if (osolve%z(i).lt.eps) then
+       osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
+       osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+       osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
+       if (osolve%x(i).le.0.5d0) osolve%u(i)=1.d0
+       osolve%kfixt(i)=1
+       osolve%temp(i)=1.d0
+     endif
+     if (osolve%z(i).gt.1.d0-eps) then
+       osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
+       osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+       osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
+       osolve%kfixt(i)=1
+       osolve%temp(i)=0.d0
+     endif
+     if (.not.vo%influid(i)) then
+       osolve%kfixt(i)=1
+       osolve%temp(i)=0.d0
+     endif
+     osolve%u(i)=osolve%u(i)*vin
+     osolve%v(i)=osolve%v(i)*vin
+     osolve%w(i)=osolve%w(i)*vin
+   enddo
+
+case ('input.jgr')
+   call define_bc_jgr (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 
+
+case ('input.sphere')
+   call define_bc_sphere (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 
+
+case ('input.3Dpunch')
+   call define_bc_3Dpunch (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 
+
+case ('input.2Dpunch')
+   call define_bc_2Dpunch (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 
+
+case ('input.2Dpunch_vert')
+   call define_bc_2Dpunch_vert (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 
+
+case ('input.folding')
+   call define_bc_folding (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 
+
+case ('input.pipo')
+   call define_bc_pipo (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 
+
+case ('input.ritske')
+   call define_bc_ritske (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 
+
+case ('input.ritske_isurf')
+   call define_bc_ritske_isurf (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 
+
+case ('input.subduction')
+   call define_bc_subduction (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 
+
+case ('input.riedel')
+   call define_bc_riedel (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 
+
+case ('input.parallipipede')
+   call define_bc_parallipipede (osolve%nnode,osolve%kfix,osolve%kfixt,osolve%x,osolve%y,osolve%z,osolve%u,osolve%v,osolve%w,osolve%temp,vo) 
+
+case default
+   if (iproc.eq.0) print *,params%infile
+   call stop_run (' pb with input file in define_bc.f90')
+
+end select
+
+! here we need to impose that the bc satisfy the bad faces constrains
+
+do ie=1,osolve%nface
+  do j=1,4
+  jp=1+mod(j,4)
+  i=osolve%iface(j+4,ie)
+  if (osolve%kfix((i-1)*3+1).eq.1) osolve%u(i)=(osolve%u(osolve%iface(j,ie))+osolve%u(osolve%iface(jp,ie)))/2.d0
+  if (osolve%kfix((i-1)*3+2).eq.1) osolve%v(i)=(osolve%v(osolve%iface(j,ie))+osolve%v(osolve%iface(jp,ie)))/2.d0
+  if (osolve%kfix((i-1)*3+3).eq.1) osolve%w(i)=(osolve%w(osolve%iface(j,ie))+osolve%w(osolve%iface(jp,ie)))/2.d0
+  if (osolve%kfixt(i).eq.1) osolve%temp(i)=(osolve%temp(osolve%iface(j,ie))+osolve%temp(osolve%iface(jp,ie)))/2.d0
+  enddo
+i=osolve%iface(9,ie)
+if (osolve%kfix((i-1)*3+1).eq.1) osolve%u(i)=sum(osolve%u(osolve%iface(1:4,ie)))/4.d0
+if (osolve%kfix((i-1)*3+2).eq.1) osolve%v(i)=sum(osolve%v(osolve%iface(1:4,ie)))/4.d0
+if (osolve%kfix((i-1)*3+3).eq.1) osolve%w(i)=sum(osolve%w(osolve%iface(1:4,ie)))/4.d0
+if (osolve%kfixt(i).eq.1) osolve%temp(i)=sum(osolve%temp(osolve%iface(1:4,ie)))/4.d0
+enddo
+
+if (params%debug==2) call output_bc (osolve)
+
+return
+end
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|