Interpolate Solution from Different Mesh

This will be a general page for how to interpolate solutions from different meshes onto a new mesh. Those meshes are assumed to be of the same domain.

The generic terms for the two meshes are "source" and "target", where source has the desired solution data and target is the mesh that will be receiving.

Laminar, Incompressible, Semi-structured Mesh

This section will assume that the source mesh is in a structured ijk form. In the future, this may be expanded to meshes created by MGEN.

Process Overview

1. A CSV of the solution must be created. This must be done on a serial running version of ParaView and must be done after a MergeBlocks filter is done.
2. A special solution file is created via Sort2StructuredGrid, which will take in the csv from step 1 and the ordered coordinates of the source file from Matlab.
   • Effectively, it creates a single solution file in the same ordering as the Matlab points. The importance of this is significant as the executable used in the next step depends on some expected ordering. If the ordering is different, a new executable will need to be used/created.
3. The interpolation is performed onto the new grid via par3DInterp3, which creates solInterp.<1-nparts> files in the solnTarget directory.
4. The interpolation is then used by PHASTA by setting Load and set 3D IC: True in the solver.inp.
   • This should only be done for 1 timestep, as it will continue to reset the IC for all proceeding timesteps. The solnTarget directory needs to be symlinked into the -procs_case directory for this to work.

1. Create CSV

• Created using ParaView
• PV must be running in Serial mode
  • Otherwise the CSV will not be in the correct order and possibly have duplicated points

1. Load in source dataset
2. Apply Mergeblocks filter
3. Save dataset as a csv with 12 digits of precision in scientific notation
   • Make sure the csv is in "pressure, u0, u1, u2, x0, x1, x2" format
   • This can be done by only loading the pressure and velocity fields into Paraview (either by editing the .phts or in the data load menu in Paraview).
4. Replace the commas with spaces
   • Can use vim or run sed -i 's/,/\ /g' test.csv
   • Though the next step looks for a .csv extension, it is a fortran formatted read and actually needs those commas replaced by spaces
5. Remove the first line of the csv file
   • Done in vi or sed (sed -i 1,1d test.csv) or tail (tail -n +2 test.csv > trimmedLine1.csv)
   • Needed for the next program
   • Better yet, we should change the next code to read past that header line and then delete this line when that is complete. We should also consider the solution in this StackOverflow answer as it shows how to make a data structure that could read the csv lines directly in the next program and avoid ALL this file manipulation with modern fortran (see HighPerformanceMark's answer).

2. Create Structured Solution File

Note: These instructions will be for the parallelSortDNSzBinJames executable, which has some highly specific requirements and command inputs.

This step will take the data from the source solution file and put it in an format/order that will make the interpolation process work much faster.

1. Symlink the source mesh's ordered coordinate file as ordered.crd
   • This may come from the files used to create the mesh (ie. for matchedNodeElementReader)
   • (Untested)This may also be created using the coordinates from the solution file
2. Rename/symlink csv to be the correct file name (in my specific case, it was dnsSolution1procLongFort.csv)
3. Create an interactive job on whatever machine you're needing to run on (ALCF Cooley in this case)
4. Load approriate MPI environment variables (soft add +mvapich2 for Cooley)
5. Run the executable via mpirun -np [nprocs] [executable path] [executable inputs]
   • This will produce source.sln.{1..nprocs} files
6. Concatenate source.sln.{1..nprocs} files in order into single source.sln file
   • The individual source.sln.{1..nprocs} files need to be concatenated into a single solution.sln file, which can be done (in zsh at least) via echo source.sln.{1..[MPIRanks]}| xargs cat > source.sln (or probably equivalent cat source.sln.{1..[MPIRanks]} > source.sln). Note these files must be concatenated in order of rank, otherwise it will be out of sequence.

Example Command: mpirun -np 24 /lus/theta-fs0/projects/PHASTA_aesp/Utilities/Sort2StructuredGrid/parallelSortDNSzBinJames 47822547 47822547 212 0.0291

• The inputs for this command are [nlines of csv] [nlines of ordered.crd] [number of element in z] [z domain width]
• Note: These inputs are specfic to this executable. Changing executable will change which is used
• Also note that the [number of elements in z] is equivalent to nsons - 1 or the number of nodes in z - 1.

3. Create Interpolated files

This step will create the solInterp.[target nprocs] files used by PHASTA to perform the interpolation.

1. Create Interpolate.../[target nprocs]-procs_case directory in the target's Chef directory and move to that directory
2. Symlink the target's POSIX geombc.[target nprocs] files (that were created by Chef) to the [target nprocs]-procs_case directory
   • The geombc.[target nprocs] files should be copied in the exact fashion that they are in the Chef created [target nprocs]-procs_case directory, including if they're "fanned out"
3. Create a directory called solnTarget
   • This may be corrected in the future, but currently if solnTarget is not present the job will fail
4. Symlink the source.sln to the directory and the ordered.crd file as source.crd
5. Run phInterp via mpirun on an interactive job.
6. This creates a series of solInterp.[target nprocs] files in the solnTarget directory

The file format for solInterp.N is quite simple. Each line corresponds to the node number in the partition and the file itself has 7 columns:

coord_x coord_y coord_z pressure velocity_x velocity_y velocity_z

Example Command: mpirun -np 64 /path/to/phInterp 16 799 281 213 0.452

• The inputs for this command are [target parts per MPIProc] [source nx] [source ny] [source nz] [z Length]
• Note that the number of processes given to mpirun times the [target parts per MPIProc] must be equal to the number of target partitions. In this case, the target partition was 1024, so 64*16 = 1024
• The way it works is that each of the 64 MPIProcs is given 16 partitions to interpolate.

4. Interpolate the solution in PHASTA

This step will take the solInterp.[target nprocs] files and load them as initial conditions.

1. Symlink the solnTarget directory into the [target nprocs]-procs_case directory.
2. Add/uncomment Load and set 3D IC: True in the solver.inp
3. Run PHASTA for a few timesteps and write out restart-dat.[target nprocs] files
4. Remove/comment out the Load and set 3D IC: True line from the solver.inp
5. The new restart files have the interpolated solution
   • Note that if you forget to remove the Load and set 3D IC: True statement from solver.inp, PHASTA will overwrite the existing solution in the restart files

Turbulent, Compressible, Unstructured Mesh

Currently, the only version of PHASTA that is set up to handle this type of solution transfer is the conrad54418/ViscousHypersonics branch (as of 6/22/22). Creating the solInterp.1 file can be automated via a script or done manually (see below).

Process Overview (scripted)

Scripts have been made for use with the compressible, but not turbulent, version of the code. An example folder with the scripts included can be found at /project/tutorials/ParaviewSolutionTransfer. The three scripts needed are called interpolateSol.py, pvCSV2customSLN_Nproc_prim.m, and parRunAll.sh. The only script of those three you need to run is parRunAll.sh, which can be found in targetFolder/solutionInterp. You then need to run PHASTA via the runPhasta.sh script, which will produce a restart.1.1 file that contains the transferred solution on the new mesh. The manual section below provides insight about what the scripts are doing.

Process Overview (manual)

1. Load existing solution into ParaView. Use the 'merge blocks' filter to convert to serial case.
2. Load the target case .pht file into ParaView
3. Use the 'Resample from dataset' filter and select the source and target blocks accordingly (ParaView’s naming is super unclear: they want source to be the coordinates where you need the solution (new mesh) and input to be the mesh that has the solution values associated with it that you want to interpolate from (which in this case is the MergeBlocks). This is confusing because it is exactly backwards to what we use for solution interpolation where we call the mesh with a solution the source and the new mesh the target).
4. Save the output as a .csv file. Write a single time step and select scientific notation to 12 decimals. NOTE: sometimes ParaView will make an error and write zeros where it can't quite find the closest point when doing the solution transfer. In this case, the .csv file needs to be manually edited to replace the zeros for pressure and temperature with realistic values.
5. The .csv file can now be reformatted and renamed with MATLAB to match the expected form of solInterp.1.
6. Advance PHASTA one step in serial, then convert to desired number of processors using Chef.

Alternative to Using MATLAB for Reformatting

For those wanting to skip the MATLAB step, awk can also do the necessary column manipulation.

Copy your file (in case you goof up):

cp PVinterp0.csv test.dat

Remove the header line:

sed -i 1,1d test.dat

Replace the commas with spaces:

sed -i 's/,/\ /g' test.dat

Rearrange the columns to be what solInterp.1 wants. It is probably a good idea to do a head -1 PVinterp0.csv to be sure as you might have more or less fields than I did but use that header to find the column numbers (starting from 1, not 0) to write x, y, z, p, u, v, w.

awk '{print $6,$7,$8,$1,$2,$3,$4}' test.dat > solInterp.1

and don’t forget to put this into a directory called solnTarget and also to turn on the flag in solver.inp Load and set 3D IC: True for the 1 step that Joe mentioned. Finally, if you are worried about that one step messing up your solution recent versions of the code can take

     iexec : 0

or

    Number of Timesteps: 0

to not take any actual steps. Note however that only very recent version of the code have the iexec conditional moved AFTER the loading of the interpolated solution but it should be pretty easy to figure out where to move that conditional. Alternatively, the second option also skips over the time stepping and writes the solution AFTER applying the boundary conditions which can be useful to see as well to confirm you have the intended BC’s set (iexec :0 won’t detect this)