zasexton
commented
9 months ago The following function is how I propose to read in time-varying solution field data across a given vtu mesh for svFSIplus. This would replicate behavior of the subroutine READVELOCITYVTU called in svFSI. However, there are some notable differences.
- the
field_namesuch as"Velocity", "Pressure", etc must be supplied by the user and will not be hard-coded into the function. - the time step of the field data is assumed to come at the end of the
vtkDataArrayname. This conforms to how such time-varying velocity and pressure data are written out to vtk files fromsvPost. However, because we do not hard-code the array name format we needed to define a robust method for extracting the time step information. -
Last, this function assumes that the
mshTypeobjects are equipped with an additionalArray3<double>(calledYs) attribute capable of storing the solution field of interest.void load_time_varying_field_vtu(const std::string file_name, const std::string field_name, mshType& mesh) { #define n_debug_load_vtu #ifdef debug_load_vtu std::cout << "[load_vtu] " << std::endl; std::cout << "[load_vtu] ===== vtk_xml_parser::load_time_varying_field_vtu ===== " << std::endl; std::cout << "[load_vtu] file_name: " << file_name << std::endl; #endif auto reader = vtkSmartPointer<vtkXMLUnstructuredGridReader>::New(); reader->SetFileName(file_name.c_str()); reader->Update(); vtkSmartPointer<vtkUnstructuredGrid> vtk_ugrid = reader->GetOutput(); vtkIdType num_nodes = vtk_ugrid->GetNumberOfPoints(); int array_count = 0; std::vector<str::pair<std::string, int>> array_names; if (num_nodes == 0) { throw std::runtime_error("Failed reading the VTK file '" + file_name + "'."); } // Store all array names for (int i = 0; i < vtk_ugrid->GetPointData()->GetNumberOfArrays(); i++) { std::string array_name = vtk_ugrid->GetPointData()->GetArrayName(i); size_t pos = array_name.find(field_name.c_str()); if (pos != std::string::npos) { auto not_digit = [](char c) { return !std::isdigit(c); }; auto it = std::find_if(array_name.rbegin(), array_name.rend(), not_digit); std::string time_step = std::string(it.base(), array_name.end()); array_count++; if (!time_step.empty()) { array_names.push_back({array_name, std::stoi(time_step)}); } else { array_names.push_back({array_name, 0}); } } } if (array_count == 0) { throw std::runtime_error("No '" + field_name + "' data found in the VTK file '" + file_name + "'."); } // Order all array names std::sort(array_names.begin(), array_names.end(), [](const std::pair<std::string, int>& a, const std::pair<std::string, int>& b) { return a.second < b.second; }); int num_components = array->GetNumberOfComponents(); mesh.Ys.resize(num_components, num_nodes, array_count); for (int i = 0; i < array_count; i++) { auto array = vtk_ugrid->GetPointData()->GetArray(array_names[i].first.c_str()); if (array == nullptr) { throw std::runtime_error("No '" + array_names[i].first + "' data found in the VTK file '" + file_name + "'."); } for (int j = 0; j < num_nodes; j++) { for (int k = 0; k < num_components; k++) { mesh.Ys(k, j, i) = array->GetComponent(j, k); } } } }
ktbolt
Problem
The general problem: We wish to use the (potentially time-varying) velocity solutions obtained from a previous simulation within a subsequent physics (namely, scalar advection-diffusion).
Solution
This has been done previously in
svFSIon side-branches where residence time calculations have been desired or where thrombosis modeling has been done. See Gutierrez et al 2021.Although it is not clear if this exact framework would work well for a more object-oriented
svFSIplusthe following scheme is generally implemented insvFSI:True/Falseflag is selected forUse pre-computed velocityIf
Truethe following subroutine is called insvFSIto read in the velocity solution vector fields across some set of time points. This subroutine can be found in theVTKXML.ffile of the attachedsvFSIsource code. Reminder this will not be found on the main branches ofsvFSIas this functionality was never formally included.DISTRIBUTEsubroutine inDISTRIBUTE.f! Distributing tagFile node info to processors flag = (ALLOCATED(tagRT)) CALL cm%bcast(flag) IF (flag) THEN
IF (cm%mas()) THEN ALLOCATE(tmpRT(gtnNo)) tmpRT = tagRT DEALLOCATE(tagRT) ELSE ALLOCATE(tmpRT(0)) END IF ALLOCATE(tagRT(tnNo)) tagRT = LOCAL(tmpRT) DEALLOCATE(tmpRT) END IF
Attached files include the source code and example input files for residence time calculations. svFSI_input.txt tagFile.txt Code.tar.gz
Additional context
No response
Code of Conduct