ElmerFoamFSI  2.0
ElmerFoamFSI is fluid-solid interaction simulation application built up from OpenFOAM CFD and Elmer CSM coupled through the IMPACT multiphysics software integration infrastructure.
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Macros Groups Pages
fsiFOAMPar.C
Go to the documentation of this file.
1 #include "fsiFOAMPar.H"
2 
3 using namespace COM;
4 
5 int fsifoam_module::Initialize(int argc,char *argv[])
6 {
7  //
8  // initConinuityErrs.H
9  // ~~~~~~~~~~~~~~~~~~~
10  sumLocalContErr = 0;
11  globalContErr = 0;
12  cumulativeContErr = 0;
13  pRefCell = 0;
14  pRefValue = 0.0;
15 
16  //
17  // setRootCase.H
18  // ~~~~~~~~~~~~
19  argsPtr = new Foam::argList(argc,argv, true, true, my_window_comm);
20  if(!argsPtr){
21  Info << "FsiFoam:Initialize: Error: Failed to create args object." << endl;
22  return(1);
23  }
24  if (!argsPtr->checkRootCase()){
25  Info << "FsiFoam:Initialize: Error: Root case setup failed." << endl;
26  return(1);
27  }
28 
29  //
30  // createTime.H
31  // ~~~~~~~~~~~~
32  Foam::Info<< "Create time\n" << Foam::endl;
33  Foam::argList &args(*argsPtr);
34  runTimePtr = new Foam::Time
35  (
36  Foam::Time::controlDictName,
37  args.rootPath(),
38  args.caseName(),
39  "system",
40  "constant",
41  !args.optionFound("noFunctionObjects")
42  );
43  if(!runTimePtr){
44  Info << "FsiFoam:Initialize: Error: Failed to setup time object." << endl;
45  return(1);
46  }
47 
48  //
49  // createDynamicFvMesh.H
50  // ~~~~~~~~~~~~
51  InitFluidMesh();
52 
53  //
54  // createFields.H
55  // ~~~~~~~~~~~~
56  InitTransportProperties();
57  CreateFluidFields();
58 
59  //
60  // createStressMesh.H
61  // ~~~~~~~~~~~~
62  InitStructuresMesh();
63 
64  //
65  // createStressFields.H
66  // ~~~~~~~~~~~~
67  CreateStructuresFields();
68 
69  //
70  // readCouplingProperties.H
71  // ~~~~~~~~~~~~
72  ReadCouplingProperties();
73 
74  //
75  // createZoneToZoneInterpolators.H
76  // ~~~~~~~~~~~~
77  CreateInterZoneInterpolators();
78 
79  // initContinuityErrs is called at the top
80 
81  //
82  // findGlobalFaceZones.H
83  // ~~~~~~~~~~~~
84  FindGlobalFaceZones();
85 
86 
87  // functions specific to IMPACT FSI
88  CreateFSISurfaceMesh();
89 
90  Foam::Info << "End of initialization of openFoamPar module." << endl;
91  return(0);
92 };
93 
94 int fsifoam_module::InitFluidMesh(){
95  // Foam::argList &args(*argsPtr);
96  Foam::Time &runTime(*runTimePtr);
97  Info << "FsiFoam:Initialize: Create dynamic mesh for time = "
98  << runTime.timeName() << endl;
99  // Info << OutStream.str() << endl;
100  // OutStream.clear();
101  // OutStream.str("");
102 
103  meshPtr = dynamicFvMesh::New(IOobject(dynamicFvMesh::defaultRegion,
104  runTime.timeName(),
105  runTime,
106  IOobject::MUST_READ));
107 
108  return(0);
109 };
110 
111 int fsifoam_module::InitTransportProperties(){
112  Foam::Time &runTime(*runTimePtr);
113  dynamicFvMesh &mesh = meshPtr();
114  Info << "FsiFoam:Initialize: Reading transportProperties" << endl;
115  transportPropertiesPtr = new IOdictionary(IOobject("transportProperties",
116  runTime.constant(),
117  mesh,
118  IOobject::MUST_READ,
119  IOobject::NO_WRITE));
120  return(0);
121 };
122 
123 int fsifoam_module::CreateFluidFields(){
124  // Foam::argList &args(*argsPtr);
125  Foam::Time &runTime(*runTimePtr);
126  dynamicFvMesh &fluidsMesh = meshPtr();
127  IOdictionary &transportProperties(*transportPropertiesPtr);
128 
129 
130  nuPtr = new dimensionedScalar(transportProperties.lookup("nu"));
131  rhoFluidPtr = new dimensionedScalar(transportProperties.lookup("rho"));
132 
133  Info << "FsiFoam:CreateFluidFields: Reading field p" << endl;
134  // Info << OutStream.str() << endl;
135  // OutStream.clear();
136  // OutStream.str("");
137 
138  pPtr = new volScalarField(IOobject("p",
139  runTime.timeName(),
140  fluidsMesh,
141  IOobject::MUST_READ,
142  IOobject::AUTO_WRITE),
143  fluidsMesh);
144 
145  volScalarField &p(*pPtr);
146 
147  Info << "FsiFoam:CreateFluidFields: Reading field U" << endl;
148  // Info << OutStream.str() << endl;
149  // OutStream.clear();
150  // OutStream.str("");
151 
152  UPtr = new volVectorField(IOobject("U",runTime.timeName(),fluidsMesh,
153  IOobject::MUST_READ,IOobject::AUTO_WRITE),fluidsMesh);
154 
155  volVectorField &U = *UPtr;
156 
157  //# include "createPhi.H"
158  Info << "FsiFoam:CreateFluidFields: Reading/calculating face flux field phi" << endl;
159 
160  phiPtr = new surfaceScalarField(IOobject("phi",runTime.timeName(),fluidsMesh,
161  IOobject::READ_IF_PRESENT,IOobject::AUTO_WRITE),
162  linearInterpolate(U) & fluidsMesh.Sf());
163 
164 
165 
166  setRefCell(p, fluidsMesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
167  return(0);
168 };
169 
170 int fsifoam_module::InitStructuresMesh(){
171  Foam::Time &runTime(*runTimePtr);
172  stressMeshPtr = new fvMesh(IOobject("solid",runTime.timeName(),runTime,IOobject::MUST_READ));
173  return(stressMeshPtr != NULL);
174 };
175 
176 int fsifoam_module::CreateStructuresFields(){
177  Foam::Time &runTime(*runTimePtr);
178  fvMesh &structMesh(*stressMeshPtr);
179 
180  Info << "FsiFoam:CreateStructuresFields: Reading incremental displacement field DU" << endl;
181  // Info << OutStream.str() << endl;
182  // OutStream.clear();
183  // OutStream.str("");
184 
185  DUPtr = new volVectorField(IOobject("DU",runTime.timeName(),structMesh,
186  IOobject::MUST_READ,IOobject::AUTO_WRITE),
187  structMesh);
188 
189  volVectorField &DURef(*DUPtr);
190 
191  // Info << "DU(0): " << DU << endl;
192 
193  gradDUPtr = new volTensorField(fvc::grad(DURef));
194  volTensorField &gradDURef(*gradDUPtr);
195 
196  UsolidPtr = new volVectorField(IOobject("Usolid",runTime.timeName(),structMesh,
197  IOobject::READ_IF_PRESENT,IOobject::AUTO_WRITE),
198  DURef);
199 
200 
201  Info << "FsiFoam:CreateStructuresFields: Reading incremental displacement field DV" << endl;
202  DVPtr = new volVectorField(IOobject("DV",runTime.timeName(),structMesh,
203  IOobject::READ_IF_PRESENT,IOobject::AUTO_WRITE),
204  fvc::ddt(DURef));
205 
206 
207  Info << "FsiFoam:CreateStructuresFields: Reading accumulated velocity field V" << endl;
208  VsPtr = new volVectorField(IOobject("Vs",runTime.timeName(),structMesh,
209  IOobject::READ_IF_PRESENT,IOobject::AUTO_WRITE),
210  structMesh,
211  dimensionedVector("zero",dimVelocity,vector::zero));
212 
213 
214  Info << "FsiFoam:CreateStructuresFields: Reading accumulated stress field sigma" << endl;
215  sigmaPtr = new volSymmTensorField(IOobject("sigma",runTime.timeName(),structMesh,
216  IOobject::READ_IF_PRESENT,IOobject::AUTO_WRITE),
217  structMesh,
218  dimensionedSymmTensor("zero",dimForce/dimArea,symmTensor::zero));
219  volSymmTensorField &sigmaref(*sigmaPtr);
220 
221 
222  Info << "FsiFoam:CreateStructuresFields: Reading incremental stress field DSigma" << endl;
223  DSigmaPtr = new volSymmTensorField(IOobject("DSigma",runTime.timeName(),structMesh,
224  IOobject::READ_IF_PRESENT,IOobject::AUTO_WRITE),
225  structMesh,
226  dimensionedSymmTensor("zero",dimForce/dimArea,symmTensor::zero));
227 
228 
229  rheologyPtr = new constitutiveModel(sigmaref,DURef);
230 
231  constitutiveModel &rheology(*rheologyPtr);
232 
233  // volScalarField &rhoTest(rheology.rho());
234 
235  rhoPtr = new volScalarField(rheology.rho());
236  muPtr = new volScalarField(rheology.mu());
237  lambdaPtr = new volScalarField(rheology.lambda());
238 
239  FPtr = new volTensorField(I + gradDURef.T());
240  volTensorField &FRef(*FPtr);
241  DFPtr = new volTensorField(FRef - I);
242  JPtr = new volScalarField(det(FRef));
243 
244  solidDdtSchemePtr = new word(structMesh.schemesDict().ddtScheme("ddt("+DURef.name()+')'));
245  return(0);
246 };
247 int fsifoam_module::ReadCouplingProperties(){
248  Foam::argList &args(*argsPtr);
249  Foam::Time &runTime(*runTimePtr);
250  dynamicFvMesh &fluidMesh(meshPtr());
251  fvMesh &structMesh(*stressMeshPtr);
252 
253  Info << "FsiFoam:ReadCouplingProperties: Reading coupling properties" << endl;
254 
255  couplingPropertiesPtr = new IOdictionary(IOobject("couplingProperties",runTime.constant(),
256  fluidMesh,IOobject::MUST_READ,IOobject::NO_WRITE));
257  IOdictionary &couplingProperties(*couplingPropertiesPtr);
258 
259  solidPatchName = word(couplingProperties.lookup("solidPatch"));
260  solidPatchID = structMesh.boundaryMesh().findPatchID(solidPatchName);
261 
262  if (solidPatchID < 0)
263  {
264  FatalErrorIn(args.executable())
265  << "FsiFoam:ReadCouplingProperties: Problem with finding solid patch"
266  << abort(FatalError);
267  }
268 
269  solidZoneName = word(couplingProperties.lookup("solidZone"));
270 
271 
272  solidZoneID =
273  structMesh.faceZones().findZoneID(solidZoneName);
274 
275  if (solidZoneID < 0)
276  {
277  FatalErrorIn(args.executable())
278  << "FsiFoam:ReadCouplingProperties: Problem with finding solid zone"
279  << abort(FatalError);
280  }
281 
282  fluidPatchName = word(couplingProperties.lookup("fluidPatch"));
283  fluidPatchID = fluidMesh.boundaryMesh().findPatchID(fluidPatchName);
284 
285 
286  if (fluidPatchID < 0)
287  {
288  FatalErrorIn(args.executable())
289  << "FsiFoam:ReadCouplingProperties: Problem with finding fluid patch"
290  << abort(FatalError);
291  }
292 
293  fluidZoneName = word(couplingProperties.lookup("fluidZone"));
294  fluidZoneID = fluidMesh.faceZones().findZoneID(fluidZoneName);
295  if (fluidZoneID < 0)
296  {
297  FatalErrorIn(args.executable())
298  << "FsiFoam:ReadCouplingProperties: Problem with finding fluid zone"
299  << abort(FatalError);
300  }
301 
302  feMotionSolver = fluidMesh.objectRegistry::foundObject<tetPointVectorField>("motionU");
303  fvMotionSolver = fluidMesh.objectRegistry::foundObject<pointVectorField>("pointMotionU");
304 
305  // Info << OutStream.str() << endl;
306  // OutStream.clear();
307  // OutStream.str("");
308 
309 
310  // Grab solid patch field
311  Info << "FsiFoam:ReadCouplingProperties:Solid patch ID: " << solidPatchID << endl;
312  volVectorField &DURef(*DUPtr);
313  if
314  (
315  DURef.boundaryField()[solidPatchID].type()
316  != solidTractionFvPatchVectorField::typeName
318  )
319  {
320  FatalErrorIn(args.executable())
321  << "FsiFoam:ReadCouplingProperties: Boundary condition on " << DURef.name()
322  << " is "
323  << DURef.boundaryField()[solidPatchID].type()
324  << "for fluid -solid interface patch, instead "
325  << solidTractionFvPatchVectorField::typeName
326  //<< tractionDisplacementIncrementFvPatchVectorField::typeName
327  << abort(FatalError);
328  }
329 
330  //tractionDisplacementIncrementFvPatchVectorField& tForce =
331  // refCast<tractionDisplacementIncrementFvPatchVectorField>
332  // tForcePtr = new solidTractionFvPatchVectorField(refCast<solidTractionFvPatchVectorField>(DU.boundaryField()[solidPatchID]));
333  // solidTractionFvPatchVectorField& tForce = *tForcePtr;
334 
335 
336  accumulatedFluidInterfaceDisplacementHeaderPtr = new IOobject("accumulatedFluidInterfaceDisplacement",
337  runTime.timeName(),
338  fluidMesh,
339  IOobject::MUST_READ);
340  IOobject &accumulatedFluidInterfaceDisplacementHeader(*accumulatedFluidInterfaceDisplacementHeaderPtr);
341 
342 
343 
344  if(accumulatedFluidInterfaceDisplacementHeader.headerOk())
345  {
346  Pout << "FsiFoam:ReadCouplingProperties: Reading accumulated fluid interface displacement" << endl;
347 
348  accumulatedFluidInterfaceDisplacementPtr =
349  new vectorIOField
350  (
351  IOobject
352  (
353  "accumulatedFluidInterfaceDisplacement",
354  runTime.timeName(),
355  fluidMesh,
356  IOobject::MUST_READ,
357  IOobject::AUTO_WRITE
358  )
359  );
360  }
361  else
362  {
363  accumulatedFluidInterfaceDisplacementPtr =
364  new vectorIOField
365  (
366  IOobject
367  (
368  "accumulatedFluidInterfaceDisplacement",
369  runTime.timeName(),
370  fluidMesh,
371  IOobject::NO_READ,
372  IOobject::AUTO_WRITE
373  ),
374  vectorField
375 
376  (
377  fluidMesh.boundaryMesh()[fluidPatchID].nPoints(),
378  vector::zero
379  )
380  );
381  }
382 
383 
384  // sumLocalContErr = 0;
385  // globalContErr = 0;
386  // cumulativeContErr = 0;
387 
388 
389 
390  return(0);
391 };
392 
393 int fsifoam_module::CreateInterZoneInterpolators()
394 {
395 
396  dynamicFvMesh &fluidMesh = meshPtr();
397  fvMesh &structMesh(*stressMeshPtr);
398 
399  if(!interpolatorFluidSolidPtr || !interpolatorSolidFluidPtr)
400  {
401  deleteDemandDrivenData(interpolatorFluidSolidPtr);
402  deleteDemandDrivenData(interpolatorSolidFluidPtr);
403 
404  Info << "FsiFoam:CreateInterZoneInterpolators: Create fluid-to-solid and solid-to-fluid interpolators" << endl;
405 
406  interpolatorFluidSolidPtr = new zoneToZoneInterpolation
407  (
408  fluidMesh.faceZones()[fluidZoneID](),
409  structMesh.faceZones()[solidZoneID](),
410  intersection::VISIBLE
411  );
412 
413  interpolatorSolidFluidPtr = new zoneToZoneInterpolation
414  (
415  structMesh.faceZones()[solidZoneID](),
416  fluidMesh.faceZones()[fluidZoneID](),
417  intersection::VISIBLE
418  );
419 
420  Info << "FsiFoam:CreateInterZoneInterpolators: Check fluid-to-solid and solid-to-fluid interpolators" << endl;
421 
422  {
423  vectorField fluidPatchFaceCentres =
424  vectorField(fluidMesh.boundaryMesh()[fluidPatchID].faceCentres());
425 
426  vectorField fluidZoneFaceCentres
427  (
428  fluidMesh.faceZones()[fluidZoneID].size(),
429  vector::zero
430  );
431 
432 
433  const label fluidPatchStart =
434  fluidMesh.boundaryMesh()[fluidPatchID].start();
435 
436  forAll (fluidPatchFaceCentres, i)
437  {
438  fluidZoneFaceCentres
439  [
440  fluidMesh.faceZones()[fluidZoneID].whichFace(fluidPatchStart + i)
441  ] =
442  fluidPatchFaceCentres[i];
443  }
444 
445  // Parallel data exchange: collect faceCentres field on all processors
446  reduce(fluidZoneFaceCentres, sumOp<vectorField>());
447 
448  vectorField solidZoneFaceCentres =
449  interpolatorFluidSolidPtr->faceInterpolate
450  (
451  fluidZoneFaceCentres
452  );
453 
454  vectorField solidPatchFaceCentres
455  (
456  structMesh.boundaryMesh()[solidPatchID].size(),
457  vector::zero
458  );
459 
460  const label solidPatchStart =
461  structMesh.boundaryMesh()[solidPatchID].start();
462 
463  forAll(solidPatchFaceCentres, i)
464  {
465  solidPatchFaceCentres[i] =
466  solidZoneFaceCentres
467  [
468  structMesh.faceZones()[solidZoneID]
469  .whichFace(solidPatchStart + i)
470  ];
471  }
472 
473  scalar maxDist = gMax
474  (
475  mag
476  (
477  solidPatchFaceCentres
478  - structMesh.boundaryMesh()[solidPatchID].faceCentres()
479  )
480  );
481 
482  Info << "FsiFoam:CreateInterZoneInterpolators: Fluid-to-solid face interpolation error: " << maxDist
483  << endl;
484  }
485 
486 
487 
488  {
489  vectorField solidPatchFaceCentres =
490  vectorField(structMesh.boundaryMesh()[solidPatchID].faceCentres());
491 
492  vectorField solidZoneFaceCentres
493  (
494  structMesh.faceZones()[solidZoneID].size(),
495  vector::zero
496  );
497 
498  const label solidPatchStart =
499  structMesh.boundaryMesh()[solidPatchID].start();
500 
501  forAll (solidPatchFaceCentres, i)
502  {
503  solidZoneFaceCentres
504  [
505  structMesh.faceZones()[solidZoneID]
506  .whichFace(solidPatchStart + i)
507  ] =
508  solidPatchFaceCentres[i];
509  }
510 
511  // Parallel data exchange: collect faceCentres field on all processors
512  reduce(solidZoneFaceCentres, sumOp<vectorField>());
513 
514  vectorField fluidZoneFaceCentres =
515  interpolatorSolidFluidPtr->faceInterpolate
516  (
517  solidZoneFaceCentres
518  );
519 
520  vectorField fluidPatchFaceCentres
521  (
522  fluidMesh.boundaryMesh()[fluidPatchID].size(),
523  vector::zero
524  );
525 
526  const label fluidPatchStart =
527  fluidMesh.boundaryMesh()[fluidPatchID].start();
528 
529  forAll(fluidPatchFaceCentres, i)
530  {
531  fluidPatchFaceCentres[i] =
532  fluidZoneFaceCentres
533  [
534  fluidMesh.faceZones()[fluidZoneID]
535  .whichFace(fluidPatchStart + i)
536  ];
537  }
538 
539  scalar maxDist = gMax
540  (
541  mag
542  (
543  fluidPatchFaceCentres
544  - fluidMesh.boundaryMesh()[fluidPatchID].faceCentres()
545  )
546  );
547 
548  Info << "FsiFoam:CreateInterZoneInterpolators: Solid-to-fluid face interpolation error: " << maxDist
549  << endl;
550  }
551 
552  }
553 
554  return(0);
555 };
556 
557 int fsifoam_module::FindGlobalFaceZones()
558 {
559 
560  fvMesh &structMesh(*stressMeshPtr);
561  Foam::argList &args = *argsPtr;
562 
563  {
564  SLList<label> globalFaceZonesSet;
565 
566  const faceZoneMesh& faceZones = structMesh.faceZones();
567 
568  forAll(faceZones, zoneI)
569  {
570  const faceZone& curFaceZone = faceZones[zoneI];
571 
572  forAll(curFaceZone, faceI)
573  {
574  // if unused face exist
575  if (curFaceZone[faceI] >= structMesh.nFaces())
576  {
577  globalFaceZonesSet.insert(zoneI);
578  break;
579  }
580  }
581  }
582 
583  globalFaceZones = labelList(globalFaceZonesSet);
584  }
585  globalToLocalFaceZonePointMap.resize(globalFaceZones.size());
586 
587  forAll(globalFaceZones, zoneI)
588  {
589  label curZoneID = globalFaceZones[zoneI];
590  labelList curMap(structMesh.faceZones()[curZoneID]().nPoints(), -1);
591  vectorField fzGlobalPoints =
592  structMesh.faceZones()[curZoneID]().localPoints();
593 
594  //- set all slave points to zero because only the master order is used
595  if(!Pstream::master())
596  {
597  fzGlobalPoints *= 0.0;
598  }
599 
600  //- pass points to all procs
601  reduce(fzGlobalPoints, sumOp<vectorField>());
602 
603  //- now every proc has the master's list of FZ points
604  //- every proc must now find the mapping from their local FZ points to
605  //- the global FZ points
606  const vectorField& fzLocalPoints =
607  structMesh.faceZones()[curZoneID]().localPoints();
608  const edgeList& fzLocalEdges =
609  structMesh.faceZones()[curZoneID]().edges();
610  const labelListList& fzPointEdges =
611  structMesh.faceZones()[curZoneID]().pointEdges();
612  scalarField minEdgeLength(fzLocalPoints.size(), GREAT);
613 
614  forAll(minEdgeLength, pI)
615  {
616  const labelList& curPointEdges = fzPointEdges[pI];
617  forAll(curPointEdges, eI)
618  {
619  scalar Le = fzLocalEdges[curPointEdges[eI]].mag(fzLocalPoints);
620  if (Le < minEdgeLength[pI])
621  {
622  minEdgeLength[pI] = Le;
623  }
624  }
625  }
626 
627  forAll(fzGlobalPoints, globalPointI)
628  {
629  // scalar minDist = GREAT;
630 
631  forAll(fzLocalPoints, procPointI)
632  {
633  scalar curDist =
634  mag
635  (
636  fzLocalPoints[procPointI]
637  - fzGlobalPoints[globalPointI]
638  );
639 
640  // if (curDist < minDist)
641  // {
642  // minDist = curDist;
643  // }
644 
645  if (curDist < 1e-4*minEdgeLength[procPointI])
646  {
647  curMap[globalPointI] = procPointI;
648  break;
649  }
650  }
651 
652  // if (curMap[globalPointI] == -1)
653  // {
654  // Pout << "minDist: " << minDist << endl;
655  // }
656  }
657 
658  forAll(curMap, globalPointI)
659  {
660  if (curMap[globalPointI] == -1)
661  {
662  FatalErrorIn(args.executable())
663  << "FsiFoam:FindGlobalFaceZones: local to global face zone point map is not correct"
664  << abort(FatalError);
665  }
666  }
667 
668  globalToLocalFaceZonePointMap[zoneI] = curMap;
669  }
670  return(0);
671 };
672 int fsifoam_module::ReadPISOControls()
673 {
674  dynamicFvMesh &fluidMesh(meshPtr());
675  dictionary &piso(fluidMesh.solutionDict().subDict("PISO"));
676 
677  nCorr = readInt(piso.lookup("nCorrectors"));
678  nNonOrthCorr = piso.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
679  momentumPredictor = piso.lookupOrDefault<Switch>("momentumPredictor", true);
680  transonic = piso.lookupOrDefault<Switch>("transonic", false);
681  nOuterCorr = piso.lookupOrDefault<int>("nOuterCorrectors", 1);
682 
683  return(0);
684 };
685 int fsifoam_module::ReadFSIControls()
686 {
687  IOdictionary &couplingProperties(*couplingPropertiesPtr);
688  // dynamicFvMesh &fluidMesh(meshPtr());
689  // Foam::Time &runTime(*runTimePtr);
690 
691  if (couplingProperties.found("couplingScheme"))
692  {
693  couplingScheme = word(couplingProperties.lookup("couplingScheme"));
694  if
695  (
696  (couplingScheme == "IQN-ILS")
697  || (couplingScheme == "Aitken")
698  || (couplingScheme == "FixedRelaxation")
699  || (couplingScheme == "NonIterative")
700  )
701  {
702  Info<< "FsiFoam:ReadFSIControl: Selecting coupling scheme " << couplingScheme << endl;
703  }
704  else
705  {
706  FatalErrorIn
707  (
708  "readFsiProperties"
709  ) << "couplingScheme: " << couplingScheme
710  << " is not a valid choice. "
711  << "Options are: IQN-ILS, Aitken, FixedRelaxation"
712  << abort(FatalError);
713  }
714  }
715 
716  interfaceDeformationLimit =
717  scalar(readScalar(couplingProperties.lookup("interfaceDeformationLimit")));
718 
719 
720  // if(dynamicMeshDictPtr)
721  // delete dynamicMeshDictPtr;
722  // dynamicMeshDictPtr = new IOdictionary(IOobject
723  // (
724  // "dynamicMeshDict",
725  // runTime.constant(),
726  // fluidMesh,
727  // IOobject::MUST_READ,
728  // IOobject::NO_WRITE,
729  // false
730  // ));
731 
732 
733 
734  fsiRelaxationFactorMin = scalar(readScalar(couplingProperties.lookup("fsiRelaxationFactor")));
735  fsiRelaxationFactor = fsiRelaxationFactorMin;
736  outerCorrTolerance = scalar(readScalar(couplingProperties.lookup("outerCorrTolerance")));
737  Info << "FsiFoam:ReadFSIControl: outerCorrTolerance = " << outerCorrTolerance << endl;
738  fsi = Switch(couplingProperties.lookup("fsi"));
739  if(fsi)
740  Info << "FsiFoam:ReadFSIControl: FSI IS ENABLED." << endl;
741  else
742  Info << "FsiFoam:ReadFSIControl: FSI IS DISABLED." << endl;
743  return(0);
744 };
745 
746 int fsifoam_module::Step(){
747  Foam::argList &args(this->ArgList());
748  Foam::Time &runTime(this->RunTime());
749  dynamicFvMesh &fluidsMesh(this->FluidMesh());
750  // IOdictionary &transportProperties(this->TransportProperties());
751  dimensionedScalar &nu(this->nu());
752  dimensionedScalar &rhoFluid(this->rhoFluid());
753  volScalarField &p(this->p());
754  volVectorField &U(this->U());
755  surfaceScalarField &phi(this->phi());
756  // label pRefCell = 0;
757  // scalar pRefValue = 0.0;
758  // setRefCell(p, fluidsMesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
759 
760  fvMesh &structuresMesh(this->StructuresMesh());
761 
762 
763  volVectorField &DU(this->DU());
764  volTensorField &gradDU(this->gradDU());
765  volVectorField &Usolid(this->Usolid());
766  volVectorField &DV(this->DV());
767  volVectorField &Vs(this->Vs());
768  volSymmTensorField &sigma(this->sigma());
769  volSymmTensorField &DSigma(this->DSigma());
770  // constitutiveModel &Rheology(){return(*rheologyPtr);};
771  volScalarField &rho(this->rhoSolid());
772  volScalarField &mu(this->mu());
773  volScalarField &lambda(this->lambda());
774  volTensorField &F(this->F());
775  volTensorField &DF(this->DF());
776  volScalarField &J(this->J());
777  word &solidDdtScheme(this->SolidDdtScheme());
778  // pointMesh &pStressMesh(){return(*pStressMeshPtr);};
779  // pointVectorField &PointDU(){return(*pointDUPtr);};
780 
781 
782  // IOdictionary &couplingProperties(this->CouplingProperties());
783  label solidPatchID(this->SolidPatchID());
784  // word solidPatchName(this->SolidPatchName());
785  label fluidPatchID(this->FluidPatchID());
786  // word fluidPatchName(this->FluidPatchName());
787  label fluidZoneID(this->FluidZoneID());
788  label solidZoneID(this->SolidZoneID());
789 
790  bool feMotionSolver(this->FEMotion());
791  bool fvMotionSolver(this->FVMotion());
792  vectorField &accumulatedFluidInterfaceDisplacement(this->AccumulatedFluidInterfaceDisplacements());
793  Info << "FsiFoam:Step: accumulatedFluidInterfaceDisplacement = " << endl;
794  Info << accumulatedFluidInterfaceDisplacement << endl;
795  solidTractionFvPatchVectorField &tForce(this->tForce());
796 
797  zoneToZoneInterpolation &interpolatorFluidSolid(this->interpFluidSolid());
798  zoneToZoneInterpolation &interpolatorSolidFluid(this->interpSolidFluid());
799 
800  scalar &sumLocalContErr(this->LocalContErr());
801  scalar &globalContErr(this->GlobalContErr());
802  scalar &cumulativeContErr(this->CumulativeContErr());
803 
804  labelList &globalFaceZones(this->GlobalFaceZones());
805  labelListList &globalToLocalFaceZonePointMap(this->GlobalToLocalFaceZonePointMap());
806  runTime++;
807  if(!runTime.end()){
808  // for (runTime++; !runTime.end(); runTime++)
809  // {
810  // Info << "Time = " << runTime.timeName() << nl << endl;
811 
812  //# include "readPISOControls.H"
813  //# include "readFsiControls.H"
814  //# include "createStressPointMesh.H"
815 
816  this->ReadPISOControls();
817  this->ReadFSIControls();
818  word couplingScheme(this->CouplingScheme());
819  // label &outerCorr(this->OuterCorr());
820  scalar &fsiRelaxationFactor(this->FSIRelaxationFactor());
821  // scalar &fsiRelaxationFactorMin(this->FSIRelaxationFactorMin());
822  Switch fsi(this->FSIEnabled());
823  int &nCorrPISO(this->NCorrPISO());
824  int &nNonOrthCorr(this->NNonOrthCorr());
825  scalar &interfaceDeformationLimit(this->InterfaceDeformationLimit());
826  scalar &outerCorrTolerance(this->OuterCorrTolerance());
827  int &nOuterCorr(this->NOuterCorr());
828 
829  // this->CreateStructuresPointMesh();
830 
831 
832  // Create point stress mesh for interpolation
833  pointMesh StructuresPointMesh(structuresMesh);
834 
835  pointPatchInterpolation patchPointInterpolator(structuresMesh);
836 
837 
838  // Create point displacement field (check fixedValue patches)
839 
840  wordList types
841  (
842  StructuresPointMesh.boundary().size(),
843  calculatedFvPatchVectorField::typeName
844  );
845 
846  // wordList types = DU.boundaryField().types();
847 
848  forAll(DU.boundaryField().types(), patchI)
849  {
850  if
851  (
852  DU.boundaryField().types()[patchI]
853  == fixedValueFvPatchVectorField::typeName
854  )
855  {
856  types[patchI] = fixedValueFvPatchVectorField::typeName;
857  }
858  }
859 
860  pointVectorField pointDU
861  (
862  IOobject
863  (
864  "pointDU",
865  runTime.timeName(),
866  structuresMesh
867  ),
868  StructuresPointMesh,
869  dimensionedVector("zero", dimLength, vector::zero),
870  types
871  );
872  //# include "createInterfaceFields.H"
873  vectorField fluidPatchPointsDispl
874  (
875  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
876  vector::zero
877  );
878 
879  vectorField fluidPatchPointsDisplOld
880  (
881  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
882  vector::zero
883  );
884 
885  vectorField solidPatchPointsDispl
886  (
887  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
888  vector::zero
889  );
890 
891  vectorField fsiResidual
892  (
893  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
894  vector::zero
895  );
896 
897  vectorField fsiResidualOld
898  (
899  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
900  vector::zero
901  );
902 
903  vectorField solidPatchTraction
904  (
905  structuresMesh.boundary()[solidPatchID].size(),
906  vector::zero
907  );
908 
909  scalarField solidPatchPressure
910  (
911  structuresMesh.boundary()[solidPatchID].size(),
912  0.0
913  );
914 
915  scalar initialFsiResidualNorm = 0;
916  scalar fsiResidualNorm = 0;
917  label outerCorr=0;
918 
919  do
920  {
921  outerCorr++;
922 
923  //# include "setInterfaceDisplacement.H"
924  Info << "\nTime = " << runTime.timeName()
925  << ", iteration: " << outerCorr << endl;
926 
927  if (outerCorr < 3 || couplingScheme == "FixedRelaxation")
928  {
929  Info << "FsiFoam:Step: Current fsi under-relaxation factor: "
930  << fsiRelaxationFactor << endl;
931 
932  fluidPatchPointsDisplOld = fluidPatchPointsDispl;
933 
934  fluidPatchPointsDispl += fsiRelaxationFactor*fsiResidual;
935  }
936  else
937  {
938  if (couplingScheme == "Aitken")
939  {
940  fsiRelaxationFactor =
941  -fsiRelaxationFactor
942  *(
943  gSum
944  (
945  fsiResidualOld
946  &(fsiResidual - fsiResidualOld)
947  )
948  /(
949  gSum
950  (
951  (fsiResidual - fsiResidualOld)
952  &(fsiResidual - fsiResidualOld)
953  )
954  )
955  );
956 
957  fsiRelaxationFactor = mag(fsiRelaxationFactor);
958 
959  Info << "FsiFoam:Step: Current fsi under-relaxation factor (Aitken): "
960  << fsiRelaxationFactor << endl;
961 
962  fluidPatchPointsDisplOld = fluidPatchPointsDispl;
963 
964  fluidPatchPointsDispl +=
965  fsiRelaxationFactor*fsiResidual;
966  }
967  // else if (couplingScheme == "IQN-ILS")
968  // {
969 
970  // }
971  }
972 
973  //# include "moveFluidMesh.H"
974  {
975  // Move fluid mesh
976  const vectorField& n =
977  fluidsMesh.boundaryMesh()[fluidPatchID].pointNormals();
978 
979  primitivePatchInterpolation patchInterpolator
980  (
981  fluidsMesh.boundaryMesh()[fluidPatchID]
982  );
983 
984  scalarField pointDeltaCoeffs =
985  patchInterpolator.faceToPointInterpolate
986  (
987  fluidsMesh.boundary()[fluidPatchID].deltaCoeffs()
988  );
989 
990  scalar delta =
991  gMax
992  (
993  mag
994  (
995  n
996  & (
997  accumulatedFluidInterfaceDisplacement
998  + fluidPatchPointsDispl
999  - fluidPatchPointsDisplOld
1000  )
1001  )
1002  *pointDeltaCoeffs
1003  );
1004 
1005  Info << "FsiFoam:Step: Maximal accumulated displacement of interface points: "
1006  << delta << endl;
1007 
1008  if(delta < interfaceDeformationLimit)
1009  {
1010  // Move only interface points
1011  // Masoud
1012  Info << "FsiFoam:Step: Moving only interface points..." << endl;
1013  // Masoud: End
1014  pointField newPoints = fluidsMesh.allPoints();
1015 
1016  const labelList& meshPoints =
1017  fluidsMesh.boundaryMesh()[fluidPatchID].meshPoints();
1018 
1019  forAll (fluidPatchPointsDispl, pointI)
1020  {
1021  //Info << "1. meshPoints[" << pointI << "] = " << meshPoints[pointI]
1022  //<< endl;
1023  newPoints[meshPoints[pointI]] +=
1024  fluidPatchPointsDispl[pointI]
1025  - fluidPatchPointsDisplOld[pointI];
1026  }
1027 
1028  twoDPointCorrector twoDCorrector(fluidsMesh);
1029 
1030  twoDCorrector.correctPoints(newPoints);
1031 
1032  fluidsMesh.movePoints(newPoints);
1033 
1034  // Accumulate interface points displacement
1035  accumulatedFluidInterfaceDisplacement +=
1036  fluidPatchPointsDispl
1037  - fluidPatchPointsDisplOld;
1038  }
1039  else
1040  {
1041  // Move whole fluid mesh
1042  Info << "FsiFoam:Step: Moving the whole mesh .... " << endl;
1043  pointField newPoints = fluidsMesh.allPoints();
1044 
1045  const labelList& meshPoints =
1046  fluidsMesh.boundaryMesh()[fluidPatchID].meshPoints();
1047 
1048  forAll (accumulatedFluidInterfaceDisplacement, pointI)
1049  {
1050  //Info << "2. meshPoints[" << pointI << "] = " << meshPoints[pointI]
1051  //<< endl;
1052  newPoints[meshPoints[pointI]] -=
1053  accumulatedFluidInterfaceDisplacement[pointI];
1054  }
1055 
1056  twoDPointCorrector twoDCorrector(fluidsMesh);
1057 
1058  twoDCorrector.correctPoints(newPoints);
1059 
1060  fluidsMesh.movePoints(newPoints);
1061 
1062  accumulatedFluidInterfaceDisplacement +=
1063  fluidPatchPointsDispl
1064  - fluidPatchPointsDisplOld;
1065 
1066 
1067  if (feMotionSolver)
1068  {
1069  tetPointVectorField& motionU =
1070  const_cast<tetPointVectorField&>
1071  (
1072  fluidsMesh.objectRegistry:: lookupObject<tetPointVectorField>
1073  (
1074  "motionU"
1075  )
1076  );
1077 
1078  fixedValueTetPolyPatchVectorField& motionUFluidPatch =
1079  refCast<fixedValueTetPolyPatchVectorField>
1080  (
1081  motionU.boundaryField()[fluidPatchID]
1082  );
1083 
1084  tetPolyPatchInterpolation tppi
1085  (
1086  refCast<const faceTetPolyPatch>(motionUFluidPatch.patch())
1087  );
1088 
1089  motionUFluidPatch ==
1090  tppi.pointToPointInterpolate
1091  (
1092  accumulatedFluidInterfaceDisplacement
1093  /runTime.deltaT().value()
1094  );
1095  }
1096  else if (fvMotionSolver)
1097  {
1098  pointVectorField& motionU =
1099  const_cast<pointVectorField&>
1100  (
1101  fluidsMesh.objectRegistry:: lookupObject<pointVectorField>
1102  (
1103  "pointMotionU"
1104  )
1105  );
1106 
1107  fixedValuePointPatchVectorField& motionUFluidPatch =
1108  refCast<fixedValuePointPatchVectorField>
1109  (
1110  motionU.boundaryField()[fluidPatchID]
1111  );
1112 
1113  motionUFluidPatch ==
1114  accumulatedFluidInterfaceDisplacement
1115  /runTime.deltaT().value();
1116  }
1117  else
1118  {
1119  FatalErrorIn(args.executable())
1120  << "FsiFoam:Step: Problem with mesh motion solver selection"
1121  << abort(FatalError);
1122  }
1123 
1124  fluidsMesh.update();
1125 
1126  accumulatedFluidInterfaceDisplacement = vector::zero;
1127  }
1128  }
1129 
1130  // Info << "DU after move fluid: " << DU << endl;
1131 
1132  //# include "solveFluid.H"
1133  if(fluidsMesh.moving())
1134  {
1135  // Make the fluxes relative
1136  phi -= fvc::meshPhi(U);
1137  }
1138 
1139 
1140  //# include "CourantNo.H"
1141  scalar CoNum = 0.0;
1142  scalar meanCoNum = 0.0;
1143  scalar velMag = 0.0;
1144 
1145  if (fluidsMesh.nInternalFaces())
1146  {
1147  surfaceScalarField magPhi = mag(phi);
1148 
1149  surfaceScalarField SfUfbyDelta =
1150  fluidsMesh.surfaceInterpolation::deltaCoeffs()*magPhi;
1151 
1152  const scalar deltaT = runTime.deltaT().value();
1153 
1154  CoNum = max(SfUfbyDelta/fluidsMesh.magSf()).value()*deltaT;
1155 
1156  meanCoNum = (sum(SfUfbyDelta)/sum(fluidsMesh.magSf())).value()*deltaT;
1157 
1158  velMag = max(magPhi/fluidsMesh.magSf()).value();
1159  }
1160 
1161  Info<< "FsiFoam:Step: Courant Number mean: " << meanCoNum
1162  << " max: " << CoNum
1163  << " velocity magnitude: " << velMag
1164  << endl;
1165 
1166  fvVectorMatrix UEqn
1167  (
1168  fvm::ddt(U)
1169  + fvm::div(phi, U)
1170  - fvm::laplacian(nu, U)
1171  );
1172 
1173  solve(UEqn == -fvc::grad(p));
1174 
1175  // --- PISO loop
1176  volScalarField rUA = 1.0/UEqn.A();
1177 
1178  for (int corr=0; corr<nCorrPISO; corr++)
1179  {
1180  U = rUA*UEqn.H();
1181  phi = (fvc::interpolate(U) & fluidsMesh.Sf());
1182 
1183  adjustPhi(phi, U, p);
1184 
1185  for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
1186  {
1187  fvScalarMatrix pEqn
1188  (
1189  fvm::laplacian(rUA, p)
1190  == fvc::div(phi)
1191  );
1192 
1193  pEqn.setReference(pRefCell, pRefValue);
1194  pEqn.solve();
1195 
1196  if (nonOrth == nNonOrthCorr)
1197  {
1198  phi -= pEqn.flux();
1199  }
1200  }
1201 
1202  //# include "continuityErrs.H"
1203  {
1204  volScalarField contErr = fvc::div(phi);
1205 
1206  sumLocalContErr = runTime.deltaT().value()*
1207  mag(contErr)().weightedAverage(fluidsMesh.V()).value();
1208 
1209  globalContErr = runTime.deltaT().value()*
1210  contErr.weightedAverage(fluidsMesh.V()).value();
1211 
1212  cumulativeContErr += globalContErr;
1213 
1214  Info<< "FsiFoam:Step: time step continuity errors : sum local = " << sumLocalContErr
1215  << ", global = " << globalContErr
1216  << ", cumulative = " << cumulativeContErr
1217  << endl;
1218  }
1219 
1220  U -= rUA*fvc::grad(p);
1221  U.correctBoundaryConditions();
1222  }
1223 
1224  // Info << "DU after solve fluid: " << DU << endl;
1225 
1226  //# include "setInterfaceForce.H"
1227  {
1228  Info << "FsiFoam:Step: Setting traction on solid patch" << endl;
1229 
1230  // vectorField fluidPatchTraction =
1231  // -rhoFluid.value()*nu.value()
1232  // *U.boundaryField()[fluidPatchID].snGrad()
1233  // + rhoFluid.value()*p.boundaryField()[fluidPatchID]
1234  // *fluidsMesh.boundary()[fluidPatchID].nf();
1235 
1236  vectorField fluidPatchTraction =
1237  -rhoFluid.value()*nu.value()
1238  *U.boundaryField()[fluidPatchID].snGrad();
1239 
1240  // Info << "Fluid patch traction: " << fluidPatchTraction << endl;
1241 
1242  scalarField fluidPatchPressure =
1243  rhoFluid.value()*p.boundaryField()[fluidPatchID];
1244 
1245  // Info << "Fluid patch pressure: " << fluidPatchPressure << endl;
1246 
1247  vectorField fluidZoneTraction
1248  (
1249  fluidsMesh.faceZones()[fluidZoneID].size(),
1250  vector::zero
1251  );
1252 
1253  // Info << "Fluid zone traction: " << fluidZoneTraction << endl;
1254 
1255  const label fluidPatchStart =
1256  fluidsMesh.boundaryMesh()[fluidPatchID].start();
1257 
1258  forAll(fluidPatchTraction, i)
1259  {
1260  fluidZoneTraction
1261  [
1262  fluidsMesh.faceZones()[fluidZoneID].whichFace(fluidPatchStart + i)
1263  ] =
1264  fluidPatchTraction[i];
1265  }
1266 
1267  // Parallel data exchange: collect pressure field on all processors
1268  reduce(fluidZoneTraction, sumOp<vectorField>());
1269 
1270 
1271  scalarField fluidZonePressure
1272  (
1273  fluidsMesh.faceZones()[fluidZoneID].size(),
1274  0.0
1275  );
1276 
1277  forAll(fluidPatchPressure, i)
1278  {
1279  fluidZonePressure
1280  [
1281  fluidsMesh.faceZones()[fluidZoneID].whichFace(fluidPatchStart + i)
1282  ] =
1283  fluidPatchPressure[i];
1284  }
1285 
1286  // Parallel data exchange: collect pressure field on all processors
1287  reduce(fluidZonePressure, sumOp<scalarField>());
1288 
1289  vectorField solidZoneTraction =
1290  interpolatorFluidSolid.faceInterpolate
1291  (
1292  fluidZoneTraction
1293  );
1294 
1295  // Info << "Solid Zone Traction: " << solidZoneTraction << endl;
1296 
1297  scalarField solidZonePressure =
1298  interpolatorFluidSolid.faceInterpolate
1299  (
1300  fluidZonePressure
1301  );
1302 
1303  // Info << "Solid Zone Pressure: " << solidZonePressure << endl;
1304 
1305  const label solidPatchStart =
1306  structuresMesh.boundaryMesh()[solidPatchID].start();
1307 
1308  forAll(solidPatchTraction, i)
1309  {
1310  solidPatchTraction[i] =
1311  solidZoneTraction
1312  [
1313  structuresMesh.faceZones()[solidZoneID]
1314  .whichFace(solidPatchStart + i)
1315  ];
1316  }
1317 
1318  forAll(solidPatchPressure, i)
1319  {
1320  solidPatchPressure[i] =
1321  solidZonePressure
1322  [
1323  structuresMesh.faceZones()[solidZoneID]
1324  .whichFace(solidPatchStart + i)
1325  ];
1326  }
1327 
1328  if (fsi)
1329  {
1330  tForce.traction() = solidPatchTraction;
1331  tForce.pressure() = solidPatchPressure;
1332  }
1333 
1334  vector totalTractionForce =
1335  gSum
1336  (
1337  solidPatchTraction
1338  *structuresMesh.magSf().boundaryField()[solidPatchID]
1339  );
1340 
1341  Info << "FsiFoam:Step: Total traction force = "
1342  << totalTractionForce << endl;
1343  }
1344 
1345  // Info << "DU after setting interface force: " << DU << endl;
1346 
1347  //# include "solveSolid.H"
1348  if (solidDdtScheme == fv::EulerDdtScheme<vector>::typeName)
1349  {
1350  Info << "FsiFoam:Step: Solve Solid: Euler" << endl;
1351  //# include "solveSolidEuler.H"
1352  {
1353  const dictionary& stressControl =
1354  structuresMesh.solutionDict().subDict("solidMechanics");
1355 
1356  int nCorrStruct(readInt(stressControl.lookup("nCorrectors")));
1357  scalar convergenceTolerance(readScalar(stressControl.lookup("DU")));
1358  //# include "readSolidMechanicsControls.H"
1359 
1360  int iCorr = 0;
1361  lduMatrix::solverPerformance solverPerf;
1362  scalar initialResidual = 0;
1363 
1364  lduMatrix::debug = 0;
1365 
1366  //# include "EulerCoeffs.H"
1367  dimensionedScalar Cn("Cn", dimless/dimTime, 1.0/runTime.deltaT().value());
1368  dimensionedScalar Co("Co", dimless/dimTime, 1.0/runTime.deltaT().value());
1369 
1370  do
1371  {
1372  DU.storePrevIter();
1373 
1374  fvVectorMatrix DUEqn
1375  (
1376  Cn*rho*fvm::ddt(DU)
1377  - Co*rho*DV.oldTime()
1378  ==
1379  fvm::laplacian(2*mu + lambda, DU, "laplacian(DDU,DU)")
1380  - fvc::laplacian(mu + lambda, DU, "laplacian(DDU,DU)")
1381  + fvc::div
1382  (
1383  mu*gradDU.T()
1384  + lambda*(I*tr(gradDU))
1385  + mu*(gradDU&gradDU.T())
1386  + 0.5*lambda*(I*tr(gradDU & gradDU.T()))
1387  + (sigma & DF.T())
1388  + (DSigma & DF.T()),
1389  "div(sigma)"
1390  )
1391  );
1392 
1393  solverPerf = DUEqn.solve();
1394 
1395  DU.relax();
1396 
1397  if(iCorr == 0)
1398  {
1399  initialResidual = solverPerf.initialResidual();
1400  }
1401 
1402  gradDU = fvc::grad(DU);
1403 
1404  DF = gradDU.T();
1405 
1406  //# include "calculateDSigma.H"
1407  {
1408  volSymmTensorField Depsilon = symm(gradDU)
1409  + 0.5*symm(gradDU & gradDU.T());
1410 
1411  DSigma = 2*mu*Depsilon + I*(lambda*tr(Depsilon));
1412  }
1413  }
1414  while
1415  (
1416  solverPerf.initialResidual() > convergenceTolerance
1417  && ++iCorr < nCorrStruct
1418  );
1419 
1420  Info << "FsiFoam:Step: Solving for " << DU.name()
1421  << ", Initial residual = " << initialResidual
1422  << ", Final residual = " << solverPerf.initialResidual()
1423  << ", No outer iterations " << iCorr << endl;
1424 
1425  DV = fvc::ddt(DU);
1426 
1427  lduMatrix::debug = 1;
1428  }
1429  }
1430  else if (solidDdtScheme == fv::backwardDdtScheme<vector>::typeName)
1431  {
1432  Info << "FsiFoam:Step: Solve Solid: Backward" << endl;
1433  //# include "solveSolidBackward.H"
1434 
1435  {
1436  //# include "readSolidMechanicsControls.H"
1437  const dictionary& stressControl =
1438  structuresMesh.solutionDict().subDict("solidMechanics");
1439 
1440  int nCorrStruct(readInt(stressControl.lookup("nCorrectors")));
1441  scalar convergenceTolerance(readScalar(stressControl.lookup("DU")));
1442 
1443  int iCorr = 0;
1444  lduMatrix::solverPerformance solverPerf;
1445  scalar initialResidual = 0;
1446 
1447  lduMatrix::debug = 0;
1448 
1449  //# include "backwardCoeffs.H"
1450  dimensionedScalar Cn("Cn", dimless/dimTime, 0.0);
1451  dimensionedScalar Co("Co", dimless/dimTime, 0.0);
1452  dimensionedScalar Coo("Coo", dimless/dimTime, 0.0);
1453 
1454  scalar deltaT = runTime.deltaT().value();
1455  scalar deltaT0 = runTime.deltaT0().value();
1456 
1457  Cn.value() = 1 + deltaT/(deltaT + deltaT0);
1458  Coo.value() = deltaT*deltaT/(deltaT0*(deltaT + deltaT0));
1459  Co.value() = Cn.value() + Coo.value();
1460 
1461  if(runTime.timeIndex() == 1)
1462  {
1463  Cn.value() = 1.0;
1464  Co.value() = 1.0;
1465  Coo.value() = 0.0;
1466  }
1467 
1468  Cn.value() /= deltaT;
1469  Co.value() /= deltaT;
1470  Coo.value() /= deltaT;
1471 
1472  // Info << "ncorr = " << nCorr << endl;
1473  // Info << "convergence Tol: " << convergenceTolerance << endl;
1474  // Info << "DU:" << DU << endl;
1475  // Info << Cn << "," << Co << "," << Coo << endl;
1476  // Info << "RHO: " << rho << endl;
1477  // Info << "MU: " << mu << endl;
1478  // Info << "LAMBDA: " << lambda << endl;
1479  // Info << "gradDU:" << gradDU << endl;
1480  // Info << "gradDU.T(): " << gradDU.T() << endl;
1481  // Info << "DF: " << DF << endl;
1482  // Info << "DF.T(): " << DF.T() << endl;
1483  // Info << "SIGMA: " << sigma << endl;
1484  // Info << "DSigma: " << DSigma << endl;
1485  // Info << "DV: " << DV << endl;
1486  // Info << "DV.oldTime(): " << DV.oldTime() << endl;
1487  // Info << "DV.oldTime().oldTime(): " << DV.oldTime().oldTime() << endl;
1488 
1489  do
1490  {
1491  DU.storePrevIter();
1492 
1493  fvVectorMatrix DUEqn
1494  (
1495  Cn*rho*fvm::ddt(DU)
1496  - Co*rho*DV.oldTime()
1497  + Coo*rho*DV.oldTime().oldTime()
1498  ==
1499  fvm::laplacian(2*mu + lambda, DU, "laplacian(DDU,DU)")
1500  - fvc::laplacian(mu + lambda, DU, "laplacian(DDU,DU)")
1501  + fvc::div
1502  (
1503  mu*gradDU.T()
1504  + lambda*(I*tr(gradDU))
1505  + mu*(gradDU&gradDU.T())
1506  + 0.5*lambda*(I*tr(gradDU & gradDU.T()))
1507  + (sigma & DF.T())
1508  + (DSigma & DF.T()),
1509  "div(sigma)"
1510  )
1511  );
1512 
1513  solverPerf = DUEqn.solve();
1514 
1515  DU.relax();
1516 
1517  if(iCorr == 0)
1518  {
1519  initialResidual = solverPerf.initialResidual();
1520  }
1521 
1522  gradDU = fvc::grad(DU);
1523 
1524  DF = gradDU.T();
1525 
1526  //# include "calculateDSigma.H"
1527  {
1528  volSymmTensorField Depsilon = symm(gradDU)
1529  + 0.5*symm(gradDU & gradDU.T());
1530 
1531  DSigma = 2*mu*Depsilon + I*(lambda*tr(Depsilon));
1532  }
1533  // Info << "Solving for " << DU.name()
1534  // << ", Initial residual = " << initialResidual
1535  // << ", Final residual = " << solverPerf.initialResidual()
1536  // << ", No outer iterations " << iCorr << endl;
1537 
1538  }
1539  while
1540  (
1541  solverPerf.initialResidual() > convergenceTolerance
1542  && ++iCorr < nCorrStruct
1543  );
1544 
1545  Info << "FsiFoam:Step: Solving for " << DU.name()
1546  << ", Initial residual = " << initialResidual
1547  << ", Final residual = " << solverPerf.initialResidual()
1548  << ", No outer iterations " << iCorr << endl;
1549 
1550  DV = fvc::ddt(DU);
1551 
1552  lduMatrix::debug = 1;
1553  }
1554  }
1555  else
1556  {
1557  FatalErrorIn(args.executable())
1558  << "FsiFoam:Step: Wrong temporal (ddt) scheme for solid solver. "
1559  << "Possible schemes are: "
1560  << fv::EulerDdtScheme<vector>::typeName << " and "
1561  << fv::backwardDdtScheme<vector>::typeName
1562  << abort(FatalError);
1563  }
1564 
1565  //# include "calcFsiResidual.H"
1566  {
1567  const vectorField& solidPatchDisplacement =
1568  DU.boundaryField()[solidPatchID];
1569 
1570  vectorField solidZoneDisplacement
1571  (
1572  structuresMesh.faceZones()[solidZoneID]().size(),
1573  vector::zero
1574  );
1575 
1576  const label solidPatchStart =
1577  structuresMesh.boundaryMesh()[solidPatchID].start();
1578 
1579  forAll(solidPatchDisplacement, i)
1580  {
1581  solidZoneDisplacement
1582  [
1583  structuresMesh.faceZones()[solidZoneID]
1584  .whichFace(solidPatchStart + i)
1585  ] =
1586  solidPatchDisplacement[i];
1587  }
1588 
1589  // Parallel data exchange: collect displacement field on all processors
1590  reduce(solidZoneDisplacement, sumOp<vectorField>());
1591 
1592  vectorField fluidZoneDisplacement =
1593  interpolatorSolidFluid.faceInterpolate
1594  (
1595  solidZoneDisplacement
1596  );
1597 
1598  vectorField fluidPatchDisplacement
1599  (
1600  fluidsMesh.boundary()[fluidPatchID].size(),
1601  vector::zero
1602  );
1603 
1604  const label fluidPatchStart =
1605  fluidsMesh.boundaryMesh()[fluidPatchID].start();
1606 
1607  forAll(fluidPatchDisplacement, i)
1608  {
1609  fluidPatchDisplacement[i] =
1610  fluidZoneDisplacement
1611  [
1612  fluidsMesh.faceZones()[fluidZoneID].whichFace(fluidPatchStart + i)
1613  ];
1614  }
1615 
1616  primitivePatchInterpolation fluidPatchInterpolator
1617  (
1618  fluidsMesh.boundaryMesh()[fluidPatchID]
1619  );
1620 
1621  solidPatchPointsDispl =
1622  fluidPatchInterpolator.faceToPointInterpolate
1623  (
1624  fluidPatchDisplacement
1625  );
1626 
1627  fsiResidualOld = fsiResidual;
1628 
1629  fsiResidual = solidPatchPointsDispl - fluidPatchPointsDispl;
1630 
1631  /*
1632  Info << "FsiFoam:Step: solidPatchPointsDispl = " << endl;
1633  Info << solidPatchPointsDispl << endl;
1634  Info << "FsiFoam:Step: fluidPatchPointsDispl = " << endl;
1635  Info << fluidPatchPointsDispl << endl;
1636  */
1637 
1638  // maxFsiResidual =
1639  // gMax
1640  // (
1641  // mag(fsiResidual)
1642  // /(mag(solidPatchPointsDispl) + SMALL)
1643  // );
1644 
1645  // Info << "Maximal fsi residual: " << maxFsiResidual << endl;
1646 
1647  fsiResidualNorm = ::sqrt(gSum(magSqr(fsiResidual)));
1648 
1649  if (outerCorr == 1)
1650  {
1651  initialFsiResidualNorm = fsiResidualNorm;
1652  }
1653 
1654  fsiResidualNorm /= initialFsiResidualNorm + SMALL;
1655 
1656  Info << "FsiFoam:Step: Current fsi residual norm: " << fsiResidualNorm << endl;
1657  }
1658  }
1659  while
1660  (
1661  (fsiResidualNorm > outerCorrTolerance)
1662  && (outerCorr < nOuterCorr)
1663  );
1664 
1665  Vs += DV;
1666 
1667  //# include "rotateSolidFields.H"
1668  {
1669  Info << "FsiFoam:Step: Rotating fields" << endl;
1670 
1671  F = I + DF;
1672 
1673  //U += DU;
1674  Usolid += DU;
1675 
1676  //epsilon += DEpsilon;
1677 
1678  sigma += DSigma;
1679 
1680  //volTensorField Finv = inv(F);
1681 
1682  J = det(F);
1683 
1684  rho = rho/J;
1685 
1686  //epsilon = symm(Finv.T() & epsilon & Finv);
1687 
1688  sigma = 1/J * symm(F & sigma & F.T());
1689  }
1690  //# include "moveSolidMeshLeastSquares.H"
1691  //--------------------------------------------------//
1692  //- move mesh
1693  //--------------------------------------------------//
1694  if(min(J.internalField()) > 0)
1695  {
1696  Info << "FsiFoam:Step: Moving mesh using least squares interpolation" << endl;
1697 
1698  leastSquaresVolPointInterpolation pointInterpolation(structuresMesh);
1699 
1700  // Create point mesh
1701  pointMesh pMesh(structuresMesh);
1702 
1703  wordList types
1704  (
1705  pMesh.boundary().size(),
1706  calculatedFvPatchVectorField::typeName
1707  );
1708 
1709  pointVectorField pointDU
1710  (
1711  IOobject
1712  (
1713  "pointDU",
1714  runTime.timeName(),
1715  structuresMesh
1716  ),
1717  pMesh,
1718  dimensionedVector("zero", dimLength, vector::zero),
1719  types
1720  );
1721 
1722  pointInterpolation.interpolate(DU, pointDU);
1723 
1724  const vectorField& pointDUI =
1725  pointDU.internalField();
1726 
1727  //- Move mesh
1728  vectorField newPoints = structuresMesh.allPoints();
1729 
1730  forAll (pointDUI, pointI)
1731  {
1732  newPoints[pointI] += pointDUI[pointI];
1733  }
1734 
1735  // Correct symmetryPlane points
1736 
1737  forAll(structuresMesh.boundaryMesh(), patchI)
1738  {
1739  if (isA<symmetryPolyPatch>(structuresMesh.boundaryMesh()[patchI]))
1740  {
1741  const labelList& meshPoints =
1742  structuresMesh.boundaryMesh()[patchI].meshPoints();
1743 
1744  vector avgN =
1745  gAverage(structuresMesh.boundaryMesh()[patchI].pointNormals());
1746 
1747  vector i(1, 0, 0);
1748  vector j(0, 1, 0);
1749  vector k(0, 0, 1);
1750 
1751  if (mag(avgN&i) > 0.95)
1752  {
1753  forAll(meshPoints, pI)
1754  {
1755  newPoints[meshPoints[pI]].x() = 0;
1756  }
1757  }
1758  else if (mag(avgN&j) > 0.95)
1759  {
1760  forAll(meshPoints, pI)
1761  {
1762  newPoints[meshPoints[pI]].y() = 0;
1763  }
1764  }
1765  else if (mag(avgN&k) > 0.95)
1766  {
1767  forAll(meshPoints, pI)
1768  {
1769  newPoints[meshPoints[pI]].z() = 0;
1770  }
1771  }
1772  }
1773  }
1774 
1775  //# include "calcUnusedNewPoints.H"
1776  forAll(globalFaceZones, zoneI)
1777  {
1778  const label curZoneID = globalFaceZones[zoneI];
1779 
1780  const labelList& curMap =
1781  globalToLocalFaceZonePointMap[zoneI];
1782 
1783  const labelList& curZoneMeshPoints =
1784  structuresMesh.faceZones()[curZoneID]().meshPoints();
1785 
1786  vectorField curGlobalZonePointDispl
1787  (
1788  curZoneMeshPoints.size(),
1789  vector::zero
1790  );
1791 
1792  //-Inter-proc points are shared by multiple procs
1793  // pointNumProc is the number of procs which a point lies on
1794  scalarField pointNumProcs(curZoneMeshPoints.size(), 0);
1795 
1796  forAll(curGlobalZonePointDispl, globalPointI)
1797  {
1798  label localPoint = curMap[globalPointI];
1799 
1800  if(curZoneMeshPoints[localPoint] < structuresMesh.nPoints())
1801  {
1802  label procPoint = curZoneMeshPoints[localPoint];
1803 
1804  curGlobalZonePointDispl[globalPointI] = pointDUI[procPoint];
1805 
1806  pointNumProcs[globalPointI] = 1;
1807  }
1808  }
1809 
1810  if (Pstream::parRun())
1811  {
1812  reduce(curGlobalZonePointDispl, sumOp<vectorField>());
1813  reduce(pointNumProcs, sumOp<scalarField>());
1814 
1815  //- now average the displacement between all procs
1816  curGlobalZonePointDispl /= pointNumProcs;
1817  }
1818 
1819  //- The curZonePointsDisplGlobal now contains the correct face zone
1820  // displacement in a global master processor order, now convert them
1821  // back into the local proc order
1822 
1823  vectorField curZonePointDispl
1824  (
1825  curZoneMeshPoints.size(),
1826  vector::zero
1827  );
1828 
1829  forAll(curGlobalZonePointDispl, globalPointI)
1830  {
1831  label localPoint = curMap[globalPointI];
1832 
1833  curZonePointDispl[localPoint] =
1834  curGlobalZonePointDispl[globalPointI];
1835  }
1836 
1837  forAll(curZonePointDispl, pointI)
1838  {
1839  // unused points
1840  if (curZoneMeshPoints[pointI] >= structuresMesh.nPoints())
1841  {
1842  newPoints[curZoneMeshPoints[pointI]] +=
1843  curZonePointDispl[pointI];
1844  }
1845  }
1846  }
1847 
1848  twoDPointCorrector twoDCorrector(structuresMesh);
1849  twoDCorrector.correctPoints(newPoints);
1850  structuresMesh.movePoints(newPoints);
1851  structuresMesh.V00();
1852  structuresMesh.moving(false);
1853  }
1854  else
1855  {
1856  FatalErrorIn(args.executable())
1857  << "FsiFoam:Step: Negative Jacobian"
1858  << exit(FatalError);
1859  }
1860  //# include "calculateStress.H"
1861  if (runTime.outputTime())
1862  {
1863  volScalarField sigmaEq
1864  (
1865  IOobject
1866  (
1867  "sigmaEq",
1868  runTime.timeName(),
1869  structuresMesh,
1870  IOobject::NO_READ,
1871  IOobject::AUTO_WRITE
1872  ),
1873  sqrt((3.0/2.0)*magSqr(dev(sigma)))
1874  );
1875 
1876  Info<< "FsiFoam:Step: Max sigmaEq = " << max(sigmaEq).value()
1877  << endl;
1878 
1879  runTime.write();
1880  }
1881 
1882  //# include "calculateLiftAndDrag.H"
1883 
1884  // Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
1885  // << " ClockTime = " << runTime.elapsedClockTime() << " s"
1886  // << endl << endl;
1887  }
1888  return 0;
1889 };
1890 
1891 int fsifoam_module::StepFluidAlone(){
1892 
1893  Foam::argList &args(this->ArgList());
1894  Foam::Time &runTime(this->RunTime());
1895  dynamicFvMesh &fluidsMesh(this->FluidMesh());
1896  // IOdictionary &transportProperties(this->TransportProperties());
1897  dimensionedScalar &nu(this->nu());
1898  dimensionedScalar &rhoFluid(this->rhoFluid());
1899  volScalarField &p(this->p());
1900  volVectorField &U(this->U());
1901  surfaceScalarField &phi(this->phi());
1902 
1903 
1904  // // IOdictionary &couplingProperties(this->CouplingProperties());
1905 
1906  label fluidPatchID(this->FluidPatchID());
1907  label fluidZoneID(this->FluidZoneID());
1908  label solidZoneID(this->SolidZoneID());
1909 
1910  bool feMotionSolver(this->FEMotion());
1911  bool fvMotionSolver(this->FVMotion());
1912  vectorField &accumulatedFluidInterfaceDisplacement(this->AccumulatedFluidInterfaceDisplacements());
1913  // Masoud : Checking AccumulatedFluidInterfaceDisplacement
1914  Info << "FsiFoam:StepFluidAlone: Checking AccumulatedFluidInterfaceDisplacementi = " << endl;
1915  Info << accumulatedFluidInterfaceDisplacement;
1916  // Masoud: End
1917  solidTractionFvPatchVectorField &tForce(this->tForce());
1918 
1919  zoneToZoneInterpolation &interpolatorFluidSolid(this->interpFluidSolid());
1920  zoneToZoneInterpolation &interpolatorSolidFluid(this->interpSolidFluid());
1921 
1922  scalar &sumLocalContErr(this->LocalContErr());
1923  scalar &globalContErr(this->GlobalContErr());
1924  scalar &cumulativeContErr(this->CumulativeContErr());
1925 
1926  labelList &globalFaceZones(this->GlobalFaceZones());
1927  labelListList &globalToLocalFaceZonePointMap(this->GlobalToLocalFaceZonePointMap());
1928  runTime++;
1929  if(!runTime.end()){
1930 
1931  this->ReadPISOControls();
1932  this->ReadFSIControls();
1933  word couplingScheme(this->CouplingScheme());
1934  // label &outerCorr(this->OuterCorr());
1935  scalar &fsiRelaxationFactor(this->FSIRelaxationFactor());
1936  // scalar &fsiRelaxationFactorMin(this->FSIRelaxationFactorMin());
1937  Switch fsi(this->FSIEnabled());
1938  int &nCorrPISO(this->NCorrPISO());
1939  int &nNonOrthCorr(this->NNonOrthCorr());
1940  scalar &interfaceDeformationLimit(this->InterfaceDeformationLimit());
1941  scalar &outerCorrTolerance(this->OuterCorrTolerance());
1942  int &nOuterCorr(this->NOuterCorr());
1943 
1944  vectorField fluidPatchPointsDispl
1945  (
1946  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
1947  vector::zero
1948  );
1949 
1950  vectorField fluidPatchPointsDisplOld
1951  (
1952  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
1953  vector::zero
1954  );
1955 
1956  vectorField solidPatchPointsDispl
1957  (
1958  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
1959  vector::zero
1960  );
1961 
1962  vectorField fsiResidual
1963  (
1964  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
1965  vector::zero
1966  );
1967 
1968  vectorField fsiResidualOld
1969  (
1970  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
1971  vector::zero
1972  );
1973 
1974 
1975  scalar initialFsiResidualNorm = 0;
1976  scalar fsiResidualNorm = 0;
1977  label outerCorr=0;
1978 
1979  do
1980  {
1981  outerCorr++;
1982  Info << "FsiFoam:StepFluidAlone: outerCorr = " << outerCorr << endl;
1983 
1984  //# include "setInterfaceDisplacement.H"
1985  Info << "\nFsiFoam:StepFluidAlone: Time = " << runTime.timeName()
1986  << ", iteration: " << outerCorr << endl;
1987 
1988  if (outerCorr < 3 || couplingScheme == "FixedRelaxation")
1989  {
1990  Info << "FsiFoam:StepFluidAlone: Current fsi under-relaxation factor: "
1991  << fsiRelaxationFactor << endl;
1992 
1993  fluidPatchPointsDisplOld = fluidPatchPointsDispl;
1994 
1995  fluidPatchPointsDispl += fsiRelaxationFactor*fsiResidual;
1996  }
1997  else
1998  {
1999  if (couplingScheme == "Aitken")
2000  {
2001  fsiRelaxationFactor =
2002  -fsiRelaxationFactor
2003  *(
2004  gSum
2005  (
2006  fsiResidualOld
2007  &(fsiResidual - fsiResidualOld)
2008  )
2009  /(
2010  gSum
2011  (
2012  (fsiResidual - fsiResidualOld)
2013  &(fsiResidual - fsiResidualOld)
2014  )
2015  )
2016  );
2017 
2018  fsiRelaxationFactor = mag(fsiRelaxationFactor);
2019 
2020  Info << "FsiFoam:StepFluidAlone: Current fsi under-relaxation factor (Aitken): "
2021  << fsiRelaxationFactor << endl;
2022 
2023  fluidPatchPointsDisplOld = fluidPatchPointsDispl;
2024 
2025  fluidPatchPointsDispl +=
2026  fsiRelaxationFactor*fsiResidual;
2027  }
2028  }
2029 
2030  {
2031 
2032  const vectorField& n =
2033  fluidsMesh.boundaryMesh()[fluidPatchID].pointNormals();
2034 
2035  primitivePatchInterpolation patchInterpolator
2036  (
2037  fluidsMesh.boundaryMesh()[fluidPatchID]
2038  );
2039 
2040  scalarField pointDeltaCoeffs =
2041  patchInterpolator.faceToPointInterpolate
2042  (
2043  fluidsMesh.boundary()[fluidPatchID].deltaCoeffs()
2044  );
2045 
2046  scalar delta =
2047  gMax
2048  (
2049  mag
2050  (
2051  n
2052  & (
2053  accumulatedFluidInterfaceDisplacement
2054  + fluidPatchPointsDispl
2055  - fluidPatchPointsDisplOld
2056  )
2057  )
2058  *pointDeltaCoeffs
2059  );
2060 
2061  Info << "FsiFoam:StepFluidAlone: Maximal accumulated displacement of interface points: "
2062  << delta << endl;
2063 
2064  if(delta < interfaceDeformationLimit)
2065  {
2066  // Move only interface points
2067  // Masoud
2068  Info << "FsiFoam:StepFluidAlone: Moving only interface...";
2069  // Masoud : End
2070  pointField newPoints = fluidsMesh.allPoints();
2071 
2072  const labelList& meshPoints =
2073  fluidsMesh.boundaryMesh()[fluidPatchID].meshPoints();
2074 
2075  forAll (fluidPatchPointsDispl, pointI)
2076  {
2077  //Info << "1. meshPoints[" << pointI << "] = " << meshPoints[pointI]
2078  //<< endl;
2079 
2080  newPoints[meshPoints[pointI]] +=
2081  fluidPatchPointsDispl[pointI]
2082  - fluidPatchPointsDisplOld[pointI];
2083  }
2084 
2085  twoDPointCorrector twoDCorrector(fluidsMesh);
2086 
2087  twoDCorrector.correctPoints(newPoints);
2088 
2089  fluidsMesh.movePoints(newPoints);
2090 
2091  // Accumulate interface points displacement
2092  accumulatedFluidInterfaceDisplacement +=
2093  fluidPatchPointsDispl
2094  - fluidPatchPointsDisplOld;
2095  }
2096  else
2097  {
2098  // Move whole fluid mesh
2099  // Masoud
2100  Info << "FsiFoam:StepFluidAlone: Moving the whole mesh...";
2101  // Masoud: End
2102  pointField newPoints = fluidsMesh.allPoints();
2103 
2104  const labelList& meshPoints =
2105  fluidsMesh.boundaryMesh()[fluidPatchID].meshPoints();
2106 
2107  forAll (accumulatedFluidInterfaceDisplacement, pointI)
2108  {
2109  //Info << "2. meshPoints[" << pointI << "] = " << meshPoints[pointI]
2110  //<< endl;
2111  newPoints[meshPoints[pointI]] -=
2112  accumulatedFluidInterfaceDisplacement[pointI];
2113  }
2114 
2115  twoDPointCorrector twoDCorrector(fluidsMesh);
2116 
2117  twoDCorrector.correctPoints(newPoints);
2118 
2119  fluidsMesh.movePoints(newPoints);
2120 
2121  accumulatedFluidInterfaceDisplacement +=
2122  fluidPatchPointsDispl
2123  - fluidPatchPointsDisplOld;
2124 
2125 
2126  if (feMotionSolver)
2127  {
2128  tetPointVectorField& motionU =
2129  const_cast<tetPointVectorField&>
2130  (
2131  fluidsMesh.objectRegistry:: lookupObject<tetPointVectorField>
2132  (
2133  "motionU"
2134  )
2135  );
2136 
2137  fixedValueTetPolyPatchVectorField& motionUFluidPatch =
2138  refCast<fixedValueTetPolyPatchVectorField>
2139  (
2140  motionU.boundaryField()[fluidPatchID]
2141  );
2142 
2143  tetPolyPatchInterpolation tppi
2144  (
2145  refCast<const faceTetPolyPatch>(motionUFluidPatch.patch())
2146  );
2147 
2148  motionUFluidPatch ==
2149  tppi.pointToPointInterpolate
2150  (
2151  accumulatedFluidInterfaceDisplacement
2152  /runTime.deltaT().value()
2153  );
2154  }
2155  else if (fvMotionSolver)
2156  {
2157  pointVectorField& motionU =
2158  const_cast<pointVectorField&>
2159  (
2160  fluidsMesh.objectRegistry:: lookupObject<pointVectorField>
2161  (
2162  "pointMotionU"
2163  )
2164  );
2165 
2166  fixedValuePointPatchVectorField& motionUFluidPatch =
2167  refCast<fixedValuePointPatchVectorField>
2168  (
2169  motionU.boundaryField()[fluidPatchID]
2170  );
2171 
2172  motionUFluidPatch ==
2173  accumulatedFluidInterfaceDisplacement
2174  /runTime.deltaT().value();
2175  }
2176  else
2177  {
2178  FatalErrorIn(args.executable())
2179  << "FsiFoam:StepFluidAlone: Problem with mesh motion solver selection"
2180  << abort(FatalError);
2181  }
2182 
2183  fluidsMesh.update();
2184 
2185  accumulatedFluidInterfaceDisplacement = vector::zero;
2186  }
2187  }
2188 
2189  if(fluidsMesh.moving())
2190  {
2191  // Make the fluxes relative
2192  phi -= fvc::meshPhi(U);
2193  }
2194 
2195 
2196  scalar CoNum = 0.0;
2197  scalar meanCoNum = 0.0;
2198  scalar velMag = 0.0;
2199 
2200  if (fluidsMesh.nInternalFaces())
2201  {
2202  surfaceScalarField magPhi = mag(phi);
2203 
2204  surfaceScalarField SfUfbyDelta =
2205  fluidsMesh.surfaceInterpolation::deltaCoeffs()*magPhi;
2206 
2207  const scalar deltaT = runTime.deltaT().value();
2208 
2209  CoNum = max(SfUfbyDelta/fluidsMesh.magSf()).value()*deltaT;
2210 
2211  meanCoNum = (sum(SfUfbyDelta)/sum(fluidsMesh.magSf())).value()*deltaT;
2212 
2213  velMag = max(magPhi/fluidsMesh.magSf()).value();
2214  }
2215 
2216  Info<< "FsiFoam:StepFluidAlone: Courant Number mean: " << meanCoNum
2217  << " max: " << CoNum
2218  << " velocity magnitude: " << velMag
2219  << endl;
2220 
2221  fvVectorMatrix UEqn
2222  (
2223  fvm::ddt(U)
2224  + fvm::div(phi, U)
2225  - fvm::laplacian(nu, U)
2226  );
2227 
2228  solve(UEqn == -fvc::grad(p));
2229 
2230  // --- PISO loop
2231  volScalarField rUA = 1.0/UEqn.A();
2232  // Masoud
2233  Info << "FsiFoam:StepFluidAlone: Performing " << nCorrPISO << " Pressure corrections." << endl;
2234  // Masoud End
2235 
2236  for (int corr=0; corr<nCorrPISO; corr++)
2237  {
2238  U = rUA*UEqn.H();
2239  phi = (fvc::interpolate(U) & fluidsMesh.Sf());
2240 
2241  adjustPhi(phi, U, p);
2242 
2243  for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
2244  {
2245  fvScalarMatrix pEqn
2246  (
2247  fvm::laplacian(rUA, p)
2248  == fvc::div(phi)
2249  );
2250 
2251  pEqn.setReference(pRefCell, pRefValue);
2252  pEqn.solve();
2253 
2254  // Added by Masoud
2255  //Info << "p = " << (this-> p()) << endl;
2256 
2257  if (nonOrth == nNonOrthCorr)
2258  {
2259  phi -= pEqn.flux();
2260  }
2261  }
2262 
2263  {
2264  volScalarField contErr = fvc::div(phi);
2265 
2266  sumLocalContErr = runTime.deltaT().value()*
2267  mag(contErr)().weightedAverage(fluidsMesh.V()).value();
2268 
2269  globalContErr = runTime.deltaT().value()*
2270  contErr.weightedAverage(fluidsMesh.V()).value();
2271 
2272  cumulativeContErr += globalContErr;
2273 
2274  Info<< "FsiFoam:StepFluidAlone: time step continuity errors : sum local = " << sumLocalContErr
2275  << ", global = " << globalContErr
2276  << ", cumulative = " << cumulativeContErr
2277  << endl;
2278  }
2279 
2280  U -= rUA*fvc::grad(p);
2281  U.correctBoundaryConditions();
2282  }
2283 
2284  {
2285  Info << "FsiFoam:StepFluidAlone: Not Setting traction on solid patch" << endl;
2286 
2287  {
2288 
2289  fsiResidualOld = fsiResidual;
2290  // Original: a fixed displacement delta will be used for the solid
2291  //this->UpdateFSISurface(solidPatchPointsDispl);
2292  // Oridinal : End
2293 
2294  // Masoud : in the initial step displacement delta will be added
2295  // and for the rest of the process, it will be consider-
2296  // as zeros vector.
2297  if (outerCorr == 1)
2298  {
2299  this->UpdateFSISurface(solidPatchPointsDispl);
2300  } else {
2301  solidPatchPointsDispl = vector::zero;
2302  //fluidPatchPointsDispl = vector::zero;
2303  }
2304  // Masoud : End
2305 
2306  // Masoud
2307  Info << "FsiFoam:StepFluidAlone: solidPatchPointsDispl = " << solidPatchPointsDispl << endl;
2308  //Masoud end
2309 
2310  // Original (solidPatchPointsDispl won't change by iteration
2311  // therefore huge pressure build-up happens)
2312  fsiResidual = (solidPatchPointsDispl - fluidPatchPointsDispl);
2313  // Original end
2314 
2315 
2316  fsiResidualNorm = ::sqrt(gSum(magSqr(fsiResidual)));
2317 
2318  if (outerCorr == 1)
2319  {
2320  initialFsiResidualNorm = fsiResidualNorm;
2321  }
2322 
2323  fsiResidualNorm /= initialFsiResidualNorm + SMALL;
2324 
2325  Info << "FsiFoam:StepFluidAlone: Current fsi residual norm: " << fsiResidualNorm << endl;
2326  }
2327  }
2328  }
2329  while
2330  (
2331  (fsiResidualNorm > outerCorrTolerance)
2332  && (outerCorr < nOuterCorr) //Masoud changed from outerCorr < nOuterCorr
2333  );
2334 
2335 
2336  if (runTime.outputTime())
2337  {
2338  runTime.write();
2339  }
2340  }
2341  return(0);
2342 }
2343 
2344 // Masoud: This is simple iteration-less scheme
2345 int fsifoam_module::StepFluidNonItr(){
2346 
2347  Foam::argList &args(this->ArgList());
2348  Foam::Time &runTime(this->RunTime());
2349  dynamicFvMesh &fluidsMesh(this->FluidMesh());
2350  // IOdictionary &transportProperties(this->TransportProperties());
2351  dimensionedScalar &nu(this->nu());
2352  dimensionedScalar &rhoFluid(this->rhoFluid());
2353  volScalarField &p(this->p());
2354  volVectorField &U(this->U());
2355  surfaceScalarField &phi(this->phi());
2356 
2357 
2358  // // IOdictionary &couplingProperties(this->CouplingProperties());
2359 
2360  label fluidPatchID(this->FluidPatchID());
2361  label fluidZoneID(this->FluidZoneID());
2362  label solidZoneID(this->SolidZoneID());
2363 
2364  bool feMotionSolver(this->FEMotion());
2365  bool fvMotionSolver(this->FVMotion());
2366  vectorField &accumulatedFluidInterfaceDisplacement(this->AccumulatedFluidInterfaceDisplacements());
2367  //Info << "Checking AccumulatedFluidInterfaceDisplacementi = " << endl;
2368  //Info << accumulatedFluidInterfaceDisplacement;
2369  solidTractionFvPatchVectorField &tForce(this->tForce());
2370 
2371  zoneToZoneInterpolation &interpolatorFluidSolid(this->interpFluidSolid());
2372  zoneToZoneInterpolation &interpolatorSolidFluid(this->interpSolidFluid());
2373 
2374  scalar &sumLocalContErr(this->LocalContErr());
2375  scalar &globalContErr(this->GlobalContErr());
2376  scalar &cumulativeContErr(this->CumulativeContErr());
2377 
2378  labelList &globalFaceZones(this->GlobalFaceZones());
2379  labelListList &globalToLocalFaceZonePointMap(this->GlobalToLocalFaceZonePointMap());
2380  runTime++;
2381  if(!runTime.end()){
2382 
2383  this->ReadPISOControls();
2384  this->ReadFSIControls();
2385  word couplingScheme(this->CouplingScheme());
2386  // label &outerCorr(this->OuterCorr());
2387  scalar &fsiRelaxationFactor(this->FSIRelaxationFactor());
2388  // scalar &fsiRelaxationFactorMin(this->FSIRelaxationFactorMin());
2389  Switch fsi(this->FSIEnabled());
2390  int &nCorrPISO(this->NCorrPISO());
2391  int &nNonOrthCorr(this->NNonOrthCorr());
2392  scalar &interfaceDeformationLimit(this->InterfaceDeformationLimit());
2393  scalar &outerCorrTolerance(this->OuterCorrTolerance());
2394  int &nOuterCorr(this->NOuterCorr());
2395 
2396  vectorField fluidPatchPointsDispl
2397  (
2398  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
2399  vector::zero
2400  );
2401 
2402  vectorField fluidPatchPointsDisplOld
2403  (
2404  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
2405  vector::zero
2406  );
2407 
2408  vectorField solidPatchPointsDispl
2409  (
2410  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
2411  vector::zero
2412  );
2413 
2414  vectorField fsiResidual
2415  (
2416  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
2417  vector::zero
2418  );
2419 
2420  vectorField fsiResidualOld
2421  (
2422  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
2423  vector::zero
2424  );
2425 
2426 
2427  scalar initialFsiResidualNorm = 0;
2428  scalar fsiResidualNorm = 0;
2429  label outerCorr=0;
2430 
2431  do
2432  {
2433  outerCorr++;
2434  Info << "FsiFoam:StepFluidNonItr: outerCorr = " << outerCorr << endl;
2435 
2436  //# include "setInterfaceDisplacement.H"
2437  Info << "FsiFoam:StepFluidNonItr: Time = " << runTime.timeName()
2438  << ", iteration: " << outerCorr << endl;
2439  // Updating fluid patch displacement
2440  this->UpdateFSISurface(fluidPatchPointsDispl);
2441  //Info << " Updating fluid interface coordinates with : " << endl;
2442  //Info << fluidPatchPointsDispl << endl;
2443  {
2444 
2445  const vectorField& n =
2446  fluidsMesh.boundaryMesh()[fluidPatchID].pointNormals();
2447 
2448  primitivePatchInterpolation patchInterpolator
2449  (
2450  fluidsMesh.boundaryMesh()[fluidPatchID]
2451  );
2452 
2453  scalarField pointDeltaCoeffs =
2454  patchInterpolator.faceToPointInterpolate
2455  (
2456  fluidsMesh.boundary()[fluidPatchID].deltaCoeffs()
2457  );
2458 
2459  scalar delta =
2460  gMax
2461  (
2462  mag
2463  (
2464  n
2465  & (
2466  accumulatedFluidInterfaceDisplacement
2467  + fluidPatchPointsDispl
2468  )
2469  )
2470  *pointDeltaCoeffs
2471  );
2472 
2473  Info << "FsiFoam:StepFluidNonItr: Maximal accumulated displacement of interface points: "
2474  << delta << endl;
2475 
2476  if(delta < interfaceDeformationLimit)
2477  {
2478  // Move only interface points
2479  Info << "FsiFoam:StepFluidNonItr: Moving only interface...";
2480  pointField newPoints = fluidsMesh.allPoints();
2481 
2482  const labelList& meshPoints =
2483  fluidsMesh.boundaryMesh()[fluidPatchID].meshPoints();
2484 
2485  forAll (fluidPatchPointsDispl, pointI)
2486  {
2487  //Info << "1. meshPoints[" << pointI << "] = " << meshPoints[pointI]
2488  //<< endl;
2489 
2490  newPoints[meshPoints[pointI]] +=
2491  fluidPatchPointsDispl[pointI];
2492  }
2493 
2494  twoDPointCorrector twoDCorrector(fluidsMesh);
2495 
2496  twoDCorrector.correctPoints(newPoints);
2497 
2498  fluidsMesh.movePoints(newPoints);
2499 
2500  // Accumulate interface points displacement
2501  accumulatedFluidInterfaceDisplacement +=
2502  fluidPatchPointsDispl;
2503  }
2504  else
2505  {
2506  // Move whole fluid mesh
2507  Info << "FsiFoam:StepFluidNonItr: Moving the whole mesh...";
2508  pointField newPoints = fluidsMesh.allPoints();
2509 
2510  const labelList& meshPoints =
2511  fluidsMesh.boundaryMesh()[fluidPatchID].meshPoints();
2512 
2513  forAll (accumulatedFluidInterfaceDisplacement, pointI)
2514  {
2515  //Info << "2. meshPoints[" << pointI << "] = " << meshPoints[pointI]
2516  //<< endl;
2517  newPoints[meshPoints[pointI]] -=
2518  accumulatedFluidInterfaceDisplacement[pointI];
2519  }
2520 
2521  twoDPointCorrector twoDCorrector(fluidsMesh);
2522 
2523  twoDCorrector.correctPoints(newPoints);
2524 
2525  fluidsMesh.movePoints(newPoints);
2526 
2527  accumulatedFluidInterfaceDisplacement +=
2528  fluidPatchPointsDispl;
2529 
2530  if (feMotionSolver)
2531  {
2532  tetPointVectorField& motionU =
2533  const_cast<tetPointVectorField&>
2534  (
2535  fluidsMesh.objectRegistry:: lookupObject<tetPointVectorField>
2536  (
2537  "motionU"
2538  )
2539  );
2540 
2541  fixedValueTetPolyPatchVectorField& motionUFluidPatch =
2542  refCast<fixedValueTetPolyPatchVectorField>
2543  (
2544  motionU.boundaryField()[fluidPatchID]
2545  );
2546 
2547  tetPolyPatchInterpolation tppi
2548  (
2549  refCast<const faceTetPolyPatch>(motionUFluidPatch.patch())
2550  );
2551 
2552  motionUFluidPatch ==
2553  tppi.pointToPointInterpolate
2554  (
2555  accumulatedFluidInterfaceDisplacement
2556  /runTime.deltaT().value()
2557  );
2558  }
2559  else if (fvMotionSolver)
2560  {
2561  pointVectorField& motionU =
2562  const_cast<pointVectorField&>
2563  (
2564  fluidsMesh.objectRegistry:: lookupObject<pointVectorField>
2565  (
2566  "pointMotionU"
2567  )
2568  );
2569 
2570  fixedValuePointPatchVectorField& motionUFluidPatch =
2571  refCast<fixedValuePointPatchVectorField>
2572  (
2573  motionU.boundaryField()[fluidPatchID]
2574  );
2575 
2576  motionUFluidPatch ==
2577  accumulatedFluidInterfaceDisplacement
2578  /runTime.deltaT().value();
2579  }
2580  else
2581  {
2582  FatalErrorIn(args.executable())
2583  << "FsiFoam:StepFluidNonItr: Problem with mesh motion solver selection"
2584  << abort(FatalError);
2585  }
2586 
2587 
2588  Info << "FsiFoam:StepFluidNonItr: runTime.deltaT() = "
2589  << runTime.deltaT().value() << endl;
2590 
2591 
2592  fluidsMesh.update();
2593 
2594  accumulatedFluidInterfaceDisplacement = vector::zero;
2595  }
2596  }
2597 
2598  if(fluidsMesh.moving())
2599  {
2600  // Make the fluxes relative
2601  Info << "FsiFoam:StepFluidNonItr: Compensating for mesh movment !" << endl;
2602  phi -= fvc::meshPhi(U);
2603  }
2604 
2605 
2606  scalar CoNum = 0.0;
2607  scalar meanCoNum = 0.0;
2608  scalar velMag = 0.0;
2609 
2610  if (fluidsMesh.nInternalFaces())
2611  {
2612  surfaceScalarField magPhi = mag(phi);
2613 
2614  surfaceScalarField SfUfbyDelta =
2615  fluidsMesh.surfaceInterpolation::deltaCoeffs()*magPhi;
2616 
2617  const scalar deltaT = runTime.deltaT().value();
2618 
2619  CoNum = max(SfUfbyDelta/fluidsMesh.magSf()).value()*deltaT;
2620 
2621  meanCoNum = (sum(SfUfbyDelta)/sum(fluidsMesh.magSf())).value()*deltaT;
2622 
2623  velMag = max(magPhi/fluidsMesh.magSf()).value();
2624  }
2625 
2626  Info<< "FsiFoam:StepFluidNonItr: Courant Number mean: " << meanCoNum
2627  << " max: " << CoNum
2628  << " velocity magnitude: " << velMag
2629  << endl;
2630 
2631  fvVectorMatrix UEqn
2632  (
2633  fvm::ddt(U)
2634  + fvm::div(phi, U)
2635  - fvm::laplacian(nu, U)
2636  );
2637 
2638  solve(UEqn == -fvc::grad(p));
2639 
2640  // --- PISO loop
2641  volScalarField rUA = 1.0/UEqn.A();
2642  //Info << "nCorrPISO = " << nCorrPISO << endl;
2643  //Info << "nNonOrthCorr = " << nNonOrthCorr << endl;
2644  for (int corr=0; corr<nCorrPISO; corr++)
2645  {
2646  U = rUA*UEqn.H();
2647  phi = (fvc::interpolate(U) & fluidsMesh.Sf());
2648 
2649  adjustPhi(phi, U, p);
2650 
2651  for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
2652  {
2653  fvScalarMatrix pEqn
2654  (
2655  fvm::laplacian(rUA, p)
2656  == fvc::div(phi)
2657  );
2658 
2659  pEqn.setReference(pRefCell, pRefValue);
2660  pEqn.solve();
2661 
2662  //Info << "p = " << p << endl;
2663 
2664  if (nonOrth == nNonOrthCorr)
2665  {
2666  phi -= pEqn.flux();
2667  }
2668  }
2669 
2670  {
2671  volScalarField contErr = fvc::div(phi);
2672 
2673  sumLocalContErr = runTime.deltaT().value()*
2674  mag(contErr)().weightedAverage(fluidsMesh.V()).value();
2675 
2676  globalContErr = runTime.deltaT().value()*
2677  contErr.weightedAverage(fluidsMesh.V()).value();
2678 
2679  cumulativeContErr += globalContErr;
2680 
2681  Info<< "FsiFoam:StepFluidNonItr: time step continuity errors : sum local = " << sumLocalContErr
2682  << ", global = " << globalContErr
2683  << ", cumulative = " << cumulativeContErr
2684  << endl;
2685  }
2686 
2687  U -= rUA*fvc::grad(p);
2688  U.correctBoundaryConditions();
2689  }
2690 
2691  }
2692  while
2693  ( outerCorr < 1 );
2694 
2695 
2696  if (runTime.outputTime())
2697  {
2698  runTime.write();
2699  }
2700  }
2701  return(0);
2702 }
2703 //Masoud End (StepFluidNonItr)
2704 
2705 
2706 //Masoud: StepFluidItr
2707 int fsifoam_module::StepFluidItr(){
2708 
2709  Foam::argList &args(this->ArgList());
2710  Foam::Time &runTime(this->RunTime());
2711  dynamicFvMesh &fluidsMesh(this->FluidMesh());
2712  // IOdictionary &transportProperties(this->TransportProperties());
2713  dimensionedScalar &nu(this->nu());
2714  dimensionedScalar &rhoFluid(this->rhoFluid());
2715  volScalarField &p(this->p());
2716  volVectorField &U(this->U());
2717  surfaceScalarField &phi(this->phi());
2718 
2719 
2720  // // IOdictionary &couplingProperties(this->CouplingProperties());
2721 
2722  label fluidPatchID(this->FluidPatchID());
2723  label fluidZoneID(this->FluidZoneID());
2724  label solidZoneID(this->SolidZoneID());
2725 
2726  bool feMotionSolver(this->FEMotion());
2727  bool fvMotionSolver(this->FVMotion());
2728  vectorField &accumulatedFluidInterfaceDisplacement(this->AccumulatedFluidInterfaceDisplacements());
2729  // Masoud : Checking AccumulatedFluidInterfaceDisplacement
2730  Info << "FsiFoam:StepFluidItr: Checking AccumulatedFluidInterfaceDisplacementi = " << endl;
2731  Info << accumulatedFluidInterfaceDisplacement;
2732  // Masoud: End
2733  solidTractionFvPatchVectorField &tForce(this->tForce());
2734 
2735  zoneToZoneInterpolation &interpolatorFluidSolid(this->interpFluidSolid());
2736  zoneToZoneInterpolation &interpolatorSolidFluid(this->interpSolidFluid());
2737 
2738  scalar &sumLocalContErr(this->LocalContErr());
2739  scalar &globalContErr(this->GlobalContErr());
2740  scalar &cumulativeContErr(this->CumulativeContErr());
2741 
2742  labelList &globalFaceZones(this->GlobalFaceZones());
2743  labelListList &globalToLocalFaceZonePointMap(this->GlobalToLocalFaceZonePointMap());
2744  runTime++;
2745  if(!runTime.end()){
2746 
2747  this->ReadPISOControls();
2748  this->ReadFSIControls();
2749  word couplingScheme(this->CouplingScheme());
2750  // label &outerCorr(this->OuterCorr());
2751  scalar &fsiRelaxationFactor(this->FSIRelaxationFactor());
2752  // scalar &fsiRelaxationFactorMin(this->FSIRelaxationFactorMin());
2753  Switch fsi(this->FSIEnabled());
2754  int &nCorrPISO(this->NCorrPISO());
2755  int &nNonOrthCorr(this->NNonOrthCorr());
2756  scalar &interfaceDeformationLimit(this->InterfaceDeformationLimit());
2757  scalar &outerCorrTolerance(this->OuterCorrTolerance());
2758  int &nOuterCorr(this->NOuterCorr());
2759 
2760  vectorField fluidPatchPointsDispl
2761  (
2762  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
2763  vector::zero
2764  );
2765 
2766  vectorField fluidPatchPointsDisplOld
2767  (
2768  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
2769  vector::zero
2770  );
2771 
2772  vectorField solidPatchPointsDispl
2773  (
2774  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
2775  vector::zero
2776  );
2777 
2778  vectorField fsiResidual
2779  (
2780  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
2781  vector::zero
2782  );
2783 
2784  vectorField fsiResidualOld
2785  (
2786  fluidsMesh.boundaryMesh()[fluidPatchID].nPoints(),
2787  vector::zero
2788  );
2789 
2790 
2791  scalar initialFsiResidualNorm = 0;
2792  scalar fsiResidualNorm = 0;
2793  label outerCorr=0;
2794 
2795  do
2796  {
2797  outerCorr++;
2798  Info << "FsiFoam:StepFluidItr: outerCorr = " << outerCorr << endl;
2799 
2800  //# include "setInterfaceDisplacement.H"
2801  Info << "\nFsiFoam:StepFluidItr: Time = " << runTime.timeName()
2802  << ", iteration: " << outerCorr << endl;
2803 
2804  if (outerCorr < 3 || couplingScheme == "FixedRelaxation")
2805  {
2806  Info << "FsiFoam:StepFluidItr: Current fsi under-relaxation factor: "
2807  << fsiRelaxationFactor << endl;
2808 
2809  fluidPatchPointsDisplOld = fluidPatchPointsDispl;
2810 
2811  fluidPatchPointsDispl += fsiRelaxationFactor*fsiResidual;
2812  }
2813  else
2814  {
2815  if (couplingScheme == "Aitken")
2816  {
2817  fsiRelaxationFactor =
2818  -fsiRelaxationFactor
2819  *(
2820  gSum
2821  (
2822  fsiResidualOld
2823  &(fsiResidual - fsiResidualOld)
2824  )
2825  /(
2826  gSum
2827  (
2828  (fsiResidual - fsiResidualOld)
2829  &(fsiResidual - fsiResidualOld)
2830  )
2831  )
2832  );
2833 
2834  fsiRelaxationFactor = mag(fsiRelaxationFactor);
2835 
2836  Info << "FsiFoam:StepFluidItr: Current fsi under-relaxation factor (Aitken): "
2837  << fsiRelaxationFactor << endl;
2838 
2839  fluidPatchPointsDisplOld = fluidPatchPointsDispl;
2840 
2841  fluidPatchPointsDispl +=
2842  fsiRelaxationFactor*fsiResidual;
2843  }
2844  }
2845 
2846  {
2847 
2848  const vectorField& n =
2849  fluidsMesh.boundaryMesh()[fluidPatchID].pointNormals();
2850 
2851  primitivePatchInterpolation patchInterpolator
2852  (
2853  fluidsMesh.boundaryMesh()[fluidPatchID]
2854  );
2855 
2856  scalarField pointDeltaCoeffs =
2857  patchInterpolator.faceToPointInterpolate
2858  (
2859  fluidsMesh.boundary()[fluidPatchID].deltaCoeffs()
2860  );
2861 
2862  scalar delta =
2863  gMax
2864  (
2865  mag
2866  (
2867  n
2868  & (
2869  accumulatedFluidInterfaceDisplacement
2870  + fluidPatchPointsDispl
2871  - fluidPatchPointsDisplOld
2872  )
2873  )
2874  *pointDeltaCoeffs
2875  );
2876 
2877  Info << "FsiFoam:StepFluidItr: Maximal accumulated displacement of interface points: "
2878  << delta << endl;
2879 
2880  if(delta < interfaceDeformationLimit)
2881  {
2882  // Move only interface points
2883  // Masoud
2884  Info << "FsiFoam:StepFluidItr: Moving only interface...";
2885  // Masoud : End
2886  pointField newPoints = fluidsMesh.allPoints();
2887 
2888  const labelList& meshPoints =
2889  fluidsMesh.boundaryMesh()[fluidPatchID].meshPoints();
2890 
2891  forAll (fluidPatchPointsDispl, pointI)
2892  {
2893  //Info << "1. meshPoints[" << pointI << "] = " << meshPoints[pointI]
2894  //<< endl;
2895 
2896  newPoints[meshPoints[pointI]] +=
2897  fluidPatchPointsDispl[pointI]
2898  - fluidPatchPointsDisplOld[pointI];
2899  }
2900 
2901  twoDPointCorrector twoDCorrector(fluidsMesh);
2902 
2903  twoDCorrector.correctPoints(newPoints);
2904 
2905  fluidsMesh.movePoints(newPoints);
2906 
2907  // Accumulate interface points displacement
2908  accumulatedFluidInterfaceDisplacement +=
2909  fluidPatchPointsDispl
2910  - fluidPatchPointsDisplOld;
2911  }
2912  else
2913  {
2914  // Move whole fluid mesh
2915  // Masoud
2916  Info << "FsiFoam:StepFluidItr: Moving the whole mesh...";
2917  // Masoud: End
2918  pointField newPoints = fluidsMesh.allPoints();
2919 
2920  const labelList& meshPoints =
2921  fluidsMesh.boundaryMesh()[fluidPatchID].meshPoints();
2922 
2923  forAll (accumulatedFluidInterfaceDisplacement, pointI)
2924  {
2925  //Info << "2. meshPoints[" << pointI << "] = " << meshPoints[pointI]
2926  //<< endl;
2927  newPoints[meshPoints[pointI]] -=
2928  accumulatedFluidInterfaceDisplacement[pointI];
2929  }
2930 
2931  twoDPointCorrector twoDCorrector(fluidsMesh);
2932 
2933  twoDCorrector.correctPoints(newPoints);
2934 
2935  fluidsMesh.movePoints(newPoints);
2936 
2937  accumulatedFluidInterfaceDisplacement +=
2938  fluidPatchPointsDispl
2939  - fluidPatchPointsDisplOld;
2940 
2941 
2942  if (feMotionSolver)
2943  {
2944  tetPointVectorField& motionU =
2945  const_cast<tetPointVectorField&>
2946  (
2947  fluidsMesh.objectRegistry:: lookupObject<tetPointVectorField>
2948  (
2949  "motionU"
2950  )
2951  );
2952 
2953  fixedValueTetPolyPatchVectorField& motionUFluidPatch =
2954  refCast<fixedValueTetPolyPatchVectorField>
2955  (
2956  motionU.boundaryField()[fluidPatchID]
2957  );
2958 
2959  tetPolyPatchInterpolation tppi
2960  (
2961  refCast<const faceTetPolyPatch>(motionUFluidPatch.patch())
2962  );
2963 
2964  motionUFluidPatch ==
2965  tppi.pointToPointInterpolate
2966  (
2967  accumulatedFluidInterfaceDisplacement
2968  /runTime.deltaT().value()
2969  );
2970  }
2971  else if (fvMotionSolver)
2972  {
2973  pointVectorField& motionU =
2974  const_cast<pointVectorField&>
2975  (
2976  fluidsMesh.objectRegistry:: lookupObject<pointVectorField>
2977  (
2978  "pointMotionU"
2979  )
2980  );
2981 
2982  fixedValuePointPatchVectorField& motionUFluidPatch =
2983  refCast<fixedValuePointPatchVectorField>
2984  (
2985  motionU.boundaryField()[fluidPatchID]
2986  );
2987 
2988  motionUFluidPatch ==
2989  accumulatedFluidInterfaceDisplacement
2990  /runTime.deltaT().value();
2991  }
2992  else
2993  {
2994  FatalErrorIn(args.executable())
2995  << "FsiFoam:StepFluidItr: Problem with mesh motion solver selection"
2996  << abort(FatalError);
2997  }
2998 
2999  fluidsMesh.update();
3000 
3001  accumulatedFluidInterfaceDisplacement = vector::zero;
3002  }
3003  }
3004 
3005  if(fluidsMesh.moving())
3006  {
3007  // Make the fluxes relative
3008  phi -= fvc::meshPhi(U);
3009  }
3010 
3011 
3012  scalar CoNum = 0.0;
3013  scalar meanCoNum = 0.0;
3014  scalar velMag = 0.0;
3015 
3016  if (fluidsMesh.nInternalFaces())
3017  {
3018  surfaceScalarField magPhi = mag(phi);
3019 
3020  surfaceScalarField SfUfbyDelta =
3021  fluidsMesh.surfaceInterpolation::deltaCoeffs()*magPhi;
3022 
3023  const scalar deltaT = runTime.deltaT().value();
3024 
3025  CoNum = max(SfUfbyDelta/fluidsMesh.magSf()).value()*deltaT;
3026 
3027  meanCoNum = (sum(SfUfbyDelta)/sum(fluidsMesh.magSf())).value()*deltaT;
3028 
3029  velMag = max(magPhi/fluidsMesh.magSf()).value();
3030  }
3031 
3032  Info<< "FsiFoam:StepFluidItr: Courant Number mean: " << meanCoNum
3033  << " max: " << CoNum
3034  << " velocity magnitude: " << velMag
3035  << endl;
3036 
3037  fvVectorMatrix UEqn
3038  (
3039  fvm::ddt(U)
3040  + fvm::div(phi, U)
3041  - fvm::laplacian(nu, U)
3042  );
3043 
3044  solve(UEqn == -fvc::grad(p));
3045 
3046  // --- PISO loop
3047  volScalarField rUA = 1.0/UEqn.A();
3048  // Masoud
3049  Info << "FsiFoam:StepFluidItr: Performing " << nCorrPISO << " Pressure corrections." << endl;
3050  // Masoud End
3051 
3052  for (int corr=0; corr<nCorrPISO; corr++)
3053  {
3054  U = rUA*UEqn.H();
3055  phi = (fvc::interpolate(U) & fluidsMesh.Sf());
3056 
3057  adjustPhi(phi, U, p);
3058 
3059  for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
3060  {
3061  fvScalarMatrix pEqn
3062  (
3063  fvm::laplacian(rUA, p)
3064  == fvc::div(phi)
3065  );
3066 
3067  pEqn.setReference(pRefCell, pRefValue);
3068  pEqn.solve();
3069 
3070  // Added by Masoud
3071  //Info << "p = " << (this-> p()) << endl;
3072 
3073  if (nonOrth == nNonOrthCorr)
3074  {
3075  phi -= pEqn.flux();
3076  }
3077  }
3078 
3079  {
3080  volScalarField contErr = fvc::div(phi);
3081 
3082  sumLocalContErr = runTime.deltaT().value()*
3083  mag(contErr)().weightedAverage(fluidsMesh.V()).value();
3084 
3085  globalContErr = runTime.deltaT().value()*
3086  contErr.weightedAverage(fluidsMesh.V()).value();
3087 
3088  cumulativeContErr += globalContErr;
3089 
3090  Info<< "FsiFoam:StepFluidItr: time step continuity errors : sum local = " << sumLocalContErr
3091  << ", global = " << globalContErr
3092  << ", cumulative = " << cumulativeContErr
3093  << endl;
3094  }
3095 
3096  U -= rUA*fvc::grad(p);
3097  U.correctBoundaryConditions();
3098  }
3099 
3100  {
3101  Info << "FsiFoam:StepFluidItr: Setting tractionis on solid patch" << endl;
3102 
3103  {
3104 
3105  fsiResidualOld = fsiResidual;
3106  this->UpdateFSISurface(solidPatchPointsDispl);
3107  Info << "FsiFoam:StepFluidItr: solidPatchPointsDispl = " << solidPatchPointsDispl << endl;
3108 
3109  fsiResidual = (solidPatchPointsDispl - fluidPatchPointsDispl);
3110  fsiResidualNorm = ::sqrt(gSum(magSqr(fsiResidual)));
3111 
3112  if (outerCorr == 1)
3113  {
3114  initialFsiResidualNorm = fsiResidualNorm;
3115  }
3116 
3117  fsiResidualNorm /= initialFsiResidualNorm + SMALL;
3118 
3119  Info << "FsiFoam:StepFluidItr: Current fsi residual norm: " << fsiResidualNorm << endl;
3120  }
3121  }
3122  }
3123  while
3124  (
3125  (fsiResidualNorm > outerCorrTolerance)
3126  && (outerCorr < nOuterCorr) //Masoud changed from outerCorr < nOuterCorr
3127  );
3128 
3129 
3130  if (runTime.outputTime())
3131  {
3132  runTime.write();
3133  }
3134  }
3135  return(0);
3136 }
3137 //Masoud End StepFluidItr
3138 
3139 int fsifoam_module::Dump(){
3140  Foam::Time &runTime(this->RunTime());
3141  fvMesh &structuresMesh(this->StructuresMesh());
3142  volSymmTensorField &sigma(this->sigma());
3143  volScalarField sigmaEq
3144  (
3145  IOobject
3146  (
3147  "sigmaEq",
3148  runTime.timeName(),
3149  structuresMesh,
3150  IOobject::NO_READ,
3151  IOobject::AUTO_WRITE
3152  ),
3153  sqrt((3.0/2.0)*magSqr(dev(sigma)))
3154  );
3155 
3156  Info<< "FsiFoam:Dump: Max sigmaEq = " << max(sigmaEq).value()
3157  << endl;
3158 
3159  runTime.write();
3160  return 0;
3161 };
3162 
3167 //static void fsifoam_module::Load(const std::string &name){
3168 void fsifoam_module::Load(const std::string &name){
3169 
3170  // anouncing default communicator
3171  MPI_Comm inComm;
3172  inComm = COM_get_default_communicator();
3173 
3174  int rank, size;
3175  MPI_Comm_rank(inComm, &rank);
3176  MPI_Comm_size(inComm, &size);
3177  Info << "FsiFoam:Load:Rank #"
3178  << rank
3179  << " on communicator "
3180  << inComm
3181  << " with "
3182  << size
3183  << " processes.\n";
3184  Info << "FsiFoam:Load:Rank #"
3185  << rank
3186  << " Loading FsiFoamModule with name "
3187  << name
3188  << "\n.";
3189 
3190 
3191 
3193  fsifoam_module *module_pointer = new fsifoam_module();
3194  COM_new_window(name, MPI_COMM_NULL);
3195  module_pointer->my_window_name = name;
3196  MPI_Comm_dup(inComm, &(module_pointer->my_window_comm));
3197  MPI_Comm_rank(module_pointer->my_window_comm, &(module_pointer->rank));
3198  MPI_Comm_size(module_pointer->my_window_comm, &(module_pointer->nproc));
3199  std::string global_name(name+".global");
3200  COM_new_dataitem(global_name.c_str(),'w',COM_VOID,1,"");
3201  COM_set_object(global_name.c_str(),0,module_pointer);
3202 
3203 
3205  std::vector<COM_Type> types(13,COM_INT);
3206 
3207  types[0] = COM_RAWDATA;
3208  types[2] = COM_VOID;
3209  COM_set_member_function( (name+".InitFoam").c_str(),
3210  (Member_func_ptr)(&fsifoam_module::InitFoam),
3211  global_name.c_str(), "biii", &types[0]);
3212 
3213  COM_set_member_function( (name+".RunFoam").c_str(),
3214  (Member_func_ptr)(&fsifoam_module::RunFoam),
3215  global_name.c_str(), "b", &types[0]);
3216 
3217 
3218  COM_set_member_function( (name+".StepFoam").c_str(),
3219  (Member_func_ptr)(&fsifoam_module::StepFoam),
3220  global_name.c_str(), "b", &types[0]);
3221 
3222  COM_set_member_function( (name+".StepFluid").c_str(),
3223  (Member_func_ptr)(&fsifoam_module::StepFluid),
3224  global_name.c_str(), "b", &types[0]);
3225 
3226  // registering communicator for this module to COM
3227  COM_new_dataitem( (name+".nproc").c_str(), 'w', COM_INT, 1, "");
3228  COM_set_size( (name+".nproc").c_str(), 0, 1);
3229  COM_set_array( (name+".nproc").c_str(), 0, &(module_pointer->nproc));
3230 /*
3231  types[1] = COM_DOUBLE;
3232  COM_set_member_function( (name+".ModifyEndTime").c_str(),
3233  (Member_func_ptr)(&fsifoam_module::ModifyEndTime),
3234  global_name.c_str(), "bi", &types[0]);
3235 */
3236 
3237 
3238  COM_window_init_done(name);
3239 }
3240 
3241 
3242 void fsifoam_module::InitFoam(int *pargc, void **pargv, int *verbIn)
3243 {
3244  int argc = *pargc;
3245  char** argv = (char**)(pargv);
3246  verbosity = *verbIn;
3247  //std::cout << "OF Module verbosity = " << verbosity << std::endl;
3248  //std::cout << "rank = " << rank << " size = " << nproc << std::endl;
3249 
3250  /*
3251  // register function status variables
3252  std::string initStatusName(my_window_name+".initStatus");
3253  COM_new_dataitem(initStatusName.c_str(),'w', COM_INT, 1, "");
3254  COM_set_array(initStatusName.c_str(), 11, &initStatus);
3255  */
3256 
3257  Initialize(argc, argv);
3258 
3260  // set up solution meta data
3261  Solution().Meta().AddField("time", 's', 1, 8, "s");
3262  Solution().Meta().AddField("endTime", 's', 1, 8, "s");
3263  Solution().Meta().AddField("initStatus", 's', 1, 4, "");
3264  Solution().Meta().AddField("runStatus", 's', 1, 4, "");
3265  Solution().Meta().AddField("dt", 's', 1, 8, "s");
3266  Solution().Meta().AddField("solidDisplacement", 'n', 3, 8, "m");
3267  Solution().Meta().AddField("pressure", 'c', 1, 8, "Pa");
3268  Solution().Meta().AddField("traction", 'c', 3, 8, "N");
3269  //Solution().Meta().AddField("nElmSurf", 'c', 1, 1, "");
3270 
3271  Mesh().nc.init(numPointsSurface, &surfaceCoordinates[0]);
3272 
3273  // this assumes that our elements are quads...we can generalize this
3274  Mesh().con.AddElements(numElementsSurface, 4, surfaceConnectivity);
3275 
3276  // create buffers for the actual data
3277  // only call this method if the process has a share
3278  // from the interface mesh
3279  CreateSoln();
3280 
3281  if (numPointsSurface==0 && numElementsSurface==0)
3282  std::cout << "Process " << Rank()
3283  << " has no share on interface data ... "
3284  << std::endl;
3285 
3286  //unsigned int nnodes = Mesh().nc.NNodes();
3287  //unsigned int nelem = Mesh().con.Nelem();
3288 
3289  // allocate our own local buffers for the data
3290 
3291  time.resize(1, -1.0);
3292  endTime.resize(1, -1.0);
3293  UpdateTime();
3294 
3295  initStatus.resize(1, -1000);
3296  runStatus.resize(1, -1000);
3297 
3298  surfacePressure.resize(numElementsSurface, -1);
3299  surfaceTraction.resize(3*numElementsSurface, -1);
3300  solidDisplacement.resize(3*numPointsSurface, 0.);
3301 
3302 
3303  // reset the buffers to be our own local buffers
3304  Solution().SetFieldBuffer("initStatus", initStatus);
3305  Solution().SetFieldBuffer("runStatus", runStatus);
3306  Solution().SetFieldBuffer("time", time);
3307  Solution().SetFieldBuffer("endTime", endTime);
3308  Solution().SetFieldBuffer("pressure", surfacePressure);
3309  Solution().SetFieldBuffer("traction", surfaceTraction);
3310  Solution().SetFieldBuffer("solidDisplacement", solidDisplacement);
3311  //Solution().SetFieldBuffer("nElmSurf", &numElementsSurface);
3312 
3313  // This always registers given solver to Pane ID = 101
3314  //SolverUtils::RegisterSolverInto(my_window_name.c_str(), *this);
3315 
3316 
3317  // This registers given solver to given Pane ID
3318  int paneId = rank + 200;
3319  // following two lines replaces
3320  // SolverUtils::RegisterSolver(my_window_name.c_str(), *this, paneId);
3321  // this function creates a duplicate window which is not acceptable
3322  SolverUtils::CreateDataItemsFromSolution(my_window_name.c_str(),Solution());
3323  if (numPointsSurface!=0 && numElementsSurface!=0)
3324  SolverUtils::AgentToPane(my_window_name.c_str(),paneId,*this);
3325 
3326 
3327  COM_window_init_done(my_window_name.c_str());
3328 
3329  initStatus[0] = 0;
3330  std::cout << "Rank " << Rank() << " finished " << std::endl;
3331 
3332  return;
3333 }
3334 
3335 
3340 //static void fsifoam_module::Unload(const std::string &name){
3341 void fsifoam_module::Unload(const std::string &name){
3342  std::cout << "FsiFoam:Unload: Unloading FsiFoamModule with name " << name
3343  << "." << std::endl;
3344  fsifoam_module *module_pointer = NULL;
3345  std::string global_name(name+".global");
3346  COM_get_object(global_name.c_str(),0,&module_pointer);
3347  COM_assertion_msg( module_pointer->validate_object()==0, "Invalid object");
3348  delete module_pointer;
3349  COM_delete_window(std::string(name));
3350 }
3351 
3355 void fsifoam_module::RunFoam(){
3356 
3357 
3358  Foam::Time &runTime(RunTime());
3359 
3360  Info << "\nFsiFoam:RunFoam: Starting time loop\n" << endl;
3361 
3362  while(!runTime.end()){
3363  Info << "FsiFoam:RunFoam: Time = " << runTime.timeName() << nl << endl;
3364  Step();
3365  Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
3366  << " ClockTime = " << runTime.elapsedClockTime() << " s"
3367  << endl << endl;
3368  if (runTime.outputTime())
3369  Dump();
3370  }
3371 
3372  Info<< "FsiFoam:RunFoam: End\n" << endl;
3373  runStatus[0] = 0;
3374 
3375  return;
3376 }
3377 
3381 void fsifoam_module::StepFoam(){
3382 
3383 
3384  Foam::Time &runTime(RunTime());
3385 
3386  Info << "\nFsiFoam:StepFoam: Stepping time loop\n" << endl;
3387 
3388  // while(!runTime.end()){
3389  if(!runTime.end()){
3390  Info << "FsiFoam:StepFoam: Time = " << runTime.timeName() << nl << endl;
3391  Step();
3392  Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
3393  << " ClockTime = " << runTime.elapsedClockTime() << " s"
3394  << endl << endl;
3395  if (runTime.outputTime())
3396  Dump();
3397 
3398  // update IMPACT time
3399  UpdateTime();
3400  UpdateFSISurfaceData();
3401  UpdateFSISurfaceMesh();
3402  Info << "FsiFoam:StepFoam: End of time step." << endl;
3403  /*
3404  Info << "FsiFoam:StepFoam: numPointSurface = " << numPointsSurface << endl;
3405  Info << "FsiFoam:StepFoam: numElementSurface = " << numElementsSurface << endl;
3406  Info << "FsiFoam:StepFoam: SurfaceCoordinates = " << endl;
3407  for (int i = 0; i < numPointsSurface; i++)
3408  {
3409  Info << surfaceCoordinates[i*3] << " "
3410  << surfaceCoordinates[i*3+1] << " "
3411  << surfaceCoordinates[i*3+2] << endl;
3412  }
3413  */
3414  }
3415  if(runTime.end()){
3416  Info<< "End\n" << endl;
3417  }
3418  runStatus[0] = 0;
3419  return;
3420 }
3421 
3425 void fsifoam_module::StepFluid(){
3426 
3427  // Probe output file and preps
3428  std::stringstream ss;
3429  std::string tmpFname;
3430  ofstream surfDmpFile;
3431  dynamicFvMesh &fluidsMesh(this->FluidMesh());
3432  pointField foamCoords = fluidsMesh.allPoints();
3433  // check if needs to probe
3434  int *setsProb = NULL;
3435  int setsProb_handle;
3436  int *probProcId, *probNdeId;
3437  bool recordProb = false;
3438  std::string dataItemName(my_window_name+".setsProb");
3439  setsProb_handle = COM_get_dataitem_handle(dataItemName.c_str());
3440  if (setsProb_handle) {
3441  COM_get_array(dataItemName.c_str(), 0, &setsProb);
3442  if (*setsProb) {
3443  dataItemName = my_window_name+".probProcId";
3444  COM_get_array(dataItemName.c_str(), 0, &probProcId);
3445  dataItemName = my_window_name+".probNdeId";
3446  COM_get_array(dataItemName.c_str(), 0, &probNdeId);
3447  recordProb = (*probProcId==rank) && (*probNdeId>=0 && *probNdeId<numPointsSurface);
3448  std::cout << "****** Rank " << rank << " probe state = " << recordProb << std::endl;
3449  }
3450  }
3451 
3452  Foam::Time &runTime(RunTime());
3453  Info << "FsiFoam:StepFluid: Stepping fluid solver..." << endl;
3454 
3455  // while(!runTime.end()){
3456  if(!runTime.end()){
3457  Info << "FsiFoam:StepFluid: Time = " << runTime.timeName() << endl;
3458 
3459  // Non-iterative stepping
3460  StepFluidNonItr();
3461 
3462  Info<< "FsiFoam:StepFluid: ExecutionTime = " << runTime.elapsedCpuTime() << " s"
3463  << " ClockTime = " << runTime.elapsedClockTime() << " s"
3464  << endl;
3465  if (runTime.outputTime())
3466  Dump();
3467 
3468  // update IMPACT data
3469  UpdateTime();
3470  UpdateFSISurfaceData(); // updates tractions and pressures
3471  UpdateFSISurfaceMesh(); // updates surface mesh
3472 
3473  // update probe
3474  if(recordProb){
3475  ss.clear();
3476  ss.str("");
3477  ss << "probe_proc_"
3478  << rank
3479  << "_nde_"
3480  << *probNdeId
3481  << ".dat";
3482  tmpFname = ss.str();
3483  surfDmpFile.open(tmpFname.c_str(), std::ios::app);
3484  //surfDmpFile << "numPointSurface = " << numPointsSurface << "\n";
3485  //surfDmpFile << "numElementSurface = " << numElementsSurface << "\n";
3486  //surfDmpFile << "SurfaceCoordinates = " << "\n";
3487  int i = *probNdeId;
3488  //for (int i = 0; i < numPointsSurface; i++)
3489  //{
3490  surfDmpFile << runTime.timeName() << " "
3491  << surfaceCoordinates[i*3] << " "
3492  << surfaceCoordinates[i*3+1] << " "
3493  << surfaceCoordinates[i*3+2] << "\n";
3494  //}
3495  surfDmpFile.close();
3496  }
3497 
3498  }
3499  if(runTime.end()){
3500  Info<< "FsiFoam:StepFluid: End\n" << endl;
3501  }
3502  runStatus[0] = 0;
3503  return;
3504 }
3505 
3506 /*
3507 void fsifoam_module::ModifyEndTime(const double &endTime){
3508  RunTime().setEndTime(endTime);
3509 }
3510 */
3511 
3512 void fsifoam_module::CreateFSISurfaceMesh(){
3513 
3514  dynamicFvMesh &fluidsMesh(this->FluidMesh());
3515  numPointsSurface = fluidsMesh.boundaryMesh()[fluidPatchID].nPoints();
3516  //numElementsSurface = fluidsMesh.boundaryMesh()[fluidPatchID].faceCentres().nPoints();
3517 
3518  // face storage
3519  const faceList& faces = fluidsMesh.faces();
3520  const labelList& owner = fluidsMesh.faceOwner();
3521  const polyBoundaryMesh& boundaryMesh = fluidsMesh.boundaryMesh();
3522  const label start = boundaryMesh[fluidPatchID].start();
3523  const label end = start + boundaryMesh[fluidPatchID].size();
3524  numElementsSurface = end - start;
3525 
3526  // point storage
3527  const pointField points = fluidsMesh.points();
3528 
3529  // node 1 and element 1 live at index zero
3530  int nodeNumber = 1;
3531  for (label faceI = start; faceI < end; ++faceI) {
3532  const labelList& vrtList = faces[faceI];
3533  forAll(vrtList, i) {
3534  const label& vertex = vrtList[i];
3535 
3536  // create forward and reverse mapping for this vertex if not already in map
3537  if (!surfaceNodeMap[vertex]) {
3538  surfaceNodeMap[vertex] = nodeNumber;
3539  interfaceToFoamNodeMap[nodeNumber] = vertex;
3540  nodeNumber++;
3541  // add to node list
3542  surfaceCoordinates.push_back(points[vertex].x());
3543  surfaceCoordinates.push_back(points[vertex].y());
3544  surfaceCoordinates.push_back(points[vertex].z());
3545  }
3546  // build local element connectivity
3547  surfaceConnectivity.push_back(vertex);
3548  }
3549  }
3550 
3551  // build map for OpenFoam global vertices to OpenFoam patch vertices
3552  const labelList& meshPoints =
3553  fluidsMesh.boundaryMesh()[fluidPatchID].meshPoints();
3554  for (int i=0; i<numPointsSurface; ++i) {
3555  foamGlobalToPatchNodeMap[meshPoints[i]] = i;
3556  }
3557 
3558  // translate OpenFOAM global node numbers to local surface node numbers
3559  // node 1 and element 1 live at index zero
3560  for (unsigned int i=0; i<surfaceConnectivity.size(); ++i) {
3561  surfaceConnectivity[i] = surfaceNodeMap[surfaceConnectivity[i]];
3562  }
3563 }
3564 
3565 void fsifoam_module::UpdateFSISurfaceMesh(){
3566 
3567  dynamicFvMesh &fluidsMesh(this->FluidMesh());
3568  pointField foamCoords = fluidsMesh.allPoints();
3569 
3570  for (int i=0; i<numPointsSurface; ++i) {
3571  int vertex = interfaceToFoamNodeMap[i+1];
3572  surfaceCoordinates[3*i] = foamCoords[vertex].x();
3573  surfaceCoordinates[3*i+1] = foamCoords[vertex].y();
3574  surfaceCoordinates[3*i+2] = foamCoords[vertex].z();
3575  }
3576 }
3577 
3578 void fsifoam_module::UpdateTime(){
3579  Foam::Time &runTime(this->RunTime());
3580  time[0] = runTime.value();
3581  endTime[0] = runTime.endTime().value();
3582 }
3583 
3584 int fsifoam_module::UpdateFSISurface(Foam::vectorField &solidDispl){
3585  //Masoud : Testing displacement
3586  Info << "FsiFoam:UpdateFSISurface: Reporting solid displacements to the fluid solver" << endl;
3587  //std::cout<< "Displacements passed to me : " << std::endl;
3588  // Masoud : End
3590  // final displacement of the fluid surface
3592  for(int i=0; i<numPointsSurface; ++i) {
3593  solidDispl[foamGlobalToPatchNodeMap[interfaceToFoamNodeMap[i+1]]].x() =
3594  solidDisplacement[3*i];
3595  solidDispl[foamGlobalToPatchNodeMap[interfaceToFoamNodeMap[i+1]]].y() =
3596  solidDisplacement[3*i+1];
3597  solidDispl[foamGlobalToPatchNodeMap[interfaceToFoamNodeMap[i+1]]].z() =
3598  solidDisplacement[3*i+2];
3599  //Masoud : Testing displacement
3600  //std::cout<< solidDisplacement[3*i] << " "
3601  // << solidDisplacement[3*i+1] << " "
3602  // << solidDisplacement[3*i+2] << std::endl;
3603  //Masoud: End
3604  }
3605  return(0);
3606 }
3607 
3608 void fsifoam_module::UpdateFSISurfaceData(){
3609 
3610  //Masoud
3611  double currTime((this->RunTime()).value());
3612  Info << "FsiFoam:UpdateFSISurface: Reporting surface loads to the solid solver " << endl;
3613 
3615  // pressure on the surface
3617  dimensionedScalar &rhoFluid(this->rhoFluid());
3618  volScalarField &p(this->p());
3619 
3620  scalarField foamSurfacePressure = rhoFluid.value()*p.boundaryField()[fluidPatchID];
3621 
3622  for(int i=0; i<numElementsSurface; ++i) {
3623  surfacePressure[i] = foamSurfacePressure[i];
3624  }
3625 
3627  // traction on the surface
3629  dimensionedScalar &nu(this->nu());
3630  volVectorField &U(this->U());
3631  // Original
3632  //vectorField foamSurfaceTraction = -rhoFluid.value()*nu.value()*
3633  // U.boundaryField()[fluidPatchID].snGrad();
3634  // Original end
3635 
3636  // Masoud change: based on following taken from StepFoam()
3637  // vectorField fluidPatchTraction =
3638  // -rhoFluid.value()*nu.value()
3639  // *U.boundaryField()[fluidPatchID].snGrad()
3640  // + rhoFluid.value()*p.boundaryField()[fluidPatchID]
3641  // *fluidsMesh.boundary()[fluidPatchID].nf();
3642 
3643  dynamicFvMesh &fluidsMesh(this->FluidMesh());
3644  vectorField foamSurfaceTraction = -rhoFluid.value()*nu.value()
3645  *U.boundaryField()[fluidPatchID].snGrad()
3646  + rhoFluid.value()*p.boundaryField()[fluidPatchID]
3647  *fluidsMesh.boundary()[fluidPatchID].nf();
3648  //Masoud: end
3649 
3650  //Masoud: sending tractions after some initial time
3651  if (currTime > -0.1) {
3652  for(int i=0; i<numElementsSurface; ++i) {
3653  surfaceTraction[3*i] = foamSurfaceTraction[i].x();
3654  surfaceTraction[3*i+1] = foamSurfaceTraction[i].y();
3655  surfaceTraction[3*i+2] = foamSurfaceTraction[i].z();
3656  }
3657  } else {
3658  for(int i=0; i<numElementsSurface; ++i) {
3659  surfaceTraction[3*i] = 0.0;
3660  surfaceTraction[3*i+1] = 0.0;
3661  surfaceTraction[3*i+2] = 0.0;
3662  }
3663  }
3664  //Masoud: end
3665 }
3666 
3670 extern "C" void OpenFoamFSIPar_load_module( const char *name) {
3671  fsifoam_module::Load(name);
3672 }
3676 extern "C" void OpenFoamFSIPar_unload_module( const char *name){
3677  fsifoam_module::Unload(name);
3678 }
3679 
fsifoam_module::my_window_comm
MPI_Comm my_window_comm
Tracks this window name.
Definition: fsiFOAMPar.H:42
fsifoam_module::UpdateFSISurfaceData
void UpdateFSISurfaceData()
Update the data registered on the FSI surface.
Definition: fsiFOAM.C:3514
fsifoam_module::StepFoam
void StepFoam()
the OpenFOAM stepping
Definition: fsiFOAM.C:3297
fsifoam_module::StepFluidItr
int StepFluidItr()
Definition: fsiFOAM.C:2674
OpenFoamFSIPar_unload_module
void OpenFoamFSIPar_unload_module(const char *name)
C/C++ bindings to unload IcoFoamModule.
Definition: fsiFOAMPar.C:3684
fsifoam_module::InitFoam
void InitFoam(int *pargc, void **pargv, int *verbIn)
function to register through IMPACT
Definition: fsiFOAM.C:3184
fsifoam_module
Definition: fsiFOAM.H:38
fsifoam_module::InitFluidMesh
int InitFluidMesh()
Definition: fsiFOAM.C:52
fsifoam_module::UpdateFSISurface
int UpdateFSISurface(Foam::vectorField &solidPatchPointsDispl)
Update the nodal coordinates of the IMAPCT and OpenFoam FSI surfaces from the IMPACT displacement dat...
Definition: fsiFOAM.C:3489
fsifoam_module::UpdateTime
void UpdateTime()
Update the time control data registered with OpenFOAM i.e.
Definition: fsiFOAM.C:3483
fsifoam_module::Load
static void Load(const std::string &name)
&quot;Loads&quot; IcoFoamModule
Definition: fsiFOAM.C:3137
fsifoam_module::Unload
static void Unload(const std::string &name)
Unloads the IcoFoamModule.
Definition: fsiFOAM.C:3257
fsifoam_module::ReadPISOControls
int ReadPISOControls()
Definition: fsiFOAM.C:629
fsifoam_module::Dump
int Dump()
Definition: fsiFOAM.C:3108
fsifoam_module::RunFoam
void RunFoam()
the OpenFOAM main
Definition: fsiFOAM.C:3271
fsifoam_module::nproc
int nproc
Definition: fsiFOAMPar.H:47
fsifoam_module::Initialize
int Initialize(int argc, char *argv[])
Definition: fsiFOAM.C:5
fsifoam_module::rank
int rank
Tracks window&#39;s communicator.
Definition: fsiFOAMPar.H:46
fsifoam_module::my_window_name
std::string my_window_name
Definition: fsiFOAM.H:41
fsifoam_module::ReadFSIControls
int ReadFSIControls()
Definition: fsiFOAM.C:642
fsifoam_module::StepFluid
void StepFluid()
the OpenFOAM stepping fluid alone
Definition: fsiFOAM.C:3344
fsifoam_module::UpdateFSISurfaceMesh
void UpdateFSISurfaceMesh()
Update of the surface mesh data coordinates.
Definition: fsiFOAM.C:3470
fsifoam_module::CreateInterZoneInterpolators
int CreateInterZoneInterpolators()
Definition: fsiFOAM.C:350
fsifoam_module::StepFluidNonItr
int StepFluidNonItr()
Definition: fsiFOAM.C:2310
fsifoam_module::CreateFSISurfaceMesh
void CreateFSISurfaceMesh()
Initialization of the surface mesh data structures from the OpenFOAM mesh data structures.
Definition: fsiFOAM.C:3417
fsifoam_module::InitTransportProperties
int InitTransportProperties()
Definition: fsiFOAM.C:68
fsifoam_module::ReadCouplingProperties
int ReadCouplingProperties()
Definition: fsiFOAM.C:204
fsifoam_module::Step
int Step()
Definition: fsiFOAM.C:703
fsifoam_module::InitStructuresMesh
int InitStructuresMesh()
Definition: fsiFOAM.C:127
fsifoam_module::CreateFluidFields
int CreateFluidFields()
Definition: fsiFOAM.C:80
fsifoam_module::CreateStructuresFields
int CreateStructuresFields()
Definition: fsiFOAM.C:133
fsifoam_module::FindGlobalFaceZones
int FindGlobalFaceZones()
Definition: fsiFOAM.C:514
fsifoam_module::StepFluidAlone
int StepFluidAlone()
Definition: fsiFOAM.C:1854
OpenFoamFSIPar_load_module
void OpenFoamFSIPar_load_module(const char *name)
C/C++ bindings to load IcoFoamModule.
Definition: fsiFOAMPar.C:3678