diff --git a/partitioned-heat-conduction/README.md b/partitioned-heat-conduction/README.md index 5afb3b4df..fa4195e86 100644 --- a/partitioned-heat-conduction/README.md +++ b/partitioned-heat-conduction/README.md @@ -1,7 +1,7 @@ --- title: Partitioned heat conduction permalink: tutorials-partitioned-heat-conduction.html -keywords: FEniCS, Nutils, Heat conduction +keywords: FEniCS, Nutils, G+Smo, Heat conduction summary: We solve a simple heat equation. The domain is partitioned and the coupling is established in a Dirichlet-Neumann fashion. --- @@ -37,6 +37,8 @@ You can either couple a solver with itself or different solvers with each other. * OpenFOAM. This case uses the custom [heatTransfer](https://github.com/precice/tutorials/blob/master/partitioned-heat-conduction/solver-openfoam/heatTransfer.C) solver (find it in `solver-openfoam` and build with `wmake`). Read more details in the [OpenFOAM adapter](https://precice.org/adapter-openfoam-overview.html). +* G+Smo. Inatsll [G+Smo](https://github.com/gismo/gismo). + ## Running the simulation You can find the corresponding `run.sh` script for running the case in the folders corresponding to the participant you want to use: @@ -73,6 +75,38 @@ If you want to use Nutils or OpenFOAM, use `cd dirichlet/neumann-nutils`, respec mpirun -n heat.py -d ``` +The G+Smo-based version of the tutorial offers IsoGeometric Analysis discretization method. To run the example you need to follow the following steps: + +- Download G+Smo and Create a Build Folder +``` +git clone git@github.com:gismo/gismo.git +cd gismo +mkdir build +``` + - Configure G+Smo +``` +cmake .. -DGISMO_OPTIONAL=";gsPreCICE" +``` +- Build the Example +``` +make partitioned-heat-conduction -j +``` +- Link the compiled executable to the gismo-executable folder within the tutorial directory +``` +cd /partitioned-heat-conduction/gismo-executable +ln -sf /bin/partitioned-heat-conduction ./gismo_executable` +``` +- Open two terminals and run +``` + cd dirichlet-gismo + ./run.sh +``` + +``` +cd neumann-gismo +./run.sh +``` + ### Note on the combination of Nutils & FEniCS You can mix the Nutils and FEniCS solver, if you like. Note that the error for a pure FEniCS simulation is lower than for a mixed one. We did not yet study the origin of this error, but assume that this is due to the fact that Nutils uses Gauss points as coupling mesh and therefore entails extrapolation in the data mapping at the top and bottom corners. @@ -85,6 +119,8 @@ For FEniCS you can visualize the content with paraview by opening the `*.pvd` fi For Nutils, please use the files `Dirichlet-*.vtk` or `Neumann-*.vtk`. Please note that these files contain the temperature as well as the reference solution. +For G+Smo, please use the file `solution.pvd` in both dirichlet-gismo and neumann-gismo directories. + ![Animation of the partitioned heat equation](images/tutorials-partitioned-heat-conduction-FEniCS-movie.gif) Visualization in paraview for `x_c = 1.5`. diff --git a/partitioned-heat-conduction/dirichlet-gismo/clean.sh b/partitioned-heat-conduction/dirichlet-gismo/clean.sh new file mode 100755 index 000000000..ad9a90c93 --- /dev/null +++ b/partitioned-heat-conduction/dirichlet-gismo/clean.sh @@ -0,0 +1,27 @@ +#!/bin/bash + +# Cleaning script for partitioned-heat-conduction directory + +echo "Cleaning unnecessary files..." + +# Remove precice-profiling directory +if [ -d "precice-profiling" ]; then + rm -rf precice-profiling + echo "Deleted 'precice-profiling' folder." +fi + +# Remove precice-run directory +if [ -d "../precice-run" ]; then + rm -rf ../precice-run + echo "Deleted 'precice-run' folder." +fi + + +# Remove files ending with .pvd, .vts, .vtp, .log, and .txt +for ext in pvd vts vtp log txt; do + find . -type f -name "*.$ext" -exec rm -f {} \; + echo "Deleted all *.$ext files." +done + +echo "Cleaning completed!" + diff --git a/partitioned-heat-conduction/dirichlet-gismo/run.sh b/partitioned-heat-conduction/dirichlet-gismo/run.sh new file mode 100755 index 000000000..2612a8cd0 --- /dev/null +++ b/partitioned-heat-conduction/dirichlet-gismo/run.sh @@ -0,0 +1,4 @@ +#!/bin/bash +set -e -u + +../gismo-executable/gismo-executable -c ../precice-config.xml --plot -s 0 diff --git a/partitioned-heat-conduction/gismo-executable/gismo-executable b/partitioned-heat-conduction/gismo-executable/gismo-executable new file mode 120000 index 000000000..e63cd46fd --- /dev/null +++ b/partitioned-heat-conduction/gismo-executable/gismo-executable @@ -0,0 +1 @@ +/Users/jli51/MyCodes/gismo_fsi/build/bin/partitioned-heat-conduction \ No newline at end of file diff --git a/partitioned-heat-conduction/gismo-executable/partitioned-heat-conduction.cpp b/partitioned-heat-conduction/gismo-executable/partitioned-heat-conduction.cpp new file mode 100644 index 000000000..54a222c32 --- /dev/null +++ b/partitioned-heat-conduction/gismo-executable/partitioned-heat-conduction.cpp @@ -0,0 +1,382 @@ +/** @file heat-equation-coupling.cpp + + @brief Heat equation participant for a double coupled heat equation + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): H.M. Verhelst (2019-2024, TUDelft), J. Li (2023-..., TUDelft) +*/ + +#include +#include +#include + +using namespace gismo; + +int main(int argc, char *argv[]) +{ + //! [Parse command line] + bool plot = false; + index_t plotmod = 1; + index_t numRefine = 2; // Number of uniform refinement (increase the number also requires to modify the rbf radius in the precice_config.xml) + index_t numElevate = 0; + short_t side = 0; + std::string precice_config("../precice_config.xml"); + + gsCmdLine cmd("Coupled heat equation using PreCICE."); + cmd.addInt( "e", "degreeElevation", + "Number of degree elevation steps to perform before solving (0: equalize degree in all directions)", numElevate ); + cmd.addInt( "r", "uniformRefine", "Number of Uniform h-refinement loops", numRefine ); + cmd.addString( "c", "config", "PreCICE config file", precice_config ); + cmd.addSwitch("plot", "Create a ParaView visualization file with the solution", plot); + cmd.addInt("m","plotmod", "Modulo for plotting, i.e. if plotmod==1, plots every timestep", plotmod); + cmd.addInt("s","side", "Patchside of interface: 0 (Dirichlet), 1 (Neumann)", side); + try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } + + //! [Read input file] + + gsMultiPatch<> patches; + if (side==0) //left + patches.addPatch(gsNurbsCreator<>::BSplineRectangle(0.0,0.0,1.0,1.0)); + else if (side==1) //right + patches.addPatch(gsNurbsCreator<>::BSplineRectangle(1.0,0.0,2.0,1.0)); + else + GISMO_ERROR("Side unknown"); + + real_t alpha = 3; + real_t beta = 1.3; + real_t time = 0; + real_t k_temp = 1; + + // Set external heat-flux to zero + gsConstantFunction<> f(beta-2-2*alpha,2); + gsFunctionExpr<> u_ex("1+x^2+" + std::to_string(alpha) + "*y^2 + " + std::to_string(beta) + "*" + std::to_string(time),2); + + gsMultiBasis<> bases(patches);//true: poly-splines (not NURBS) + + bases.setDegree( bases.maxCwiseDegree() + numElevate); + + // h-refine each basis + for (int r =0; r < numRefine; ++r) + bases.uniformRefine(); + numRefine = 0; + + gsInfo << "Patches: "<< patches.nPatches() <<", degree: "<< bases.minCwiseDegree() <<"\n"; + +// ---------------------------------------------------------------------------------------------- + typedef gsExprAssembler<>::geometryMap geometryMap; + typedef gsExprAssembler<>::space space; + typedef gsExprAssembler<>::solution solution; + + gsExprAssembler<> A(1,1); + + gsInfo<<"Active options:\n"<< A.options() <<"\n"; + + A.setIntegrationElements(bases); + + gsExprEvaluator<> ev(A); + +// ---------------------------------------------------------------------------------------------- + boxSide couplingSide; + if (side==0) //left + couplingSide = boxSide(2); + else if (side==1) //right + couplingSide = boxSide(1); + else + GISMO_ERROR("Side unknown"); + + patchSide couplingInterface(0,couplingSide); + typename gsBasis::domainIter domIt = bases.basis(couplingInterface.patch).makeDomainIterator(couplingInterface.side()); + index_t rows = patches.targetDim(); + gsMatrix<> nodes; + // Start iteration over elements + gsVector<> tmp; + index_t k=0; + + gsOptionList quadOptions = A.options(); + // NEED A DIFFERENT RULE FOR dirichlet::interpolation --> see: gsGaussRule bdQuRule(basis, 1.0, 1, iter->side().direction()); + /* + quadOptions.addInt("quRule","Quadrature rule [1:GaussLegendre,2:GaussLobatto]",1); + quadOptions.addReal("quA","Number of quadrature points: quA*deg + quB",1.0); + quadOptions.addInt("quB","Number of quadrature points: quA*deg + quB",1); + */ + + index_t quadSize = 0; + typename gsQuadRule::uPtr QuRule; // Quadrature rule ---->OUT + for (; domIt->good(); domIt->next(), k++ ) + { + QuRule = gsQuadrature::getPtr(bases.basis(couplingInterface.patch), quadOptions,couplingInterface.side().direction()); + quadSize+=QuRule->numNodes(); + } + gsMatrix<> uv(rows,quadSize); // Coordinates of the quadrature points in parameter space + gsMatrix<> xy(rows,quadSize); // Coordinates of the quadrature points in physical space + + index_t offset = 0; + + for (domIt->reset(); domIt->good(); domIt->next(), k++ ) + { + QuRule = gsQuadrature::getPtr(bases.basis(couplingInterface.patch), quadOptions,couplingInterface.side().direction()); + // Map the Quadrature rule to the element + QuRule->mapTo( domIt->lowerCorner(), domIt->upperCorner(), + nodes, tmp); + uv.block(0,offset,rows,QuRule->numNodes()) = nodes; + + gsMatrix<> tmp2; + patches.patch(couplingInterface.patch).eval_into(nodes,tmp2); + xy.block(0,offset,rows,QuRule->numNodes()) = patches.patch(couplingInterface.patch).eval(nodes); + offset += QuRule->numNodes(); + } + + std::string membername; + if (side==0) //left + membername = "Dirichlet"; + else if (side==1) //right + membername = "Neumann"; + else + GISMO_ERROR("Side unknown"); + + gsPreCICE interface(membername, precice_config); + std::string meshName = membername + "-Mesh"; + interface.addMesh(meshName,xy); + + interface.requiresInitialData(); + + real_t precice_dt = interface.initialize(); + + std::string tempName = "Temperature"; + std::string fluxName = "Heat-Flux"; + +// ---------------------------------------------------------------------------------------------- + + gsBoundaryConditions<> bcInfo; + gsPreCICEFunction g_CD(&interface,meshName,(side==0 ? tempName : fluxName),patches,1); + gsPreCICEFunction g_CN(&interface,meshName,(side==0 ? tempName : fluxName),patches,1); + gsFunction<> * g_C = &u_ex; + + bcInfo.addCondition(0, boundary::south, condition_type::dirichlet, &u_ex, 0, false, 0); + bcInfo.addCondition(0, boundary::north, condition_type::dirichlet, &u_ex, 0, false, 0); + if (side==0) + { + bcInfo.addCondition(0, couplingSide,condition_type::dirichlet , g_C, 0, false, 0); + bcInfo.addCondition(0, couplingSide.opposite(),condition_type::dirichlet, &u_ex, 0, false, 0); + } + else + { + bcInfo.addCondition(0, couplingSide,condition_type::neumann , &g_CN); + bcInfo.addCondition(0, couplingSide.opposite(),condition_type::dirichlet, &u_ex, 0, false, 0); + } + + bcInfo.setGeoMap(patches); + // gsDebugVar(bcInfo); + +// ---------------------------------------------------------------------------------------------- + + // Generate system matrix and load vector + // gsInfo << "Assembling mass and stiffness...\n"; + + // Set the geometry map + geometryMap G = A.getMap(patches); + + // Set the discretization space + space u = A.getSpace(bases); + + // Set the source term + auto ff = A.getCoeff(f, G); + + // Set the solution + gsMatrix<> solVector, solVector_ini; + solution u_sol = A.getSolution(u, solVector); + solution u_ini = A.getSolution(u, solVector); + + u.setup(bcInfo, dirichlet::homogeneous, 0); + A.initSystem(); + gsDebugVar(A.numDofs()); + A.assemble( u * u.tr() * meas(G)); + gsSparseMatrix<> M = A.matrix(); + + + // A Conjugate Gradient linear solver with a diagonal (Jacobi) preconditionner + gsSparseSolver<>::CGDiagonal solver; + + real_t dt = 0.01; + + // Project u_wall as initial condition (violates Dirichlet side on precice interface) + auto uex = A.getCoeff(u_ex, G); + // RHS of the projection + u.setup(bcInfo, dirichlet::l2Projection, 0); + A.initSystem(); + A.assemble( u * u.tr() * meas(G), u * uex * meas(G) ); + solver.compute(A.matrix()); + solVector_ini = solVector = solver.solve(A.rhs()); + + gsMatrix<> result(1,uv.cols()), tmp2; + // Write initial data + if (interface.requiresWritingCheckpoint()) + { + for (index_t k=0; k!=uv.cols(); k++) + { + if (side==0) + { + gsWarn<<"Write the flux here!!!\n"; + tmp2 = ev.eval( - jac(u_sol) * nv(G).normalized(),uv.col(k)); + result(0,k) = tmp2.at(0); + } + else + { + tmp2 = ev.eval(u_sol,uv.col(k)); + result(0,k) = tmp2.at(0); + } + } + gsDebugVar(result); + interface.writeData(meshName,side==0 ? fluxName : tempName,xy,result); + } + // interface.initialize_data(); + + // interface.readBlockScalarData(meshName,side==0 ? tempName : fluxName,xy,result); + // gsDebugVar(result); + + // Initialize the RHS for assembly + if (side==0) + g_C = &g_CD; + gsDebugVar(bcInfo); + A.initSystem(); + A.assemble( k_temp * igrad(u, G) * igrad(u, G).tr() * meas(G), u * uex * meas(G) ); + gsSparseMatrix<> K = A.matrix(); + + // Assemble the RHS + gsVector<> F = dt*A.rhs() + M*solVector; + gsVector<> F0 = F; + gsVector<> F_checkpoint = F; + gsMatrix<> solVector_checkpoint = solVector; + + gsParaviewCollection collection("solution"); + gsParaviewCollection exact_collection("exact_solution"); + + index_t timestep = 0; + index_t timestep_checkpoint = 0; + real_t t_checkpoint = 0; + if (plot) + { + std::string fileName = "solution_" + util::to_string(timestep); + ev.options().setSwitch("plot.elements", true); + ev.options().setInt("plot.npts", 1000); + ev.writeParaview( u_sol , G, fileName); + for (size_t p=0; p!=patches.nPatches(); p++) + { + fileName = "solution_" + util::to_string(timestep) + std::to_string(p); + collection.addTimestep(fileName,time,".vts"); + } + + // fileName = "exact_solution_" + util::to_string(timestep); + // ev.writeParaview( uex , G, fileName); + // for (size_t p=0; p!=patches.nPatches(); p++) + // { + // fileName = "exact_solution_" + util::to_string(timestep) + std::to_string(p); + // exact_collection.addTimestep(fileName,time,".vts"); + // } + + ev.writeParaview( u_sol , G, "initial_solution"); + + } + + // time += precice_dt; + while (interface.isCouplingOngoing()) + { + // read temperature from interface + u_ex = gsFunctionExpr<>("1+x^2+" + std::to_string(alpha) + "*y^2 + " + std::to_string(beta) + "*" + std::to_string(time),2); + u.setup(bcInfo, dirichlet::l2Projection, 0); // NOTE: dirichlet::l2Projection is used to project the Dirichlet BCs + + A.initSystem(); + A.assemble( k_temp * igrad(u, G) * igrad(u, G).tr() * meas(G), u * ff * meas(G) ); + K = A.matrix(); + auto g_Neumann = A.getBdrFunction(G); + A.assembleBdr(bcInfo.get("Neumann"), u * g_Neumann.val() * nv(G).norm() ); + F = dt*A.rhs() + M*solVector; + + // save checkpoint + if (interface.requiresWritingCheckpoint()) + { + gsDebugVar("Write checkpoint"); + F_checkpoint = F0; + t_checkpoint = time; + timestep_checkpoint = timestep; + solVector_checkpoint = solVector; + } + + // potentially adjust non-matching timestep sizes + dt = std::min(dt,precice_dt); + + // solve gismo timestep + gsInfo << "Solving timestep " << timestep*dt << "..."; + solVector = solver.compute(M + dt*K).solve(F); + gsInfo<<"Finished\n"; + // write heat fluxes to interface + gsMatrix<> result(1,uv.cols()), tmp; + for (index_t k=0; k!=uv.cols(); k++) + { + // gsDebugVar(ev.eval(nv(G),uv.col(k))); + // tmp = ev.eval(k_temp * igrad(u_sol,G),uv.col(k)); + // Only exchange y component + // result(0,k) = -tmp.at(1); + // + // + if (side==0) + { + gsWarn<<"Write the flux here!!!\n"; + tmp = ev.eval(jac(u_sol) * nv(G).normalized(),uv.col(k)); + result(0,k) = tmp.at(0); + } + else + { + tmp = ev.eval(u_sol,uv.col(k)); + result(0,k) = tmp.at(0); + } + } + interface.writeData(meshName,side==0 ? fluxName : tempName,xy,result); + + // do the coupling + precice_dt = interface.advance(dt); + + // advance variables + time += dt; + timestep += 1; + F0 = F; + + if (interface.requiresReadingCheckpoint()) + { + gsDebugVar("Read checkpoint"); + F0 = F_checkpoint; + time = t_checkpoint; + timestep = timestep_checkpoint; + solVector = solVector_checkpoint; + } + else + { + if (timestep % plotmod==0 && plot) + { + std::string fileName = "solution_" + util::to_string(timestep); + ev.options().setSwitch("plot.elements", true); + ev.options().setInt("plot.npts", 1000); + ev.writeParaview( u_sol , G, fileName); + for (size_t p=0; p!=patches.nPatches(); p++) + { + fileName = "solution_" + util::to_string(timestep) + std::to_string(p); + collection.addTimestep(fileName,time,".vts"); + } + } + } + } + + if (plot) + { + collection.save(); + exact_collection.save(); + } + + + return EXIT_SUCCESS; +} diff --git a/partitioned-heat-conduction/images/tutorial-partitioned-heat-conduction-gismo.gif b/partitioned-heat-conduction/images/tutorial-partitioned-heat-conduction-gismo.gif new file mode 100644 index 000000000..b6456897b Binary files /dev/null and b/partitioned-heat-conduction/images/tutorial-partitioned-heat-conduction-gismo.gif differ diff --git a/partitioned-heat-conduction/neumann-gismo/clean.sh b/partitioned-heat-conduction/neumann-gismo/clean.sh new file mode 100755 index 000000000..195480415 --- /dev/null +++ b/partitioned-heat-conduction/neumann-gismo/clean.sh @@ -0,0 +1,19 @@ +#!/bin/bash + +# Cleaning script for partitioned-heat-conduction directory + +echo "Cleaning unnecessary files..." + +# Remove precice-profiling directory +if [ -d "precice-profiling" ]; then + rm -rf precice-profiling + echo "Deleted 'precice-profiling' folder." +fi + +# Remove files ending with .pvd, .vts, .vtp, .log, and .txt +for ext in pvd vts vtp log txt; do + find . -type f -name "*.$ext" -exec rm -f {} \; + echo "Deleted all *.$ext files." +done + +echo "Cleaning completed!" diff --git a/partitioned-heat-conduction/neumann-gismo/run.sh b/partitioned-heat-conduction/neumann-gismo/run.sh new file mode 100755 index 000000000..6e6bddd53 --- /dev/null +++ b/partitioned-heat-conduction/neumann-gismo/run.sh @@ -0,0 +1,4 @@ +#!/bin/bash +set -e -u + +../gismo-executable/gismo-executable -c ../precice-config.xml --plot -s 1