Checkpointing does not work for subcycling

[derekrisseeuw]

For implicit coupling the writeCheckpoint is called by the function:
Precice_WriteIterationCheckPoint

This function takes the values of the last timestep as checkpoint values. However, for implicit coupled simulations where the calculix is subcycling, the last calculix timestep is not the same as the last coupling timestep. Therefore, the wrong checkpoint is written.

[uekerman]

Thanks for pointing us to this. I guess that this has simply never been implemented.

[derekrisseeuw]

Thank you for your reaction Benjamin.

This problem can be easily overcome in the calculix input file by specifying:
*DYNAMIC,DIRECT

Where the 'direct' keyword forces calculix to use the timestep specified (this should be the precice timestep).

[nkr0]

Doesn't putting write/read checkpoints inside if *required prevent this? My understanding is that the solver has to be at one of the ends of a coupling step for these to return true.

If that is the case, CalculiX sub-cycles by storing states in vini and vold without touching simulationData->coupling_init_v, which is only written/read at the beginning/end of coupling steps. precice resets vold to the beginning of the coupling time step, when Precice_IsReadCheckpointRequired returns true, i.e., only when the coupling step is not converged at the end of a coupling step (not at the end of CalculiX sub-cycles). Similarly, Precice_IsWriteCheckpointRequired returns true only at the beginning of a coupling time step to store a checkpoint in simulationData->coupling_init_v.

calculix-adapter/nonlingeo_precice.c

Lines 1115 to 1121 in b01641e

	memcpy(&vini[0],&vold[0],sizeof(double)*mt**nk);
	if( Precice_IsWriteCheckpointRequired() )
	{
	Precice_WriteIterationCheckpoint( &simulationData, vini );
	Precice_FulfilledWriteCheckpoint();
	}

calculix-adapter/nonlingeo_precice.c

Lines 2664 to 2678 in b01641e

	memcpy(&vold[0],&v[0],sizeof(double)*mt**nk);
	Precice_WriteCouplingData( &simulationData );
	/* Adapter: Advance the coupling */
	Precice_Advance( &simulationData );
	/* Adapter: If the coupling does not converge, read the checkpoint */
	if( Precice_IsReadCheckpointRequired() )
	{
	if( *nmethod == 4 )
	{
	Precice_ReadIterationCheckpoint( &simulationData, vold );
	icutb++;
	}
	Precice_FulfilledReadCheckpoint();
	}

Read/WriteCouplingData, AdjustSolverTimestep, and Advance executes every sub-cycle, to

calculix-adapter/nonlingeo_precice.c

Lines 1097 to 1099 in b01641e

	Precice_AdjustSolverTimestep( &simulationData );
	/* Adapter read coupling data if available */
	Precice_ReadCouplingData( &simulationData );

calculix-adapter/nonlingeo_precice.c

Line 2666 in b01641e

Precice_WriteCouplingData( &simulationData );

read/write simulationData->preciceInterfaces,

calculix-adapter/adapter/PreciceInterface.c

Line 162 in b01641e

PreciceInterface ** interfaces = sim->preciceInterfaces;

calculix-adapter/adapter/PreciceInterface.c

Line 229 in b01641e

PreciceInterface ** interfaces = sim->preciceInterfaces;

set simulationData->dtheta (to keep calculix within the coupling step),

calculix-adapter/adapter/PreciceInterface.c

Line 84 in b01641e

*sim->dtheta = fmin( sim->precice_dt / *sim->tper, *sim->dtheta );

and increments simulationData->precice_dt (to track progress within a coupling step), respectively.

calculix-adapter/adapter/PreciceInterface.c

Line 96 in b01641e

sim->precice_dt = precicec_advance( sim->solver_dt );

[nkr0]

Furthermore, CalculiX's check for sub-cycle convergence is icutb==0.

calculix-adapter/nonlingeo_precice.c

Lines 2683 to 2725 in b01641e

	/* printing the energies (only for dynamic calculations) */
	if((icutb==0)&&(*nmethod==4)&&(*ithermal<2)){
	printf(" initial energy (at start of step) = %e\n\n",energyref);
	printf(" since start of the step: \n");
	printf(" external work = %e\n",allwk);
	printf(" work performed by the damping forces = %e\n",dampwk);
	printf(" netto work = %e\n\n",allwk+dampwk);
	printf(" actual energy: \n");
	printf(" internal energy = %e\n",energy[0]);
	printf(" kinetic energy = %e\n",energy[1]);
	printf(" elastic contact energy = %e\n",energy[2]);
	printf(" energy lost due to friction = %e\n",energy[3]);
	printf(" total energy = %e\n\n",energy[0]+energy[1]+energy[2]+energy[3]);
	printf(" energy increase = %e\n\n",energy[0]+energy[1]+energy[2]+energy[3]-energyref);
	printf(" energy balance (absolute) = %e \n",energy[0]+energy[1]+energy[2]+energy[3]-energyref-allwk-dampwk);
	/* Belytschko criterion */
	denergymax=energy[0];
	if(denergymax<energy[1]) denergymax=energy[1];
	if(denergymax<fabs(allwk)) denergymax=fabs(allwk);
	if(denergymax>1.e-30){
	printf(" energy balance (relative) = %f %% \n\n",
	fabs((energy[0]+energy[1]+energy[2]+energy[3]-energyref-allwk-dampwk)/
	denergymax*100.));
	}else{
	printf(" energy balance (relative) =0 %% \n\n");
	}
	// # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
	// MPADD start
	// printf(" work done by the damping forces = %e\n", dampwk);
	// neini=*ne;
	// printf(" contact elements end of increment = %"ITGFORMAT"\n\n", *ne - ne0);
	// MPADD end
	// # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
	}

This controls memcpy between vini and vold (vini is Calculix's sub-cycle "checkpoint") and precice checkpoints are also under this check, making sure precice checkpoints are only read/set when calculix sub-cycles are converged.

calculix-adapter/nonlingeo_precice.c

Lines 1101 to 1121 in b01641e

	if(icutb==0){
	/* previous increment converged: update the initial values */
	iinc++;
	jprint++;
	/* store number of elements (important for implicit dynamic
	contact */
	neini=*ne;
	/* vold is copied into vini */
	memcpy(&vini[0],&vold[0],sizeof(double)*mt**nk);
	if( Precice_IsWriteCheckpointRequired() )
	{
	Precice_WriteIterationCheckpoint( &simulationData, vini );
	Precice_FulfilledWriteCheckpoint();
	}

calculix-adapter/nonlingeo_precice.c

Lines 2628 to 2681 in b01641e

	if( icutb == 0 )
	{
	/* Adapter: Write coupling data */
	NNEW(v,double,mt**nk);
	NNEW(fn,double,mt**nk);
	NNEW(stn,double,6**nk);
	NNEW(stx,double,6*mi[0]**ne);
	NNEW(inum,ITG,*nk);
	memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
	// iout=-1 means that the displacements and temperatures are assumed to be known and used to calculate strains, stresses...., with no result output
	iout=-1;
	icmd=3;// calculate only stress (not stiffness)
	results(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
	elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
	ielorien,norien,orab,ntmat_,t0,t1act,ithermal,
	prestr,iprestr,filab,eme,emn,een,iperturb,
	f,fn,nactdof,&iout,qa,vold,b,nodeboun,
	ndirboun,xbounact,nboun,ipompc,
	nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
	&bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
	xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
	&icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,
	emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
	iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
	fmpc,nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
	&reltime,&ne0,thicke,shcon,nshcon,
	sideload,xloadact,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
	mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
	islavsurf,ielprop,prop,energyini,energy,&kscale,iponoel,
	inoel,nener,orname,network,ipobody,xbodyact,ibody,typeboun);
	simulationData.fn = fn;
	memcpy(&vold[0],&v[0],sizeof(double)*mt**nk);
	Precice_WriteCouplingData( &simulationData );
	/* Adapter: Advance the coupling */
	Precice_Advance( &simulationData );
	/* Adapter: If the coupling does not converge, read the checkpoint */
	if( Precice_IsReadCheckpointRequired() )
	{
	if( *nmethod == 4 )
	{
	Precice_ReadIterationCheckpoint( &simulationData, vold );
	icutb++;
	}
	Precice_FulfilledReadCheckpoint();
	}
	SFREE(v);SFREE(stn);SFREE(stx);SFREE(fn);SFREE(inum);
	}

You can see in the blob above that when precice reads a checkpoint, we increment icutb to unset CalculiX's converged state.

[MakisH]

I assume that the branch https://github.com/precice/calculix-adapter/tree/i9 refers to this issue. Should we still keep it around? Is this issue still relevant?

[MakisH]

This issue is still valid with the adapter 2.19.0. Easy to reproduce: get the perpendicular-flap tutorial and, in flap.inp, replace:

*DYNAMIC, ALPHA=0.0, DIRECT
- 1.E-2, 5.0
+ 3.E-3, 5.0

This makes the OpenFOAM side crash due to a mesh failure:

We definitely need to at least document this.

[boris-martin]

I've been toying around to fix this, but this opens new questions.
Roughly speaking, in terms of data, it should be sufficient to have checkpoints "go back deep enough". I have prototypes who seem to go into the correct direction.
However I'm not sure how to handle outputs: by default, every sub-step is recorded. If we subcycle, we might get steps rewritten may times.
The only way I see it would be to add a condition like "if implicit coupling is used, output frequency must match exactly the time window size (or every N windows could work too)". If we do that, should we override manual frequency control ? Just warn ?
This would probably mean that switching between explicit and implicit will mess the output frequency, which is a bit weird.

[MakisH]

Roughly speaking, in terms of data, it should be sufficient to have checkpoints "go back deep enough".

Important to ensure here, is that checkpoints should only be written on every coupling time window (and every time a coupling time window is repeated), not every solver time step. When reading a checkpoint, CalculiX should go back to the state of the last coupling time window.

However I'm not sure how to handle outputs: by default, every sub-step is recorded. If we subcycle, we might get steps rewritten may times.

This is something I am not completely sure how it works also in OpenFOAM at the moment. If I remember correctly, results get overwritten, and this is fine.

The only way I see it would be to add a condition like "if implicit coupling is used, output frequency must match exactly the time window size (or every N windows could work too)". If we do that, should we override manual frequency control ? Just warn ?

We could always give a warning (or error) till we fix the actual problem.

This would probably mean that switching between explicit and implicit will mess the output frequency, which is a bit weird.

I am not sure I understand this point.

[boris-martin]

Roughly speaking, in terms of data, it should be sufficient to have checkpoints "go back deep enough".

Important to ensure here, is that checkpoints should only be written on every coupling time window (and every time a coupling time window is repeated), not every solver time step. When reading a checkpoint, CalculiX should go back to the state of the last coupling time window.

Yes. From what I understood the issue lies mostly in the fact that not sufficient data is saved but since we check that writing checkpoints is required, this should be fine. We mostly need to upgrade to reset correctly the internal time, step number etc, as currently only the DOFs are saved.

However I'm not sure how to handle outputs: by default, every sub-step is recorded. If we subcycle, we might get steps rewritten may times.

This is something I am not completely sure how it works also in OpenFOAM at the moment. If I remember correctly, results get overwritten, and this is fine.

I'll have to check if CalculiX can overwrite. Since it's a unique file, I'm not so sure, it's not like you can write many times the file "output_step_37" many times.

This would probably mean that switching between explicit and implicit will mess the output frequency, which is a bit weird.

I am not sure I understand this point.

Imagine a time window size of 1.0 and an internal step of 0.4.

• In explicit coupling we'd have writes at 0, 0.4, 0.8, 1.0, 1.4, 1.8, 2.0, etc. (Which is fine, except for the fact that in the output VTK the difference in step is not visible)
• In implicit coupling, if we modify to only write at the end of a coupling window because we can't solve the overwriting issue etc, we'd have writes at 0, 1.0, 2.0, ... That's probably not the worst option, but having a behavior dependent on the coupling kind is not ideal.